Actual source code: ex20.c

  1: static char help[] = "Nonlinear Radiative Transport PDE with multigrid in 3d.\n\
  2: Uses 3-dimensional distributed arrays.\n\
  3: A 3-dim simplified Radiative Transport test problem is used, with analytic Jacobian. \n\
  4: \n\
  5:   Solves the linear systems via multilevel methods \n\
  6: \n\
  7: The command line\n\
  8: options are:\n\
  9:   -tleft <tl>, where <tl> indicates the left Diriclet BC \n\
 10:   -tright <tr>, where <tr> indicates the right Diriclet BC \n\
 11:   -beta <beta>, where <beta> indicates the exponent in T \n\n";

 13: /*

 15:     This example models the partial differential equation

 17:          - Div(alpha* T^beta (GRAD T)) = 0.

 19:     where beta = 2.5 and alpha = 1.0

 21:     BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = dT/dn_up = dT/dn_down = 0.

 23:     in the unit square, which is uniformly discretized in each of x and
 24:     y in this simple encoding.  The degrees of freedom are cell centered.

 26:     A finite volume approximation with the usual 7-point stencil
 27:     is used to discretize the boundary value problem to obtain a
 28:     nonlinear system of equations.

 30:     This code was contributed by Nickolas Jovanovic based on ex18.c

 32: */

 34: #include <petscsnes.h>
 35: #include <petscdm.h>
 36: #include <petscdmda.h>

 38: /* User-defined application context */

 40: typedef struct {
 41:   PetscReal tleft, tright;   /* Dirichlet boundary conditions */
 42:   PetscReal beta, bm1, coef; /* nonlinear diffusivity parameterizations */
 43:   PetscBool converged;       /* For user convergence test, whether to mark as converged */
 44: } AppCtx;

 46: #define POWFLOP 5 /* assume a pow() takes five flops */

 48: extern PetscErrorCode FormInitialGuess(SNES, Vec, void *);
 49: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
 50: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
 51: extern PetscErrorCode TestConvergence(SNES, PetscInt, PetscReal, PetscReal, PetscReal, SNESConvergedReason *, void *);

 53: int main(int argc, char **argv)
 54: {
 55:   SNES      snes;
 56:   AppCtx    user;
 57:   PetscInt  its, lits;
 58:   PetscReal litspit = 0; /* avoid uninitialized warning */
 59:   DM        da;
 60:   PetscBool use_convergence_test = PETSC_FALSE;

 62:   PetscFunctionBeginUser;
 63:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 64:   /* set problem parameters */
 65:   user.tleft  = 1.0;
 66:   user.tright = 0.1;
 67:   user.beta   = 2.5;
 68:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-tleft", &user.tleft, NULL));
 69:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-tright", &user.tright, NULL));
 70:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-beta", &user.beta, NULL));
 71:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_convergence_test", &use_convergence_test, NULL));
 72:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-mark_converged", &user.converged, NULL));
 73:   user.bm1  = user.beta - 1.0;
 74:   user.coef = user.beta / 2.0;

 76:   /*
 77:       Set the DMDA (grid structure) for the grids.
 78:   */
 79:   PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 5, 5, 5, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, 0, &da));
 80:   PetscCall(DMSetFromOptions(da));
 81:   PetscCall(DMSetUp(da));
 82:   PetscCall(DMSetApplicationContext(da, &user));

 84:   /*
 85:      Create the nonlinear solver
 86:   */
 87:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
 88:   PetscCall(SNESSetDM(snes, da));
 89:   PetscCall(SNESSetFunction(snes, NULL, FormFunction, &user));
 90:   PetscCall(SNESSetJacobian(snes, NULL, NULL, FormJacobian, &user));
 91:   if (use_convergence_test) PetscCall(SNESSetConvergenceTest(snes, TestConvergence, &user, NULL));
 92:   PetscCall(SNESSetFromOptions(snes));
 93:   PetscCall(SNESSetComputeInitialGuess(snes, FormInitialGuess, NULL));

 95:   PetscCall(SNESSolve(snes, NULL, NULL));
 96:   PetscCall(SNESGetIterationNumber(snes, &its));
 97:   PetscCall(SNESGetLinearSolveIterations(snes, &lits));
 98:   if (its) litspit = ((PetscReal)lits) / ((PetscReal)its);
 99:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
100:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of Linear iterations = %" PetscInt_FMT "\n", lits));
101:   if (its) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Average Linear its / SNES = %e\n", (double)litspit));

