1 static const char help[] = "Tries to solve u`` + u^{2} = f for an easy case and an impossible case.\n\n"; 2 3 /* 4 This example was contributed by Peter Graf to show how SNES fails when given a nonlinear problem with no solution. 5 6 Run with -n 14 to see fail to converge and -n 15 to see convergence 7 8 The option -second_order uses a different discretization of the Neumann boundary condition and always converges 9 10 */ 11 12 #include <petscsnes.h> 13 14 PetscBool second_order = PETSC_FALSE; 15 #define X0DOT -2.0 16 #define X1 5.0 17 #define KPOW 2.0 18 const PetscScalar sperturb = 1.1; 19 20 /* 21 User-defined routines 22 */ 23 PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *); 24 PetscErrorCode FormFunction(SNES, Vec, Vec, void *); 25 26 int main(int argc, char **argv) 27 { 28 SNES snes; /* SNES context */ 29 Vec x, r, F; /* vectors */ 30 Mat J; /* Jacobian */ 31 PetscInt it, n = 11, i; 32 PetscReal h, xp = 0.0; 33 PetscScalar v; 34 const PetscScalar a = X0DOT; 35 const PetscScalar b = X1; 36 const PetscScalar k = KPOW; 37 PetscScalar v2; 38 PetscScalar *xx; 39 40 PetscFunctionBeginUser; 41 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 42 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 43 PetscCall(PetscOptionsGetBool(NULL, NULL, "-second_order", &second_order, NULL)); 44 h = 1.0 / (n - 1); 45 46 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 47 Create nonlinear solver context 48 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 49 50 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 51 52 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 53 Create vector data structures; set function evaluation routine 54 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 55 56 PetscCall(VecCreate(PETSC_COMM_SELF, &x)); 57 PetscCall(VecSetSizes(x, PETSC_DECIDE, n)); 58 PetscCall(VecSetFromOptions(x)); 59 PetscCall(VecDuplicate(x, &r)); 60 PetscCall(VecDuplicate(x, &F)); 61 62 PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)F)); 63 64 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 65 Create matrix data structures; set Jacobian evaluation routine 66 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 67 68 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, &J)); 69 70 /* 71 Note that in this case we create separate matrices for the Jacobian 72 and preconditioner matrix. Both of these are computed in the 73 routine FormJacobian() 74 */ 75 /* PetscCall(SNESSetJacobian(snes,NULL,JPrec,FormJacobian,0)); */ 76 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, 0)); 77 /* PetscCall(SNESSetJacobian(snes,J,JPrec,FormJacobian,0)); */ 78 79 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 80 Customize nonlinear solver; set runtime options 81 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 82 83 PetscCall(SNESSetFromOptions(snes)); 84 85 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 86 Initialize application: 87 Store right-hand side of PDE and exact solution 88 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 89 90 /* set right-hand side and initial guess to be exact solution of continuum problem */ 91 #define SQR(x) ((x) * (x)) 92 xp = 0.0; 93 for (i = 0; i < n; i++) { 94 v = k * (k - 1.) * (b - a) * PetscPowScalar(xp, k - 2.) + SQR(a * xp) + SQR(b - a) * PetscPowScalar(xp, 2. * k) + 2. * a * (b - a) * PetscPowScalar(xp, k + 1.); 95 PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES)); 96 v2 = a * xp + (b - a) * PetscPowScalar(xp, k); 97 PetscCall(VecSetValues(x, 1, &i, &v2, INSERT_VALUES)); 98 xp += h; 99 } 100 101 /* perturb initial guess */ 102 PetscCall(VecGetArray(x, &xx)); 103 for (i = 0; i < n; i++) { 104 v2 = xx[i] * sperturb; 105 PetscCall(VecSetValues(x, 1, &i, &v2, INSERT_VALUES)); 106 } 107 PetscCall(VecRestoreArray(x, &xx)); 108 109 PetscCall(SNESSolve(snes, NULL, x)); 110 PetscCall(SNESGetIterationNumber(snes, &it)); 111 PetscCall(PetscPrintf(PETSC_COMM_SELF, "SNES iterations = %" PetscInt_FMT "\n\n", it)); 112 113 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 114 Free work space. All PETSc objects should be destroyed when they 115 are no longer needed. 116 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 117 118 PetscCall(VecDestroy(&x)); 119 PetscCall(VecDestroy(&r)); 120 PetscCall(VecDestroy(&F)); 121 PetscCall(MatDestroy(&J)); 122 PetscCall(SNESDestroy(&snes)); 123 PetscCall(PetscFinalize()); 124 return 0; 125 } 126 127 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *dummy) 128 { 129 const PetscScalar *xx; 130 PetscScalar *ff, *FF, d, d2; 131 PetscInt i, n; 132 133 PetscFunctionBeginUser; 134 PetscCall(VecGetArrayRead(x, &xx)); 135 PetscCall(VecGetArray(f, &ff)); 136 PetscCall(VecGetArray((Vec)dummy, &FF)); 137 PetscCall(VecGetSize(x, &n)); 138 d = (PetscReal)(n - 1); 139 d2 = d * d; 140 141 if (second_order) ff[0] = d * (0.5 * d * (-xx[2] + 4. * xx[1] - 3. * xx[0]) - X0DOT); 142 else ff[0] = d * (d * (xx[1] - xx[0]) - X0DOT); 143 144 for (i = 1; i < n - 1; i++) ff[i] = d2 * (xx[i - 1] - 2. * xx[i] + xx[i + 1]) + xx[i] * xx[i] - FF[i]; 145 146 ff[n - 1] = d * d * (xx[n - 1] - X1); 147 PetscCall(VecRestoreArrayRead(x, &xx)); 148 PetscCall(VecRestoreArray(f, &ff)); 149 PetscCall(VecRestoreArray((Vec)dummy, &FF)); 150 PetscFunctionReturn(PETSC_SUCCESS); 151 } 152 153 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat prejac, void *dummy) 154 { 155 const PetscScalar *xx; 156 PetscScalar A[3], d, d2; 157 PetscInt i, n, j[3]; 158 159 PetscFunctionBeginUser; 160 PetscCall(VecGetSize(x, &n)); 161 PetscCall(VecGetArrayRead(x, &xx)); 162 d = (PetscReal)(n - 1); 163 d2 = d * d; 164 165 i = 0; 166 if (second_order) { 167 j[0] = 0; 168 j[1] = 1; 169 j[2] = 2; 170 A[0] = -3. * d * d * 0.5; 171 A[1] = 4. * d * d * 0.5; 172 A[2] = -1. * d * d * 0.5; 173 PetscCall(MatSetValues(prejac, 1, &i, 3, j, A, INSERT_VALUES)); 174 } else { 175 j[0] = 0; 176 j[1] = 1; 177 A[0] = -d * d; 178 A[1] = d * d; 179 PetscCall(MatSetValues(prejac, 1, &i, 2, j, A, INSERT_VALUES)); 180 } 181 for (i = 1; i < n - 1; i++) { 182 j[0] = i - 1; 183 j[1] = i; 184 j[2] = i + 1; 185 A[0] = d2; 186 A[1] = -2. * d2 + 2. * xx[i]; 187 A[2] = d2; 188 PetscCall(MatSetValues(prejac, 1, &i, 3, j, A, INSERT_VALUES)); 189 } 190 191 i = n - 1; 192 A[0] = d * d; 193 PetscCall(MatSetValues(prejac, 1, &i, 1, &i, &A[0], INSERT_VALUES)); 194 195 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 196 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 197 PetscCall(MatAssemblyBegin(prejac, MAT_FINAL_ASSEMBLY)); 198 PetscCall(MatAssemblyEnd(prejac, MAT_FINAL_ASSEMBLY)); 199 200 PetscCall(VecRestoreArrayRead(x, &xx)); 201 PetscFunctionReturn(PETSC_SUCCESS); 202 } 203 204 /*TEST 205 206 test: 207 args: -n 14 -snes_monitor_short -snes_converged_reason 208 requires: !single 209 210 test: 211 suffix: 2 212 args: -n 15 -snes_monitor_short -snes_converged_reason 213 requires: !single 214 215 test: 216 suffix: 3 217 args: -n 14 -second_order -snes_monitor_short -snes_converged_reason 218 requires: !single 219 220 TEST*/ 221