Actual source code: ex15.c

  1: static const char help[] = "p-Bratu nonlinear PDE in 2d.\n\
2: We solve the  p-Laplacian (nonlinear diffusion) combined with\n\
3: the Bratu (solid fuel ignition) nonlinearity in a 2D rectangular\n\
4: domain, using distributed arrays (DAs) to partition the parallel grid.\n\
5: The command line options include:\n\
6:   -p <2>: p' in p-Laplacian term\n\
7:   -epsilon <1e-05>: Strain-regularization in p-Laplacian\n\
8:   -lambda <6>: Bratu parameter\n\
9:   -blocks <bx,by>: number of coefficient interfaces in x and y direction\n\
10:   -kappa <1e-3>: diffusivity in odd regions\n\
11: \n";

13: /*F
14:     The $p$-Bratu problem is a combination of the $p$-Laplacian (nonlinear diffusion) and the Brutu solid fuel ignition problem.
15:     This problem is modeled by the partial differential equation

17: \begin{equation*}
18:         -\nabla\cdot (\eta \nabla u) - \lambda \exp(u) = 0
19: \end{equation*}

21:     on $\Omega = (-1,1)^2$ with closure

23: \begin{align*}
24:         \eta(\gamma) &= (\epsilon^2 + \gamma)^{(p-2)/2} & \gamma &= \frac 1 2 |\nabla u|^2
25: \end{align*}

27:     and boundary conditions $u = 0$ for $(x,y) \in \partial \Omega$

29:     A 9-point finite difference stencil is used to discretize
30:     the boundary value problem to obtain a nonlinear system of equations.
31:     This would be a 5-point stencil if not for the $p$-Laplacian's nonlinearity.
32: F*/

34: /*
35:       mpiexec -n 2 ./ex15 -snes_monitor -ksp_monitor log_view
36: */

38: /*
39:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
40:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
41:    file automatically includes:
42:      petsc.h       - base PETSc routines   petscvec.h - vectors
43:      petscsys.h    - system routines       petscmat.h - matrices
44:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
45:      petscviewer.h - viewers               petscpc.h  - preconditioners
46:      petscksp.h   - linear solvers
47: */

49: #include <petscdm.h>
50: #include <petscdmda.h>
51: #include <petscsnes.h>

53: typedef enum {
54:   JAC_BRATU,
55:   JAC_PICARD,
56:   JAC_STAR,
57:   JAC_NEWTON
58: } JacType;
59: static const char *const JacTypes[] = {"BRATU", "PICARD", "STAR", "NEWTON", "JacType", "JAC_", 0};

61: /*
62:    User-defined application context - contains data needed by the
63:    application-provided call-back routines, FormJacobianLocal() and
64:    FormFunctionLocal().
65: */
66: typedef struct {
67:   PetscReal lambda;  /* Bratu parameter */
68:   PetscReal p;       /* Exponent in p-Laplacian */
69:   PetscReal epsilon; /* Regularization */
70:   PetscReal source;  /* Source term */
71:   JacType   jtype;   /* What type of Jacobian to assemble */
72:   PetscBool picard;
73:   PetscInt  blocks[2];
74:   PetscReal kappa;
75:   PetscInt  initial; /* initial conditions type */
76: } AppCtx;

78: /*
79:    User-defined routines
80: */
81: static PetscErrorCode FormRHS(AppCtx *, DM, Vec);
82: static PetscErrorCode FormInitialGuess(AppCtx *, DM, Vec);
83: static PetscErrorCode FormFunctionLocal(DMDALocalInfo *, PetscScalar **, PetscScalar **, AppCtx *);
84: static PetscErrorCode FormFunctionPicardLocal(DMDALocalInfo *, PetscScalar **, PetscScalar **, AppCtx *);
85: static PetscErrorCode FormJacobianLocal(DMDALocalInfo *, PetscScalar **, Mat, Mat, AppCtx *);
86: static PetscErrorCode NonlinearGS(SNES, Vec, Vec, void *);

88: typedef struct _n_PreCheck *PreCheck;
89: struct _n_PreCheck {
90:   MPI_Comm    comm;
91:   PetscReal   angle;
92:   Vec         Ylast;
93:   PetscViewer monitor;
94: };
95: PetscErrorCode PreCheckCreate(MPI_Comm, PreCheck *);
96: PetscErrorCode PreCheckDestroy(PreCheck *);
97: PetscErrorCode PreCheckFunction(SNESLineSearch, Vec, Vec, PetscBool *, void *);
98: PetscErrorCode PreCheckSetFromOptions(PreCheck);

