Actual source code: ex20.c


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

 14: /*

 16:     This example models the partial differential equation

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

 20:     where beta = 2.5 and alpha = 1.0

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

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

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

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

 33: */

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

 39: /* User-defined application context */

 41: typedef struct {
 42:   PetscReal tleft, tright;   /* Dirichlet boundary conditions */
 43:   PetscReal beta, bm1, coef; /* nonlinear diffusivity parameterizations */
 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 *);

 52: int main(int argc, char **argv)
 53: {
 54:   SNES      snes;
 55:   AppCtx    user;
 56:   PetscInt  its, lits;
 57:   PetscReal litspit;
 58:   DM        da;

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

 72:   /*
 73:       Set the DMDA (grid structure) for the grids.
 74:   */
 75:   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));
 76:   PetscCall(DMSetFromOptions(da));
 77:   PetscCall(DMSetUp(da));
 78:   PetscCall(DMSetApplicationContext(da, &user));

 80:   /*
 81:      Create the nonlinear solver
 82:   */
 83:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
 84:   PetscCall(SNESSetDM(snes, da));
 85:   PetscCall(SNESSetFunction(snes, NULL, FormFunction, &user));
 86:   PetscCall(SNESSetJacobian(snes, NULL, NULL, FormJacobian, &user));
 87:   PetscCall(SNESSetFromOptions(snes));
 88:   PetscCall(SNESSetComputeInitialGuess(snes, FormInitialGuess, NULL));

 90:   PetscCall(SNESSolve(snes, NULL, NULL));
 91:   PetscCall(SNESGetIterationNumber(snes, &its));
 92:   PetscCall(SNESGetLinearSolveIterations(snes, &lits));
 93:   litspit = ((PetscReal)lits) / ((PetscReal)its);
 94:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
 95:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of Linear iterations = %" PetscInt_FMT "\n", lits));
 96:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Average Linear its / SNES = %e\n", (double)litspit));

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

111:   PetscFunctionBeginUser;
112:   PetscCall(SNESGetDM(snes, &da));
113:   PetscCall(DMGetApplicationContext(da, &user));
114:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
115:   PetscCall(DMDAVecGetArray(da, X, &x));