103:   PetscCall(SNESDestroy(&snes));
104:   PetscCall(DMDestroy(&da));
105:   PetscCall(PetscFinalize());
106:   return 0;
107: }
108: /* --------------------  Form initial approximation ----------------- */
109: PetscErrorCode FormInitialGuess(SNES snes, Vec X, void *ctx)
110: {
111:   AppCtx        *user;
112:   PetscInt       i, j, k, xs, ys, xm, ym, zs, zm;
113:   PetscScalar ***x;
114:   DM             da;

116:   PetscFunctionBeginUser;
117:   PetscCall(SNESGetDM(snes, &da));
118:   PetscCall(DMGetApplicationContext(da, &user));
119:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
120:   PetscCall(DMDAVecGetArray(da, X, &x));

122:   /* Compute initial guess */
123:   for (k = zs; k < zs + zm; k++) {
124:     for (j = ys; j < ys + ym; j++) {
125:       for (i = xs; i < xs + xm; i++) x[k][j][i] = user->tleft;
126:     }
127:   }
128:   PetscCall(DMDAVecRestoreArray(da, X, &x));
129:   PetscFunctionReturn(PETSC_SUCCESS);
130: }
131: /* --------------------  Evaluate Function F(x) --------------------- */
132: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr)
133: {
134:   AppCtx        *user = (AppCtx *)ptr;
135:   PetscInt       i, j, k, mx, my, mz, xs, ys, zs, xm, ym, zm;
136:   PetscScalar    zero = 0.0, one = 1.0;
137:   PetscScalar    hx, hy, hz, hxhydhz, hyhzdhx, hzhxdhy;
138:   PetscScalar    t0, tn, ts, te, tw, an, as, ae, aw, dn, ds, de, dw, fn = 0.0, fs = 0.0, fe = 0.0, fw = 0.0;
139:   PetscScalar    tleft, tright, beta, td, ad, dd, fd = 0.0, tu, au, du = 0.0, fu = 0.0;
140:   PetscScalar ***x, ***f;
141:   Vec            localX;
142:   DM             da;

144:   PetscFunctionBeginUser;
145:   PetscCall(SNESGetDM(snes, &da));
146:   PetscCall(DMGetLocalVector(da, &localX));
147:   PetscCall(DMDAGetInfo(da, NULL, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
148:   hx      = one / (PetscReal)(mx - 1);
149:   hy      = one / (PetscReal)(my - 1);
150:   hz      = one / (PetscReal)(mz - 1);
151:   hxhydhz = hx * hy / hz;
152:   hyhzdhx = hy * hz / hx;
153:   hzhxdhy = hz * hx / hy;
154:   tleft   = user->tleft;
155:   tright  = user->tright;
156:   beta    = user->beta;

158:   /* Get ghost points */
159:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
160:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
161:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
162:   PetscCall(DMDAVecGetArray(da, localX, &x));
163:   PetscCall(DMDAVecGetArray(da, F, &f));

165:   /* Evaluate function */
166:   for (k = zs; k < zs + zm; k++) {
167:     for (j = ys; j < ys + ym; j++) {
168:       for (i = xs; i < xs + xm; i++) {
169:         t0 = x[k][j][i];

171:         if (i > 0 && i < mx - 1 && j > 0 && j < my - 1 && k > 0 && k < mz - 1) {
172:           /* general interior volume */

174:           tw = x[k][j][i - 1];
175:           aw = 0.5 * (t0 + tw);
176:           dw = PetscPowScalar(aw, beta);
177:           fw = dw * (t0 - tw);

179:           te = x[k][j][i + 1];
180:           ae = 0.5 * (t0 + te);
181:           de = PetscPowScalar(ae, beta);
182:           fe = de * (te - t0);

184:           ts = x[k][j - 1][i];
185:           as = 0.5 * (t0 + ts);
186:           ds = PetscPowScalar(as, beta);
187:           fs = ds * (t0 - ts);

189:           tn = x[k][j + 1][i];
190:           an = 0.5 * (t0 + tn);
191:           dn = PetscPowScalar(an, beta);
192:           fn = dn * (tn - t0);

194:           td = x[k - 1][j][i];
195:           ad = 0.5 * (t0 + td);
196:           dd = PetscPowScalar(ad, beta);
197:           fd = dd * (t0 - td);

199:           tu = x[k + 1][j][i];
200:           au = 0.5 * (t0 + tu);
201:           du = PetscPowScalar(au, beta);
202:           fu = du * (tu - t0);

204:         } else if (i == 0) {
205:           /* left-hand (west) boundary */
206:           tw = tleft;
207:           aw = 0.5 * (t0 + tw);
208:           dw = PetscPowScalar(aw, beta);
209:           fw = dw * (t0 - tw);

211:           te = x[k][j][i + 1];
212:           ae = 0.5 * (t0 + te);
213:           de = PetscPowScalar(ae, beta);
214:           fe = de * (te - t0);

216:           if (j > 0) {
217:             ts = x[k][j - 1][i];
218:             as = 0.5 * (t0 + ts);
219:             ds = PetscPowScalar(as, beta);
220:             fs = ds * (t0 - ts);
221:           } else {
222:             fs = zero;
223:           }

225:           if (j < my - 1) {
226:             tn = x[k][j + 1][i];
227:             an = 0.5 * (t0 + tn);
228:             dn = PetscPowScalar(an, beta);
229:             fn = dn * (tn - t0);
230:           } else {
231:             fn = zero;
232:           }

234:           if (k > 0) {
235:             td = x[k - 1][j][i];
236:             ad = 0.5 * (t0 + td);
237:             dd = PetscPowScalar(ad, beta);
238:             fd = dd * (t0 - td);
239:           } else {
240:             fd = zero;
241:           }

243:           if (k < mz - 1) {
244:             tu = x[k + 1][j][i];
245:             au = 0.5 * (t0 + tu);
246:             du = PetscPowScalar(au, beta);
247:             fu = du * (tu - t0);
248:           } else {
249:             fu = zero;
250:           }

252:         } else if (i == mx - 1) {
253:           /* right-hand (east) boundary */
254:           tw = x[k][j][i - 1];
255:           aw = 0.5 * (t0 + tw);
256:           dw = PetscPowScalar(aw, beta);
257:           fw = dw * (t0 - tw);

259:           te = tright;
260:           ae = 0.5 * (t0 + te);
261:           de = PetscPowScalar(ae, beta);
262:           fe = de * (te - t0);

264:           if (j > 0) {
265:             ts = x[k][j - 1][i];
266:             as = 0.5 * (t0 + ts);
267:             ds = PetscPowScalar(as, beta);
268:             fs = ds * (t0 - ts);
269:           } else {
270:             fs = zero;
271:           }

273:           if (j < my - 1) {
274:             tn = x[k][j + 1][i];
275:             an = 0.5 * (t0 + tn);
276:             dn = PetscPowScalar(an, beta);
277:             fn = dn * (tn - t0);
278:           } else {
279:             fn = zero;
280:           }

282:           if (k > 0) {
283:             td = x[k - 1][j][i];
284:             ad = 0.5 * (t0 + td);
285:             dd = PetscPowScalar(ad, beta);
286:             fd = dd * (t0 - td);
287:           } else {
288:             fd = zero;
289:           }

291:           if (k < mz - 1) {
292:             tu = x[k + 1][j][i];
293:             au = 0.5 * (t0 + tu);
294:             du = PetscPowScalar(au, beta);
295:             fu = du * (tu - t0);
296:           } else {
297:             fu = zero;
298:           }

300:         } else if (j == 0) {
301:           /* bottom (south) boundary, and i <> 0 or mx-1 */
302:           tw = x[k][j][i - 1];
303:           aw = 0.5 * (t0 + tw);
304:           dw = PetscPowScalar(aw, beta);
305:           fw = dw * (t0 - tw);

307:           te = x[k][j][i + 1];
308:           ae = 0.5 * (t0 + te);
309:           de = PetscPowScalar(ae, beta);
310:           fe = de * (te - t0);

312:           fs = zero;

314:           tn = x[k][j + 1][i];
315:           an = 0.5 * (t0 + tn);
316:           dn = PetscPowScalar(an, beta);
317:           fn = dn * (tn - t0);

319:           if (k > 0) {
320:             td = x[k - 1][j][i];
321:             ad = 0.5 * (t0 + td);
322:             dd = PetscPowScalar(ad, beta);
323:             fd = dd * (t0 - td);
324:           } else {
325:             fd = zero;
326:           }

328:           if (k < mz - 1) {
329:             tu = x[k + 1][j][i];
330:             au = 0.5 * (t0 + tu);
331:             du = PetscPowScalar(au, beta);
332:             fu = du * (tu - t0);
333:           } else {
334:             fu = zero;
335:           }

337:         } else if (j == my - 1) {
338:           /* top (north) boundary, and i <> 0 or mx-1 */
339:           tw = x[k][j][i - 1];
340:           aw = 0.5 * (t0 + tw);
341:           dw = PetscPowScalar(aw, beta);
342:           fw = dw * (t0 - tw);

344:           te = x[k][j][i + 1];
345:           ae = 0.5 * (t0 + te);
346:           de = PetscPowScalar(ae, beta);
347:           fe = de * (te - t0);

349:           ts = x[k][j - 1][i];
350:           as = 0.5 * (t0 + ts);
351:           ds = PetscPowScalar(as, beta);
352:           fs = ds * (t0 - ts);

354:           fn = zero;

356:           if (k > 0) {
357:             td = x[k - 1][j][i];
358:             ad = 0.5 * (t0 + td);
359:             dd = PetscPowScalar(ad, beta);
360:             fd = dd * (t0 - td);
361:           } else {
362:             fd = zero;
363:           }

365:           if (k < mz - 1) {
366:             tu = x[k + 1][j][i];
367:             au = 0.5 * (t0 + tu);
368:             du = PetscPowScalar(au, beta);
369:             fu = du * (tu - t0);
370:           } else {
371:             fu = zero;
372:           }

374:         } else if (k == 0) {
375:           /* down boundary (interior only) */
376:           tw = x[k][j][i - 1];
377:           aw = 0.5 * (t0 + tw);
378:           dw = PetscPowScalar(aw, beta);
379:           fw = dw * (t0 - tw);

381:           te = x[k][j][i + 1];
382:           ae = 0.5 * (t0 + te);
383:           de = PetscPowScalar(ae, beta);
384:           fe = de * (te - t0);

386:           ts = x[k][j - 1][i];
387:           as = 0.5 * (t0 + ts);
388:           ds = PetscPowScalar(as, beta);
389:           fs = ds * (t0 - ts);

391:           tn = x[k][j + 1][i];
392:           an = 0.5 * (t0 + tn);
393:           dn = PetscPowScalar(an, beta);
394:           fn = dn * (tn - t0);

396:           fd = zero;

398:           tu = x[k + 1][j][i];
399:           au = 0.5 * (t0 + tu);
400:           du = PetscPowScalar(au, beta);
401:           fu = du * (tu - t0);

403:         } else if (k == mz - 1) {
404:           /* up boundary (interior only) */
405:           tw = x[k][j][i - 1];
406:           aw = 0.5 * (t0 + tw);
407:           dw = PetscPowScalar(aw, beta);
408:           fw = dw * (t0 - tw);

410:           te = x[k][j][i + 1];
411:           ae = 0.5 * (t0 + te);
412:           de = PetscPowScalar(ae, beta);
413:           fe = de * (te - t0);

415:           ts = x[k][j - 1][i];
416:           as = 0.5 * (t0 + ts);
417:           ds = PetscPowScalar(as, beta);
418:           fs = ds * (t0 - ts);

420:           tn = x[k][j + 1][i];
421:           an = 0.5 * (t0 + tn);
422:           dn = PetscPowScalar(an, beta);
423:           fn = dn * (tn - t0);

425:           td = x[k - 1][j][i];
426:           ad = 0.5 * (t0 + td);
427:           dd = PetscPowScalar(ad, beta);
428:           fd = dd * (t0 - td);

430:           fu = zero;
431:         }

433:         f[k][j][i] = -hyhzdhx * (fe - fw) - hzhxdhy * (fn - fs) - hxhydhz * (fu - fd);
434:       }
435:     }
436:   }
437:   PetscCall(DMDAVecRestoreArray(da, localX, &x));
438:   PetscCall(DMDAVecRestoreArray(da, F, &f));
439:   PetscCall(DMRestoreLocalVector(da, &localX));
440:   PetscCall(PetscLogFlops((22.0 + 4.0 * POWFLOP) * ym * xm));
441:   PetscFunctionReturn(PETSC_SUCCESS);
442: }
443: /* --------------------  Evaluate Jacobian F(x) --------------------- */
444: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr)
445: {
446:   AppCtx        *user = (AppCtx *)ptr;
447:   PetscInt       i, j, k, mx, my, mz, xs, ys, zs, xm, ym, zm;
448:   PetscScalar    one = 1.0;
449:   PetscScalar    hx, hy, hz, hxhydhz, hyhzdhx, hzhxdhy;
450:   PetscScalar    t0, tn, ts, te, tw, an, as, ae, aw, dn, ds, de, dw;
451:   PetscScalar    tleft, tright, beta, td, ad, dd, tu, au, du, v[7], bm1, coef;
452:   PetscScalar ***x, bn, bs, be, bw, bu, bd, gn, gs, ge, gw, gu, gd;
453:   Vec            localX;
454:   MatStencil     c[7], row;
455:   DM             da;

457:   PetscFunctionBeginUser;
458:   PetscCall(SNESGetDM(snes, &da));
459:   PetscCall(DMGetLocalVector(da, &localX));
460:   PetscCall(DMDAGetInfo(da, NULL, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
461:   hx      = one / (PetscReal)(mx - 1);
462:   hy      = one / (PetscReal)(my - 1);
463:   hz      = one / (PetscReal)(mz - 1);
464:   hxhydhz = hx * hy / hz;
465:   hyhzdhx = hy * hz / hx;
466:   hzhxdhy = hz * hx / hy;
467:   tleft   = user->tleft;
468:   tright  = user->tright;
469:   beta    = user->beta;
470:   bm1     = user->bm1;
471:   coef    = user->coef;

473:   /* Get ghost points */
474:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
475:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
476:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
477:   PetscCall(DMDAVecGetArray(da, localX, &x));

479:   /* Evaluate Jacobian of function */
480:   for (k = zs; k < zs + zm; k++) {
481:     for (j = ys; j < ys + ym; j++) {
482:       for (i = xs; i < xs + xm; i++) {
483:         t0    = x[k][j][i];
484:         row.k = k;
485:         row.j = j;
486:         row.i = i;
487:         if (i > 0 && i < mx - 1 && j > 0 && j < my - 1 && k > 0 && k < mz - 1) {
488:           /* general interior volume */

490:           tw = x[k][j][i - 1];
491:           aw = 0.5 * (t0 + tw);
492:           bw = PetscPowScalar(aw, bm1);
493:           /* dw = bw * aw */
494:           dw = PetscPowScalar(aw, beta);
495:           gw = coef * bw * (t0 - tw);

497:           te = x[k][j][i + 1];
498:           ae = 0.5 * (t0 + te);
499:           be = PetscPowScalar(ae, bm1);
500:           /* de = be * ae; */
501:           de = PetscPowScalar(ae, beta);
502:           ge = coef * be * (te - t0);

504:           ts = x[k][j - 1][i];
505:           as = 0.5 * (t0 + ts);
506:           bs = PetscPowScalar(as, bm1);
507:           /* ds = bs * as; */
508:           ds = PetscPowScalar(as, beta);
509:           gs = coef * bs * (t0 - ts);

511:           tn = x[k][j + 1][i];
512:           an = 0.5 * (t0 + tn);
513:           bn = PetscPowScalar(an, bm1);
514:           /* dn = bn * an; */
515:           dn = PetscPowScalar(an, beta);
516:           gn = coef * bn * (tn - t0);

518:           td = x[k - 1][j][i];
519:           ad = 0.5 * (t0 + td);
520:           bd = PetscPowScalar(ad, bm1);
521:           /* dd = bd * ad; */
522:           dd = PetscPowScalar(ad, beta);
523:           gd = coef * bd * (t0 - td);

525:           tu = x[k + 1][j][i];
526:           au = 0.5 * (t0 + tu);
527:           bu = PetscPowScalar(au, bm1);
528:           /* du = bu * au; */
529:           du = PetscPowScalar(au, beta);
530:           gu = coef * bu * (tu - t0);

532:           c[0].k = k - 1;
533:           c[0].j = j;
534:           c[0].i = i;
535:           v[0]   = -hxhydhz * (dd - gd);
536:           c[1].k = k;
537:           c[1].j = j - 1;
538:           c[1].i = i;
539:           v[1]   = -hzhxdhy * (ds - gs);
540:           c[2].k = k;
541:           c[2].j = j;
542:           c[2].i = i - 1;
543:           v[2]   = -hyhzdhx * (dw - gw);
544:           c[3].k = k;
545:           c[3].j = j;
546:           c[3].i = i;
547:           v[3]   = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
548:           c[4].k = k;
549:           c[4].j = j;
550:           c[4].i = i + 1;
551:           v[4]   = -hyhzdhx * (de + ge);
552:           c[5].k = k;
553:           c[5].j = j + 1;
554:           c[5].i = i;
555:           v[5]   = -hzhxdhy * (dn + gn);
556:           c[6].k = k + 1;
557:           c[6].j = j;
558:           c[6].i = i;
559:           v[6]   = -hxhydhz * (du + gu);
560:           PetscCall(MatSetValuesStencil(jac, 1, &row, 7, c, v, INSERT_VALUES));

562:         } else if (i == 0) {
563:           /* left-hand plane boundary */
564:           tw = tleft;
565:           aw = 0.5 * (t0 + tw);
566:           bw = PetscPowScalar(aw, bm1);
567:           /* dw = bw * aw */
568:           dw = PetscPowScalar(aw, beta);
569:           gw = coef * bw * (t0 - tw);

571:           te = x[k][j][i + 1];
572:           ae = 0.5 * (t0 + te);
573:           be = PetscPowScalar(ae, bm1);
574:           /* de = be * ae; */
575:           de = PetscPowScalar(ae, beta);
576:           ge = coef * be * (te - t0);

578:           /* left-hand bottom edge */
579:           if (j == 0) {
580:             tn = x[k][j + 1][i];
581:             an = 0.5 * (t0 + tn);
582:             bn = PetscPowScalar(an, bm1);
583:             /* dn = bn * an; */
584:             dn = PetscPowScalar(an, beta);
585:             gn = coef * bn * (tn - t0);

587:             /* left-hand bottom down corner */
588:             if (k == 0) {
589:               tu = x[k + 1][j][i];
590:               au = 0.5 * (t0 + tu);
591:               bu = PetscPowScalar(au, bm1);
592:               /* du = bu * au; */
593:               du = PetscPowScalar(au, beta);
594:               gu = coef * bu * (tu - t0);

596:               c[0].k = k;
597:               c[0].j = j;
598:               c[0].i = i;
599:               v[0]   = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
600:               c[1].k = k;
601:               c[1].j = j;
602:               c[1].i = i + 1;
603:               v[1]   = -hyhzdhx * (de + ge);
604:               c[2].k = k;
605:               c[2].j = j + 1;
606:               c[2].i = i;
607:               v[2]   = -hzhxdhy * (dn + gn);
608:               c[3].k = k + 1;
609:               c[3].j = j;
610:               c[3].i = i;
611:               v[3]   = -hxhydhz * (du + gu);
612:               PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));

614:               /* left-hand bottom interior edge */
615:             } else if (k < mz - 1) {
616:               tu = x[k + 1][j][i];
617:               au = 0.5 * (t0 + tu);
618:               bu = PetscPowScalar(au, bm1);
619:               /* du = bu * au; */
620:               du = PetscPowScalar(au, beta);
621:               gu = coef * bu * (tu - t0);

623:               td = x[k - 1][j][i];
624:               ad = 0.5 * (t0 + td);
625:               bd = PetscPowScalar(ad, bm1);
626:               /* dd = bd * ad; */
627:               dd = PetscPowScalar(ad, beta);
628:               gd = coef * bd * (td - t0);

630:               c[0].k = k - 1;
631:               c[0].j = j;
632:               c[0].i = i;
633:               v[0]   = -hxhydhz * (dd - gd);
634:               c[1].k = k;
635:               c[1].j = j;
636:               c[1].i = i;
637:               v[1]   = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
638:               c[2].k = k;
639:               c[2].j = j;
640:               c[2].i = i + 1;
641:               v[2]   = -hyhzdhx * (de + ge);
642:               c[3].k = k;
643:               c[3].j = j + 1;
644:               c[3].i = i;
645:               v[3]   = -hzhxdhy * (dn + gn);
646:               c[4].k = k + 1;
647:               c[4].j = j;
648:               c[4].i = i;
649:               v[4]   = -hxhydhz * (du + gu);
650:               PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));

652:               /* left-hand bottom up corner */
653:             } else {
654:               td = x[k - 1][j][i];
655:               ad = 0.5 * (t0 + td);
656:               bd = PetscPowScalar(ad, bm1);
657:               /* dd = bd * ad; */
658:               dd = PetscPowScalar(ad, beta);
659:               gd = coef * bd * (td - t0);

661:               c[0].k = k - 1;
662:               c[0].j = j;
663:               c[0].i = i;
664:               v[0]   = -hxhydhz * (dd - gd);
665:               c[1].k = k;
666:               c[1].j = j;
667:               c[1].i = i;
668:               v[1]   = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
669:               c[2].k = k;
670:               c[2].j = j;
671:               c[2].i = i + 1;
672:               v[2]   = -hyhzdhx * (de + ge);
673:               c[3].k = k;
674:               c[3].j = j + 1;
675:               c[3].i = i;
676:               v[3]   = -hzhxdhy * (dn + gn);
677:               PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
678:             }

680:             /* left-hand top edge */
681:           } else if (j == my - 1) {
682:             ts = x[k][j - 1][i];
683:             as = 0.5 * (t0 + ts);
684:             bs = PetscPowScalar(as, bm1);
685:             /* ds = bs * as; */
686:             ds = PetscPowScalar(as, beta);
687:             gs = coef * bs * (ts - t0);

689:             /* left-hand top down corner */
690:             if (k == 0) {
691:               tu = x[k + 1][j][i];
692:               au = 0.5 * (t0 + tu);
693:               bu = PetscPowScalar(au, bm1);
694:               /* du = bu * au; */
695:               du = PetscPowScalar(au, beta);
696:               gu = coef * bu * (tu - t0);

698:               c[0].k = k;
699:               c[0].j = j - 1;
700:               c[0].i = i;
701:               v[0]   = -hzhxdhy * (ds - gs);
702:               c[1].k = k;
703:               c[1].j = j;
704:               c[1].i = i;
705:               v[1]   = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
706:               c[2].k = k;
707:               c[2].j = j;
708:               c[2].i = i + 1;
709:               v[2]   = -hyhzdhx * (de + ge);
710:               c[3].k = k + 1;
711:               c[3].j = j;
712:               c[3].i = i;
713:               v[3]   = -hxhydhz * (du + gu);
714:               PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));

716:               /* left-hand top interior edge */
717:             } else if (k < mz - 1) {
718:               tu = x[k + 1][j][i];
719:               au = 0.5 * (t0 + tu);
720:               bu = PetscPowScalar(au, bm1);
721:               /* du = bu * au; */
722:               du = PetscPowScalar(au, beta);
723:               gu = coef * bu * (tu - t0);

725:               td = x[k - 1][j][i];
726:               ad = 0.5 * (t0 + td);
727:               bd = PetscPowScalar(ad, bm1);
728:               /* dd = bd * ad; */
729:               dd = PetscPowScalar(ad, beta);
730:               gd = coef * bd * (td - t0);

732:               c[0].k = k - 1;
733:               c[0].j = j;
734:               c[0].i = i;
735:               v[0]   = -hxhydhz * (dd - gd);
736:               c[1].k = k;
737:               c[1].j = j - 1;
738:               c[1].i = i;
739:               v[1]   = -hzhxdhy * (ds - gs);
740:               c[2].k = k;
741:               c[2].j = j;
742:               c[2].i = i;
743:               v[2]   = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
744:               c[3].k = k;
745:               c[3].j = j;
746:               c[3].i = i + 1;
747:               v[3]   = -hyhzdhx * (de + ge);
748:               c[4].k = k + 1;
749:               c[4].j = j;
750:               c[4].i = i;
751:               v[4]   = -hxhydhz * (du + gu);
752:               PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));

754:               /* left-hand top up corner */
755:             } else {
756:               td = x[k - 1][j][i];
757:               ad = 0.5 * (t0 + td);
758:               bd = PetscPowScalar(ad, bm1);
759:               /* dd = bd * ad; */
760:               dd = PetscPowScalar(ad, beta);
761:               gd = coef * bd * (td - t0);

763:               c[0].k = k - 1;
764:               c[0].j = j;
765:               c[0].i = i;
766:               v[0]   = -hxhydhz * (dd - gd);
767:               c[1].k = k;
768:               c[1].j = j - 1;
769:               c[1].i = i;
770:               v[1]   = -hzhxdhy * (ds - gs);
771:               c[2].k = k;
772:               c[2].j = j;
773:               c[2].i = i;
774:               v[2]   = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
775:               c[3].k = k;
776:               c[3].j = j;
777:               c[3].i = i + 1;
778:               v[3]   = -hyhzdhx * (de + ge);
779:               PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
780:             }

782:           } else {
783:             ts = x[k][j - 1][i];
784:             as = 0.5 * (t0 + ts);
785:             bs = PetscPowScalar(as, bm1);
786:             /* ds = bs * as; */
787:             ds = PetscPowScalar(as, beta);
788:             gs = coef * bs * (t0 - ts);

790:             tn = x[k][j + 1][i];
791:             an = 0.5 * (t0 + tn);
792:             bn = PetscPowScalar(an, bm1);
793:             /* dn = bn * an; */
794:             dn = PetscPowScalar(an, beta);
795:             gn = coef * bn * (tn - t0);

797:             /* left-hand down interior edge */
798:             if (k == 0) {
799:               tu = x[k + 1][j][i];
800:               au = 0.5 * (t0 + tu);
801:               bu = PetscPowScalar(au, bm1);
802:               /* du = bu * au; */
803:               du = PetscPowScalar(au, beta);
804:               gu = coef * bu * (tu - t0);

806:               c[0].k = k;
807:               c[0].j = j - 1;
808:               c[0].i = i;
809:               v[0]   = -hzhxdhy * (ds - gs);
810:               c[1].k = k;
811:               c[1].j = j;
812:               c[1].i = i;
813:               v[1]   = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
814:               c[2].k = k;
815:               c[2].j = j;
816:               c[2].i = i + 1;
817:               v[2]   = -hyhzdhx * (de + ge);
818:               c[3].k = k;
819:               c[3].j = j + 1;
820:               c[3].i = i;
821:               v[3]   = -hzhxdhy * (dn + gn);
822:               c[4].k = k + 1;
823:               c[4].j = j;
824:               c[4].i = i;
825:               v[4]   = -hxhydhz * (du + gu);
826:               PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));