100: int main(int argc, char **argv)
101: {
102:   SNES                snes;       /* nonlinear solver */
103:   Vec                 x, r, b;    /* solution, residual, rhs vectors */
104:   AppCtx              user;       /* user-defined work context */
105:   PetscInt            its;        /* iterations for convergence */
106:   SNESConvergedReason reason;     /* Check convergence */
107:   PetscBool           alloc_star; /* Only allocate for the STAR stencil  */
108:   PetscBool           write_output;
109:   char                filename[PETSC_MAX_PATH_LEN] = "ex15.vts";
110:   PetscReal           bratu_lambda_max = 6.81, bratu_lambda_min = 0.;
111:   DM                  da;              /* distributed array data structure */
112:   PreCheck            precheck = NULL; /* precheck context for version in this file */
113:   PetscInt            use_precheck;    /* 0=none, 1=version in this file, 2=SNES-provided version */
114:   PetscReal           precheck_angle;  /* When manually setting the SNES-provided precheck function */
115:   SNESLineSearch      linesearch;

117:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118:      Initialize program
119:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

121:   PetscFunctionBeginUser;
122:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

124:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125:      Initialize problem parameters
126:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
127:   user.lambda    = 0.0;
128:   user.p         = 2.0;
129:   user.epsilon   = 1e-5;
130:   user.source    = 0.1;
131:   user.jtype     = JAC_NEWTON;
132:   user.initial   = -1;
133:   user.blocks[0] = 1;
134:   user.blocks[1] = 1;
135:   user.kappa     = 1e-3;
136:   alloc_star     = PETSC_FALSE;
137:   use_precheck   = 0;
138:   precheck_angle = 10.;
139:   user.picard    = PETSC_FALSE;
140:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "p-Bratu options", __FILE__);
141:   {
142:     PetscInt two = 2;
143:     PetscCall(PetscOptionsReal("-lambda", "Bratu parameter", "", user.lambda, &user.lambda, NULL));
144:     PetscCall(PetscOptionsReal("-p", "Exponent p' in p-Laplacian", "", user.p, &user.p, NULL));
145:     PetscCall(PetscOptionsReal("-epsilon", "Strain-regularization in p-Laplacian", "", user.epsilon, &user.epsilon, NULL));
146:     PetscCall(PetscOptionsReal("-source", "Constant source term", "", user.source, &user.source, NULL));
147:     PetscCall(PetscOptionsEnum("-jtype", "Jacobian approximation to assemble", "", JacTypes, (PetscEnum)user.jtype, (PetscEnum *)&user.jtype, NULL));
148:     PetscCall(PetscOptionsName("-picard", "Solve with defect-correction Picard iteration", "", &user.picard));
149:     if (user.picard) {
150:       user.jtype = JAC_PICARD;
151:       PetscCheck(user.p == 3, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Picard iteration is only supported for p == 3");
152:       /* the Picard linearization only requires a 5 point stencil, while the Newton linearization requires a 9 point stencil */
153:       /* hence allocating the 5 point stencil gives the same convergence as the 9 point stencil since the extra stencil points are not used */
154:       PetscCall(PetscOptionsBool("-alloc_star", "Allocate for STAR stencil (5-point)", "", alloc_star, &alloc_star, NULL));
155:     }
156:     PetscCall(PetscOptionsInt("-precheck", "Use a pre-check correction intended for use with Picard iteration 1=this version, 2=library", "", use_precheck, &use_precheck, NULL));
157:     PetscCall(PetscOptionsInt("-initial", "Initial conditions type (-1: default, 0: zero-valued, 1: peaked guess)", "", user.initial, &user.initial, NULL));
158:     if (use_precheck == 2) { /* Using library version, get the angle */
159:       PetscCall(PetscOptionsReal("-precheck_angle", "Angle in degrees between successive search directions necessary to activate step correction", "", precheck_angle, &precheck_angle, NULL));
160:     }
161:     PetscCall(PetscOptionsIntArray("-blocks", "number of coefficient interfaces in x and y direction", "", user.blocks, &two, NULL));
162:     PetscCall(PetscOptionsReal("-kappa", "diffusivity in odd regions", "", user.kappa, &user.kappa, NULL));
163:     PetscCall(PetscOptionsString("-o", "Output solution in vts format", "", filename, filename, sizeof(filename), &write_output));
164:   }
165:   PetscOptionsEnd();
166:   if (user.lambda > bratu_lambda_max || user.lambda < bratu_lambda_min) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "WARNING: lambda %g out of range for p=2\n", (double)user.lambda));

168:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169:      Create nonlinear solver context
170:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
171:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));

173:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174:      Create distributed array (DMDA) to manage parallel grid and vectors
175:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
176:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, alloc_star ? DMDA_STENCIL_STAR : DMDA_STENCIL_BOX, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
177:   PetscCall(DMSetFromOptions(da));
178:   PetscCall(DMSetUp(da));

180:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181:      Extract global vectors from DM; then duplicate for remaining
182:      vectors that are the same types
183:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
184:   PetscCall(DMCreateGlobalVector(da, &x));
185:   PetscCall(VecDuplicate(x, &r));
186:   PetscCall(VecDuplicate(x, &b));

188:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189:      User can override with:
190:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
191:                 (unless user explicitly sets preconditioner)
192:      -snes_mf_operator : form preconditioning matrix as set by the user,
193:                          but use matrix-free approx for Jacobian-vector
194:                          products within Newton-Krylov method

196:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

198:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199:      Set local function evaluation routine
200:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
201:   PetscCall(DMSetApplicationContext(da, &user));
202:   PetscCall(SNESSetDM(snes, da));
203:   if (user.picard) {
204:     /*
205:         This is not really right requiring the user to call SNESSetFunction/Jacobian but the DMDASNESSetPicardLocal() cannot access
206:         the SNES to set it
207:     */
208:     PetscCall(DMDASNESSetPicardLocal(da, INSERT_VALUES, (PetscErrorCode(*)(DMDALocalInfo *, void *, void *, void *))FormFunctionPicardLocal, (PetscErrorCode(*)(DMDALocalInfo *, void *, Mat, Mat, void *))FormJacobianLocal, &user));
209:     PetscCall(SNESSetFunction(snes, NULL, SNESPicardComputeFunction, &user));
210:     PetscCall(SNESSetJacobian(snes, NULL, NULL, SNESPicardComputeJacobian, &user));
211:   } else {
212:     PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (PetscErrorCode(*)(DMDALocalInfo *, void *, void *, void *))FormFunctionLocal, &user));
213:     PetscCall(DMDASNESSetJacobianLocal(da, (PetscErrorCode(*)(DMDALocalInfo *, void *, Mat, Mat, void *))FormJacobianLocal, &user));
214:   }

216:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217:      Customize nonlinear solver; set runtime options
218:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
219:   PetscCall(SNESSetFromOptions(snes));
220:   PetscCall(SNESSetNGS(snes, NonlinearGS, &user));
221:   PetscCall(SNESGetLineSearch(snes, &linesearch));
222:   /* Set up the precheck context if requested */
223:   if (use_precheck == 1) { /* Use the precheck routines in this file */
224:     PetscCall(PreCheckCreate(PETSC_COMM_WORLD, &precheck));
225:     PetscCall(PreCheckSetFromOptions(precheck));
226:     PetscCall(SNESLineSearchSetPreCheck(linesearch, PreCheckFunction, precheck));
227:   } else if (use_precheck == 2) { /* Use the version provided by the library */
228:     PetscCall(SNESLineSearchSetPreCheck(linesearch, SNESLineSearchPreCheckPicard, &precheck_angle));
229:   }

231:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
232:      Evaluate initial guess
233:      Note: The user should initialize the vector, x, with the initial guess
234:      for the nonlinear solver prior to calling SNESSolve().  In particular,
235:      to employ an initial guess of zero, the user should explicitly set
236:      this vector to zero by calling VecSet().
237:   */

239:   PetscCall(FormInitialGuess(&user, da, x));
240:   PetscCall(FormRHS(&user, da, b));

242:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
243:      Solve nonlinear system
244:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
245:   PetscCall(SNESSolve(snes, b, x));
246:   PetscCall(SNESGetIterationNumber(snes, &its));
247:   PetscCall(SNESGetConvergedReason(snes, &reason));

249:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s Number of nonlinear iterations = %" PetscInt_FMT "\n", SNESConvergedReasons[reason], its));

251:   if (write_output) {
252:     PetscViewer viewer;
253:     PetscCall(PetscViewerVTKOpen(PETSC_COMM_WORLD, filename, FILE_MODE_WRITE, &viewer));
254:     PetscCall(VecView(x, viewer));
255:     PetscCall(PetscViewerDestroy(&viewer));
256:   }

258:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
259:      Free work space.  All PETSc objects should be destroyed when they
260:      are no longer needed.
261:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

263:   PetscCall(VecDestroy(&x));
264:   PetscCall(VecDestroy(&r));
265:   PetscCall(VecDestroy(&b));
266:   PetscCall(SNESDestroy(&snes));
267:   PetscCall(DMDestroy(&da));
268:   PetscCall(PreCheckDestroy(&precheck));
269:   PetscCall(PetscFinalize());
270:   return 0;
271: }

273: /* ------------------------------------------------------------------- */
274: /*
275:    FormInitialGuess - Forms initial approximation.

277:    Input Parameters:
278:    user - user-defined application context
279:    X - vector

281:    Output Parameter:
282:    X - vector
283:  */
284: static PetscErrorCode FormInitialGuess(AppCtx *user, DM da, Vec X)
285: {
286:   PetscInt      i, j, Mx, My, xs, ys, xm, ym;
287:   PetscReal     temp1, temp, hx, hy;
288:   PetscScalar **x;

290:   PetscFunctionBeginUser;
291:   PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));

293:   hx    = 1.0 / (PetscReal)(Mx - 1);
294:   hy    = 1.0 / (PetscReal)(My - 1);
295:   temp1 = user->lambda / (user->lambda + 1.);

297:   /*
298:      Get a pointer to vector data.
299:        - For default PETSc vectors, VecGetArray() returns a pointer to
300:          the data array.  Otherwise, the routine is implementation dependent.
301:        - You MUST call VecRestoreArray() when you no longer need access to
302:          the array.
303:   */
304:   PetscCall(DMDAVecGetArray(da, X, &x));

306:   /*
307:      Get local grid boundaries (for 2-dimensional DA):
308:        xs, ys   - starting grid indices (no ghost points)
309:        xm, ym   - widths of local grid (no ghost points)

311:   */
312:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));

314:   /*
315:      Compute initial guess over the locally owned part of the grid
316:   */
317:   for (j = ys; j < ys + ym; j++) {
318:     temp = (PetscReal)(PetscMin(j, My - j - 1)) * hy;
319:     for (i = xs; i < xs + xm; i++) {
320:       if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) {
321:         /* boundary conditions are all zero Dirichlet */
322:         x[j][i] = 0.0;
323:       } else {
324:         if (user->initial == -1) {
325:           if (user->lambda != 0) {
326:             x[j][i] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, Mx - i - 1)) * hx, temp));
327:           } else {
328:             /* The solution above is an exact solution for lambda=0, this avoids "accidentally" starting
329:              * with an exact solution. */
330:             const PetscReal xx = 2 * (PetscReal)i / (Mx - 1) - 1, yy = 2 * (PetscReal)j / (My - 1) - 1;
331:             x[j][i] = (1 - xx * xx) * (1 - yy * yy) * xx * yy;
332:           }
333:         } else if (user->initial == 0) {
334:           x[j][i] = 0.;
335:         } else if (user->initial == 1) {
336:           const PetscReal xx = 2 * (PetscReal)i / (Mx - 1) - 1, yy = 2 * (PetscReal)j / (My - 1) - 1;
337:           x[j][i] = (1 - xx * xx) * (1 - yy * yy) * xx * yy;
338:         } else {
339:           if (user->lambda != 0) {
340:             x[j][i] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, Mx - i - 1)) * hx, temp));
341:           } else {
342:             x[j][i] = 0.5 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, Mx - i - 1)) * hx, temp));
343:           }
344:         }
345:       }
346:     }
347:   }
348:   /*
349:      Restore vector
350:   */
351:   PetscCall(DMDAVecRestoreArray(da, X, &x));
352:   PetscFunctionReturn(PETSC_SUCCESS);
353: }

355: /*
356:    FormRHS - Forms constant RHS for the problem.

358:    Input Parameters:
359:    user - user-defined application context
360:    B - RHS vector

362:    Output Parameter:
363:    B - vector
364:  */
365: static PetscErrorCode FormRHS(AppCtx *user, DM da, Vec B)
366: {
367:   PetscInt      i, j, Mx, My, xs, ys, xm, ym;
368:   PetscReal     hx, hy;
369:   PetscScalar **b;

371:   PetscFunctionBeginUser;
372:   PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));

