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: } AppCtx;

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

 47: extern PetscErrorCode FormInitialGuess(SNES, Vec, void *);
 48: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
 49: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

306:           fs = zero;

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

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

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

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

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

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

348:           fn = zero;

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

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

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

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

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

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

390:           fd = zero;

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

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

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

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

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

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

424:           fu = zero;
425:         }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1635: /*TEST

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

1642: TEST*/