Actual source code: ex14.c


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

  9: /* ------------------------------------------------------------------------

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

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

 16:     with boundary conditions

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

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

 24:   ------------------------------------------------------------------------- */

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

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

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

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

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

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

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

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

 88:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 89:      Create distributed array (DMDA) to manage parallel grid and vectors
 90:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 91:   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));
 92:   PetscCall(DMSetFromOptions(user.da));
 93:   PetscCall(DMSetUp(user.da));
 94:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 95:      Extract global vectors from DMDA; then duplicate for remaining
 96:      vectors that are the same types
 97:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 98:   PetscCall(DMCreateGlobalVector(user.da, &x));
 99:   PetscCall(VecDuplicate(x, &r));

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

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

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

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

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

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

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

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

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

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

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

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

209:   PetscFunctionBeginUser;
210:   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));

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

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

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

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

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

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

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

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

279:   PetscFunctionBeginUser;
280:   PetscCall(SNESGetDM(snes, &da));
281:   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));

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

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

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

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

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

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

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

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

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

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

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

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

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

394:   PetscFunctionBeginUser;
395:   PetscCall(SNESGetDM(snes, &da));
396:   PetscCall(DMGetLocalVector(da, &localX));
397:   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));

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

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

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

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

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

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

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

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

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

511: /*TEST

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

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

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

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

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

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

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

550: TEST*/