374:   hx = 1.0 / (PetscReal)(Mx - 1);
375:   hy = 1.0 / (PetscReal)(My - 1);
376:   PetscCall(DMDAVecGetArray(da, B, &b));
377:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
378:   for (j = ys; j < ys + ym; j++) {
379:     for (i = xs; i < xs + xm; i++) {
380:       if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) {
381:         b[j][i] = 0.0;
382:       } else {
383:         b[j][i] = hx * hy * user->source;
384:       }
385:     }
386:   }
387:   PetscCall(DMDAVecRestoreArray(da, B, &b));
388:   PetscFunctionReturn(PETSC_SUCCESS);
389: }

391: static inline PetscReal kappa(const AppCtx *ctx, PetscReal x, PetscReal y)
392: {
393:   return (((PetscInt)(x * ctx->blocks[0])) + ((PetscInt)(y * ctx->blocks[1]))) % 2 ? ctx->kappa : 1.0;
394: }
395: /* p-Laplacian diffusivity */
396: static inline PetscScalar eta(const AppCtx *ctx, PetscReal x, PetscReal y, PetscScalar ux, PetscScalar uy)
397: {
398:   return kappa(ctx, x, y) * PetscPowScalar(PetscSqr(ctx->epsilon) + 0.5 * (ux * ux + uy * uy), 0.5 * (ctx->p - 2.));
399: }
400: static inline PetscScalar deta(const AppCtx *ctx, PetscReal x, PetscReal y, PetscScalar ux, PetscScalar uy)
401: {
402:   return (ctx->p == 2) ? 0 : kappa(ctx, x, y) * PetscPowScalar(PetscSqr(ctx->epsilon) + 0.5 * (ux * ux + uy * uy), 0.5 * (ctx->p - 4)) * 0.5 * (ctx->p - 2.);
403: }

405: /* ------------------------------------------------------------------- */
406: /*
407:    FormFunctionLocal - Evaluates nonlinear function, F(x).
408:  */
409: static PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, PetscScalar **x, PetscScalar **f, AppCtx *user)
410: {
411:   PetscReal   hx, hy, dhx, dhy, sc;
412:   PetscInt    i, j;
413:   PetscScalar eu;

415:   PetscFunctionBeginUser;
416:   hx  = 1.0 / (PetscReal)(info->mx - 1);
417:   hy  = 1.0 / (PetscReal)(info->my - 1);
418:   sc  = hx * hy * user->lambda;
419:   dhx = 1 / hx;
420:   dhy = 1 / hy;
421:   /*
422:      Compute function over the locally owned part of the grid
423:   */
424:   for (j = info->ys; j < info->ys + info->ym; j++) {
425:     for (i = info->xs; i < info->xs + info->xm; i++) {
426:       PetscReal xx = i * hx, yy = j * hy;
427:       if (i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1) {
428:         f[j][i] = x[j][i];
429:       } else {
430:         const PetscScalar u = x[j][i], ux_E = dhx * (x[j][i + 1] - x[j][i]), uy_E = 0.25 * dhy * (x[j + 1][i] + x[j + 1][i + 1] - x[j - 1][i] - x[j - 1][i + 1]), ux_W = dhx * (x[j][i] - x[j][i - 1]), uy_W = 0.25 * dhy * (x[j + 1][i - 1] + x[j + 1][i] - x[j - 1][i - 1] - x[j - 1][i]), ux_N = 0.25 * dhx * (x[j][i + 1] + x[j + 1][i + 1] - x[j][i - 1] - x[j + 1][i - 1]), uy_N = dhy * (x[j + 1][i] - x[j][i]), ux_S = 0.25 * dhx * (x[j - 1][i + 1] + x[j][i + 1] - x[j - 1][i - 1] - x[j][i - 1]), uy_S = dhy * (x[j][i] - x[j - 1][i]), e_E = eta(user, xx, yy, ux_E, uy_E), e_W = eta(user, xx, yy, ux_W, uy_W), e_N = eta(user, xx, yy, ux_N, uy_N), e_S = eta(user, xx, yy, ux_S, uy_S), uxx = -hy * (e_E * ux_E - e_W * ux_W), uyy = -hx * (e_N * uy_N - e_S * uy_S);
431:         if (sc) eu = PetscExpScalar(u);
432:         else eu = 0.;
433:         /* For p=2, these terms decay to:
434:          uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx
435:          uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy
436:         */
437:         f[j][i] = uxx + uyy - sc * eu;
438:       }
439:     }
440:   }
441:   PetscCall(PetscLogFlops(info->xm * info->ym * (72.0)));
442:   PetscFunctionReturn(PETSC_SUCCESS);
443: }

445: /*
446:     This is the opposite sign of the part of FormFunctionLocal that excludes the A(x) x part of the operation,
447:     that is FormFunction applies A(x) x - b(x) while this applies b(x) because for Picard we think of it as solving A(x) x = b(x)

449: */
450: static PetscErrorCode FormFunctionPicardLocal(DMDALocalInfo *info, PetscScalar **x, PetscScalar **f, AppCtx *user)
451: {
452:   PetscReal hx, hy, sc;
453:   PetscInt  i, j;

455:   PetscFunctionBeginUser;
456:   hx = 1.0 / (PetscReal)(info->mx - 1);
457:   hy = 1.0 / (PetscReal)(info->my - 1);
458:   sc = hx * hy * user->lambda;
459:   /*
460:      Compute function over the locally owned part of the grid
461:   */
462:   for (j = info->ys; j < info->ys + info->ym; j++) {
463:     for (i = info->xs; i < info->xs + info->xm; i++) {
464:       if (!(i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1)) {
465:         const PetscScalar u = x[j][i];
466:         f[j][i]             = sc * PetscExpScalar(u);
467:       } else {
468:         f[j][i] = 0.0; /* this is zero because the A(x) x term forces the x to be zero on the boundary */
469:       }
470:     }
471:   }
472:   PetscCall(PetscLogFlops(info->xm * info->ym * 2.0));
473:   PetscFunctionReturn(PETSC_SUCCESS);
474: }

