Actual source code: ex20.c
1: static char help[] = "Nonlinear Radiative Transport PDE with multigrid in 3d.\n\
2: Uses 3-dimensional distributed arrays.\n\
3: A 3-dim simplified Radiative Transport test problem is used, with analytic Jacobian. \n\
4: \n\
5: Solves the linear systems via multilevel methods \n\
6: \n\
7: The command line\n\
8: options are:\n\
9: -tleft <tl>, where <tl> indicates the left Diriclet BC \n\
10: -tright <tr>, where <tr> indicates the right Diriclet BC \n\
11: -beta <beta>, where <beta> indicates the exponent in T \n\n";
13: /*
15: This example models the partial differential equation
17: - Div(alpha* T^beta (GRAD T)) = 0.
19: where beta = 2.5 and alpha = 1.0
21: BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = dT/dn_up = dT/dn_down = 0.
23: in the unit square, which is uniformly discretized in each of x and
24: y in this simple encoding. The degrees of freedom are cell centered.
26: A finite volume approximation with the usual 7-point stencil
27: is used to discretize the boundary value problem to obtain a
28: nonlinear system of equations.
30: This code was contributed by Nickolas Jovanovic based on ex18.c
32: */
34: #include <petscsnes.h>
35: #include <petscdm.h>
36: #include <petscdmda.h>
38: /* User-defined application context */
40: typedef struct {
41: PetscReal tleft, tright; /* Dirichlet boundary conditions */
42: PetscReal beta, bm1, coef; /* nonlinear diffusivity parameterizations */
43: PetscBool converged; /* For user convergence test, whether to mark as converged */
44: } AppCtx;
46: #define POWFLOP 5 /* assume a pow() takes five flops */
48: extern PetscErrorCode FormInitialGuess(SNES, Vec, void *);
49: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
50: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
51: extern PetscErrorCode TestConvergence(SNES, PetscInt, PetscReal, PetscReal, PetscReal, SNESConvergedReason *, void *);
53: int main(int argc, char **argv)
54: {
55: SNES snes;
56: AppCtx user;
57: PetscInt its, lits;
58: PetscReal litspit = 0; /* avoid uninitialized warning */
59: DM da;
60: PetscBool use_convergence_test = PETSC_FALSE;
62: PetscFunctionBeginUser;
63: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
64: /* set problem parameters */
65: user.tleft = 1.0;
66: user.tright = 0.1;
67: user.beta = 2.5;
68: PetscCall(PetscOptionsGetReal(NULL, NULL, "-tleft", &user.tleft, NULL));
69: PetscCall(PetscOptionsGetReal(NULL, NULL, "-tright", &user.tright, NULL));
70: PetscCall(PetscOptionsGetReal(NULL, NULL, "-beta", &user.beta, NULL));
71: PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_convergence_test", &use_convergence_test, NULL));
72: PetscCall(PetscOptionsGetBool(NULL, NULL, "-mark_converged", &user.converged, NULL));
73: user.bm1 = user.beta - 1.0;
74: user.coef = user.beta / 2.0;
76: /*
77: Set the DMDA (grid structure) for the grids.
78: */
79: PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 5, 5, 5, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, 0, &da));
80: PetscCall(DMSetFromOptions(da));
81: PetscCall(DMSetUp(da));
82: PetscCall(DMSetApplicationContext(da, &user));
84: /*
85: Create the nonlinear solver
86: */
87: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
88: PetscCall(SNESSetDM(snes, da));
89: PetscCall(SNESSetFunction(snes, NULL, FormFunction, &user));
90: PetscCall(SNESSetJacobian(snes, NULL, NULL, FormJacobian, &user));
91: if (use_convergence_test) PetscCall(SNESSetConvergenceTest(snes, TestConvergence, &user, NULL));
92: PetscCall(SNESSetFromOptions(snes));
93: PetscCall(SNESSetComputeInitialGuess(snes, FormInitialGuess, NULL));
95: PetscCall(SNESSolve(snes, NULL, NULL));
96: PetscCall(SNESGetIterationNumber(snes, &its));
97: PetscCall(SNESGetLinearSolveIterations(snes, &lits));
98: if (its) litspit = ((PetscReal)lits) / ((PetscReal)its);
99: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
100: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of Linear iterations = %" PetscInt_FMT "\n", lits));
101: if (its) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Average Linear its / SNES = %e\n", (double)litspit));
103: PetscCall(SNESDestroy(&snes));
104: PetscCall(DMDestroy(&da));
105: PetscCall(PetscFinalize());
106: return 0;
107: }
108: /* -------------------- Form initial approximation ----------------- */
109: PetscErrorCode FormInitialGuess(SNES snes, Vec X, void *ctx)
110: {
111: AppCtx *user;
112: PetscInt i, j, k, xs, ys, xm, ym, zs, zm;
113: PetscScalar ***x;
114: DM da;
116: PetscFunctionBeginUser;
117: PetscCall(SNESGetDM(snes, &da));
118: PetscCall(DMGetApplicationContext(da, &user));
119: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
120: PetscCall(DMDAVecGetArray(da, X, &x));
122: /* Compute initial guess */
123: for (k = zs; k < zs + zm; k++) {
124: for (j = ys; j < ys + ym; j++) {
125: for (i = xs; i < xs + xm; i++) x[k][j][i] = user->tleft;
126: }
127: }
128: PetscCall(DMDAVecRestoreArray(da, X, &x));
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
131: /* -------------------- Evaluate Function F(x) --------------------- */
132: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr)
133: {
134: AppCtx *user = (AppCtx *)ptr;
135: PetscInt i, j, k, mx, my, mz, xs, ys, zs, xm, ym, zm;
136: PetscScalar zero = 0.0, one = 1.0;
137: PetscScalar hx, hy, hz, hxhydhz, hyhzdhx, hzhxdhy;
138: PetscScalar t0, tn, ts, te, tw, an, as, ae, aw, dn, ds, de, dw, fn = 0.0, fs = 0.0, fe = 0.0, fw = 0.0;
139: PetscScalar tleft, tright, beta, td, ad, dd, fd = 0.0, tu, au, du = 0.0, fu = 0.0;
140: PetscScalar ***x, ***f;
141: Vec localX;
142: DM da;
144: PetscFunctionBeginUser;
145: PetscCall(SNESGetDM(snes, &da));
146: PetscCall(DMGetLocalVector(da, &localX));
147: PetscCall(DMDAGetInfo(da, NULL, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
148: hx = one / (PetscReal)(mx - 1);
149: hy = one / (PetscReal)(my - 1);
150: hz = one / (PetscReal)(mz - 1);
151: hxhydhz = hx * hy / hz;
152: hyhzdhx = hy * hz / hx;
153: hzhxdhy = hz * hx / hy;
154: tleft = user->tleft;
155: tright = user->tright;
156: beta = user->beta;
158: /* Get ghost points */
159: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
160: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
161: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
162: PetscCall(DMDAVecGetArray(da, localX, &x));
163: PetscCall(DMDAVecGetArray(da, F, &f));
165: /* Evaluate function */
166: for (k = zs; k < zs + zm; k++) {
167: for (j = ys; j < ys + ym; j++) {
168: for (i = xs; i < xs + xm; i++) {
169: t0 = x[k][j][i];
171: if (i > 0 && i < mx - 1 && j > 0 && j < my - 1 && k > 0 && k < mz - 1) {
172: /* general interior volume */
174: tw = x[k][j][i - 1];
175: aw = 0.5 * (t0 + tw);
176: dw = PetscPowScalar(aw, beta);
177: fw = dw * (t0 - tw);
179: te = x[k][j][i + 1];
180: ae = 0.5 * (t0 + te);
181: de = PetscPowScalar(ae, beta);
182: fe = de * (te - t0);
184: ts = x[k][j - 1][i];
185: as = 0.5 * (t0 + ts);
186: ds = PetscPowScalar(as, beta);
187: fs = ds * (t0 - ts);
189: tn = x[k][j + 1][i];
190: an = 0.5 * (t0 + tn);
191: dn = PetscPowScalar(an, beta);
192: fn = dn * (tn - t0);
194: td = x[k - 1][j][i];
195: ad = 0.5 * (t0 + td);
196: dd = PetscPowScalar(ad, beta);
197: fd = dd * (t0 - td);
199: tu = x[k + 1][j][i];
200: au = 0.5 * (t0 + tu);
201: du = PetscPowScalar(au, beta);
202: fu = du * (tu - t0);
204: } else if (i == 0) {
205: /* left-hand (west) boundary */
206: tw = tleft;
207: aw = 0.5 * (t0 + tw);
208: dw = PetscPowScalar(aw, beta);
209: fw = dw * (t0 - tw);
211: te = x[k][j][i + 1];
212: ae = 0.5 * (t0 + te);
213: de = PetscPowScalar(ae, beta);
214: fe = de * (te - t0);
216: if (j > 0) {
217: ts = x[k][j - 1][i];
218: as = 0.5 * (t0 + ts);
219: ds = PetscPowScalar(as, beta);
220: fs = ds * (t0 - ts);
221: } else {
222: fs = zero;
223: }
225: if (j < my - 1) {
226: tn = x[k][j + 1][i];
227: an = 0.5 * (t0 + tn);
228: dn = PetscPowScalar(an, beta);
229: fn = dn * (tn - t0);
230: } else {
231: fn = zero;
232: }
234: if (k > 0) {
235: td = x[k - 1][j][i];
236: ad = 0.5 * (t0 + td);
237: dd = PetscPowScalar(ad, beta);
238: fd = dd * (t0 - td);
239: } else {
240: fd = zero;
241: }
243: if (k < mz - 1) {
244: tu = x[k + 1][j][i];
245: au = 0.5 * (t0 + tu);
246: du = PetscPowScalar(au, beta);
247: fu = du * (tu - t0);
248: } else {
249: fu = zero;
250: }
252: } else if (i == mx - 1) {
253: /* right-hand (east) boundary */
254: tw = x[k][j][i - 1];
255: aw = 0.5 * (t0 + tw);
256: dw = PetscPowScalar(aw, beta);
257: fw = dw * (t0 - tw);
259: te = tright;
260: ae = 0.5 * (t0 + te);
261: de = PetscPowScalar(ae, beta);
262: fe = de * (te - t0);
264: if (j > 0) {
265: ts = x[k][j - 1][i];
266: as = 0.5 * (t0 + ts);
267: ds = PetscPowScalar(as, beta);
268: fs = ds * (t0 - ts);
269: } else {
270: fs = zero;
271: }
273: if (j < my - 1) {
274: tn = x[k][j + 1][i];
275: an = 0.5 * (t0 + tn);
276: dn = PetscPowScalar(an, beta);
277: fn = dn * (tn - t0);
278: } else {
279: fn = zero;
280: }
282: if (k > 0) {
283: td = x[k - 1][j][i];
284: ad = 0.5 * (t0 + td);
285: dd = PetscPowScalar(ad, beta);
286: fd = dd * (t0 - td);
287: } else {
288: fd = zero;
289: }
291: if (k < mz - 1) {
292: tu = x[k + 1][j][i];
293: au = 0.5 * (t0 + tu);
294: du = PetscPowScalar(au, beta);
295: fu = du * (tu - t0);
296: } else {
297: fu = zero;
298: }
300: } else if (j == 0) {
301: /* bottom (south) boundary, and i <> 0 or mx-1 */
302: tw = x[k][j][i - 1];
303: aw = 0.5 * (t0 + tw);
304: dw = PetscPowScalar(aw, beta);
305: fw = dw * (t0 - tw);
307: te = x[k][j][i + 1];
308: ae = 0.5 * (t0 + te);
309: de = PetscPowScalar(ae, beta);
310: fe = de * (te - t0);
312: fs = zero;
314: tn = x[k][j + 1][i];
315: an = 0.5 * (t0 + tn);
316: dn = PetscPowScalar(an, beta);
317: fn = dn * (tn - t0);
319: if (k > 0) {
320: td = x[k - 1][j][i];
321: ad = 0.5 * (t0 + td);
322: dd = PetscPowScalar(ad, beta);
323: fd = dd * (t0 - td);
324: } else {
325: fd = zero;
326: }
328: if (k < mz - 1) {
329: tu = x[k + 1][j][i];
330: au = 0.5 * (t0 + tu);
331: du = PetscPowScalar(au, beta);
332: fu = du * (tu - t0);
333: } else {
334: fu = zero;
335: }
337: } else if (j == my - 1) {
338: /* top (north) boundary, and i <> 0 or mx-1 */
339: tw = x[k][j][i - 1];
340: aw = 0.5 * (t0 + tw);
341: dw = PetscPowScalar(aw, beta);
342: fw = dw * (t0 - tw);
344: te = x[k][j][i + 1];
345: ae = 0.5 * (t0 + te);
346: de = PetscPowScalar(ae, beta);
347: fe = de * (te - t0);
349: ts = x[k][j - 1][i];
350: as = 0.5 * (t0 + ts);
351: ds = PetscPowScalar(as, beta);
352: fs = ds * (t0 - ts);
354: fn = zero;
356: if (k > 0) {
357: td = x[k - 1][j][i];
358: ad = 0.5 * (t0 + td);
359: dd = PetscPowScalar(ad, beta);
360: fd = dd * (t0 - td);
361: } else {
362: fd = zero;
363: }
365: if (k < mz - 1) {
366: tu = x[k + 1][j][i];
367: au = 0.5 * (t0 + tu);
368: du = PetscPowScalar(au, beta);
369: fu = du * (tu - t0);
370: } else {
371: fu = zero;
372: }
374: } else if (k == 0) {
375: /* down boundary (interior only) */
376: tw = x[k][j][i - 1];
377: aw = 0.5 * (t0 + tw);
378: dw = PetscPowScalar(aw, beta);
379: fw = dw * (t0 - tw);
381: te = x[k][j][i + 1];
382: ae = 0.5 * (t0 + te);
383: de = PetscPowScalar(ae, beta);
384: fe = de * (te - t0);
386: ts = x[k][j - 1][i];
387: as = 0.5 * (t0 + ts);
388: ds = PetscPowScalar(as, beta);
389: fs = ds * (t0 - ts);
391: tn = x[k][j + 1][i];
392: an = 0.5 * (t0 + tn);
393: dn = PetscPowScalar(an, beta);
394: fn = dn * (tn - t0);
396: fd = zero;
398: tu = x[k + 1][j][i];
399: au = 0.5 * (t0 + tu);
400: du = PetscPowScalar(au, beta);
401: fu = du * (tu - t0);
403: } else if (k == mz - 1) {
404: /* up boundary (interior only) */
405: tw = x[k][j][i - 1];
406: aw = 0.5 * (t0 + tw);
407: dw = PetscPowScalar(aw, beta);
408: fw = dw * (t0 - tw);
410: te = x[k][j][i + 1];
411: ae = 0.5 * (t0 + te);
412: de = PetscPowScalar(ae, beta);
413: fe = de * (te - t0);
415: ts = x[k][j - 1][i];
416: as = 0.5 * (t0 + ts);
417: ds = PetscPowScalar(as, beta);
418: fs = ds * (t0 - ts);
420: tn = x[k][j + 1][i];
421: an = 0.5 * (t0 + tn);
422: dn = PetscPowScalar(an, beta);
423: fn = dn * (tn - t0);
425: td = x[k - 1][j][i];
426: ad = 0.5 * (t0 + td);
427: dd = PetscPowScalar(ad, beta);
428: fd = dd * (t0 - td);
430: fu = zero;
431: }
433: f[k][j][i] = -hyhzdhx * (fe - fw) - hzhxdhy * (fn - fs) - hxhydhz * (fu - fd);
434: }
435: }
436: }
437: PetscCall(DMDAVecRestoreArray(da, localX, &x));
438: PetscCall(DMDAVecRestoreArray(da, F, &f));
439: PetscCall(DMRestoreLocalVector(da, &localX));
440: PetscCall(PetscLogFlops((22.0 + 4.0 * POWFLOP) * ym * xm));
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
443: /* -------------------- Evaluate Jacobian F(x) --------------------- */
444: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr)
445: {
446: AppCtx *user = (AppCtx *)ptr;
447: PetscInt i, j, k, mx, my, mz, xs, ys, zs, xm, ym, zm;
448: PetscScalar one = 1.0;
449: PetscScalar hx, hy, hz, hxhydhz, hyhzdhx, hzhxdhy;
450: PetscScalar t0, tn, ts, te, tw, an, as, ae, aw, dn, ds, de, dw;
451: PetscScalar tleft, tright, beta, td, ad, dd, tu, au, du, v[7], bm1, coef;
452: PetscScalar ***x, bn, bs, be, bw, bu, bd, gn, gs, ge, gw, gu, gd;
453: Vec localX;
454: MatStencil c[7], row;
455: DM da;
457: PetscFunctionBeginUser;
458: PetscCall(SNESGetDM(snes, &da));
459: PetscCall(DMGetLocalVector(da, &localX));
460: PetscCall(DMDAGetInfo(da, NULL, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
461: hx = one / (PetscReal)(mx - 1);
462: hy = one / (PetscReal)(my - 1);
463: hz = one / (PetscReal)(mz - 1);
464: hxhydhz = hx * hy / hz;
465: hyhzdhx = hy * hz / hx;
466: hzhxdhy = hz * hx / hy;
467: tleft = user->tleft;
468: tright = user->tright;
469: beta = user->beta;
470: bm1 = user->bm1;
471: coef = user->coef;
473: /* Get ghost points */
474: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
475: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
476: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
477: PetscCall(DMDAVecGetArray(da, localX, &x));
479: /* Evaluate Jacobian of function */
480: for (k = zs; k < zs + zm; k++) {
481: for (j = ys; j < ys + ym; j++) {
482: for (i = xs; i < xs + xm; i++) {
483: t0 = x[k][j][i];
484: row.k = k;
485: row.j = j;
486: row.i = i;
487: if (i > 0 && i < mx - 1 && j > 0 && j < my - 1 && k > 0 && k < mz - 1) {
488: /* general interior volume */
490: tw = x[k][j][i - 1];
491: aw = 0.5 * (t0 + tw);
492: bw = PetscPowScalar(aw, bm1);
493: /* dw = bw * aw */
494: dw = PetscPowScalar(aw, beta);
495: gw = coef * bw * (t0 - tw);
497: te = x[k][j][i + 1];
498: ae = 0.5 * (t0 + te);
499: be = PetscPowScalar(ae, bm1);
500: /* de = be * ae; */
501: de = PetscPowScalar(ae, beta);
502: ge = coef * be * (te - t0);
504: ts = x[k][j - 1][i];
505: as = 0.5 * (t0 + ts);
506: bs = PetscPowScalar(as, bm1);
507: /* ds = bs * as; */
508: ds = PetscPowScalar(as, beta);
509: gs = coef * bs * (t0 - ts);
511: tn = x[k][j + 1][i];
512: an = 0.5 * (t0 + tn);
513: bn = PetscPowScalar(an, bm1);
514: /* dn = bn * an; */
515: dn = PetscPowScalar(an, beta);
516: gn = coef * bn * (tn - t0);
518: td = x[k - 1][j][i];
519: ad = 0.5 * (t0 + td);
520: bd = PetscPowScalar(ad, bm1);
521: /* dd = bd * ad; */
522: dd = PetscPowScalar(ad, beta);
523: gd = coef * bd * (t0 - td);
525: tu = x[k + 1][j][i];
526: au = 0.5 * (t0 + tu);
527: bu = PetscPowScalar(au, bm1);
528: /* du = bu * au; */
529: du = PetscPowScalar(au, beta);
530: gu = coef * bu * (tu - t0);
532: c[0].k = k - 1;
533: c[0].j = j;
534: c[0].i = i;
535: v[0] = -hxhydhz * (dd - gd);
536: c[1].k = k;
537: c[1].j = j - 1;
538: c[1].i = i;
539: v[1] = -hzhxdhy * (ds - gs);
540: c[2].k = k;
541: c[2].j = j;
542: c[2].i = i - 1;
543: v[2] = -hyhzdhx * (dw - gw);
544: c[3].k = k;
545: c[3].j = j;
546: c[3].i = i;
547: v[3] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
548: c[4].k = k;
549: c[4].j = j;
550: c[4].i = i + 1;
551: v[4] = -hyhzdhx * (de + ge);
552: c[5].k = k;
553: c[5].j = j + 1;
554: c[5].i = i;
555: v[5] = -hzhxdhy * (dn + gn);
556: c[6].k = k + 1;
557: c[6].j = j;
558: c[6].i = i;
559: v[6] = -hxhydhz * (du + gu);
560: PetscCall(MatSetValuesStencil(jac, 1, &row, 7, c, v, INSERT_VALUES));
562: } else if (i == 0) {
563: /* left-hand plane boundary */
564: tw = tleft;
565: aw = 0.5 * (t0 + tw);
566: bw = PetscPowScalar(aw, bm1);
567: /* dw = bw * aw */
568: dw = PetscPowScalar(aw, beta);
569: gw = coef * bw * (t0 - tw);
571: te = x[k][j][i + 1];
572: ae = 0.5 * (t0 + te);
573: be = PetscPowScalar(ae, bm1);
574: /* de = be * ae; */
575: de = PetscPowScalar(ae, beta);
576: ge = coef * be * (te - t0);
578: /* left-hand bottom edge */
579: if (j == 0) {
580: tn = x[k][j + 1][i];
581: an = 0.5 * (t0 + tn);
582: bn = PetscPowScalar(an, bm1);
583: /* dn = bn * an; */
584: dn = PetscPowScalar(an, beta);
585: gn = coef * bn * (tn - t0);
587: /* left-hand bottom down corner */
588: if (k == 0) {
589: tu = x[k + 1][j][i];
590: au = 0.5 * (t0 + tu);
591: bu = PetscPowScalar(au, bm1);
592: /* du = bu * au; */
593: du = PetscPowScalar(au, beta);
594: gu = coef * bu * (tu - t0);
596: c[0].k = k;
597: c[0].j = j;
598: c[0].i = i;
599: v[0] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
600: c[1].k = k;
601: c[1].j = j;
602: c[1].i = i + 1;
603: v[1] = -hyhzdhx * (de + ge);
604: c[2].k = k;
605: c[2].j = j + 1;
606: c[2].i = i;
607: v[2] = -hzhxdhy * (dn + gn);
608: c[3].k = k + 1;
609: c[3].j = j;
610: c[3].i = i;
611: v[3] = -hxhydhz * (du + gu);
612: PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
614: /* left-hand bottom interior edge */
615: } else if (k < mz - 1) {
616: tu = x[k + 1][j][i];
617: au = 0.5 * (t0 + tu);
618: bu = PetscPowScalar(au, bm1);
619: /* du = bu * au; */
620: du = PetscPowScalar(au, beta);
621: gu = coef * bu * (tu - t0);
623: td = x[k - 1][j][i];
624: ad = 0.5 * (t0 + td);
625: bd = PetscPowScalar(ad, bm1);
626: /* dd = bd * ad; */
627: dd = PetscPowScalar(ad, beta);
628: gd = coef * bd * (td - t0);
630: c[0].k = k - 1;
631: c[0].j = j;
632: c[0].i = i;
633: v[0] = -hxhydhz * (dd - gd);
634: c[1].k = k;
635: c[1].j = j;
636: c[1].i = i;
637: v[1] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
638: c[2].k = k;
639: c[2].j = j;
640: c[2].i = i + 1;
641: v[2] = -hyhzdhx * (de + ge);
642: c[3].k = k;
643: c[3].j = j + 1;
644: c[3].i = i;
645: v[3] = -hzhxdhy * (dn + gn);
646: c[4].k = k + 1;
647: c[4].j = j;
648: c[4].i = i;
649: v[4] = -hxhydhz * (du + gu);
650: PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
652: /* left-hand bottom up corner */
653: } else {
654: td = x[k - 1][j][i];
655: ad = 0.5 * (t0 + td);
656: bd = PetscPowScalar(ad, bm1);
657: /* dd = bd * ad; */
658: dd = PetscPowScalar(ad, beta);
659: gd = coef * bd * (td - t0);
661: c[0].k = k - 1;
662: c[0].j = j;
663: c[0].i = i;
664: v[0] = -hxhydhz * (dd - gd);
665: c[1].k = k;
666: c[1].j = j;
667: c[1].i = i;
668: v[1] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
669: c[2].k = k;
670: c[2].j = j;
671: c[2].i = i + 1;
672: v[2] = -hyhzdhx * (de + ge);
673: c[3].k = k;
674: c[3].j = j + 1;
675: c[3].i = i;
676: v[3] = -hzhxdhy * (dn + gn);
677: PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
678: }
680: /* left-hand top edge */
681: } else if (j == my - 1) {
682: ts = x[k][j - 1][i];
683: as = 0.5 * (t0 + ts);
684: bs = PetscPowScalar(as, bm1);
685: /* ds = bs * as; */
686: ds = PetscPowScalar(as, beta);
687: gs = coef * bs * (ts - t0);
689: /* left-hand top down corner */
690: if (k == 0) {
691: tu = x[k + 1][j][i];
692: au = 0.5 * (t0 + tu);
693: bu = PetscPowScalar(au, bm1);
694: /* du = bu * au; */
695: du = PetscPowScalar(au, beta);
696: gu = coef * bu * (tu - t0);
698: c[0].k = k;
699: c[0].j = j - 1;
700: c[0].i = i;
701: v[0] = -hzhxdhy * (ds - gs);
702: c[1].k = k;
703: c[1].j = j;
704: c[1].i = i;
705: v[1] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
706: c[2].k = k;
707: c[2].j = j;
708: c[2].i = i + 1;
709: v[2] = -hyhzdhx * (de + ge);
710: c[3].k = k + 1;
711: c[3].j = j;
712: c[3].i = i;
713: v[3] = -hxhydhz * (du + gu);
714: PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
716: /* left-hand top interior edge */
717: } else if (k < mz - 1) {
718: tu = x[k + 1][j][i];
719: au = 0.5 * (t0 + tu);
720: bu = PetscPowScalar(au, bm1);
721: /* du = bu * au; */
722: du = PetscPowScalar(au, beta);
723: gu = coef * bu * (tu - t0);
725: td = x[k - 1][j][i];
726: ad = 0.5 * (t0 + td);
727: bd = PetscPowScalar(ad, bm1);
728: /* dd = bd * ad; */
729: dd = PetscPowScalar(ad, beta);
730: gd = coef * bd * (td - t0);
732: c[0].k = k - 1;
733: c[0].j = j;
734: c[0].i = i;
735: v[0] = -hxhydhz * (dd - gd);
736: c[1].k = k;
737: c[1].j = j - 1;
738: c[1].i = i;
739: v[1] = -hzhxdhy * (ds - gs);
740: c[2].k = k;
741: c[2].j = j;
742: c[2].i = i;
743: v[2] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
744: c[3].k = k;
745: c[3].j = j;
746: c[3].i = i + 1;
747: v[3] = -hyhzdhx * (de + ge);
748: c[4].k = k + 1;
749: c[4].j = j;
750: c[4].i = i;
751: v[4] = -hxhydhz * (du + gu);
752: PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
754: /* left-hand top up corner */
755: } else {
756: td = x[k - 1][j][i];
757: ad = 0.5 * (t0 + td);
758: bd = PetscPowScalar(ad, bm1);
759: /* dd = bd * ad; */
760: dd = PetscPowScalar(ad, beta);
761: gd = coef * bd * (td - t0);
763: c[0].k = k - 1;
764: c[0].j = j;
765: c[0].i = i;
766: v[0] = -hxhydhz * (dd - gd);
767: c[1].k = k;
768: c[1].j = j - 1;
769: c[1].i = i;
770: v[1] = -hzhxdhy * (ds - gs);
771: c[2].k = k;
772: c[2].j = j;
773: c[2].i = i;
774: v[2] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
775: c[3].k = k;
776: c[3].j = j;
777: c[3].i = i + 1;
778: v[3] = -hyhzdhx * (de + ge);
779: PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
780: }
782: } else {
783: ts = x[k][j - 1][i];
784: as = 0.5 * (t0 + ts);
785: bs = PetscPowScalar(as, bm1);
786: /* ds = bs * as; */
787: ds = PetscPowScalar(as, beta);
788: gs = coef * bs * (t0 - ts);
790: tn = x[k][j + 1][i];
791: an = 0.5 * (t0 + tn);
792: bn = PetscPowScalar(an, bm1);
793: /* dn = bn * an; */
794: dn = PetscPowScalar(an, beta);
795: gn = coef * bn * (tn - t0);
797: /* left-hand down interior edge */
798: if (k == 0) {
799: tu = x[k + 1][j][i];
800: au = 0.5 * (t0 + tu);
801: bu = PetscPowScalar(au, bm1);
802: /* du = bu * au; */
803: du = PetscPowScalar(au, beta);
804: gu = coef * bu * (tu - t0);
806: c[0].k = k;
807: c[0].j = j - 1;
808: c[0].i = i;
809: v[0] = -hzhxdhy * (ds - gs);
810: c[1].k = k;
811: c[1].j = j;
812: c[1].i = i;
813: v[1] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
814: c[2].k = k;
815: c[2].j = j;
816: c[2].i = i + 1;
817: v[2] = -hyhzdhx * (de + ge);
818: c[3].k = k;
819: c[3].j = j + 1;
820: c[3].i = i;
821: v[3] = -hzhxdhy * (dn + gn);
822: c[4].k = k + 1;
823: c[4].j = j;
824: c[4].i = i;
825: v[4] = -hxhydhz * (du + gu);
826: PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
828: } else if (k == mz - 1) { /* left-hand up interior edge */
830: td = x[k - 1][j][i];
831: ad = 0.5 * (t0 + td);
832: bd = PetscPowScalar(ad, bm1);
833: /* dd = bd * ad; */
834: dd = PetscPowScalar(ad, beta);
835: gd = coef * bd * (t0 - td);
837: c[0].k = k - 1;
838: c[0].j = j;
839: c[0].i = i;
840: v[0] = -hxhydhz * (dd - gd);
841: c[1].k = k;
842: c[1].j = j - 1;
843: c[1].i = i;
844: v[1] = -hzhxdhy * (ds - gs);
845: c[2].k = k;
846: c[2].j = j;
847: c[2].i = i;
848: v[2] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
849: c[3].k = k;
850: c[3].j = j;
851: c[3].i = i + 1;
852: v[3] = -hyhzdhx * (de + ge);
853: c[4].k = k;
854: c[4].j = j + 1;
855: c[4].i = i;
856: v[4] = -hzhxdhy * (dn + gn);
857: PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
858: } else { /* left-hand interior plane */
860: td = x[k - 1][j][i];
861: ad = 0.5 * (t0 + td);
862: bd = PetscPowScalar(ad, bm1);
863: /* dd = bd * ad; */
864: dd = PetscPowScalar(ad, beta);
865: gd = coef * bd * (t0 - td);
867: tu = x[k + 1][j][i];
868: au = 0.5 * (t0 + tu);
869: bu = PetscPowScalar(au, bm1);
870: /* du = bu * au; */
871: du = PetscPowScalar(au, beta);
872: gu = coef * bu * (tu - t0);
874: c[0].k = k - 1;
875: c[0].j = j;
876: c[0].i = i;
877: v[0] = -hxhydhz * (dd - gd);
878: c[1].k = k;
879: c[1].j = j - 1;
880: c[1].i = i;
881: v[1] = -hzhxdhy * (ds - gs);
882: c[2].k = k;
883: c[2].j = j;
884: c[2].i = i;
885: v[2] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
886: c[3].k = k;
887: c[3].j = j;
888: c[3].i = i + 1;
889: v[3] = -hyhzdhx * (de + ge);
890: c[4].k = k;
891: c[4].j = j + 1;
892: c[4].i = i;
893: v[4] = -hzhxdhy * (dn + gn);
894: c[5].k = k + 1;
895: c[5].j = j;
896: c[5].i = i;
897: v[5] = -hxhydhz * (du + gu);
898: PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
899: }
900: }
902: } else if (i == mx - 1) {
903: /* right-hand plane boundary */
904: tw = x[k][j][i - 1];
905: aw = 0.5 * (t0 + tw);
906: bw = PetscPowScalar(aw, bm1);
907: /* dw = bw * aw */
908: dw = PetscPowScalar(aw, beta);
909: gw = coef * bw * (t0 - tw);
911: te = tright;
912: ae = 0.5 * (t0 + te);
913: be = PetscPowScalar(ae, bm1);
914: /* de = be * ae; */
915: de = PetscPowScalar(ae, beta);
916: ge = coef * be * (te - t0);
918: /* right-hand bottom edge */
919: if (j == 0) {
920: tn = x[k][j + 1][i];
921: an = 0.5 * (t0 + tn);
922: bn = PetscPowScalar(an, bm1);
923: /* dn = bn * an; */
924: dn = PetscPowScalar(an, beta);
925: gn = coef * bn * (tn - t0);
927: /* right-hand bottom down corner */
928: if (k == 0) {
929: tu = x[k + 1][j][i];
930: au = 0.5 * (t0 + tu);
931: bu = PetscPowScalar(au, bm1);
932: /* du = bu * au; */
933: du = PetscPowScalar(au, beta);
934: gu = coef * bu * (tu - t0);
936: c[0].k = k;
937: c[0].j = j;
938: c[0].i = i - 1;
939: v[0] = -hyhzdhx * (dw - gw);
940: c[1].k = k;
941: c[1].j = j;
942: c[1].i = i;
943: v[1] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
944: c[2].k = k;
945: c[2].j = j + 1;
946: c[2].i = i;
947: v[2] = -hzhxdhy * (dn + gn);
948: c[3].k = k + 1;
949: c[3].j = j;
950: c[3].i = i;
951: v[3] = -hxhydhz * (du + gu);
952: PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
954: /* right-hand bottom interior edge */
955: } else if (k < mz - 1) {
956: tu = x[k + 1][j][i];
957: au = 0.5 * (t0 + tu);
958: bu = PetscPowScalar(au, bm1);
959: /* du = bu * au; */
960: du = PetscPowScalar(au, beta);
961: gu = coef * bu * (tu - t0);
963: td = x[k - 1][j][i];
964: ad = 0.5 * (t0 + td);
965: bd = PetscPowScalar(ad, bm1);
966: /* dd = bd * ad; */
967: dd = PetscPowScalar(ad, beta);
968: gd = coef * bd * (td - t0);
970: c[0].k = k - 1;
971: c[0].j = j;
972: c[0].i = i;
973: v[0] = -hxhydhz * (dd - gd);
974: c[1].k = k;
975: c[1].j = j;
976: c[1].i = i - 1;
977: v[1] = -hyhzdhx * (dw - gw);
978: c[2].k = k;
979: c[2].j = j;
980: c[2].i = i;
981: v[2] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
982: c[3].k = k;
983: c[3].j = j + 1;
984: c[3].i = i;
985: v[3] = -hzhxdhy * (dn + gn);
986: c[4].k = k + 1;
987: c[4].j = j;
988: c[4].i = i;
989: v[4] = -hxhydhz * (du + gu);
990: PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
992: /* right-hand bottom up corner */
993: } else {
994: td = x[k - 1][j][i];
995: ad = 0.5 * (t0 + td);
996: bd = PetscPowScalar(ad, bm1);
997: /* dd = bd * ad; */
998: dd = PetscPowScalar(ad, beta);
999: gd = coef * bd * (td - t0);
1001: c[0].k = k - 1;
1002: c[0].j = j;
1003: c[0].i = i;
1004: v[0] = -hxhydhz * (dd - gd);
1005: c[1].k = k;
1006: c[1].j = j;
1007: c[1].i = i - 1;
1008: v[1] = -hyhzdhx * (dw - gw);
1009: c[2].k = k;
1010: c[2].j = j;
1011: c[2].i = i;
1012: v[2] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1013: c[3].k = k;
1014: c[3].j = j + 1;
1015: c[3].i = i;
1016: v[3] = -hzhxdhy * (dn + gn);
1017: PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
1018: }
1020: /* right-hand top edge */
1021: } else if (j == my - 1) {
1022: ts = x[k][j - 1][i];
1023: as = 0.5 * (t0 + ts);
1024: bs = PetscPowScalar(as, bm1);
1025: /* ds = bs * as; */
1026: ds = PetscPowScalar(as, beta);
1027: gs = coef * bs * (ts - t0);
1029: /* right-hand top down corner */
1030: if (k == 0) {
1031: tu = x[k + 1][j][i];
1032: au = 0.5 * (t0 + tu);
1033: bu = PetscPowScalar(au, bm1);
1034: /* du = bu * au; */
1035: du = PetscPowScalar(au, beta);
1036: gu = coef * bu * (tu - t0);
1038: c[0].k = k;
1039: c[0].j = j - 1;
1040: c[0].i = i;
1041: v[0] = -hzhxdhy * (ds - gs);
1042: c[1].k = k;
1043: c[1].j = j;
1044: c[1].i = i - 1;
1045: v[1] = -hyhzdhx * (dw - gw);
1046: c[2].k = k;
1047: c[2].j = j;
1048: c[2].i = i;
1049: v[2] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
1050: c[3].k = k + 1;
1051: c[3].j = j;
1052: c[3].i = i;
1053: v[3] = -hxhydhz * (du + gu);
1054: PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
1056: /* right-hand top interior edge */
1057: } else if (k < mz - 1) {
1058: tu = x[k + 1][j][i];
1059: au = 0.5 * (t0 + tu);
1060: bu = PetscPowScalar(au, bm1);
1061: /* du = bu * au; */
1062: du = PetscPowScalar(au, beta);
1063: gu = coef * bu * (tu - t0);
1065: td = x[k - 1][j][i];
1066: ad = 0.5 * (t0 + td);
1067: bd = PetscPowScalar(ad, bm1);
1068: /* dd = bd * ad; */
1069: dd = PetscPowScalar(ad, beta);
1070: gd = coef * bd * (td - t0);
1072: c[0].k = k - 1;
1073: c[0].j = j;
1074: c[0].i = i;
1075: v[0] = -hxhydhz * (dd - gd);
1076: c[1].k = k;
1077: c[1].j = j - 1;
1078: c[1].i = i;
1079: v[1] = -hzhxdhy * (ds - gs);
1080: c[2].k = k;
1081: c[2].j = j;
1082: c[2].i = i - 1;
1083: v[2] = -hyhzdhx * (dw - gw);
1084: c[3].k = k;
1085: c[3].j = j;
1086: c[3].i = i;
1087: v[3] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
1088: c[4].k = k + 1;
1089: c[4].j = j;
1090: c[4].i = i;
1091: v[4] = -hxhydhz * (du + gu);
1092: PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
1094: /* right-hand top up corner */
1095: } else {
1096: td = x[k - 1][j][i];
1097: ad = 0.5 * (t0 + td);
1098: bd = PetscPowScalar(ad, bm1);
1099: /* dd = bd * ad; */
1100: dd = PetscPowScalar(ad, beta);
1101: gd = coef * bd * (td - t0);
1103: c[0].k = k - 1;
1104: c[0].j = j;
1105: c[0].i = i;
1106: v[0] = -hxhydhz * (dd - gd);
1107: c[1].k = k;
1108: c[1].j = j - 1;
1109: c[1].i = i;
1110: v[1] = -hzhxdhy * (ds - gs);
1111: c[2].k = k;
1112: c[2].j = j;
1113: c[2].i = i - 1;
1114: v[2] = -hyhzdhx * (dw - gw);
1115: c[3].k = k;
1116: c[3].j = j;
1117: c[3].i = i;
1118: v[3] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1119: PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
1120: }
1122: } else {
1123: ts = x[k][j - 1][i];
1124: as = 0.5 * (t0 + ts);
1125: bs = PetscPowScalar(as, bm1);
1126: /* ds = bs * as; */
1127: ds = PetscPowScalar(as, beta);
1128: gs = coef * bs * (t0 - ts);
1130: tn = x[k][j + 1][i];
1131: an = 0.5 * (t0 + tn);
1132: bn = PetscPowScalar(an, bm1);
1133: /* dn = bn * an; */
1134: dn = PetscPowScalar(an, beta);
1135: gn = coef * bn * (tn - t0);
1137: /* right-hand down interior edge */
1138: if (k == 0) {
1139: tu = x[k + 1][j][i];
1140: au = 0.5 * (t0 + tu);
1141: bu = PetscPowScalar(au, bm1);
1142: /* du = bu * au; */
1143: du = PetscPowScalar(au, beta);
1144: gu = coef * bu * (tu - t0);
1146: c[0].k = k;
1147: c[0].j = j - 1;
1148: c[0].i = i;
1149: v[0] = -hzhxdhy * (ds - gs);
1150: c[1].k = k;
1151: c[1].j = j;
1152: c[1].i = i - 1;
1153: v[1] = -hyhzdhx * (dw - gw);
1154: c[2].k = k;
1155: c[2].j = j;
1156: c[2].i = i;
1157: v[2] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
1158: c[3].k = k;
1159: c[3].j = j + 1;
1160: c[3].i = i;
1161: v[3] = -hzhxdhy * (dn + gn);
1162: c[4].k = k + 1;
1163: c[4].j = j;
1164: c[4].i = i;
1165: v[4] = -hxhydhz * (du + gu);
1166: PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
1168: } else if (k == mz - 1) { /* right-hand up interior edge */
1170: td = x[k - 1][j][i];
1171: ad = 0.5 * (t0 + td);
1172: bd = PetscPowScalar(ad, bm1);
1173: /* dd = bd * ad; */
1174: dd = PetscPowScalar(ad, beta);
1175: gd = coef * bd * (t0 - td);
1177: c[0].k = k - 1;
1178: c[0].j = j;
1179: c[0].i = i;
1180: v[0] = -hxhydhz * (dd - gd);
1181: c[1].k = k;
1182: c[1].j = j - 1;
1183: c[1].i = i;
1184: v[1] = -hzhxdhy * (ds - gs);
1185: c[2].k = k;
1186: c[2].j = j;
1187: c[2].i = i - 1;
1188: v[2] = -hyhzdhx * (dw - gw);
1189: c[3].k = k;
1190: c[3].j = j;
1191: c[3].i = i;
1192: v[3] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1193: c[4].k = k;
1194: c[4].j = j + 1;
1195: c[4].i = i;
1196: v[4] = -hzhxdhy * (dn + gn);
1197: PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
1199: } else { /* right-hand interior plane */
1201: td = x[k - 1][j][i];
1202: ad = 0.5 * (t0 + td);
1203: bd = PetscPowScalar(ad, bm1);
1204: /* dd = bd * ad; */
1205: dd = PetscPowScalar(ad, beta);
1206: gd = coef * bd * (t0 - td);
1208: tu = x[k + 1][j][i];
1209: au = 0.5 * (t0 + tu);
1210: bu = PetscPowScalar(au, bm1);
1211: /* du = bu * au; */
1212: du = PetscPowScalar(au, beta);
1213: gu = coef * bu * (tu - t0);
1215: c[0].k = k - 1;
1216: c[0].j = j;
1217: c[0].i = i;
1218: v[0] = -hxhydhz * (dd - gd);
1219: c[1].k = k;
1220: c[1].j = j - 1;
1221: c[1].i = i;
1222: v[1] = -hzhxdhy * (ds - gs);
1223: c[2].k = k;
1224: c[2].j = j;
1225: c[2].i = i - 1;
1226: v[2] = -hyhzdhx * (dw - gw);
1227: c[3].k = k;
1228: c[3].j = j;
1229: c[3].i = i;
1230: v[3] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
1231: c[4].k = k;
1232: c[4].j = j + 1;
1233: c[4].i = i;
1234: v[4] = -hzhxdhy * (dn + gn);
1235: c[5].k = k + 1;
1236: c[5].j = j;
1237: c[5].i = i;
1238: v[5] = -hxhydhz * (du + gu);
1239: PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
1240: }
1241: }
1243: } else if (j == 0) {
1244: tw = x[k][j][i - 1];
1245: aw = 0.5 * (t0 + tw);
1246: bw = PetscPowScalar(aw, bm1);
1247: /* dw = bw * aw */
1248: dw = PetscPowScalar(aw, beta);
1249: gw = coef * bw * (t0 - tw);
1251: te = x[k][j][i + 1];
1252: ae = 0.5 * (t0 + te);
1253: be = PetscPowScalar(ae, bm1);
1254: /* de = be * ae; */
1255: de = PetscPowScalar(ae, beta);
1256: ge = coef * be * (te - t0);
1258: tn = x[k][j + 1][i];
1259: an = 0.5 * (t0 + tn);
1260: bn = PetscPowScalar(an, bm1);
1261: /* dn = bn * an; */
1262: dn = PetscPowScalar(an, beta);
1263: gn = coef * bn * (tn - t0);
1265: /* bottom down interior edge */
1266: if (k == 0) {
1267: tu = x[k + 1][j][i];
1268: au = 0.5 * (t0 + tu);
1269: bu = PetscPowScalar(au, bm1);
1270: /* du = bu * au; */
1271: du = PetscPowScalar(au, beta);
1272: gu = coef * bu * (tu - t0);
1274: c[0].k = k;
1275: c[0].j = j;
1276: c[0].i = i - 1;
1277: v[0] = -hyhzdhx * (dw - gw);
1278: c[1].k = k;
1279: c[1].j = j;
1280: c[1].i = i;
1281: v[1] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
1282: c[2].k = k;
1283: c[2].j = j;
1284: c[2].i = i + 1;
1285: v[2] = -hyhzdhx * (de + ge);
1286: c[3].k = k;
1287: c[3].j = j + 1;
1288: c[3].i = i;
1289: v[3] = -hzhxdhy * (dn + gn);
1290: c[4].k = k + 1;
1291: c[4].j = j;
1292: c[4].i = i;
1293: v[4] = -hxhydhz * (du + gu);
1294: PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
1296: } else if (k == mz - 1) { /* bottom up interior edge */
1298: td = x[k - 1][j][i];
1299: ad = 0.5 * (t0 + td);
1300: bd = PetscPowScalar(ad, bm1);
1301: /* dd = bd * ad; */
1302: dd = PetscPowScalar(ad, beta);
1303: gd = coef * bd * (td - t0);
1305: c[0].k = k - 1;
1306: c[0].j = j;
1307: c[0].i = i;
1308: v[0] = -hxhydhz * (dd - gd);
1309: c[1].k = k;
1310: c[1].j = j;
1311: c[1].i = i - 1;
1312: v[1] = -hyhzdhx * (dw - gw);
1313: c[2].k = k;
1314: c[2].j = j;
1315: c[2].i = i;
1316: v[2] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1317: c[3].k = k;
1318: c[3].j = j;
1319: c[3].i = i + 1;
1320: v[3] = -hyhzdhx * (de + ge);
1321: c[4].k = k;
1322: c[4].j = j + 1;
1323: c[4].i = i;
1324: v[4] = -hzhxdhy * (dn + gn);
1325: PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
1327: } else { /* bottom interior plane */
1329: tu = x[k + 1][j][i];
1330: au = 0.5 * (t0 + tu);
1331: bu = PetscPowScalar(au, bm1);
1332: /* du = bu * au; */
1333: du = PetscPowScalar(au, beta);
1334: gu = coef * bu * (tu - t0);
1336: td = x[k - 1][j][i];
1337: ad = 0.5 * (t0 + td);
1338: bd = PetscPowScalar(ad, bm1);
1339: /* dd = bd * ad; */
1340: dd = PetscPowScalar(ad, beta);
1341: gd = coef * bd * (td - t0);
1343: c[0].k = k - 1;
1344: c[0].j = j;
1345: c[0].i = i;
1346: v[0] = -hxhydhz * (dd - gd);
1347: c[1].k = k;
1348: c[1].j = j;
1349: c[1].i = i - 1;
1350: v[1] = -hyhzdhx * (dw - gw);
1351: c[2].k = k;
1352: c[2].j = j;
1353: c[2].i = i;
1354: v[2] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
1355: c[3].k = k;
1356: c[3].j = j;
1357: c[3].i = i + 1;
1358: v[3] = -hyhzdhx * (de + ge);
1359: c[4].k = k;
1360: c[4].j = j + 1;
1361: c[4].i = i;
1362: v[4] = -hzhxdhy * (dn + gn);
1363: c[5].k = k + 1;
1364: c[5].j = j;
1365: c[5].i = i;
1366: v[5] = -hxhydhz * (du + gu);
1367: PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
1368: }
1370: } else if (j == my - 1) {
1371: tw = x[k][j][i - 1];
1372: aw = 0.5 * (t0 + tw);
1373: bw = PetscPowScalar(aw, bm1);
1374: /* dw = bw * aw */
1375: dw = PetscPowScalar(aw, beta);
1376: gw = coef * bw * (t0 - tw);
1378: te = x[k][j][i + 1];
1379: ae = 0.5 * (t0 + te);
1380: be = PetscPowScalar(ae, bm1);
1381: /* de = be * ae; */
1382: de = PetscPowScalar(ae, beta);
1383: ge = coef * be * (te - t0);
1385: ts = x[k][j - 1][i];
1386: as = 0.5 * (t0 + ts);
1387: bs = PetscPowScalar(as, bm1);
1388: /* ds = bs * as; */
1389: ds = PetscPowScalar(as, beta);
1390: gs = coef * bs * (t0 - ts);
1392: /* top down interior edge */
1393: if (k == 0) {
1394: tu = x[k + 1][j][i];
1395: au = 0.5 * (t0 + tu);
1396: bu = PetscPowScalar(au, bm1);
1397: /* du = bu * au; */
1398: du = PetscPowScalar(au, beta);
1399: gu = coef * bu * (tu - t0);
1401: c[0].k = k;
1402: c[0].j = j - 1;
1403: c[0].i = i;
1404: v[0] = -hzhxdhy * (ds - gs);
1405: c[1].k = k;
1406: c[1].j = j;
1407: c[1].i = i - 1;
1408: v[1] = -hyhzdhx * (dw - gw);
1409: c[2].k = k;
1410: c[2].j = j;
1411: c[2].i = i;
1412: v[2] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
1413: c[3].k = k;
1414: c[3].j = j;
1415: c[3].i = i + 1;
1416: v[3] = -hyhzdhx * (de + ge);
1417: c[4].k = k + 1;
1418: c[4].j = j;
1419: c[4].i = i;
1420: v[4] = -hxhydhz * (du + gu);
1421: PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
1423: } else if (k == mz - 1) { /* top up interior edge */
1425: td = x[k - 1][j][i];
1426: ad = 0.5 * (t0 + td);
1427: bd = PetscPowScalar(ad, bm1);
1428: /* dd = bd * ad; */
1429: dd = PetscPowScalar(ad, beta);
1430: gd = coef * bd * (td - t0);
1432: c[0].k = k - 1;
1433: c[0].j = j;
1434: c[0].i = i;
1435: v[0] = -hxhydhz * (dd - gd);
1436: c[1].k = k;
1437: c[1].j = j - 1;
1438: c[1].i = i;
1439: v[1] = -hzhxdhy * (ds - gs);
1440: c[2].k = k;
1441: c[2].j = j;
1442: c[2].i = i - 1;
1443: v[2] = -hyhzdhx * (dw - gw);
1444: c[3].k = k;
1445: c[3].j = j;
1446: c[3].i = i;
1447: v[3] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1448: c[4].k = k;
1449: c[4].j = j;
1450: c[4].i = i + 1;
1451: v[4] = -hyhzdhx * (de + ge);
1452: PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
1454: } else { /* top interior plane */
1456: tu = x[k + 1][j][i];
1457: au = 0.5 * (t0 + tu);
1458: bu = PetscPowScalar(au, bm1);
1459: /* du = bu * au; */
1460: du = PetscPowScalar(au, beta);
1461: gu = coef * bu * (tu - t0);
1463: td = x[k - 1][j][i];
1464: ad = 0.5 * (t0 + td);
1465: bd = PetscPowScalar(ad, bm1);
1466: /* dd = bd * ad; */
1467: dd = PetscPowScalar(ad, beta);
1468: gd = coef * bd * (td - t0);
1470: c[0].k = k - 1;
1471: c[0].j = j;
1472: c[0].i = i;
1473: v[0] = -hxhydhz * (dd - gd);
1474: c[1].k = k;
1475: c[1].j = j - 1;
1476: c[1].i = i;
1477: v[1] = -hzhxdhy * (ds - gs);
1478: c[2].k = k;
1479: c[2].j = j;
1480: c[2].i = i - 1;
1481: v[2] = -hyhzdhx * (dw - gw);
1482: c[3].k = k;
1483: c[3].j = j;
1484: c[3].i = i;
1485: v[3] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
1486: c[4].k = k;
1487: c[4].j = j;
1488: c[4].i = i + 1;
1489: v[4] = -hyhzdhx * (de + ge);
1490: c[5].k = k + 1;
1491: c[5].j = j;
1492: c[5].i = i;
1493: v[5] = -hxhydhz * (du + gu);
1494: PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
1495: }
1497: } else if (k == 0) {
1498: /* down interior plane */
1500: tw = x[k][j][i - 1];
1501: aw = 0.5 * (t0 + tw);
1502: bw = PetscPowScalar(aw, bm1);
1503: /* dw = bw * aw */
1504: dw = PetscPowScalar(aw, beta);
1505: gw = coef * bw * (t0 - tw);
1507: te = x[k][j][i + 1];
1508: ae = 0.5 * (t0 + te);
1509: be = PetscPowScalar(ae, bm1);
1510: /* de = be * ae; */
1511: de = PetscPowScalar(ae, beta);
1512: ge = coef * be * (te - t0);
1514: ts = x[k][j - 1][i];
1515: as = 0.5 * (t0 + ts);
1516: bs = PetscPowScalar(as, bm1);
1517: /* ds = bs * as; */
1518: ds = PetscPowScalar(as, beta);
1519: gs = coef * bs * (t0 - ts);
1521: tn = x[k][j + 1][i];
1522: an = 0.5 * (t0 + tn);
1523: bn = PetscPowScalar(an, bm1);
1524: /* dn = bn * an; */
1525: dn = PetscPowScalar(an, beta);
1526: gn = coef * bn * (tn - t0);
1528: tu = x[k + 1][j][i];
1529: au = 0.5 * (t0 + tu);
1530: bu = PetscPowScalar(au, bm1);
1531: /* du = bu * au; */
1532: du = PetscPowScalar(au, beta);
1533: gu = coef * bu * (tu - t0);
1535: c[0].k = k;
1536: c[0].j = j - 1;
1537: c[0].i = i;
1538: v[0] = -hzhxdhy * (ds - gs);
1539: c[1].k = k;
1540: c[1].j = j;
1541: c[1].i = i - 1;
1542: v[1] = -hyhzdhx * (dw - gw);
1543: c[2].k = k;
1544: c[2].j = j;
1545: c[2].i = i;
1546: v[2] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
1547: c[3].k = k;
1548: c[3].j = j;
1549: c[3].i = i + 1;
1550: v[3] = -hyhzdhx * (de + ge);
1551: c[4].k = k;
1552: c[4].j = j + 1;
1553: c[4].i = i;
1554: v[4] = -hzhxdhy * (dn + gn);
1555: c[5].k = k + 1;
1556: c[5].j = j;
1557: c[5].i = i;
1558: v[5] = -hxhydhz * (du + gu);
1559: PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
1561: } else if (k == mz - 1) {
1562: /* up interior plane */
1564: tw = x[k][j][i - 1];
1565: aw = 0.5 * (t0 + tw);
1566: bw = PetscPowScalar(aw, bm1);
1567: /* dw = bw * aw */
1568: dw = PetscPowScalar(aw, beta);
1569: gw = coef * bw * (t0 - tw);
1571: te = x[k][j][i + 1];
1572: ae = 0.5 * (t0 + te);
1573: be = PetscPowScalar(ae, bm1);
1574: /* de = be * ae; */
1575: de = PetscPowScalar(ae, beta);
1576: ge = coef * be * (te - t0);
1578: ts = x[k][j - 1][i];
1579: as = 0.5 * (t0 + ts);
1580: bs = PetscPowScalar(as, bm1);
1581: /* ds = bs * as; */
1582: ds = PetscPowScalar(as, beta);
1583: gs = coef * bs * (t0 - ts);
1585: tn = x[k][j + 1][i];
1586: an = 0.5 * (t0 + tn);
1587: bn = PetscPowScalar(an, bm1);
1588: /* dn = bn * an; */
1589: dn = PetscPowScalar(an, beta);
1590: gn = coef * bn * (tn - t0);
1592: td = x[k - 1][j][i];
1593: ad = 0.5 * (t0 + td);
1594: bd = PetscPowScalar(ad, bm1);
1595: /* dd = bd * ad; */
1596: dd = PetscPowScalar(ad, beta);
1597: gd = coef * bd * (t0 - td);
1599: c[0].k = k - 1;
1600: c[0].j = j;
1601: c[0].i = i;
1602: v[0] = -hxhydhz * (dd - gd);
1603: c[1].k = k;
1604: c[1].j = j - 1;
1605: c[1].i = i;
1606: v[1] = -hzhxdhy * (ds - gs);
1607: c[2].k = k;
1608: c[2].j = j;
1609: c[2].i = i - 1;
1610: v[2] = -hyhzdhx * (dw - gw);
1611: c[3].k = k;
1612: c[3].j = j;
1613: c[3].i = i;
1614: v[3] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1615: c[4].k = k;
1616: c[4].j = j;
1617: c[4].i = i + 1;
1618: v[4] = -hyhzdhx * (de + ge);
1619: c[5].k = k;
1620: c[5].j = j + 1;
1621: c[5].i = i;
1622: v[5] = -hzhxdhy * (dn + gn);
1623: PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
1624: }
1625: }
1626: }
1627: }
1628: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
1629: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
1630: if (jac != J) {
1631: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1632: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1633: }
1634: PetscCall(DMDAVecRestoreArray(da, localX, &x));
1635: PetscCall(DMRestoreLocalVector(da, &localX));
1637: PetscCall(PetscLogFlops((41.0 + 8.0 * POWFLOP) * xm * ym));
1638: PetscFunctionReturn(PETSC_SUCCESS);
1639: }
1641: PetscErrorCode TestConvergence(SNES snes, PETSC_UNUSED PetscInt it, PETSC_UNUSED PetscReal xnorm, PETSC_UNUSED PetscReal snorm, PETSC_UNUSED PetscReal fnorm, SNESConvergedReason *reason, void *ctx)
1642: {
1643: AppCtx *user = (AppCtx *)ctx;
1645: PetscFunctionBeginUser;
1646: if (user->converged) *reason = SNES_CONVERGED_USER;
1647: else *reason = SNES_DIVERGED_USER;
1648: PetscFunctionReturn(PETSC_SUCCESS);
1649: }
1651: /*TEST
1653: test:
1654: nsize: 4
1655: args: -snes_monitor_short -pc_mg_type full -ksp_type fgmres -pc_type mg -snes_view -pc_mg_levels 2 -pc_mg_galerkin pmat
1656: requires: !single
1658: test:
1659: suffix: 2
1660: args: -snes_converged_reason -use_convergence_test -mark_converged
1662: test:
1663: suffix: 3
1664: args: -snes_converged_reason -use_convergence_test -mark_converged 0
1666: TEST*/