Actual source code: ex14.c

  1: static char help[] = "Bratu nonlinear PDE in 3d.\n\
  2: We solve the  Bratu (SFI - solid fuel ignition) problem in a 3D rectangular\n\
  3: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\
  4: The command line options include:\n\
  5:   -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
  6:      problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)\n\n";

  8: /* ------------------------------------------------------------------------

 10:     Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 11:     the partial differential equation

 13:             -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,

 15:     with boundary conditions

 17:              u = 0  for  x = 0, x = 1, y = 0, y = 1, z = 0, z = 1

 19:     A finite difference approximation with the usual 7-point stencil
 20:     is used to discretize the boundary value problem to obtain a nonlinear
 21:     system of equations.

 23:   ------------------------------------------------------------------------- */

 25: /*
 26:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 27:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 28:    file automatically includes:
 29:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 30:      petscmat.h - matrices
 31:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 32:      petscviewer.h - viewers               petscpc.h  - preconditioners
 33:      petscksp.h   - linear solvers
 34: */
 35: #include <petscdm.h>
 36: #include <petscdmda.h>
 37: #include <petscsnes.h>

 39: /*
 40:    User-defined application context - contains data needed by the
 41:    application-provided call-back routines, FormJacobian() and
 42:    FormFunction().
 43: */
 44: typedef struct {
 45:   PetscReal param; /* test problem parameter */
 46:   DM        da;    /* distributed array data structure */
 47: } AppCtx;

 49: /*
 50:    User-defined routines
 51: */
 52: extern PetscErrorCode FormFunctionLocal(SNES, Vec, Vec, void *);
 53: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
 54: extern PetscErrorCode FormInitialGuess(AppCtx *, Vec);
 55: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);

 57: int main(int argc, char **argv)
 58: {
 59:   SNES          snes;     /* nonlinear solver */
 60:   Vec           x, r;     /* solution, residual vectors */
 61:   Mat           J = NULL; /* Jacobian matrix */
 62:   AppCtx        user;     /* user-defined work context */
 63:   PetscInt      its;      /* iterations for convergence */
 64:   MatFDColoring matfdcoloring = NULL;
 65:   PetscBool     matrix_free = PETSC_FALSE, coloring = PETSC_FALSE, coloring_ds = PETSC_FALSE, local_coloring = PETSC_FALSE;
 66:   PetscReal     bratu_lambda_max = 6.81, bratu_lambda_min = 0., fnorm;

 68:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 69:      Initialize program
 70:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 72:   PetscFunctionBeginUser;
 73:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

 75:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 76:      Initialize problem parameters
 77:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 78:   user.param = 6.0;
 79:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-par", &user.param, NULL));
 80:   PetscCheck(user.param < bratu_lambda_max && user.param > bratu_lambda_min, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Lambda is out of range");

 82:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 83:      Create nonlinear solver context
 84:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 85:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));

 87:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 88:      Create distributed array (DMDA) to manage parallel grid and vectors
 89:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 90:   PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, 4, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, NULL, &user.da));
 91:   PetscCall(DMSetFromOptions(user.da));
 92:   PetscCall(DMSetUp(user.da));
 93:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 94:      Extract global vectors from DMDA; then duplicate for remaining
 95:      vectors that are the same types
 96:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 97:   PetscCall(DMCreateGlobalVector(user.da, &x));
 98:   PetscCall(VecDuplicate(x, &r));

100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:      Set function evaluation routine and vector
102:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103:   PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)&user));

105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106:      Create matrix data structure; set Jacobian evaluation routine

108:      Set Jacobian matrix data structure and default Jacobian evaluation
109:      routine. User can override with:
110:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
111:                 (unless user explicitly sets preconditioner)
112:      -snes_mf_operator : form preconditioning matrix as set by the user,
113:                          but use matrix-free approx for Jacobian-vector
114:                          products within Newton-Krylov method
115:      -fdcoloring : using finite differences with coloring to compute the Jacobian