476: /*
477:    FormJacobianLocal - Evaluates Jacobian matrix.
478: */
479: static PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, PetscScalar **x, Mat J, Mat B, AppCtx *user)
480: {
481:   PetscInt    i, j;
482:   MatStencil  col[9], row;
483:   PetscScalar v[9];
484:   PetscReal   hx, hy, hxdhy, hydhx, dhx, dhy, sc;

486:   PetscFunctionBeginUser;
487:   PetscCall(MatZeroEntries(B));
488:   hx    = 1.0 / (PetscReal)(info->mx - 1);
489:   hy    = 1.0 / (PetscReal)(info->my - 1);
490:   sc    = hx * hy * user->lambda;
491:   dhx   = 1 / hx;
492:   dhy   = 1 / hy;
493:   hxdhy = hx / hy;
494:   hydhx = hy / hx;

496:   /*
497:      Compute entries for the locally owned part of the Jacobian.
498:       - PETSc parallel matrix formats are partitioned by
499:         contiguous chunks of rows across the processors.
500:       - Each processor needs to insert only elements that it owns
501:         locally (but any non-local elements will be sent to the
502:         appropriate processor during matrix assembly).
503:       - Here, we set all entries for a particular row at once.
504:   */
505:   for (j = info->ys; j < info->ys + info->ym; j++) {
506:     for (i = info->xs; i < info->xs + info->xm; i++) {
507:       PetscReal xx = i * hx, yy = j * hy;
508:       row.j = j;
509:       row.i = i;
510:       /* boundary points */
511:       if (i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1) {
512:         v[0] = 1.0;
513:         PetscCall(MatSetValuesStencil(B, 1, &row, 1, &row, v, INSERT_VALUES));
514:       } else {
515:         /* interior grid points */
516:         const PetscScalar ux_E = dhx * (x[j][i + 1] - x[j][i]), uy_E = 0.25 * dhy * (x[j + 1][i] + x[j + 1][i + 1] - x[j - 1][i] - x[j - 1][i + 1]), ux_W = dhx * (x[j][i] - x[j][i - 1]), uy_W = 0.25 * dhy * (x[j + 1][i - 1] + x[j + 1][i] - x[j - 1][i - 1] - x[j - 1][i]), ux_N = 0.25 * dhx * (x[j][i + 1] + x[j + 1][i + 1] - x[j][i - 1] - x[j + 1][i - 1]), uy_N = dhy * (x[j + 1][i] - x[j][i]), ux_S = 0.25 * dhx * (x[j - 1][i + 1] + x[j][i + 1] - x[j - 1][i - 1] - x[j][i - 1]), uy_S = dhy * (x[j][i] - x[j - 1][i]), u = x[j][i], e_E = eta(user, xx, yy, ux_E, uy_E), e_W = eta(user, xx, yy, ux_W, uy_W), e_N = eta(user, xx, yy, ux_N, uy_N), e_S = eta(user, xx, yy, ux_S, uy_S), de_E = deta(user, xx, yy, ux_E, uy_E), de_W = deta(user, xx, yy, ux_W, uy_W), de_N = deta(user, xx, yy, ux_N, uy_N), de_S = deta(user, xx, yy, ux_S, uy_S), skew_E = de_E * ux_E * uy_E, skew_W = de_W * ux_W * uy_W, skew_N = de_N * ux_N * uy_N, skew_S = de_S * ux_S * uy_S, cross_EW = 0.25 * (skew_E - skew_W), cross_NS = 0.25 * (skew_N - skew_S), newt_E = e_E + de_E * PetscSqr(ux_E), newt_W = e_W + de_W * PetscSqr(ux_W), newt_N = e_N + de_N * PetscSqr(uy_N), newt_S = e_S + de_S * PetscSqr(uy_S);
517:         /* interior grid points */
518:         switch (user->jtype) {
519:         case JAC_BRATU:
520:           /* Jacobian from p=2 */
521:           v[0]     = -hxdhy;
522:           col[0].j = j - 1;
523:           col[0].i = i;
524:           v[1]     = -hydhx;
525:           col[1].j = j;
526:           col[1].i = i - 1;
527:           v[2]     = 2.0 * (hydhx + hxdhy) - sc * PetscExpScalar(u);
528:           col[2].j = row.j;
529:           col[2].i = row.i;
530:           v[3]     = -hydhx;
531:           col[3].j = j;
532:           col[3].i = i + 1;
533:           v[4]     = -hxdhy;
534:           col[4].j = j + 1;
535:           col[4].i = i;
536:           PetscCall(MatSetValuesStencil(B, 1, &row, 5, col, v, INSERT_VALUES));
537:           break;
538:         case JAC_PICARD:
539:           /* Jacobian arising from Picard linearization */
540:           v[0]     = -hxdhy * e_S;
541:           col[0].j = j - 1;
542:           col[0].i = i;
543:           v[1]     = -hydhx * e_W;
544:           col[1].j = j;
545:           col[1].i = i - 1;
546:           v[2]     = (e_W + e_E) * hydhx + (e_S + e_N) * hxdhy;
547:           col[2].j = row.j;
548:           col[2].i = row.i;
549:           v[3]     = -hydhx * e_E;
550:           col[3].j = j;
551:           col[3].i = i + 1;
552:           v[4]     = -hxdhy * e_N;
553:           col[4].j = j + 1;
554:           col[4].i = i;
555:           PetscCall(MatSetValuesStencil(B, 1, &row, 5, col, v, INSERT_VALUES));
556:           break;
557:         case JAC_STAR:
558:           /* Full Jacobian, but only a star stencil */
559:           col[0].j = j - 1;
560:           col[0].i = i;
561:           col[1].j = j;
562:           col[1].i = i - 1;
563:           col[2].j = j;
564:           col[2].i = i;
565:           col[3].j = j;
566:           col[3].i = i + 1;
567:           col[4].j = j + 1;
568:           col[4].i = i;
569:           v[0]     = -hxdhy * newt_S + cross_EW;
570:           v[1]     = -hydhx * newt_W + cross_NS;
571:           v[2]     = hxdhy * (newt_N + newt_S) + hydhx * (newt_E + newt_W) - sc * PetscExpScalar(u);
572:           v[3]     = -hydhx * newt_E - cross_NS;
573:           v[4]     = -hxdhy * newt_N - cross_EW;
574:           PetscCall(MatSetValuesStencil(B, 1, &row, 5, col, v, INSERT_VALUES));
575:           break;
576:         case JAC_NEWTON:
577:           /* The Jacobian is

581:           */
582:           col[0].j = j - 1;
583:           col[0].i = i - 1;
584:           col[1].j = j - 1;
585:           col[1].i = i;
586:           col[2].j = j - 1;
587:           col[2].i = i + 1;
588:           col[3].j = j;
589:           col[3].i = i - 1;
590:           col[4].j = j;
591:           col[4].i = i;
592:           col[5].j = j;
593:           col[5].i = i + 1;
594:           col[6].j = j + 1;
595:           col[6].i = i - 1;
596:           col[7].j = j + 1;
597:           col[7].i = i;
598:           col[8].j = j + 1;
599:           col[8].i = i + 1;
600:           v[0]     = -0.25 * (skew_S + skew_W);
601:           v[1]     = -hxdhy * newt_S + cross_EW;
602:           v[2]     = 0.25 * (skew_S + skew_E);
603:           v[3]     = -hydhx * newt_W + cross_NS;
604:           v[4]     = hxdhy * (newt_N + newt_S) + hydhx * (newt_E + newt_W) - sc * PetscExpScalar(u);
605:           v[5]     = -hydhx * newt_E - cross_NS;
606:           v[6]     = 0.25 * (skew_N + skew_W);
607:           v[7]     = -hxdhy * newt_N - cross_EW;
608:           v[8]     = -0.25 * (skew_N + skew_E);
609:           PetscCall(MatSetValuesStencil(B, 1, &row, 9, col, v, INSERT_VALUES));
610:           break;
611:         default:
612:           SETERRQ(PetscObjectComm((PetscObject)info->da), PETSC_ERR_SUP, "Jacobian type %d not implemented", user->jtype);
613:         }
614:       }
615:     }
616:   }

618:   /*
619:      Assemble matrix, using the 2-step process:
620:        MatAssemblyBegin(), MatAssemblyEnd().
621:   */
622:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
623:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

625:   if (J != B) {
626:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
627:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
628:   }

630:   /*
631:      Tell the matrix we will never add a new nonzero location to the
632:      matrix. If we do, it will generate an error.
633:   */
634:   if (user->jtype == JAC_NEWTON) PetscCall(PetscLogFlops(info->xm * info->ym * (131.0)));
635:   PetscFunctionReturn(PETSC_SUCCESS);
636: }