828:             } else if (k == mz - 1) { /* left-hand up interior edge */

830:               td = x[k - 1][j][i];
831:               ad = 0.5 * (t0 + td);
832:               bd = PetscPowScalar(ad, bm1);
833:               /* dd = bd * ad; */
834:               dd = PetscPowScalar(ad, beta);
835:               gd = coef * bd * (t0 - td);

837:               c[0].k = k - 1;
838:               c[0].j = j;
839:               c[0].i = i;
840:               v[0]   = -hxhydhz * (dd - gd);
841:               c[1].k = k;
842:               c[1].j = j - 1;
843:               c[1].i = i;
844:               v[1]   = -hzhxdhy * (ds - gs);
845:               c[2].k = k;
846:               c[2].j = j;
847:               c[2].i = i;
848:               v[2]   = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
849:               c[3].k = k;
850:               c[3].j = j;
851:               c[3].i = i + 1;
852:               v[3]   = -hyhzdhx * (de + ge);
853:               c[4].k = k;
854:               c[4].j = j + 1;
855:               c[4].i = i;
856:               v[4]   = -hzhxdhy * (dn + gn);
857:               PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
858:             } else { /* left-hand interior plane */

860:               td = x[k - 1][j][i];
861:               ad = 0.5 * (t0 + td);
862:               bd = PetscPowScalar(ad, bm1);
863:               /* dd = bd * ad; */
864:               dd = PetscPowScalar(ad, beta);
865:               gd = coef * bd * (t0 - td);

867:               tu = x[k + 1][j][i];
868:               au = 0.5 * (t0 + tu);
869:               bu = PetscPowScalar(au, bm1);
870:               /* du = bu * au; */
871:               du = PetscPowScalar(au, beta);
872:               gu = coef * bu * (tu - t0);

874:               c[0].k = k - 1;
875:               c[0].j = j;
876:               c[0].i = i;
877:               v[0]   = -hxhydhz * (dd - gd);
878:               c[1].k = k;
879:               c[1].j = j - 1;
880:               c[1].i = i;
881:               v[1]   = -hzhxdhy * (ds - gs);
882:               c[2].k = k;
883:               c[2].j = j;
884:               c[2].i = i;
885:               v[2]   = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
886:               c[3].k = k;
887:               c[3].j = j;
888:               c[3].i = i + 1;
889:               v[3]   = -hyhzdhx * (de + ge);
890:               c[4].k = k;
891:               c[4].j = j + 1;
892:               c[4].i = i;
893:               v[4]   = -hzhxdhy * (dn + gn);
894:               c[5].k = k + 1;
895:               c[5].j = j;
896:               c[5].i = i;
897:               v[5]   = -hxhydhz * (du + gu);
898:               PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
899:             }
900:           }