117:      Note one can use -matfd_coloring wp or ds the only reason for the -fdcoloring_ds option
118:      below is to test the call to MatFDColoringSetType().
119:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf", &matrix_free, NULL));
121:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-fdcoloring", &coloring, NULL));
122:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-fdcoloring_ds", &coloring_ds, NULL));
123:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-fdcoloring_local", &local_coloring, NULL));
124:   if (!matrix_free) {
125:     PetscCall(DMSetMatType(user.da, MATAIJ));
126:     PetscCall(DMCreateMatrix(user.da, &J));
127:     if (coloring) {
128:       ISColoring iscoloring;
129:       if (!local_coloring) {
130:         PetscCall(DMCreateColoring(user.da, IS_COLORING_GLOBAL, &iscoloring));
131:         PetscCall(MatFDColoringCreate(J, iscoloring, &matfdcoloring));
132:         PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))FormFunction, &user));
133:       } else {
134:         PetscCall(DMCreateColoring(user.da, IS_COLORING_LOCAL, &iscoloring));
135:         PetscCall(MatFDColoringCreate(J, iscoloring, &matfdcoloring));
136:         PetscCall(MatFDColoringUseDM(J, matfdcoloring));
137:         PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))FormFunctionLocal, &user));
138:       }
139:       if (coloring_ds) PetscCall(MatFDColoringSetType(matfdcoloring, MATMFFD_DS));
140:       PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
141:       PetscCall(MatFDColoringSetUp(J, iscoloring, matfdcoloring));
142:       PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, matfdcoloring));
143:       PetscCall(ISColoringDestroy(&iscoloring));
144:     } else {
145:       PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, &user));
146:     }
147:   }

149:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150:      Customize nonlinear solver; set runtime options
151:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152:   PetscCall(SNESSetDM(snes, user.da));
153:   PetscCall(SNESSetFromOptions(snes));

155:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156:      Evaluate initial guess
157:      Note: The user should initialize the vector, x, with the initial guess
158:      for the nonlinear solver prior to calling SNESSolve().  In particular,
159:      to employ an initial guess of zero, the user should explicitly set
160:      this vector to zero by calling VecSet().
161:   */
162:   PetscCall(FormInitialGuess(&user, x));

164:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165:      Solve nonlinear system
166:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167:   PetscCall(SNESSolve(snes, NULL, x));
168:   PetscCall(SNESGetIterationNumber(snes, &its));

170:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171:      Explicitly check norm of the residual of the solution
172:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173:   PetscCall(FormFunction(snes, x, r, (void *)&user));
174:   PetscCall(VecNorm(r, NORM_2, &fnorm));
175:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT " fnorm %g\n", its, (double)fnorm));

177:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178:      Free work space.  All PETSc objects should be destroyed when they
179:      are no longer needed.
180:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

182:   PetscCall(MatDestroy(&J));
183:   PetscCall(VecDestroy(&x));
184:   PetscCall(VecDestroy(&r));
185:   PetscCall(SNESDestroy(&snes));
186:   PetscCall(DMDestroy(&user.da));
187:   PetscCall(MatFDColoringDestroy(&matfdcoloring));
188:   PetscCall(PetscFinalize());
189:   return 0;
190: }
191: /* ------------------------------------------------------------------- */
192: /*
193:    FormInitialGuess - Forms initial approximation.

195:    Input Parameters:
196:    user - user-defined application context
197:    X - vector

199:    Output Parameter:
200:    X - vector
201:  */
202: PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
203: {
204:   PetscInt       i, j, k, Mx, My, Mz, xs, ys, zs, xm, ym, zm;
205:   PetscReal      lambda, temp1, hx, hy, hz, tempk, tempj;
206:   PetscScalar ***x;

208:   PetscFunctionBeginUser;
209:   PetscCall(DMDAGetInfo(user->da, PETSC_IGNORE, &Mx, &My, &Mz, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));

211:   lambda = user->param;
212:   hx     = 1.0 / (PetscReal)(Mx - 1);
213:   hy     = 1.0 / (PetscReal)(My - 1);
214:   hz     = 1.0 / (PetscReal)(Mz - 1);
215:   temp1  = lambda / (lambda + 1.0);

217:   /*
218:      Get a pointer to vector data.
219:        - For default PETSc vectors, VecGetArray() returns a pointer to
220:          the data array.  Otherwise, the routine is implementation dependent.
221:        - You MUST call VecRestoreArray() when you no longer need access to
222:          the array.
223:   */
224:   PetscCall(DMDAVecGetArray(user->da, X, &x));

226:   /*
227:      Get local grid boundaries (for 3-dimensional DMDA):
228:        xs, ys, zs   - starting grid indices (no ghost points)
229:        xm, ym, zm   - widths of local grid (no ghost points)

231:   */
232:   PetscCall(DMDAGetCorners(user->da, &xs, &ys, &zs, &xm, &ym, &zm));

234:   /*
235:      Compute initial guess over the locally owned part of the grid
236:   */
237:   for (k = zs; k < zs + zm; k++) {
238:     tempk = (PetscReal)(PetscMin(k, Mz - k - 1)) * hz;
239:     for (j = ys; j < ys + ym; j++) {
240:       tempj = PetscMin((PetscReal)(PetscMin(j, My - j - 1)) * hy, tempk);
241:       for (i = xs; i < xs + xm; i++) {
242:         if (i == 0 || j == 0 || k == 0 || i == Mx - 1 || j == My - 1 || k == Mz - 1) {
243:           /* boundary conditions are all zero Dirichlet */
244:           x[k][j][i] = 0.0;
245:         } else {
246:           x[k][j][i] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, Mx - i - 1)) * hx, tempj));
247:         }
248:       }
249:     }
250:   }

