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*/