117:   /* Compute initial guess */
118:   for (k = zs; k < zs + zm; k++) {
119:     for (j = ys; j < ys + ym; j++) {
120:       for (i = xs; i < xs + xm; i++) x[k][j][i] = user->tleft;
121:     }
122:   }
123:   PetscCall(DMDAVecRestoreArray(da, X, &x));
124:   PetscFunctionReturn(PETSC_SUCCESS);
125: }
126: /* --------------------  Evaluate Function F(x) --------------------- */
127: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr)
128: {
129:   AppCtx        *user = (AppCtx *)ptr;
130:   PetscInt       i, j, k, mx, my, mz, xs, ys, zs, xm, ym, zm;
131:   PetscScalar    zero = 0.0, one = 1.0;
132:   PetscScalar    hx, hy, hz, hxhydhz, hyhzdhx, hzhxdhy;
133:   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;
134:   PetscScalar    tleft, tright, beta, td, ad, dd, fd = 0.0, tu, au, du = 0.0, fu = 0.0;
135:   PetscScalar ***x, ***f;
136:   Vec            localX;
137:   DM             da;

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

153:   /* Get ghost points */
154:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
155:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
156:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
157:   PetscCall(DMDAVecGetArray(da, localX, &x));
158:   PetscCall(DMDAVecGetArray(da, F, &f));

160:   /* Evaluate function */
161:   for (k = zs; k < zs + zm; k++) {
162:     for (j = ys; j < ys + ym; j++) {
163:       for (i = xs; i < xs + xm; i++) {
164:         t0 = x[k][j][i];

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

169:           tw = x[k][j][i - 1];
170:           aw = 0.5 * (t0 + tw);
171:           dw = PetscPowScalar(aw, beta);
172:           fw = dw * (t0 - tw);

174:           te = x[k][j][i + 1];
175:           ae = 0.5 * (t0 + te);
176:           de = PetscPowScalar(ae, beta);
177:           fe = de * (te - t0);

179:           ts = x[k][j - 1][i];
180:           as = 0.5 * (t0 + ts);
181:           ds = PetscPowScalar(as, beta);
182:           fs = ds * (t0 - ts);

184:           tn = x[k][j + 1][i];
185:           an = 0.5 * (t0 + tn);
186:           dn = PetscPowScalar(an, beta);
187:           fn = dn * (tn - t0);

189:           td = x[k - 1][j][i];
190:           ad = 0.5 * (t0 + td);
191:           dd = PetscPowScalar(ad, beta);
192:           fd = dd * (t0 - td);

194:           tu = x[k + 1][j][i];
195:           au = 0.5 * (t0 + tu);
196:           du = PetscPowScalar(au, beta);
197:           fu = du * (tu - t0);

199:         } else if (i == 0) {
200:           /* left-hand (west) boundary */
201:           tw = tleft;
202:           aw = 0.5 * (t0 + tw);
203:           dw = PetscPowScalar(aw, beta);
204:           fw = dw * (t0 - tw);

206:           te = x[k][j][i + 1];
207:           ae = 0.5 * (t0 + te);
208:           de = PetscPowScalar(ae, beta);
209:           fe = de * (te - t0);

211:           if (j > 0) {
212:             ts = x[k][j - 1][i];
213:             as = 0.5 * (t0 + ts);
214:             ds = PetscPowScalar(as, beta);
215:             fs = ds * (t0 - ts);
216:           } else {
217:             fs = zero;
218:           }

220:           if (j < my - 1) {
221:             tn = x[k][j + 1][i];
222:             an = 0.5 * (t0 + tn);
223:             dn = PetscPowScalar(an, beta);
224:             fn = dn * (tn - t0);
225:           } else {
226:             fn = zero;
227:           }

229:           if (k > 0) {
230:             td = x[k - 1][j][i];
231:             ad = 0.5 * (t0 + td);
232:             dd = PetscPowScalar(ad, beta);
233:             fd = dd * (t0 - td);
234:           } else {
235:             fd = zero;
236:           }

238:           if (k < mz - 1) {
239:             tu = x[k + 1][j][i];
240:             au = 0.5 * (t0 + tu);
241:             du = PetscPowScalar(au, beta);
242:             fu = du * (tu - t0);
243:           } else {
244:             fu = zero;
245:           }

247:         } else if (i == mx - 1) {
248:           /* right-hand (east) boundary */
249:           tw = x[k][j][i - 1];
250:           aw = 0.5 * (t0 + tw);
251:           dw = PetscPowScalar(aw, beta);
252:           fw = dw * (t0 - tw);

254:           te = tright;
255:           ae = 0.5 * (t0 + te);
256:           de = PetscPowScalar(ae, beta);
257:           fe = de * (te - t0);

259:           if (j > 0) {
260:             ts = x[k][j - 1][i];
261:             as = 0.5 * (t0 + ts);
262:             ds = PetscPowScalar(as, beta);
263:             fs = ds * (t0 - ts);
264:           } else {
265:             fs = zero;
266:           }

268:           if (j < my - 1) {
269:             tn = x[k][j + 1][i];
270:             an = 0.5 * (t0 + tn);
271:             dn = PetscPowScalar(an, beta);
272:             fn = dn * (tn - t0);
273:           } else {
274:             fn = zero;
275:           }

277:           if (k > 0) {
278:             td = x[k - 1][j][i];
279:             ad = 0.5 * (t0 + td);
280:             dd = PetscPowScalar(ad, beta);
281:             fd = dd * (t0 - td);
282:           } else {
283:             fd = zero;
284:           }

286:           if (k < mz - 1) {
287:             tu = x[k + 1][j][i];
288:             au = 0.5 * (t0 + tu);
289:             du = PetscPowScalar(au, beta);
290:             fu = du * (tu - t0);
291:           } else {
292:             fu = zero;
293:           }

295:         } else if (j == 0) {
296:           /* bottom (south) boundary, and i <> 0 or mx-1 */
297:           tw = x[k][j][i - 1];
298:           aw = 0.5 * (t0 + tw);
299:           dw = PetscPowScalar(aw, beta);
300:           fw = dw * (t0 - tw);

302:           te = x[k][j][i + 1];
303:           ae = 0.5 * (t0 + te);
304:           de = PetscPowScalar(ae, beta);
305:           fe = de * (te - t0);

307:           fs = zero;

309:           tn = x[k][j + 1][i];
310:           an = 0.5 * (t0 + tn);
311:           dn = PetscPowScalar(an, beta);
312:           fn = dn * (tn - t0);

314:           if (k > 0) {
315:             td = x[k - 1][j][i];
316:             ad = 0.5 * (t0 + td);
317:             dd = PetscPowScalar(ad, beta);
318:             fd = dd * (t0 - td);
319:           } else {
320:             fd = zero;
321:           }

323:           if (k < mz - 1) {
324:             tu = x[k + 1][j][i];
325:             au = 0.5 * (t0 + tu);
326:             du = PetscPowScalar(au, beta);
327:             fu = du * (tu - t0);
328:           } else {
329:             fu = zero;
330:           }

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

339:           te = x[k][j][i + 1];
340:           ae = 0.5 * (t0 + te);
341:           de = PetscPowScalar(ae, beta);
342:           fe = de * (te - t0);

344:           ts = x[k][j - 1][i];
345:           as = 0.5 * (t0 + ts);
346:           ds = PetscPowScalar(as, beta);
347:           fs = ds * (t0 - ts);

349:           fn = zero;

351:           if (k > 0) {
352:             td = x[k - 1][j][i];
353:             ad = 0.5 * (t0 + td);
354:             dd = PetscPowScalar(ad, beta);
355:             fd = dd * (t0 - td);
356:           } else {
357:             fd = zero;
358:           }

360:           if (k < mz - 1) {
361:             tu = x[k + 1][j][i];
362:             au = 0.5 * (t0 + tu);
363:             du = PetscPowScalar(au, beta);
364:             fu = du * (tu - t0);
365:           } else {
366:             fu = zero;
367:           }

369:         } else if (k == 0) {
370:           /* down boundary (interior only) */
371:           tw = x[k][j][i - 1];
372:           aw = 0.5 * (t0 + tw);
373:           dw = PetscPowScalar(aw, beta);
374:           fw = dw * (t0 - tw);

376:           te = x[k][j][i + 1];
377:           ae = 0.5 * (t0 + te);
378:           de = PetscPowScalar(ae, beta);
379:           fe = de * (te - t0);

381:           ts = x[k][j - 1][i];
382:           as = 0.5 * (t0 + ts);
383:           ds = PetscPowScalar(as, beta);
384:           fs = ds * (t0 - ts);

386:           tn = x[k][j + 1][i];
387:           an = 0.5 * (t0 + tn);
388:           dn = PetscPowScalar(an, beta);
389:           fn = dn * (tn - t0);

391:           fd = zero;

393:           tu = x[k + 1][j][i];
394:           au = 0.5 * (t0 + tu);
395:           du = PetscPowScalar(au, beta);
396:           fu = du * (tu - t0);

398:         } else if (k == mz - 1) {
399:           /* up boundary (interior only) */
400:           tw = x[k][j][i - 1];
401:           aw = 0.5 * (t0 + tw);
402:           dw = PetscPowScalar(aw, beta);
403:           fw = dw * (t0 - tw);

405:           te = x[k][j][i + 1];
406:           ae = 0.5 * (t0 + te);
407:           de = PetscPowScalar(ae, beta);
408:           fe = de * (te - t0);

410:           ts = x[k][j - 1][i];
411:           as = 0.5 * (t0 + ts);
412:           ds = PetscPowScalar(as, beta);
413:           fs = ds * (t0 - ts);

415:           tn = x[k][j + 1][i];
416:           an = 0.5 * (t0 + tn);
417:           dn = PetscPowScalar(an, beta);
418:           fn = dn * (tn - t0);

420:           td = x[k - 1][j][i];
421:           ad = 0.5 * (t0 + td);
422:           dd = PetscPowScalar(ad, beta);
423:           fd = dd * (t0 - td);

425:           fu = zero;
426:         }

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

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

468:   /* Get ghost points */
469:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
470:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
471:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
472:   PetscCall(DMDAVecGetArray(da, localX, &x));

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

485:           tw = x[k][j][i - 1];
486:           aw = 0.5 * (t0 + tw);
487:           bw = PetscPowScalar(aw, bm1);
488:           /* dw = bw * aw */
489:           dw = PetscPowScalar(aw, beta);
490:           gw = coef * bw * (t0 - tw);

492:           te = x[k][j][i + 1];
493:           ae = 0.5 * (t0 + te);
494:           be = PetscPowScalar(ae, bm1);
495:           /* de = be * ae; */
496:           de = PetscPowScalar(ae, beta);
497:           ge = coef * be * (te - t0);

499:           ts = x[k][j - 1][i];
500:           as = 0.5 * (t0 + ts);
501:           bs = PetscPowScalar(as, bm1);
502:           /* ds = bs * as; */
503:           ds = PetscPowScalar(as, beta);
504:           gs = coef * bs * (t0 - ts);

506:           tn = x[k][j + 1][i];
507:           an = 0.5 * (t0 + tn);
508:           bn = PetscPowScalar(an, bm1);
509:           /* dn = bn * an; */
510:           dn = PetscPowScalar(an, beta);
511:           gn = coef * bn * (tn - t0);

513:           td = x[k - 1][j][i];
514:           ad = 0.5 * (t0 + td);
515:           bd = PetscPowScalar(ad, bm1);
516:           /* dd = bd * ad; */
517:           dd = PetscPowScalar(ad, beta);
518:           gd = coef * bd * (t0 - td);

520:           tu = x[k + 1][j][i];
521:           au = 0.5 * (t0 + tu);
522:           bu = PetscPowScalar(au, bm1);
523:           /* du = bu * au; */
524:           du = PetscPowScalar(au, beta);
525:           gu = coef * bu * (tu - t0);

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

557:         } else if (i == 0) {
558:           /* left-hand plane boundary */
559:           tw = tleft;
560:           aw = 0.5 * (t0 + tw);
561:           bw = PetscPowScalar(aw, bm1);
562:           /* dw = bw * aw */
563:           dw = PetscPowScalar(aw, beta);
564:           gw = coef * bw * (t0 - tw);

566:           te = x[k][j][i + 1];
567:           ae = 0.5 * (t0 + te);
568:           be = PetscPowScalar(ae, bm1);
569:           /* de = be * ae; */
570:           de = PetscPowScalar(ae, beta);
571:           ge = coef * be * (te - t0);

573:           /* left-hand bottom edge */
574:           if (j == 0) {
575:             tn = x[k][j + 1][i];
576:             an = 0.5 * (t0 + tn);
577:             bn = PetscPowScalar(an, bm1);
578:             /* dn = bn * an; */
579:             dn = PetscPowScalar(an, beta);
580:             gn = coef * bn * (tn - t0);

582:             /* left-hand bottom down corner */
583:             if (k == 0) {
584:               tu = x[k + 1][j][i];
585:               au = 0.5 * (t0 + tu);
586:               bu = PetscPowScalar(au, bm1);
587:               /* du = bu * au; */
588:               du = PetscPowScalar(au, beta);
589:               gu = coef * bu * (tu - t0);

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

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

618:               td = x[k - 1][j][i];
619:               ad = 0.5 * (t0 + td);
620:               bd = PetscPowScalar(ad, bm1);
621:               /* dd = bd * ad; */
622:               dd = PetscPowScalar(ad, beta);
623:               gd = coef * bd * (td - t0);

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

647:               /* left-hand bottom up corner */
648:             } else {
649:               td = x[k - 1][j][i];
650:               ad = 0.5 * (t0 + td);
651:               bd = PetscPowScalar(ad, bm1);
652:               /* dd = bd * ad; */
653:               dd = PetscPowScalar(ad, beta);
654:               gd = coef * bd * (td - t0);

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

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

684:             /* left-hand top down corner */
685:             if (k == 0) {
686:               tu = x[k + 1][j][i];
687:               au = 0.5 * (t0 + tu);
688:               bu = PetscPowScalar(au, bm1);
689:               /* du = bu * au; */
690:               du = PetscPowScalar(au, beta);
691:               gu = coef * bu * (tu - t0);

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

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

720:               td = x[k - 1][j][i];
721:               ad = 0.5 * (t0 + td);
722:               bd = PetscPowScalar(ad, bm1);
723:               /* dd = bd * ad; */
724:               dd = PetscPowScalar(ad, beta);
725:               gd = coef * bd * (td - t0);

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

749:               /* left-hand top up corner */
750:             } else {
751:               td = x[k - 1][j][i];
752:               ad = 0.5 * (t0 + td);
753:               bd = PetscPowScalar(ad, bm1);
754:               /* dd = bd * ad; */
755:               dd = PetscPowScalar(ad, beta);
756:               gd = coef * bd * (td - t0);

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

777:           } else {
778:             ts = x[k][j - 1][i];
779:             as = 0.5 * (t0 + ts);
780:             bs = PetscPowScalar(as, bm1);
781:             /* ds = bs * as; */
782:             ds = PetscPowScalar(as, beta);
783:             gs = coef * bs * (t0 - ts);

785:             tn = x[k][j + 1][i];
786:             an = 0.5 * (t0 + tn);
787:             bn = PetscPowScalar(an, bm1);
788:             /* dn = bn * an; */
789:             dn = PetscPowScalar(an, beta);
790:             gn = coef * bn * (tn - t0);

792:             /* left-hand down interior edge */
793:             if (k == 0) {
794:               tu = x[k + 1][j][i];
795:               au = 0.5 * (t0 + tu);
796:               bu = PetscPowScalar(au, bm1);
797:               /* du = bu * au; */
798:               du = PetscPowScalar(au, beta);
799:               gu = coef * bu * (tu - t0);

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

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

825:               td = x[k - 1][j][i];
826:               ad = 0.5 * (t0 + td);
827:               bd = PetscPowScalar(ad, bm1);
828:               /* dd = bd * ad; */
829:               dd = PetscPowScalar(ad, beta);
830:               gd = coef * bd * (t0 - td);

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

855:               td = x[k - 1][j][i];
856:               ad = 0.5 * (t0 + td);
857:               bd = PetscPowScalar(ad, bm1);
858:               /* dd = bd * ad; */
859:               dd = PetscPowScalar(ad, beta);
860:               gd = coef * bd * (t0 - td);

862:               tu = x[k + 1][j][i];
863:               au = 0.5 * (t0 + tu);
864:               bu = PetscPowScalar(au, bm1);
865:               /* du = bu * au; */
866:               du = PetscPowScalar(au, beta);
867:               gu = coef * bu * (tu - t0);

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

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

906:           te = tright;
907:           ae = 0.5 * (t0 + te);
908:           be = PetscPowScalar(ae, bm1);
909:           /* de = be * ae; */
910:           de = PetscPowScalar(ae, beta);
911:           ge = coef * be * (te - t0);

913:           /* right-hand bottom edge */
914:           if (j == 0) {
915:             tn = x[k][j + 1][i];
916:             an = 0.5 * (t0 + tn);
917:             bn = PetscPowScalar(an, bm1);
918:             /* dn = bn * an; */
919:             dn = PetscPowScalar(an, beta);
920:             gn = coef * bn * (tn - t0);

922:             /* right-hand bottom down corner */
923:             if (k == 0) {
924:               tu = x[k + 1][j][i];
925:               au = 0.5 * (t0 + tu);
926:               bu = PetscPowScalar(au, bm1);
927:               /* du = bu * au; */
928:               du = PetscPowScalar(au, beta);
929:               gu = coef * bu * (tu - t0);

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

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

958:               td = x[k - 1][j][i];
959:               ad = 0.5 * (t0 + td);
960:               bd = PetscPowScalar(ad, bm1);
961:               /* dd = bd * ad; */
962:               dd = PetscPowScalar(ad, beta);
963:               gd = coef * bd * (td - t0);

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

987:               /* right-hand bottom up corner */
988:             } else {
989:               td = x[k - 1][j][i];
990:               ad = 0.5 * (t0 + td);
991:               bd = PetscPowScalar(ad, bm1);
992:               /* dd = bd * ad; */
993:               dd = PetscPowScalar(ad, beta);
994:               gd = coef * bd * (td - t0);

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

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

1024:             /* right-hand top down corner */
1025:             if (k == 0) {
1026:               tu = x[k + 1][j][i];
1027:               au = 0.5 * (t0 + tu);
1028:               bu = PetscPowScalar(au, bm1);
1029:               /* du = bu * au; */
1030:               du = PetscPowScalar(au, beta);
1031:               gu = coef * bu * (tu - t0);

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

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

1060:               td = x[k - 1][j][i];
1061:               ad = 0.5 * (t0 + td);
1062:               bd = PetscPowScalar(ad, bm1);
1063:               /* dd = bd * ad; */
1064:               dd = PetscPowScalar(ad, beta);
1065:               gd = coef * bd * (td - t0);

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

1089:               /* right-hand top up corner */
1090:             } else {
1091:               td = x[k - 1][j][i];
1092:               ad = 0.5 * (t0 + td);
1093:               bd = PetscPowScalar(ad, bm1);
1094:               /* dd = bd * ad; */
1095:               dd = PetscPowScalar(ad, beta);
1096:               gd = coef * bd * (td - t0);

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

1117:           } else {
1118:             ts = x[k][j - 1][i];
1119:             as = 0.5 * (t0 + ts);
1120:             bs = PetscPowScalar(as, bm1);
1121:             /* ds = bs * as; */
1122:             ds = PetscPowScalar(as, beta);
1123:             gs = coef * bs * (t0 - ts);

1125:             tn = x[k][j + 1][i];
1126:             an = 0.5 * (t0 + tn);
1127:             bn = PetscPowScalar(an, bm1);
1128:             /* dn = bn * an; */
1129:             dn = PetscPowScalar(an, beta);
1130:             gn = coef * bn * (tn - t0);

1132:             /* right-hand down interior edge */
1133:             if (k == 0) {
1134:               tu = x[k + 1][j][i];
1135:               au = 0.5 * (t0 + tu);
1136:               bu = PetscPowScalar(au, bm1);
1137:               /* du = bu * au; */
1138:               du = PetscPowScalar(au, beta);
1139:               gu = coef * bu * (tu - t0);

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

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

1165:               td = x[k - 1][j][i];
1166:               ad = 0.5 * (t0 + td);
1167:               bd = PetscPowScalar(ad, bm1);
1168:               /* dd = bd * ad; */
1169:               dd = PetscPowScalar(ad, beta);
1170:               gd = coef * bd * (t0 - td);

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

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

1196:               td = x[k - 1][j][i];
1197:               ad = 0.5 * (t0 + td);
1198:               bd = PetscPowScalar(ad, bm1);
1199:               /* dd = bd * ad; */
1200:               dd = PetscPowScalar(ad, beta);
1201:               gd = coef * bd * (t0 - td);

1203:               tu = x[k + 1][j][i];
1204:               au = 0.5 * (t0 + tu);
1205:               bu = PetscPowScalar(au, bm1);
1206:               /* du = bu * au; */
1207:               du = PetscPowScalar(au, beta);
1208:               gu = coef * bu * (tu - t0);

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

1238:         } else if (j == 0) {
1239:           tw = x[k][j][i - 1];
1240:           aw = 0.5 * (t0 + tw);
1241:           bw = PetscPowScalar(aw, bm1);
1242:           /* dw = bw * aw */
1243:           dw = PetscPowScalar(aw, beta);
1244:           gw = coef * bw * (t0 - tw);

1246:           te = x[k][j][i + 1];
1247:           ae = 0.5 * (t0 + te);
1248:           be = PetscPowScalar(ae, bm1);
1249:           /* de = be * ae; */
1250:           de = PetscPowScalar(ae, beta);
1251:           ge = coef * be * (te - t0);

1253:           tn = x[k][j + 1][i];
1254:           an = 0.5 * (t0 + tn);
1255:           bn = PetscPowScalar(an, bm1);
1256:           /* dn = bn * an; */
1257:           dn = PetscPowScalar(an, beta);
1258:           gn = coef * bn * (tn - t0);

1260:           /* bottom down interior edge */
1261:           if (k == 0) {
1262:             tu = x[k + 1][j][i];
1263:             au = 0.5 * (t0 + tu);
1264:             bu = PetscPowScalar(au, bm1);
1265:             /* du = bu * au; */
1266:             du = PetscPowScalar(au, beta);
1267:             gu = coef * bu * (tu - t0);

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

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

1293:             td = x[k - 1][j][i];
1294:             ad = 0.5 * (t0 + td);
1295:             bd = PetscPowScalar(ad, bm1);
1296:             /* dd = bd * ad; */
1297:             dd = PetscPowScalar(ad, beta);
1298:             gd = coef * bd * (td - t0);

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

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

1324:             tu = x[k + 1][j][i];
1325:             au = 0.5 * (t0 + tu);
1326:             bu = PetscPowScalar(au, bm1);
1327:             /* du = bu * au; */
1328:             du = PetscPowScalar(au, beta);
1329:             gu = coef * bu * (tu - t0);

1331:             td = x[k - 1][j][i];
1332:             ad = 0.5 * (t0 + td);
1333:             bd = PetscPowScalar(ad, bm1);
1334:             /* dd = bd * ad; */
1335:             dd = PetscPowScalar(ad, beta);
1336:             gd = coef * bd * (td - t0);

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

1365:         } else if (j == my - 1) {
1366:           tw = x[k][j][i - 1];
1367:           aw = 0.5 * (t0 + tw);
1368:           bw = PetscPowScalar(aw, bm1);
1369:           /* dw = bw * aw */
1370:           dw = PetscPowScalar(aw, beta);
1371:           gw = coef * bw * (t0 - tw);

1373:           te = x[k][j][i + 1];
1374:           ae = 0.5 * (t0 + te);
1375:           be = PetscPowScalar(ae, bm1);
1376:           /* de = be * ae; */
1377:           de = PetscPowScalar(ae, beta);
1378:           ge = coef * be * (te - t0);

1380:           ts = x[k][j - 1][i];
1381:           as = 0.5 * (t0 + ts);
1382:           bs = PetscPowScalar(as, bm1);
1383:           /* ds = bs * as; */
1384:           ds = PetscPowScalar(as, beta);
1385:           gs = coef * bs * (t0 - ts);

1387:           /* top down interior edge */
1388:           if (k == 0) {
1389:             tu = x[k + 1][j][i];
1390:             au = 0.5 * (t0 + tu);
1391:             bu = PetscPowScalar(au, bm1);
1392:             /* du = bu * au; */
1393:             du = PetscPowScalar(au, beta);
1394:             gu = coef * bu * (tu - t0);

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

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

1420:             td = x[k - 1][j][i];
1421:             ad = 0.5 * (t0 + td);
1422:             bd = PetscPowScalar(ad, bm1);
1423:             /* dd = bd * ad; */
1424:             dd = PetscPowScalar(ad, beta);
1425:             gd = coef * bd * (td - t0);

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

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

1451:             tu = x[k + 1][j][i];
1452:             au = 0.5 * (t0 + tu);
1453:             bu = PetscPowScalar(au, bm1);
1454:             /* du = bu * au; */
1455:             du = PetscPowScalar(au, beta);
1456:             gu = coef * bu * (tu - t0);

1458:             td = x[k - 1][j][i];
1459:             ad = 0.5 * (t0 + td);
1460:             bd = PetscPowScalar(ad, bm1);
1461:             /* dd = bd * ad; */
1462:             dd = PetscPowScalar(ad, beta);
1463:             gd = coef * bd * (td - t0);

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

1492:         } else if (k == 0) {
1493:           /* down interior plane */

1495:           tw = x[k][j][i - 1];
1496:           aw = 0.5 * (t0 + tw);
1497:           bw = PetscPowScalar(aw, bm1);
1498:           /* dw = bw * aw */
1499:           dw = PetscPowScalar(aw, beta);
1500:           gw = coef * bw * (t0 - tw);

1502:           te = x[k][j][i + 1];
1503:           ae = 0.5 * (t0 + te);
1504:           be = PetscPowScalar(ae, bm1);
1505:           /* de = be * ae; */
1506:           de = PetscPowScalar(ae, beta);
1507:           ge = coef * be * (te - t0);

1509:           ts = x[k][j - 1][i];
1510:           as = 0.5 * (t0 + ts);
1511:           bs = PetscPowScalar(as, bm1);
1512:           /* ds = bs * as; */
1513:           ds = PetscPowScalar(as, beta);
1514:           gs = coef * bs * (t0 - ts);

1516:           tn = x[k][j + 1][i];
1517:           an = 0.5 * (t0 + tn);
1518:           bn = PetscPowScalar(an, bm1);
1519:           /* dn = bn * an; */
1520:           dn = PetscPowScalar(an, beta);
1521:           gn = coef * bn * (tn - t0);

1523:           tu = x[k + 1][j][i];
1524:           au = 0.5 * (t0 + tu);
1525:           bu = PetscPowScalar(au, bm1);
1526:           /* du = bu * au; */
1527:           du = PetscPowScalar(au, beta);
1528:           gu = coef * bu * (tu - t0);

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

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

1559:           tw = x[k][j][i - 1];
1560:           aw = 0.5 * (t0 + tw);
1561:           bw = PetscPowScalar(aw, bm1);
1562:           /* dw = bw * aw */
1563:           dw = PetscPowScalar(aw, beta);
1564:           gw = coef * bw * (t0 - tw);

1566:           te = x[k][j][i + 1];
1567:           ae = 0.5 * (t0 + te);
1568:           be = PetscPowScalar(ae, bm1);
1569:           /* de = be * ae; */
1570:           de = PetscPowScalar(ae, beta);
1571:           ge = coef * be * (te - t0);

1573:           ts = x[k][j - 1][i];
1574:           as = 0.5 * (t0 + ts);
1575:           bs = PetscPowScalar(as, bm1);
1576:           /* ds = bs * as; */
1577:           ds = PetscPowScalar(as, beta);
1578:           gs = coef * bs * (t0 - ts);

1580:           tn = x[k][j + 1][i];
1581:           an = 0.5 * (t0 + tn);
1582:           bn = PetscPowScalar(an, bm1);
1583:           /* dn = bn * an; */
1584:           dn = PetscPowScalar(an, beta);
1585:           gn = coef * bn * (tn - t0);

1587:           td = x[k - 1][j][i];
1588:           ad = 0.5 * (t0 + td);
1589:           bd = PetscPowScalar(ad, bm1);
1590:           /* dd = bd * ad; */
1591:           dd = PetscPowScalar(ad, beta);
1592:           gd = coef * bd * (t0 - td);

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

1632:   PetscCall(PetscLogFlops((41.0 + 8.0 * POWFLOP) * xm * ym));
1633:   PetscFunctionReturn(PETSC_SUCCESS);
1634: }

1636: /*TEST

1638:    test:
1639:       nsize: 4
1640:       args: -snes_monitor_short -pc_mg_type full -ksp_type fgmres -pc_type mg -snes_view -pc_mg_levels 2 -pc_mg_galerkin pmat
1641:       requires: !single

1643: TEST*/