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 matrix used to compute the preconditioner. 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