902:         } else if (i == mx - 1) {
903:           /* right-hand plane boundary */
904:           tw = x[k][j][i - 1];
905:           aw = 0.5 * (t0 + tw);
906:           bw = PetscPowScalar(aw, bm1);
907:           /* dw = bw * aw */
908:           dw = PetscPowScalar(aw, beta);
909:           gw = coef * bw * (t0 - tw);

911:           te = tright;
912:           ae = 0.5 * (t0 + te);
913:           be = PetscPowScalar(ae, bm1);
914:           /* de = be * ae; */
915:           de = PetscPowScalar(ae, beta);
916:           ge = coef * be * (te - t0);

918:           /* right-hand bottom edge */
919:           if (j == 0) {
920:             tn = x[k][j + 1][i];
921:             an = 0.5 * (t0 + tn);
922:             bn = PetscPowScalar(an, bm1);
923:             /* dn = bn * an; */
924:             dn = PetscPowScalar(an, beta);
925:             gn = coef * bn * (tn - t0);

927:             /* right-hand bottom down corner */
928:             if (k == 0) {
929:               tu = x[k + 1][j][i];
930:               au = 0.5 * (t0 + tu);
931:               bu = PetscPowScalar(au, bm1);
932:               /* du = bu * au; */
933:               du = PetscPowScalar(au, beta);
934:               gu = coef * bu * (tu - t0);

936:               c[0].k = k;
937:               c[0].j = j;
938:               c[0].i = i - 1;
939:               v[0]   = -hyhzdhx * (dw - gw);
940:               c[1].k = k;
941:               c[1].j = j;
942:               c[1].i = i;
943:               v[1]   = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
944:               c[2].k = k;
945:               c[2].j = j + 1;
946:               c[2].i = i;
947:               v[2]   = -hzhxdhy * (dn + gn);
948:               c[3].k = k + 1;
949:               c[3].j = j;
950:               c[3].i = i;
951:               v[3]   = -hxhydhz * (du + gu);
952:               PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));

954:               /* right-hand bottom interior edge */
955:             } else if (k < mz - 1) {
956:               tu = x[k + 1][j][i];
957:               au = 0.5 * (t0 + tu);
958:               bu = PetscPowScalar(au, bm1);
959:               /* du = bu * au; */
960:               du = PetscPowScalar(au, beta);
961:               gu = coef * bu * (tu - t0);

963:               td = x[k - 1][j][i];
964:               ad = 0.5 * (t0 + td);
965:               bd = PetscPowScalar(ad, bm1);
966:               /* dd = bd * ad; */
967:               dd = PetscPowScalar(ad, beta);
968:               gd = coef * bd * (td - t0);

970:               c[0].k = k - 1;
971:               c[0].j = j;
972:               c[0].i = i;
973:               v[0]   = -hxhydhz * (dd - gd);
974:               c[1].k = k;
975:               c[1].j = j;
976:               c[1].i = i - 1;
977:               v[1]   = -hyhzdhx * (dw - gw);
978:               c[2].k = k;
979:               c[2].j = j;
980:               c[2].i = i;
981:               v[2]   = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
982:               c[3].k = k;
983:               c[3].j = j + 1;
984:               c[3].i = i;
985:               v[3]   = -hzhxdhy * (dn + gn);
986:               c[4].k = k + 1;
987:               c[4].j = j;
988:               c[4].i = i;
989:               v[4]   = -hxhydhz * (du + gu);
990:               PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));