252:   /*
253:      Restore vector
254:   */
255:   PetscCall(DMDAVecRestoreArray(user->da, X, &x));
256:   PetscFunctionReturn(PETSC_SUCCESS);
257: }
258: /* ------------------------------------------------------------------- */
259: /*
260:    FormFunctionLocal - Evaluates nonlinear function, F(x) on a ghosted patch

262:    Input Parameters:
263: .  snes - the SNES context
264: .  localX - input vector, this contains the ghosted region needed
265: .  ptr - optional user-defined context, as set by SNESSetFunction()

267:    Output Parameter:
268: .  F - function vector, this does not contain a ghosted region
269:  */
270: PetscErrorCode FormFunctionLocal(SNES snes, Vec localX, Vec F, void *ptr)
271: {
272:   AppCtx     *user = (AppCtx *)ptr;
273:   PetscInt    i, j, k, Mx, My, Mz, xs, ys, zs, xm, ym, zm;
274:   PetscReal   two = 2.0, lambda, hx, hy, hz, hxhzdhy, hyhzdhx, hxhydhz, sc;
275:   PetscScalar u_north, u_south, u_east, u_west, u_up, u_down, u, u_xx, u_yy, u_zz, ***x, ***f;
276:   DM          da;

278:   PetscFunctionBeginUser;
279:   PetscCall(SNESGetDM(snes, &da));
280:   PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, &Mz, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));

282:   lambda  = user->param;
283:   hx      = 1.0 / (PetscReal)(Mx - 1);
284:   hy      = 1.0 / (PetscReal)(My - 1);
285:   hz      = 1.0 / (PetscReal)(Mz - 1);
286:   sc      = hx * hy * hz * lambda;
287:   hxhzdhy = hx * hz / hy;
288:   hyhzdhx = hy * hz / hx;
289:   hxhydhz = hx * hy / hz;

291:   /*
292:      Get pointers to vector data
293:   */
294:   PetscCall(DMDAVecGetArrayRead(da, localX, &x));
295:   PetscCall(DMDAVecGetArray(da, F, &f));

297:   /*
298:      Get local grid boundaries
299:   */
300:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));

