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, NULL, 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

    579:            -div [ eta (grad u) + deta (grad u0 . grad u) grad u0 ] - (eE u0) u

    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) {
    812:     PetscCall(DMDAVecRestoreArrayRead(da, localB, &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*/