992:               /* right-hand bottom up corner */
993:             } else {
994:               td = x[k - 1][j][i];
995:               ad = 0.5 * (t0 + td);
996:               bd = PetscPowScalar(ad, bm1);
997:               /* dd = bd * ad; */
998:               dd = PetscPowScalar(ad, beta);
999:               gd = coef * bd * (td - t0);

1001:               c[0].k = k - 1;
1002:               c[0].j = j;
1003:               c[0].i = i;
1004:               v[0]   = -hxhydhz * (dd - gd);
1005:               c[1].k = k;
1006:               c[1].j = j;
1007:               c[1].i = i - 1;
1008:               v[1]   = -hyhzdhx * (dw - gw);
1009:               c[2].k = k;
1010:               c[2].j = j;
1011:               c[2].i = i;
1012:               v[2]   = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1013:               c[3].k = k;
1014:               c[3].j = j + 1;
1015:               c[3].i = i;
1016:               v[3]   = -hzhxdhy * (dn + gn);
1017:               PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
1018:             }

1020:             /* right-hand top edge */
1021:           } else if (j == my - 1) {
1022:             ts = x[k][j - 1][i];
1023:             as = 0.5 * (t0 + ts);
1024:             bs = PetscPowScalar(as, bm1);
1025:             /* ds = bs * as; */
1026:             ds = PetscPowScalar(as, beta);
1027:             gs = coef * bs * (ts - t0);

1029:             /* right-hand top down corner */
1030:             if (k == 0) {
1031:               tu = x[k + 1][j][i];
1032:               au = 0.5 * (t0 + tu);
1033:               bu = PetscPowScalar(au, bm1);
1034:               /* du = bu * au; */
1035:               du = PetscPowScalar(au, beta);
1036:               gu = coef * bu * (tu - t0);

1038:               c[0].k = k;
1039:               c[0].j = j - 1;
1040:               c[0].i = i;
1041:               v[0]   = -hzhxdhy * (ds - gs);
1042:               c[1].k = k;
1043:               c[1].j = j;
1044:               c[1].i = i - 1;
1045:               v[1]   = -hyhzdhx * (dw - gw);
1046:               c[2].k = k;
1047:               c[2].j = j;
1048:               c[2].i = i;
1049:               v[2]   = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
1050:               c[3].k = k + 1;
1051:               c[3].j = j;
1052:               c[3].i = i;
1053:               v[3]   = -hxhydhz * (du + gu);
1054:               PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));