302:   /*
303:      Compute function over the locally owned part of the grid
304:   */
305:   for (k = zs; k < zs + zm; k++) {
306:     for (j = ys; j < ys + ym; j++) {
307:       for (i = xs; i < xs + xm; i++) {
308:         if (i == 0 || j == 0 || k == 0 || i == Mx - 1 || j == My - 1 || k == Mz - 1) {
309:           f[k][j][i] = x[k][j][i];
310:         } else {
311:           u          = x[k][j][i];
312:           u_east     = x[k][j][i + 1];
313:           u_west     = x[k][j][i - 1];
314:           u_north    = x[k][j + 1][i];
315:           u_south    = x[k][j - 1][i];
316:           u_up       = x[k + 1][j][i];
317:           u_down     = x[k - 1][j][i];
318:           u_xx       = (-u_east + two * u - u_west) * hyhzdhx;
319:           u_yy       = (-u_north + two * u - u_south) * hxhzdhy;
320:           u_zz       = (-u_up + two * u - u_down) * hxhydhz;
321:           f[k][j][i] = u_xx + u_yy + u_zz - sc * PetscExpScalar(u);
322:         }
323:       }
324:     }
325:   }

327:   /*
328:      Restore vectors
329:   */
330:   PetscCall(DMDAVecRestoreArrayRead(da, localX, &x));
331:   PetscCall(DMDAVecRestoreArray(da, F, &f));
332:   PetscCall(PetscLogFlops(11.0 * ym * xm));
333:   PetscFunctionReturn(PETSC_SUCCESS);
334: }
335: /* ------------------------------------------------------------------- */
336: /*
337:    FormFunction - Evaluates nonlinear function, F(x) on the entire domain

339:    Input Parameters:
340: .  snes - the SNES context
341: .  X - input vector
342: .  ptr - optional user-defined context, as set by SNESSetFunction()

344:    Output Parameter:
345: .  F - function vector
346:  */
347: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr)
348: {
349:   Vec localX;
350:   DM  da;

352:   PetscFunctionBeginUser;
353:   PetscCall(SNESGetDM(snes, &da));
354:   PetscCall(DMGetLocalVector(da, &localX));

356:   /*
357:      Scatter ghost points to local vector,using the 2-step process
358:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
359:      By placing code between these two statements, computations can be
360:      done while messages are in transition.
361:   */
362:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
363:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));

365:   PetscCall(FormFunctionLocal(snes, localX, F, ptr));
366:   PetscCall(DMRestoreLocalVector(da, &localX));
367:   PetscFunctionReturn(PETSC_SUCCESS);
368: }
369: /* ------------------------------------------------------------------- */
370: /*
371:    FormJacobian - Evaluates Jacobian matrix.

373:    Input Parameters:
374: .  snes - the SNES context
375: .  x - input vector
376: .  ptr - optional user-defined context, as set by SNESSetJacobian()

378:    Output Parameters:
379: .  A - Jacobian matrix
380: .  B - optionally different preconditioning matrix

382: */
383: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr)
384: {
385:   AppCtx     *user = (AppCtx *)ptr; /* user-defined application context */
386:   Vec         localX;
387:   PetscInt    i, j, k, Mx, My, Mz;
388:   MatStencil  col[7], row;
389:   PetscInt    xs, ys, zs, xm, ym, zm;
390:   PetscScalar lambda, v[7], hx, hy, hz, hxhzdhy, hyhzdhx, hxhydhz, sc, ***x;
391:   DM          da;

393:   PetscFunctionBeginUser;
394:   PetscCall(SNESGetDM(snes, &da));
395:   PetscCall(DMGetLocalVector(da, &localX));
396:   PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, &Mz, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));

398:   lambda  = user->param;
399:   hx      = 1.0 / (PetscReal)(Mx - 1);
400:   hy      = 1.0 / (PetscReal)(My - 1);
401:   hz      = 1.0 / (PetscReal)(Mz - 1);
402:   sc      = hx * hy * hz * lambda;
403:   hxhzdhy = hx * hz / hy;
404:   hyhzdhx = hy * hz / hx;
405:   hxhydhz = hx * hy / hz;

407:   /*
408:      Scatter ghost points to local vector, using the 2-step process
409:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
410:      By placing code between these two statements, computations can be
411:      done while messages are in transition.
412:   */
413:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
414:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));

416:   /*
417:      Get pointer to vector data
418:   */
419:   PetscCall(DMDAVecGetArrayRead(da, localX, &x));

421:   /*
422:      Get local grid boundaries
423:   */
424:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));

