Actual source code: ex59.c


  2: static const char help[] = "Tries to solve u`` + u^{2} = f for an easy case and an impossible case.\n\n";

  4: /*
  5:        This example was contributed by Peter Graf to show how SNES fails when given a nonlinear problem with no solution.

  7:        Run with -n 14 to see fail to converge and -n 15 to see convergence

  9:        The option -second_order uses a different discretization of the Neumann boundary condition and always converges

 11: */

 13: #include <petscsnes.h>

 15: PetscBool second_order = PETSC_FALSE;
 16: #define X0DOT -2.0
 17: #define X1    5.0
 18: #define KPOW  2.0
 19: const PetscScalar sperturb = 1.1;

 21: /*
 22:    User-defined routines
 23: */
 24: PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
 25: PetscErrorCode FormFunction(SNES, Vec, Vec, void *);

 27: int main(int argc, char **argv)
 28: {
 29:   SNES              snes;    /* SNES context */
 30:   Vec               x, r, F; /* vectors */
 31:   Mat               J;       /* Jacobian */
 32:   PetscInt          it, n = 11, i;
 33:   PetscReal         h, xp = 0.0;
 34:   PetscScalar       v;
 35:   const PetscScalar a = X0DOT;
 36:   const PetscScalar b = X1;
 37:   const PetscScalar k = KPOW;
 38:   PetscScalar       v2;
 39:   PetscScalar      *xx;

 41:   PetscFunctionBeginUser;
 42:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 43:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 44:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-second_order", &second_order, NULL));
 45:   h = 1.0 / (n - 1);

 47:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 48:      Create nonlinear solver context
 49:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 51:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));

 53:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 54:      Create vector data structures; set function evaluation routine
 55:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 57:   PetscCall(VecCreate(PETSC_COMM_SELF, &x));
 58:   PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
 59:   PetscCall(VecSetFromOptions(x));
 60:   PetscCall(VecDuplicate(x, &r));
 61:   PetscCall(VecDuplicate(x, &F));

 63:   PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)F));

 65:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 66:      Create matrix data structures; set Jacobian evaluation routine
 67:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 69:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, &J));

 71:   /*
 72:      Note that in this case we create separate matrices for the Jacobian
 73:      and preconditioner matrix.  Both of these are computed in the
 74:      routine FormJacobian()
 75:   */
 76:   /*  PetscCall(SNESSetJacobian(snes,NULL,JPrec,FormJacobian,0)); */
 77:   PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, 0));
 78:   /*  PetscCall(SNESSetJacobian(snes,J,JPrec,FormJacobian,0)); */

 80:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 81:      Customize nonlinear solver; set runtime options
 82:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 84:   PetscCall(SNESSetFromOptions(snes));

 86:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 87:      Initialize application:
 88:      Store right-hand-side of PDE and exact solution
 89:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 91:   /* set right hand side and initial guess to be exact solution of continuim problem */
 92: #define SQR(x) ((x) * (x))
 93:   xp = 0.0;
 94:   for (i = 0; i < n; i++) {
 95:     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.);
 96:     PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES));
 97:     v2 = a * xp + (b - a) * PetscPowScalar(xp, k);
 98:     PetscCall(VecSetValues(x, 1, &i, &v2, INSERT_VALUES));
 99:     xp += h;
100:   }

102:   /* perturb initial guess */
103:   PetscCall(VecGetArray(x, &xx));
104:   for (i = 0; i < n; i++) {
105:     v2 = xx[i] * sperturb;
106:     PetscCall(VecSetValues(x, 1, &i, &v2, INSERT_VALUES));
107:   }
108:   PetscCall(VecRestoreArray(x, &xx));

110:   PetscCall(SNESSolve(snes, NULL, x));
111:   PetscCall(SNESGetIterationNumber(snes, &it));
112:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "SNES iterations = %" PetscInt_FMT "\n\n", it));

114:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115:      Free work space.  All PETSc objects should be destroyed when they
116:      are no longer needed.
117:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

