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