1056:               /* right-hand top interior edge */
1057:             } else if (k < mz - 1) {
1058:               tu = x[k + 1][j][i];
1059:               au = 0.5 * (t0 + tu);
1060:               bu = PetscPowScalar(au, bm1);
1061:               /* du = bu * au; */
1062:               du = PetscPowScalar(au, beta);
1063:               gu = coef * bu * (tu - t0);

1065:               td = x[k - 1][j][i];
1066:               ad = 0.5 * (t0 + td);
1067:               bd = PetscPowScalar(ad, bm1);
1068:               /* dd = bd * ad; */
1069:               dd = PetscPowScalar(ad, beta);
1070:               gd = coef * bd * (td - t0);

1072:               c[0].k = k - 1;
1073:               c[0].j = j;
1074:               c[0].i = i;
1075:               v[0]   = -hxhydhz * (dd - gd);
1076:               c[1].k = k;
1077:               c[1].j = j - 1;
1078:               c[1].i = i;
1079:               v[1]   = -hzhxdhy * (ds - gs);
1080:               c[2].k = k;
1081:               c[2].j = j;
1082:               c[2].i = i - 1;
1083:               v[2]   = -hyhzdhx * (dw - gw);
1084:               c[3].k = k;
1085:               c[3].j = j;
1086:               c[3].i = i;
1087:               v[3]   = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
1088:               c[4].k = k + 1;
1089:               c[4].j = j;
1090:               c[4].i = i;
1091:               v[4]   = -hxhydhz * (du + gu);
1092:               PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));

1094:               /* right-hand top up corner */
1095:             } else {
1096:               td = x[k - 1][j][i];
1097:               ad = 0.5 * (t0 + td);
1098:               bd = PetscPowScalar(ad, bm1);
1099:               /* dd = bd * ad; */
1100:               dd = PetscPowScalar(ad, beta);
1101:               gd = coef * bd * (td - t0);

1103:               c[0].k = k - 1;
1104:               c[0].j = j;
1105:               c[0].i = i;
1106:               v[0]   = -hxhydhz * (dd - gd);
1107:               c[1].k = k;
1108:               c[1].j = j - 1;
1109:               c[1].i = i;
1110:               v[1]   = -hzhxdhy * (ds - gs);
1111:               c[2].k = k;
1112:               c[2].j = j;
1113:               c[2].i = i - 1;
1114:               v[2]   = -hyhzdhx * (dw - gw);
1115:               c[3].k = k;
1116:               c[3].j = j;
1117:               c[3].i = i;
1118:               v[3]   = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1119:               PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
1120:             }

1122:           } else {
1123:             ts = x[k][j - 1][i];
1124:             as = 0.5 * (t0 + ts);
1125:             bs = PetscPowScalar(as, bm1);
1126:             /* ds = bs * as; */
1127:             ds = PetscPowScalar(as, beta);
1128:             gs = coef * bs * (t0 - ts);

1130:             tn = x[k][j + 1][i];
1131:             an = 0.5 * (t0 + tn);
1132:             bn = PetscPowScalar(an, bm1);
1133:             /* dn = bn * an; */
1134:             dn = PetscPowScalar(an, beta);
1135:             gn = coef * bn * (tn - t0);

1137:             /* right-hand down interior edge */
1138:             if (k == 0) {
1139:               tu = x[k + 1][j][i];
1140:               au = 0.5 * (t0 + tu);
1141:               bu = PetscPowScalar(au, bm1);
1142:               /* du = bu * au; */
1143:               du = PetscPowScalar(au, beta);
1144:               gu = coef * bu * (tu - t0);

1146:               c[0].k = k;
1147:               c[0].j = j - 1;
1148:               c[0].i = i;
1149:               v[0]   = -hzhxdhy * (ds - gs);
1150:               c[1].k = k;
1151:               c[1].j = j;
1152:               c[1].i = i - 1;
1153:               v[1]   = -hyhzdhx * (dw - gw);
1154:               c[2].k = k;
1155:               c[2].j = j;
1156:               c[2].i = i;
1157:               v[2]   = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
1158:               c[3].k = k;
1159:               c[3].j = j + 1;
1160:               c[3].i = i;
1161:               v[3]   = -hzhxdhy * (dn + gn);
1162:               c[4].k = k + 1;
1163:               c[4].j = j;
1164:               c[4].i = i;
1165:               v[4]   = -hxhydhz * (du + gu);
1166:               PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));

1168:             } else if (k == mz - 1) { /* right-hand up interior edge */

1170:               td = x[k - 1][j][i];
1171:               ad = 0.5 * (t0 + td);
1172:               bd = PetscPowScalar(ad, bm1);
1173:               /* dd = bd * ad; */
1174:               dd = PetscPowScalar(ad, beta);
1175:               gd = coef * bd * (t0 - td);

1177:               c[0].k = k - 1;
1178:               c[0].j = j;
1179:               c[0].i = i;
1180:               v[0]   = -hxhydhz * (dd - gd);
1181:               c[1].k = k;
1182:               c[1].j = j - 1;
1183:               c[1].i = i;
1184:               v[1]   = -hzhxdhy * (ds - gs);
1185:               c[2].k = k;
1186:               c[2].j = j;
1187:               c[2].i = i - 1;
1188:               v[2]   = -hyhzdhx * (dw - gw);
1189:               c[3].k = k;
1190:               c[3].j = j;
1191:               c[3].i = i;
1192:               v[3]   = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1193:               c[4].k = k;
1194:               c[4].j = j + 1;
1195:               c[4].i = i;
1196:               v[4]   = -hzhxdhy * (dn + gn);
1197:               PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));

1199:             } else { /* right-hand interior plane */

1201:               td = x[k - 1][j][i];
1202:               ad = 0.5 * (t0 + td);
1203:               bd = PetscPowScalar(ad, bm1);
1204:               /* dd = bd * ad; */
1205:               dd = PetscPowScalar(ad, beta);
1206:               gd = coef * bd * (t0 - td);

1208:               tu = x[k + 1][j][i];
1209:               au = 0.5 * (t0 + tu);
1210:               bu = PetscPowScalar(au, bm1);
1211:               /* du = bu * au; */
1212:               du = PetscPowScalar(au, beta);
1213:               gu = coef * bu * (tu - t0);

1215:               c[0].k = k - 1;
1216:               c[0].j = j;
1217:               c[0].i = i;
1218:               v[0]   = -hxhydhz * (dd - gd);
1219:               c[1].k = k;
1220:               c[1].j = j - 1;
1221:               c[1].i = i;
1222:               v[1]   = -hzhxdhy * (ds - gs);
1223:               c[2].k = k;
1224:               c[2].j = j;
1225:               c[2].i = i - 1;
1226:               v[2]   = -hyhzdhx * (dw - gw);
1227:               c[3].k = k;
1228:               c[3].j = j;
1229:               c[3].i = i;
1230:               v[3]   = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
1231:               c[4].k = k;
1232:               c[4].j = j + 1;
1233:               c[4].i = i;
1234:               v[4]   = -hzhxdhy * (dn + gn);
1235:               c[5].k = k + 1;
1236:               c[5].j = j;
1237:               c[5].i = i;
1238:               v[5]   = -hxhydhz * (du + gu);
1239:               PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
1240:             }
1241:           }