638: /***********************************************************
639:  * PreCheck implementation
640:  ***********************************************************/
641: PetscErrorCode PreCheckSetFromOptions(PreCheck precheck)
642: {
643:   PetscBool flg;

645:   PetscFunctionBeginUser;
646:   PetscOptionsBegin(precheck->comm, NULL, "PreCheck Options", "none");
647:   PetscCall(PetscOptionsReal("-precheck_angle", "Angle in degrees between successive search directions necessary to activate step correction", "", precheck->angle, &precheck->angle, NULL));
648:   flg = PETSC_FALSE;
649:   PetscCall(PetscOptionsBool("-precheck_monitor", "Monitor choices made by precheck routine", "", flg, &flg, NULL));
650:   if (flg) PetscCall(PetscViewerASCIIOpen(precheck->comm, "stdout", &precheck->monitor));
651:   PetscOptionsEnd();
652:   PetscFunctionReturn(PETSC_SUCCESS);
653: }

655: /*
656:   Compare the direction of the current and previous step, modify the current step accordingly
657: */
658: PetscErrorCode PreCheckFunction(SNESLineSearch linesearch, Vec X, Vec Y, PetscBool *changed, void *ctx)
659: {
660:   PreCheck    precheck;
661:   Vec         Ylast;
662:   PetscScalar dot;
663:   PetscInt    iter;
664:   PetscReal   ynorm, ylastnorm, theta, angle_radians;
665:   SNES        snes;

667:   PetscFunctionBeginUser;
668:   PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
669:   precheck = (PreCheck)ctx;
670:   if (!precheck->Ylast) PetscCall(VecDuplicate(Y, &precheck->Ylast));
671:   Ylast = precheck->Ylast;
672:   PetscCall(SNESGetIterationNumber(snes, &iter));
673:   if (iter < 1) {
674:     PetscCall(VecCopy(Y, Ylast));
675:     *changed = PETSC_FALSE;
676:     PetscFunctionReturn(PETSC_SUCCESS);
677:   }

679:   PetscCall(VecDot(Y, Ylast, &dot));
680:   PetscCall(VecNorm(Y, NORM_2, &ynorm));
681:   PetscCall(VecNorm(Ylast, NORM_2, &ylastnorm));
682:   /* Compute the angle between the vectors Y and Ylast, clip to keep inside the domain of acos() */
683:   theta         = PetscAcosReal((PetscReal)PetscClipInterval(PetscAbsScalar(dot) / (ynorm * ylastnorm), -1.0, 1.0));
684:   angle_radians = precheck->angle * PETSC_PI / 180.;
685:   if (PetscAbsReal(theta) < angle_radians || PetscAbsReal(theta - PETSC_PI) < angle_radians) {
686:     /* Modify the step Y */
687:     PetscReal alpha, ydiffnorm;
688:     PetscCall(VecAXPY(Ylast, -1.0, Y));
689:     PetscCall(VecNorm(Ylast, NORM_2, &ydiffnorm));
690:     alpha = ylastnorm / ydiffnorm;
691:     PetscCall(VecCopy(Y, Ylast));
692:     PetscCall(VecScale(Y, alpha));
693:     if (precheck->monitor) PetscCall(PetscViewerASCIIPrintf(precheck->monitor, "Angle %g degrees less than threshold %g, corrected step by alpha=%g\n", (double)(theta * 180. / PETSC_PI), (double)precheck->angle, (double)alpha));
694:   } else {
695:     PetscCall(VecCopy(Y, Ylast));
696:     *changed = PETSC_FALSE;
697:     if (precheck->monitor) PetscCall(PetscViewerASCIIPrintf(precheck->monitor, "Angle %g degrees exceeds threshold %g, no correction applied\n", (double)(theta * 180. / PETSC_PI), (double)precheck->angle));
698:   }
699:   PetscFunctionReturn(PETSC_SUCCESS);
700: }