426:   /*
427:      Compute entries for the locally owned part of the Jacobian.
428:       - Currently, all PETSc parallel matrix formats are partitioned by
429:         contiguous chunks of rows across the processors.
430:       - Each processor needs to insert only elements that it owns
431:         locally (but any non-local elements will be sent to the
432:         appropriate processor during matrix assembly).
433:       - Here, we set all entries for a particular row at once.
434:       - We can set matrix entries either using either
435:         MatSetValuesLocal() or MatSetValues(), as discussed above.
436:   */
437:   for (k = zs; k < zs + zm; k++) {
438:     for (j = ys; j < ys + ym; j++) {
439:       for (i = xs; i < xs + xm; i++) {
440:         row.k = k;
441:         row.j = j;
442:         row.i = i;
443:         /* boundary points */
444:         if (i == 0 || j == 0 || k == 0 || i == Mx - 1 || j == My - 1 || k == Mz - 1) {
445:           v[0] = 1.0;
446:           PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES));
447:         } else {
448:           /* interior grid points */
449:           v[0]     = -hxhydhz;
450:           col[0].k = k - 1;
451:           col[0].j = j;
452:           col[0].i = i;
453:           v[1]     = -hxhzdhy;
454:           col[1].k = k;
455:           col[1].j = j - 1;
456:           col[1].i = i;
457:           v[2]     = -hyhzdhx;
458:           col[2].k = k;
459:           col[2].j = j;
460:           col[2].i = i - 1;
461:           v[3]     = 2.0 * (hyhzdhx + hxhzdhy + hxhydhz) - sc * PetscExpScalar(x[k][j][i]);
462:           col[3].k = row.k;
463:           col[3].j = row.j;
464:           col[3].i = row.i;
465:           v[4]     = -hyhzdhx;
466:           col[4].k = k;
467:           col[4].j = j;
468:           col[4].i = i + 1;
469:           v[5]     = -hxhzdhy;
470:           col[5].k = k;
471:           col[5].j = j + 1;
472:           col[5].i = i;
473:           v[6]     = -hxhydhz;
474:           col[6].k = k + 1;
475:           col[6].j = j;
476:           col[6].i = i;
477:           PetscCall(MatSetValuesStencil(jac, 1, &row, 7, col, v, INSERT_VALUES));
478:         }
479:       }
480:     }
481:   }
482:   PetscCall(DMDAVecRestoreArrayRead(da, localX, &x));
483:   PetscCall(DMRestoreLocalVector(da, &localX));

485:   /*
486:      Assemble matrix, using the 2-step process:
487:        MatAssemblyBegin(), MatAssemblyEnd().
488:   */
489:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
490:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));

492:   /*
493:      Normally since the matrix has already been assembled above; this
494:      would do nothing. But in the matrix-free mode -snes_mf_operator
495:      this tells the "matrix-free" matrix that a new linear system solve
496:      is about to be done.
497:   */

499:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
500:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));

502:   /*
503:      Tell the matrix we will never add a new nonzero location to the
504:      matrix. If we do, it will generate an error.
505:   */
506:   PetscCall(MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
507:   PetscFunctionReturn(PETSC_SUCCESS);
508: }

510: /*TEST

512:    test:
513:       nsize: 4
514:       args: -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

516:    test:
517:       suffix: 2
518:       nsize: 4
519:       args: -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

521:    test:
522:       suffix: 3
523:       nsize: 4
524:       args: -fdcoloring -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

526:    test:
527:       suffix: 3_ds
528:       nsize: 4
529:       args: -fdcoloring -fdcoloring_ds -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

531:    test:
532:       suffix: 4
533:       nsize: 4
534:       args: -fdcoloring_local -fdcoloring -ksp_monitor_short -da_refine 1
535:       requires: !single

537:    test:
538:       suffix: 5
539:       nsize: 4
540:       args: -fdcoloring_local -fdcoloring -ksp_monitor_short -da_refine 1 -snes_type newtontrdc
541:       requires: !single

543:    test:
544:       suffix: 6
545:       nsize: 4
546:       args: -fdcoloring_local -fdcoloring -da_refine 1 -snes_type newtontr -snes_tr_fallback_type dogleg
547:       requires: !single

549: TEST*/