1243:         } else if (j == 0) {
1244:           tw = x[k][j][i - 1];
1245:           aw = 0.5 * (t0 + tw);
1246:           bw = PetscPowScalar(aw, bm1);
1247:           /* dw = bw * aw */
1248:           dw = PetscPowScalar(aw, beta);
1249:           gw = coef * bw * (t0 - tw);

1251:           te = x[k][j][i + 1];
1252:           ae = 0.5 * (t0 + te);
1253:           be = PetscPowScalar(ae, bm1);
1254:           /* de = be * ae; */
1255:           de = PetscPowScalar(ae, beta);
1256:           ge = coef * be * (te - t0);

1258:           tn = x[k][j + 1][i];
1259:           an = 0.5 * (t0 + tn);
1260:           bn = PetscPowScalar(an, bm1);
1261:           /* dn = bn * an; */
1262:           dn = PetscPowScalar(an, beta);
1263:           gn = coef * bn * (tn - t0);

1265:           /* bottom down interior edge */
1266:           if (k == 0) {
1267:             tu = x[k + 1][j][i];
1268:             au = 0.5 * (t0 + tu);
1269:             bu = PetscPowScalar(au, bm1);
1270:             /* du = bu * au; */
1271:             du = PetscPowScalar(au, beta);
1272:             gu = coef * bu * (tu - t0);

1274:             c[0].k = k;
1275:             c[0].j = j;
1276:             c[0].i = i - 1;
1277:             v[0]   = -hyhzdhx * (dw - gw);
1278:             c[1].k = k;
1279:             c[1].j = j;
1280:             c[1].i = i;
1281:             v[1]   = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
1282:             c[2].k = k;
1283:             c[2].j = j;
1284:             c[2].i = i + 1;
1285:             v[2]   = -hyhzdhx * (de + ge);
1286:             c[3].k = k;
1287:             c[3].j = j + 1;
1288:             c[3].i = i;
1289:             v[3]   = -hzhxdhy * (dn + gn);
1290:             c[4].k = k + 1;
1291:             c[4].j = j;
1292:             c[4].i = i;
1293:             v[4]   = -hxhydhz * (du + gu);
1294:             PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));

1296:           } else if (k == mz - 1) { /* bottom up interior edge */

1298:             td = x[k - 1][j][i];
1299:             ad = 0.5 * (t0 + td);
1300:             bd = PetscPowScalar(ad, bm1);
1301:             /* dd = bd * ad; */
1302:             dd = PetscPowScalar(ad, beta);
1303:             gd = coef * bd * (td - t0);

1305:             c[0].k = k - 1;
1306:             c[0].j = j;
1307:             c[0].i = i;
1308:             v[0]   = -hxhydhz * (dd - gd);
1309:             c[1].k = k;
1310:             c[1].j = j;
1311:             c[1].i = i - 1;
1312:             v[1]   = -hyhzdhx * (dw - gw);
1313:             c[2].k = k;
1314:             c[2].j = j;
1315:             c[2].i = i;
1316:             v[2]   = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1317:             c[3].k = k;
1318:             c[3].j = j;
1319:             c[3].i = i + 1;
1320:             v[3]   = -hyhzdhx * (de + ge);
1321:             c[4].k = k;
1322:             c[4].j = j + 1;
1323:             c[4].i = i;
1324:             v[4]   = -hzhxdhy * (dn + gn);
1325:             PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));

1327:           } else { /* bottom interior plane */

1329:             tu = x[k + 1][j][i];
1330:             au = 0.5 * (t0 + tu);
1331:             bu = PetscPowScalar(au, bm1);
1332:             /* du = bu * au; */
1333:             du = PetscPowScalar(au, beta);
1334:             gu = coef * bu * (tu - t0);

1336:             td = x[k - 1][j][i];
1337:             ad = 0.5 * (t0 + td);
1338:             bd = PetscPowScalar(ad, bm1);
1339:             /* dd = bd * ad; */
1340:             dd = PetscPowScalar(ad, beta);
1341:             gd = coef * bd * (td - t0);

1343:             c[0].k = k - 1;
1344:             c[0].j = j;
1345:             c[0].i = i;
1346:             v[0]   = -hxhydhz * (dd - gd);
1347:             c[1].k = k;
1348:             c[1].j = j;
1349:             c[1].i = i - 1;
1350:             v[1]   = -hyhzdhx * (dw - gw);
1351:             c[2].k = k;
1352:             c[2].j = j;
1353:             c[2].i = i;
1354:             v[2]   = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
1355:             c[3].k = k;
1356:             c[3].j = j;
1357:             c[3].i = i + 1;
1358:             v[3]   = -hyhzdhx * (de + ge);
1359:             c[4].k = k;
1360:             c[4].j = j + 1;
1361:             c[4].i = i;
1362:             v[4]   = -hzhxdhy * (dn + gn);
1363:             c[5].k = k + 1;
1364:             c[5].j = j;
1365:             c[5].i = i;
1366:             v[5]   = -hxhydhz * (du + gu);
1367:             PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
1368:           }

1370:         } else if (j == my - 1) {
1371:           tw = x[k][j][i - 1];
1372:           aw = 0.5 * (t0 + tw);
1373:           bw = PetscPowScalar(aw, bm1);
1374:           /* dw = bw * aw */
1375:           dw = PetscPowScalar(aw, beta);
1376:           gw = coef * bw * (t0 - tw);

1378:           te = x[k][j][i + 1];
1379:           ae = 0.5 * (t0 + te);
1380:           be = PetscPowScalar(ae, bm1);
1381:           /* de = be * ae; */
1382:           de = PetscPowScalar(ae, beta);
1383:           ge = coef * be * (te - t0);

1385:           ts = x[k][j - 1][i];
1386:           as = 0.5 * (t0 + ts);
1387:           bs = PetscPowScalar(as, bm1);
1388:           /* ds = bs * as; */
1389:           ds = PetscPowScalar(as, beta);
1390:           gs = coef * bs * (t0 - ts);

1392:           /* top down interior edge */
1393:           if (k == 0) {
1394:             tu = x[k + 1][j][i];
1395:             au = 0.5 * (t0 + tu);
1396:             bu = PetscPowScalar(au, bm1);
1397:             /* du = bu * au; */
1398:             du = PetscPowScalar(au, beta);
1399:             gu = coef * bu * (tu - t0);

1401:             c[0].k = k;
1402:             c[0].j = j - 1;
1403:             c[0].i = i;
1404:             v[0]   = -hzhxdhy * (ds - gs);
1405:             c[1].k = k;
1406:             c[1].j = j;
1407:             c[1].i = i - 1;
1408:             v[1]   = -hyhzdhx * (dw - gw);
1409:             c[2].k = k;
1410:             c[2].j = j;
1411:             c[2].i = i;
1412:             v[2]   = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
1413:             c[3].k = k;
1414:             c[3].j = j;
1415:             c[3].i = i + 1;
1416:             v[3]   = -hyhzdhx * (de + ge);
1417:             c[4].k = k + 1;
1418:             c[4].j = j;
1419:             c[4].i = i;
1420:             v[4]   = -hxhydhz * (du + gu);
1421:             PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));

1423:           } else if (k == mz - 1) { /* top up interior edge */

1425:             td = x[k - 1][j][i];
1426:             ad = 0.5 * (t0 + td);
1427:             bd = PetscPowScalar(ad, bm1);
1428:             /* dd = bd * ad; */
1429:             dd = PetscPowScalar(ad, beta);
1430:             gd = coef * bd * (td - t0);

1432:             c[0].k = k - 1;
1433:             c[0].j = j;
1434:             c[0].i = i;
1435:             v[0]   = -hxhydhz * (dd - gd);
1436:             c[1].k = k;
1437:             c[1].j = j - 1;
1438:             c[1].i = i;
1439:             v[1]   = -hzhxdhy * (ds - gs);
1440:             c[2].k = k;
1441:             c[2].j = j;
1442:             c[2].i = i - 1;
1443:             v[2]   = -hyhzdhx * (dw - gw);
1444:             c[3].k = k;
1445:             c[3].j = j;
1446:             c[3].i = i;
1447:             v[3]   = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1448:             c[4].k = k;
1449:             c[4].j = j;
1450:             c[4].i = i + 1;
1451:             v[4]   = -hyhzdhx * (de + ge);
1452:             PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));