702: PetscErrorCode PreCheckDestroy(PreCheck *precheck)
703: {
704:   PetscFunctionBeginUser;
705:   if (!*precheck) PetscFunctionReturn(PETSC_SUCCESS);
706:   PetscCall(VecDestroy(&(*precheck)->Ylast));
707:   PetscCall(PetscViewerDestroy(&(*precheck)->monitor));
708:   PetscCall(PetscFree(*precheck));
709:   PetscFunctionReturn(PETSC_SUCCESS);
710: }

712: PetscErrorCode PreCheckCreate(MPI_Comm comm, PreCheck *precheck)
713: {
714:   PetscFunctionBeginUser;
715:   PetscCall(PetscNew(precheck));

717:   (*precheck)->comm  = comm;
718:   (*precheck)->angle = 10.; /* only active if angle is less than 10 degrees */
719:   PetscFunctionReturn(PETSC_SUCCESS);
720: }

722: /*
723:       Applies some sweeps on nonlinear Gauss-Seidel on each process

725:  */
726: PetscErrorCode NonlinearGS(SNES snes, Vec X, Vec B, void *ctx)
727: {
728:   PetscInt      i, j, k, xs, ys, xm, ym, its, tot_its, sweeps, l, m;
729:   PetscReal     hx, hy, hxdhy, hydhx, dhx, dhy, sc;
730:   PetscScalar **x, **b, bij, F, F0 = 0, J, y, u, eu;
731:   PetscReal     atol, rtol, stol;
732:   DM            da;
733:   AppCtx       *user = (AppCtx *)ctx;
734:   Vec           localX, localB;
735:   DMDALocalInfo info;

737:   PetscFunctionBeginUser;
738:   PetscCall(SNESGetDM(snes, &da));
739:   PetscCall(DMDAGetLocalInfo(da, &info));

741:   hx    = 1.0 / (PetscReal)(info.mx - 1);
742:   hy    = 1.0 / (PetscReal)(info.my - 1);
743:   sc    = hx * hy * user->lambda;
744:   dhx   = 1 / hx;
745:   dhy   = 1 / hy;
746:   hxdhy = hx / hy;
747:   hydhx = hy / hx;

749:   tot_its = 0;
750:   PetscCall(SNESNGSGetSweeps(snes, &sweeps));
751:   PetscCall(SNESNGSGetTolerances(snes, &atol, &rtol, &stol, &its));
752:   PetscCall(DMGetLocalVector(da, &localX));
753:   if (B) PetscCall(DMGetLocalVector(da, &localB));
754:   if (B) {
755:     PetscCall(DMGlobalToLocalBegin(da, B, INSERT_VALUES, localB));
756:     PetscCall(DMGlobalToLocalEnd(da, B, INSERT_VALUES, localB));
757:   }
758:   if (B) PetscCall(DMDAVecGetArrayRead(da, localB, &b));
759:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
760:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
761:   PetscCall(DMDAVecGetArray(da, localX, &x));
762:   for (l = 0; l < sweeps; l++) {
763:     /*
764:      Get local grid boundaries (for 2-dimensional DMDA):
765:      xs, ys   - starting grid indices (no ghost points)
766:      xm, ym   - widths of local grid (no ghost points)
767:      */
768:     PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
769:     for (m = 0; m < 2; m++) {
770:       for (j = ys; j < ys + ym; j++) {
771:         for (i = xs + (m + j) % 2; i < xs + xm; i += 2) {
772:           PetscReal xx = i * hx, yy = j * hy;
773:           if (B) bij = b[j][i];
774:           else bij = 0.;

776:           if (i == 0 || j == 0 || i == info.mx - 1 || j == info.my - 1) {
777:             /* boundary conditions are all zero Dirichlet */
778:             x[j][i] = 0.0 + bij;
779:           } else {
780:             const PetscScalar u_E = x[j][i + 1], u_W = x[j][i - 1], u_N = x[j + 1][i], u_S = x[j - 1][i];
781:             const PetscScalar uy_E = 0.25 * dhy * (x[j + 1][i] + x[j + 1][i + 1] - x[j - 1][i] - x[j - 1][i + 1]), uy_W = 0.25 * dhy * (x[j + 1][i - 1] + x[j + 1][i] - x[j - 1][i - 1] - x[j - 1][i]), ux_N = 0.25 * dhx * (x[j][i + 1] + x[j + 1][i + 1] - x[j][i - 1] - x[j + 1][i - 1]), ux_S = 0.25 * dhx * (x[j - 1][i + 1] + x[j][i + 1] - x[j - 1][i - 1] - x[j][i - 1]);
782:             u = x[j][i];
783:             for (k = 0; k < its; k++) {
784:               const PetscScalar ux_E = dhx * (u_E - u), ux_W = dhx * (u - u_W), uy_N = dhy * (u_N - u), uy_S = dhy * (u - u_S), e_E = eta(user, xx, yy, ux_E, uy_E), e_W = eta(user, xx, yy, ux_W, uy_W), e_N = eta(user, xx, yy, ux_N, uy_N), e_S = eta(user, xx, yy, ux_S, uy_S), de_E = deta(user, xx, yy, ux_E, uy_E), de_W = deta(user, xx, yy, ux_W, uy_W), de_N = deta(user, xx, yy, ux_N, uy_N), de_S = deta(user, xx, yy, ux_S, uy_S), newt_E = e_E + de_E * PetscSqr(ux_E), newt_W = e_W + de_W * PetscSqr(ux_W), newt_N = e_N + de_N * PetscSqr(uy_N), newt_S = e_S + de_S * PetscSqr(uy_S), uxx = -hy * (e_E * ux_E - e_W * ux_W), uyy = -hx * (e_N * uy_N - e_S * uy_S);

786:               if (sc) eu = PetscExpScalar(u);
787:               else eu = 0;

789:               F = uxx + uyy - sc * eu - bij;
790:               if (k == 0) F0 = F;
791:               J = hxdhy * (newt_N + newt_S) + hydhx * (newt_E + newt_W) - sc * eu;
792:               y = F / J;
793:               u -= y;
794:               tot_its++;
795:               if (atol > PetscAbsReal(PetscRealPart(F)) || rtol * PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) || stol * PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) break;
796:             }
797:             x[j][i] = u;
798:           }
799:         }
800:       }
801:     }
802:     /*
803: x     Restore vector
804:      */
805:   }
806:   PetscCall(DMDAVecRestoreArray(da, localX, &x));
807:   PetscCall(DMLocalToGlobalBegin(da, localX, INSERT_VALUES, X));
808:   PetscCall(DMLocalToGlobalEnd(da, localX, INSERT_VALUES, X));
809:   PetscCall(PetscLogFlops(tot_its * (118.0)));
810:   PetscCall(DMRestoreLocalVector(da, &localX));
811:   if (B) {
813:     PetscCall(DMRestoreLocalVector(da, &localB));
814:   }
815:   PetscFunctionReturn(PETSC_SUCCESS);
816: }