119:   PetscCall(VecDestroy(&x));
120:   PetscCall(VecDestroy(&r));
121:   PetscCall(VecDestroy(&F));
122:   PetscCall(MatDestroy(&J));
123:   PetscCall(SNESDestroy(&snes));
124:   PetscCall(PetscFinalize());
125:   return 0;
126: }

128: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *dummy)
129: {
130:   const PetscScalar *xx;
131:   PetscScalar       *ff, *FF, d, d2;
132:   PetscInt           i, n;

134:   PetscFunctionBeginUser;
135:   PetscCall(VecGetArrayRead(x, &xx));
136:   PetscCall(VecGetArray(f, &ff));
137:   PetscCall(VecGetArray((Vec)dummy, &FF));
138:   PetscCall(VecGetSize(x, &n));
139:   d  = (PetscReal)(n - 1);
140:   d2 = d * d;

142:   if (second_order) ff[0] = d * (0.5 * d * (-xx[2] + 4. * xx[1] - 3. * xx[0]) - X0DOT);
143:   else ff[0] = d * (d * (xx[1] - xx[0]) - X0DOT);

145:   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];

147:   ff[n - 1] = d * d * (xx[n - 1] - X1);
148:   PetscCall(VecRestoreArrayRead(x, &xx));
149:   PetscCall(VecRestoreArray(f, &ff));
150:   PetscCall(VecRestoreArray((Vec)dummy, &FF));
151:   PetscFunctionReturn(PETSC_SUCCESS);
152: }

154: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat prejac, void *dummy)
155: {
156:   const PetscScalar *xx;
157:   PetscScalar        A[3], d, d2;
158:   PetscInt           i, n, j[3];

160:   PetscFunctionBeginUser;
161:   PetscCall(VecGetSize(x, &n));
162:   PetscCall(VecGetArrayRead(x, &xx));
163:   d  = (PetscReal)(n - 1);
164:   d2 = d * d;

166:   i = 0;
167:   if (second_order) {
168:     j[0] = 0;
169:     j[1] = 1;
170:     j[2] = 2;
171:     A[0] = -3. * d * d * 0.5;
172:     A[1] = 4. * d * d * 0.5;
173:     A[2] = -1. * d * d * 0.5;
174:     PetscCall(MatSetValues(prejac, 1, &i, 3, j, A, INSERT_VALUES));
175:   } else {
176:     j[0] = 0;
177:     j[1] = 1;
178:     A[0] = -d * d;
179:     A[1] = d * d;
180:     PetscCall(MatSetValues(prejac, 1, &i, 2, j, A, INSERT_VALUES));
181:   }
182:   for (i = 1; i < n - 1; i++) {
183:     j[0] = i - 1;
184:     j[1] = i;
185:     j[2] = i + 1;
186:     A[0] = d2;
187:     A[1] = -2. * d2 + 2. * xx[i];
188:     A[2] = d2;
189:     PetscCall(MatSetValues(prejac, 1, &i, 3, j, A, INSERT_VALUES));
190:   }

192:   i    = n - 1;
193:   A[0] = d * d;
194:   PetscCall(MatSetValues(prejac, 1, &i, 1, &i, &A[0], INSERT_VALUES));

196:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
197:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
198:   PetscCall(MatAssemblyBegin(prejac, MAT_FINAL_ASSEMBLY));
199:   PetscCall(MatAssemblyEnd(prejac, MAT_FINAL_ASSEMBLY));

201:   PetscCall(VecRestoreArrayRead(x, &xx));
202:   PetscFunctionReturn(PETSC_SUCCESS);
203: }

205: /*TEST

207:    test:
208:       args: -n 14 -snes_monitor_short -snes_converged_reason
209:       requires: !single

211:    test:
212:       suffix: 2
213:       args: -n 15 -snes_monitor_short -snes_converged_reason
214:       requires: !single

216:    test:
217:       suffix: 3
218:       args: -n 14 -second_order -snes_monitor_short -snes_converged_reason
219:       requires: !single

221: TEST*/