1454:           } else { /* top interior plane */

1456:             tu = x[k + 1][j][i];
1457:             au = 0.5 * (t0 + tu);
1458:             bu = PetscPowScalar(au, bm1);
1459:             /* du = bu * au; */
1460:             du = PetscPowScalar(au, beta);
1461:             gu = coef * bu * (tu - t0);

1463:             td = x[k - 1][j][i];
1464:             ad = 0.5 * (t0 + td);
1465:             bd = PetscPowScalar(ad, bm1);
1466:             /* dd = bd * ad; */
1467:             dd = PetscPowScalar(ad, beta);
1468:             gd = coef * bd * (td - t0);

1470:             c[0].k = k - 1;
1471:             c[0].j = j;
1472:             c[0].i = i;
1473:             v[0]   = -hxhydhz * (dd - gd);
1474:             c[1].k = k;
1475:             c[1].j = j - 1;
1476:             c[1].i = i;
1477:             v[1]   = -hzhxdhy * (ds - gs);
1478:             c[2].k = k;
1479:             c[2].j = j;
1480:             c[2].i = i - 1;
1481:             v[2]   = -hyhzdhx * (dw - gw);
1482:             c[3].k = k;
1483:             c[3].j = j;
1484:             c[3].i = i;
1485:             v[3]   = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
1486:             c[4].k = k;
1487:             c[4].j = j;
1488:             c[4].i = i + 1;
1489:             v[4]   = -hyhzdhx * (de + ge);
1490:             c[5].k = k + 1;
1491:             c[5].j = j;
1492:             c[5].i = i;
1493:             v[5]   = -hxhydhz * (du + gu);
1494:             PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
1495:           }

1497:         } else if (k == 0) {
1498:           /* down interior plane */

1500:           tw = x[k][j][i - 1];
1501:           aw = 0.5 * (t0 + tw);
1502:           bw = PetscPowScalar(aw, bm1);
1503:           /* dw = bw * aw */
1504:           dw = PetscPowScalar(aw, beta);
1505:           gw = coef * bw * (t0 - tw);

1507:           te = x[k][j][i + 1];
1508:           ae = 0.5 * (t0 + te);
1509:           be = PetscPowScalar(ae, bm1);
1510:           /* de = be * ae; */
1511:           de = PetscPowScalar(ae, beta);
1512:           ge = coef * be * (te - t0);

1514:           ts = x[k][j - 1][i];
1515:           as = 0.5 * (t0 + ts);
1516:           bs = PetscPowScalar(as, bm1);
1517:           /* ds = bs * as; */
1518:           ds = PetscPowScalar(as, beta);
1519:           gs = coef * bs * (t0 - ts);

1521:           tn = x[k][j + 1][i];
1522:           an = 0.5 * (t0 + tn);
1523:           bn = PetscPowScalar(an, bm1);
1524:           /* dn = bn * an; */
1525:           dn = PetscPowScalar(an, beta);
1526:           gn = coef * bn * (tn - t0);

1528:           tu = x[k + 1][j][i];
1529:           au = 0.5 * (t0 + tu);
1530:           bu = PetscPowScalar(au, bm1);
1531:           /* du = bu * au; */
1532:           du = PetscPowScalar(au, beta);
1533:           gu = coef * bu * (tu - t0);

1535:           c[0].k = k;
1536:           c[0].j = j - 1;
1537:           c[0].i = i;
1538:           v[0]   = -hzhxdhy * (ds - gs);
1539:           c[1].k = k;
1540:           c[1].j = j;
1541:           c[1].i = i - 1;
1542:           v[1]   = -hyhzdhx * (dw - gw);
1543:           c[2].k = k;
1544:           c[2].j = j;
1545:           c[2].i = i;
1546:           v[2]   = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
1547:           c[3].k = k;
1548:           c[3].j = j;
1549:           c[3].i = i + 1;
1550:           v[3]   = -hyhzdhx * (de + ge);
1551:           c[4].k = k;
1552:           c[4].j = j + 1;
1553:           c[4].i = i;
1554:           v[4]   = -hzhxdhy * (dn + gn);
1555:           c[5].k = k + 1;
1556:           c[5].j = j;
1557:           c[5].i = i;
1558:           v[5]   = -hxhydhz * (du + gu);
1559:           PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));

1561:         } else if (k == mz - 1) {
1562:           /* up interior plane */

1564:           tw = x[k][j][i - 1];
1565:           aw = 0.5 * (t0 + tw);
1566:           bw = PetscPowScalar(aw, bm1);
1567:           /* dw = bw * aw */
1568:           dw = PetscPowScalar(aw, beta);
1569:           gw = coef * bw * (t0 - tw);

1571:           te = x[k][j][i + 1];
1572:           ae = 0.5 * (t0 + te);
1573:           be = PetscPowScalar(ae, bm1);
1574:           /* de = be * ae; */
1575:           de = PetscPowScalar(ae, beta);
1576:           ge = coef * be * (te - t0);

1578:           ts = x[k][j - 1][i];
1579:           as = 0.5 * (t0 + ts);
1580:           bs = PetscPowScalar(as, bm1);
1581:           /* ds = bs * as; */
1582:           ds = PetscPowScalar(as, beta);
1583:           gs = coef * bs * (t0 - ts);

1585:           tn = x[k][j + 1][i];
1586:           an = 0.5 * (t0 + tn);
1587:           bn = PetscPowScalar(an, bm1);
1588:           /* dn = bn * an; */
1589:           dn = PetscPowScalar(an, beta);
1590:           gn = coef * bn * (tn - t0);

1592:           td = x[k - 1][j][i];
1593:           ad = 0.5 * (t0 + td);
1594:           bd = PetscPowScalar(ad, bm1);
1595:           /* dd = bd * ad; */
1596:           dd = PetscPowScalar(ad, beta);
1597:           gd = coef * bd * (t0 - td);

1599:           c[0].k = k - 1;
1600:           c[0].j = j;
1601:           c[0].i = i;
1602:           v[0]   = -hxhydhz * (dd - gd);
1603:           c[1].k = k;
1604:           c[1].j = j - 1;
1605:           c[1].i = i;
1606:           v[1]   = -hzhxdhy * (ds - gs);
1607:           c[2].k = k;
1608:           c[2].j = j;
1609:           c[2].i = i - 1;
1610:           v[2]   = -hyhzdhx * (dw - gw);
1611:           c[3].k = k;
1612:           c[3].j = j;
1613:           c[3].i = i;
1614:           v[3]   = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1615:           c[4].k = k;
1616:           c[4].j = j;
1617:           c[4].i = i + 1;
1618:           v[4]   = -hyhzdhx * (de + ge);
1619:           c[5].k = k;
1620:           c[5].j = j + 1;
1621:           c[5].i = i;
1622:           v[5]   = -hzhxdhy * (dn + gn);
1623:           PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
1624:         }
1625:       }
1626:     }
1627:   }
1628:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
1629:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
1630:   if (jac != J) {
1631:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1632:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1633:   }
1634:   PetscCall(DMDAVecRestoreArray(da, localX, &x));
1635:   PetscCall(DMRestoreLocalVector(da, &localX));

1637:   PetscCall(PetscLogFlops((41.0 + 8.0 * POWFLOP) * xm * ym));
1638:   PetscFunctionReturn(PETSC_SUCCESS);
1639: }

1641: PetscErrorCode TestConvergence(SNES snes, PETSC_UNUSED PetscInt it, PETSC_UNUSED PetscReal xnorm, PETSC_UNUSED PetscReal snorm, PETSC_UNUSED PetscReal fnorm, SNESConvergedReason *reason, void *ctx)
1642: {
1643:   AppCtx *user = (AppCtx *)ctx;

1645:   PetscFunctionBeginUser;
1646:   if (user->converged) *reason = SNES_CONVERGED_USER;
1647:   else *reason = SNES_DIVERGED_USER;
1648:   PetscFunctionReturn(PETSC_SUCCESS);
1649: }

1651: /*TEST

1653:    test:
1654:       nsize: 4
1655:       args: -snes_monitor_short -pc_mg_type full -ksp_type fgmres -pc_type mg -snes_view -pc_mg_levels 2 -pc_mg_galerkin pmat
1656:       requires: !single

1658:    test:
1659:       suffix: 2
1660:       args: -snes_converged_reason -use_convergence_test -mark_converged

1662:    test:
1663:       suffix: 3
1664:       args: -snes_converged_reason -use_convergence_test -mark_converged 0

1666: TEST*/