818: /*TEST

820:    test:
821:       nsize: 2
822:       args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype NEWTON
823:       requires: !single

825:    test:
826:       suffix: 2
827:       nsize: 2
828:       args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -precheck 1
829:       requires: !single

831:    test:
832:       suffix: 3
833:       nsize: 2
834:       args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -picard -precheck 1 -p 3
835:       requires: !single

837:    test:
838:       suffix: 3_star
839:       nsize: 2
840:       args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -picard -precheck 1 -p 3 -alloc_star
841:       output_file: output/ex15_3.out
842:       requires: !single

844:    test:
845:       suffix: 4
846:       args: -snes_monitor_short -snes_type newtonls -npc_snes_type ngs -snes_npc_side left -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -ksp_monitor_short -pc_type none
847:       requires: !single

849:    test:
850:       suffix: 5
851:       args: -snes_monitor_short -snes_type newtontr -npc_snes_type ngs -snes_npc_side right -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -ksp_converged_reason -pc_type ilu
852:       requires: !single

854:    test:
855:       suffix: lag_jac
856:       nsize: 4
857:       args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 6.0 -lambda 0 -jtype NEWTON -snes_type ngmres -npc_snes_type newtonls -npc_snes_lag_jacobian 5 -npc_pc_type asm -npc_ksp_converged_reason -npc_snes_lag_jacobian_persists
858:       requires: !single

860:    test:
861:       suffix: lag_pc
862:       nsize: 4
863:       args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 6.0 -lambda 0 -jtype NEWTON -snes_type ngmres -npc_snes_type newtonls -npc_snes_lag_preconditioner 5 -npc_pc_type asm -npc_ksp_converged_reason -npc_snes_lag_preconditioner_persists
864:       requires: !single

866:    test:
867:       suffix: nleqerr
868:       args: -snes_monitor_short -snes_type newtonls -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -snes_linesearch_monitor -pc_type lu -snes_linesearch_type nleqerr
869:       requires: !single

871:    test:
872:       suffix: mf
873:       args: -snes_monitor_short -pc_type lu -da_refine 4 -p 3 -ksp_rtol 1.e-12 -snes_mf_operator
874:       requires: !single

876:    test:
877:       suffix: mf_picard
878:       args: -snes_monitor_short -pc_type lu -da_refine 4 -p 3 -ksp_rtol 1.e-12 -snes_mf_operator -picard
879:       requires: !single
880:       output_file: output/ex15_mf.out

882:    test:
883:       suffix: fd_picard
884:       args: -snes_monitor_short -pc_type lu -da_refine 2 -p 3 -ksp_rtol 1.e-12 -snes_fd -picard
885:       requires: !single

887:    test:
888:       suffix: fd_color_picard
889:       args: -snes_monitor_short -pc_type lu -da_refine 4 -p 3 -ksp_rtol 1.e-12 -snes_fd_color -picard
890:       requires: !single
891:       output_file: output/ex15_mf.out

893: TEST*/