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: /*T
15: Concepts: SNES^solving a system of nonlinear equations
16: Concepts: DM^using distributed arrays
17: Concepts: multigrid;
18: Processors: n
19: T*/
21: /*
23: This example models the partial differential equation
25: - Div(alpha* T^beta (GRAD T)) = 0.
27: where beta = 2.5 and alpha = 1.0
29: BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = dT/dn_up = dT/dn_down = 0.
31: in the unit square, which is uniformly discretized in each of x and
32: y in this simple encoding. The degrees of freedom are cell centered.
34: A finite volume approximation with the usual 7-point stencil
35: is used to discretize the boundary value problem to obtain a
36: nonlinear system of equations.
38: This code was contributed by Nickolas Jovanovic based on ex18.c
40: */
42: #include <petscsnes.h>
43: #include <petscdm.h>
44: #include <petscdmda.h>
46: /* User-defined application context */
48: typedef struct {
49: PetscReal tleft,tright; /* Dirichlet boundary conditions */
50: PetscReal beta,bm1,coef; /* nonlinear diffusivity parameterizations */
51: } AppCtx;
53: #define POWFLOP 5 /* assume a pow() takes five flops */
55: extern PetscErrorCode FormInitialGuess(SNES,Vec,void*);
56: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
57: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
59: int main(int argc,char **argv)
60: {
61: SNES snes;
62: AppCtx user;
63: PetscInt its,lits;
64: PetscReal litspit;
65: DM da;
67: PetscInitialize(&argc,&argv,NULL,help);
68: /* set problem parameters */
69: user.tleft = 1.0;
70: user.tright = 0.1;
71: user.beta = 2.5;
72: PetscOptionsGetReal(NULL,NULL,"-tleft",&user.tleft,NULL);
73: PetscOptionsGetReal(NULL,NULL,"-tright",&user.tright,NULL);
74: PetscOptionsGetReal(NULL,NULL,"-beta",&user.beta,NULL);
75: user.bm1 = user.beta - 1.0;
76: user.coef = user.beta/2.0;
78: /*
79: Set the DMDA (grid structure) for the grids.
80: */
81: 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);
82: DMSetFromOptions(da);
83: DMSetUp(da);
84: DMSetApplicationContext(da,&user);
86: /*
87: Create the nonlinear solver
88: */
89: SNESCreate(PETSC_COMM_WORLD,&snes);
90: SNESSetDM(snes,da);
91: SNESSetFunction(snes,NULL,FormFunction,&user);
92: SNESSetJacobian(snes,NULL,NULL,FormJacobian,&user);
93: SNESSetFromOptions(snes);
94: SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);
96: SNESSolve(snes,NULL,NULL);
97: SNESGetIterationNumber(snes,&its);
98: SNESGetLinearSolveIterations(snes,&lits);
99: litspit = ((PetscReal)lits)/((PetscReal)its);
100: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
101: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
102: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",(double)litspit);
104: SNESDestroy(&snes);
105: DMDestroy(&da);
106: PetscFinalize();
107: return 0;
108: }
109: /* -------------------- Form initial approximation ----------------- */
110: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
111: {
112: AppCtx *user;
113: PetscInt i,j,k,xs,ys,xm,ym,zs,zm;
114: PetscScalar ***x;
115: DM da;
118: SNESGetDM(snes,&da);
119: DMGetApplicationContext(da,&user);
120: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
121: DMDAVecGetArray(da,X,&x);
123: /* Compute initial guess */
124: for (k=zs; k<zs+zm; k++) {
125: for (j=ys; j<ys+ym; j++) {
126: for (i=xs; i<xs+xm; i++) {
127: x[k][j][i] = user->tleft;
128: }
129: }
130: }
131: DMDAVecRestoreArray(da,X,&x);
132: return 0;
133: }
134: /* -------------------- Evaluate Function F(x) --------------------- */
135: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
136: {
137: AppCtx *user = (AppCtx*)ptr;
138: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
139: PetscScalar zero = 0.0,one = 1.0;
140: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
141: 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;
142: PetscScalar tleft,tright,beta,td,ad,dd,fd=0.0,tu,au,du=0.0,fu=0.0;
143: PetscScalar ***x,***f;
144: Vec localX;
145: DM da;
148: SNESGetDM(snes,&da);
149: DMGetLocalVector(da,&localX);
150: DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
151: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
152: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
153: tleft = user->tleft; tright = user->tright;
154: beta = user->beta;
156: /* Get ghost points */
157: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
158: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
159: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
160: DMDAVecGetArray(da,localX,&x);
161: DMDAVecGetArray(da,F,&f);
163: /* Evaluate function */
164: for (k=zs; k<zs+zm; k++) {
165: for (j=ys; j<ys+ym; j++) {
166: for (i=xs; i<xs+xm; i++) {
167: t0 = x[k][j][i];
169: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
171: /* general interior volume */
173: tw = x[k][j][i-1];
174: aw = 0.5*(t0 + tw);
175: dw = PetscPowScalar(aw,beta);
176: fw = dw*(t0 - tw);
178: te = x[k][j][i+1];
179: ae = 0.5*(t0 + te);
180: de = PetscPowScalar(ae,beta);
181: fe = de*(te - t0);
183: ts = x[k][j-1][i];
184: as = 0.5*(t0 + ts);
185: ds = PetscPowScalar(as,beta);
186: fs = ds*(t0 - ts);
188: tn = x[k][j+1][i];
189: an = 0.5*(t0 + tn);
190: dn = PetscPowScalar(an,beta);
191: fn = dn*(tn - t0);
193: td = x[k-1][j][i];
194: ad = 0.5*(t0 + td);
195: dd = PetscPowScalar(ad,beta);
196: fd = dd*(t0 - td);
198: tu = x[k+1][j][i];
199: au = 0.5*(t0 + tu);
200: du = PetscPowScalar(au,beta);
201: fu = du*(tu - t0);
203: } 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) {
254: /* right-hand (east) boundary */
255: tw = x[k][j][i-1];
256: aw = 0.5*(t0 + tw);
257: dw = PetscPowScalar(aw,beta);
258: fw = dw*(t0 - tw);
260: te = tright;
261: ae = 0.5*(t0 + te);
262: de = PetscPowScalar(ae,beta);
263: fe = de*(te - t0);
265: if (j > 0) {
266: ts = x[k][j-1][i];
267: as = 0.5*(t0 + ts);
268: ds = PetscPowScalar(as,beta);
269: fs = ds*(t0 - ts);
270: } else {
271: fs = zero;
272: }
274: if (j < my-1) {
275: tn = x[k][j+1][i];
276: an = 0.5*(t0 + tn);
277: dn = PetscPowScalar(an,beta);
278: fn = dn*(tn - t0);
279: } else {
280: fn = zero;
281: }
283: if (k > 0) {
284: td = x[k-1][j][i];
285: ad = 0.5*(t0 + td);
286: dd = PetscPowScalar(ad,beta);
287: fd = dd*(t0 - td);
288: } else {
289: fd = zero;
290: }
292: if (k < mz-1) {
293: tu = x[k+1][j][i];
294: au = 0.5*(t0 + tu);
295: du = PetscPowScalar(au,beta);
296: fu = du*(tu - t0);
297: } else {
298: fu = zero;
299: }
301: } else if (j == 0) {
303: /* bottom (south) boundary, and i <> 0 or mx-1 */
304: tw = x[k][j][i-1];
305: aw = 0.5*(t0 + tw);
306: dw = PetscPowScalar(aw,beta);
307: fw = dw*(t0 - tw);
309: te = x[k][j][i+1];
310: ae = 0.5*(t0 + te);
311: de = PetscPowScalar(ae,beta);
312: fe = de*(te - t0);
314: fs = zero;
316: tn = x[k][j+1][i];
317: an = 0.5*(t0 + tn);
318: dn = PetscPowScalar(an,beta);
319: fn = dn*(tn - t0);
321: if (k > 0) {
322: td = x[k-1][j][i];
323: ad = 0.5*(t0 + td);
324: dd = PetscPowScalar(ad,beta);
325: fd = dd*(t0 - td);
326: } else {
327: fd = zero;
328: }
330: if (k < mz-1) {
331: tu = x[k+1][j][i];
332: au = 0.5*(t0 + tu);
333: du = PetscPowScalar(au,beta);
334: fu = du*(tu - t0);
335: } else {
336: fu = zero;
337: }
339: } else if (j == my-1) {
341: /* top (north) boundary, and i <> 0 or mx-1 */
342: tw = x[k][j][i-1];
343: aw = 0.5*(t0 + tw);
344: dw = PetscPowScalar(aw,beta);
345: fw = dw*(t0 - tw);
347: te = x[k][j][i+1];
348: ae = 0.5*(t0 + te);
349: de = PetscPowScalar(ae,beta);
350: fe = de*(te - t0);
352: ts = x[k][j-1][i];
353: as = 0.5*(t0 + ts);
354: ds = PetscPowScalar(as,beta);
355: fs = ds*(t0 - ts);
357: fn = zero;
359: if (k > 0) {
360: td = x[k-1][j][i];
361: ad = 0.5*(t0 + td);
362: dd = PetscPowScalar(ad,beta);
363: fd = dd*(t0 - td);
364: } else {
365: fd = zero;
366: }
368: if (k < mz-1) {
369: tu = x[k+1][j][i];
370: au = 0.5*(t0 + tu);
371: du = PetscPowScalar(au,beta);
372: fu = du*(tu - t0);
373: } else {
374: fu = zero;
375: }
377: } else if (k == 0) {
379: /* down boundary (interior only) */
380: tw = x[k][j][i-1];
381: aw = 0.5*(t0 + tw);
382: dw = PetscPowScalar(aw,beta);
383: fw = dw*(t0 - tw);
385: te = x[k][j][i+1];
386: ae = 0.5*(t0 + te);
387: de = PetscPowScalar(ae,beta);
388: fe = de*(te - t0);
390: ts = x[k][j-1][i];
391: as = 0.5*(t0 + ts);
392: ds = PetscPowScalar(as,beta);
393: fs = ds*(t0 - ts);
395: tn = x[k][j+1][i];
396: an = 0.5*(t0 + tn);
397: dn = PetscPowScalar(an,beta);
398: fn = dn*(tn - t0);
400: fd = zero;
402: tu = x[k+1][j][i];
403: au = 0.5*(t0 + tu);
404: du = PetscPowScalar(au,beta);
405: fu = du*(tu - t0);
407: } else if (k == mz-1) {
409: /* up boundary (interior only) */
410: tw = x[k][j][i-1];
411: aw = 0.5*(t0 + tw);
412: dw = PetscPowScalar(aw,beta);
413: fw = dw*(t0 - tw);
415: te = x[k][j][i+1];
416: ae = 0.5*(t0 + te);
417: de = PetscPowScalar(ae,beta);
418: fe = de*(te - t0);
420: ts = x[k][j-1][i];
421: as = 0.5*(t0 + ts);
422: ds = PetscPowScalar(as,beta);
423: fs = ds*(t0 - ts);
425: tn = x[k][j+1][i];
426: an = 0.5*(t0 + tn);
427: dn = PetscPowScalar(an,beta);
428: fn = dn*(tn - t0);
430: td = x[k-1][j][i];
431: ad = 0.5*(t0 + td);
432: dd = PetscPowScalar(ad,beta);
433: fd = dd*(t0 - td);
435: fu = zero;
436: }
438: f[k][j][i] = -hyhzdhx*(fe-fw) - hzhxdhy*(fn-fs) - hxhydhz*(fu-fd);
439: }
440: }
441: }
442: DMDAVecRestoreArray(da,localX,&x);
443: DMDAVecRestoreArray(da,F,&f);
444: DMRestoreLocalVector(da,&localX);
445: PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
446: return 0;
447: }
448: /* -------------------- Evaluate Jacobian F(x) --------------------- */
449: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
450: {
451: AppCtx *user = (AppCtx*)ptr;
452: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
453: PetscScalar one = 1.0;
454: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
455: PetscScalar t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw;
456: PetscScalar tleft,tright,beta,td,ad,dd,tu,au,du,v[7],bm1,coef;
457: PetscScalar ***x,bn,bs,be,bw,bu,bd,gn,gs,ge,gw,gu,gd;
458: Vec localX;
459: MatStencil c[7],row;
460: DM da;
463: SNESGetDM(snes,&da);
464: DMGetLocalVector(da,&localX);
465: DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
466: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
467: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
468: tleft = user->tleft; tright = user->tright;
469: beta = user->beta; bm1 = user->bm1; coef = user->coef;
471: /* Get ghost points */
472: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
473: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
474: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
475: DMDAVecGetArray(da,localX,&x);
477: /* Evaluate Jacobian of function */
478: for (k=zs; k<zs+zm; k++) {
479: for (j=ys; j<ys+ym; j++) {
480: for (i=xs; i<xs+xm; i++) {
481: t0 = x[k][j][i];
482: row.k = k; row.j = j; row.i = i;
483: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
485: /* general interior volume */
487: tw = x[k][j][i-1];
488: aw = 0.5*(t0 + tw);
489: bw = PetscPowScalar(aw,bm1);
490: /* dw = bw * aw */
491: dw = PetscPowScalar(aw,beta);
492: gw = coef*bw*(t0 - tw);
494: te = x[k][j][i+1];
495: ae = 0.5*(t0 + te);
496: be = PetscPowScalar(ae,bm1);
497: /* de = be * ae; */
498: de = PetscPowScalar(ae,beta);
499: ge = coef*be*(te - t0);
501: ts = x[k][j-1][i];
502: as = 0.5*(t0 + ts);
503: bs = PetscPowScalar(as,bm1);
504: /* ds = bs * as; */
505: ds = PetscPowScalar(as,beta);
506: gs = coef*bs*(t0 - ts);
508: tn = x[k][j+1][i];
509: an = 0.5*(t0 + tn);
510: bn = PetscPowScalar(an,bm1);
511: /* dn = bn * an; */
512: dn = PetscPowScalar(an,beta);
513: gn = coef*bn*(tn - t0);
515: td = x[k-1][j][i];
516: ad = 0.5*(t0 + td);
517: bd = PetscPowScalar(ad,bm1);
518: /* dd = bd * ad; */
519: dd = PetscPowScalar(ad,beta);
520: gd = coef*bd*(t0 - td);
522: tu = x[k+1][j][i];
523: au = 0.5*(t0 + tu);
524: bu = PetscPowScalar(au,bm1);
525: /* du = bu * au; */
526: du = PetscPowScalar(au,beta);
527: gu = coef*bu*(tu - t0);
529: c[0].k = k-1; c[0].j = j; c[0].i = i; v[0] = -hxhydhz*(dd - gd);
530: c[1].k = k; c[1].j = j-1; c[1].i = i;
531: v[1] = -hzhxdhy*(ds - gs);
532: c[2].k = k; c[2].j = j; c[2].i = i-1;
533: v[2] = -hyhzdhx*(dw - gw);
534: c[3].k = k; c[3].j = j; c[3].i = i;
535: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
536: c[4].k = k; c[4].j = j; c[4].i = i+1;
537: v[4] = -hyhzdhx*(de + ge);
538: c[5].k = k; c[5].j = j+1; c[5].i = i;
539: v[5] = -hzhxdhy*(dn + gn);
540: c[6].k = k+1; c[6].j = j; c[6].i = i;
541: v[6] = -hxhydhz*(du + gu);
542: MatSetValuesStencil(jac,1,&row,7,c,v,INSERT_VALUES);
544: } else if (i == 0) {
546: /* left-hand plane boundary */
547: tw = tleft;
548: aw = 0.5*(t0 + tw);
549: bw = PetscPowScalar(aw,bm1);
550: /* dw = bw * aw */
551: dw = PetscPowScalar(aw,beta);
552: gw = coef*bw*(t0 - tw);
554: te = x[k][j][i+1];
555: ae = 0.5*(t0 + te);
556: be = PetscPowScalar(ae,bm1);
557: /* de = be * ae; */
558: de = PetscPowScalar(ae,beta);
559: ge = coef*be*(te - t0);
561: /* left-hand bottom edge */
562: if (j == 0) {
564: tn = x[k][j+1][i];
565: an = 0.5*(t0 + tn);
566: bn = PetscPowScalar(an,bm1);
567: /* dn = bn * an; */
568: dn = PetscPowScalar(an,beta);
569: gn = coef*bn*(tn - t0);
571: /* left-hand bottom down corner */
572: if (k == 0) {
574: tu = x[k+1][j][i];
575: au = 0.5*(t0 + tu);
576: bu = PetscPowScalar(au,bm1);
577: /* du = bu * au; */
578: du = PetscPowScalar(au,beta);
579: gu = coef*bu*(tu - t0);
581: c[0].k = k; c[0].j = j; c[0].i = i;
582: v[0] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
583: c[1].k = k; c[1].j = j; c[1].i = i+1;
584: v[1] = -hyhzdhx*(de + ge);
585: c[2].k = k; c[2].j = j+1; c[2].i = i;
586: v[2] = -hzhxdhy*(dn + gn);
587: c[3].k = k+1; c[3].j = j; c[3].i = i;
588: v[3] = -hxhydhz*(du + gu);
589: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
591: /* left-hand bottom interior edge */
592: } else if (k < mz-1) {
594: tu = x[k+1][j][i];
595: au = 0.5*(t0 + tu);
596: bu = PetscPowScalar(au,bm1);
597: /* du = bu * au; */
598: du = PetscPowScalar(au,beta);
599: gu = coef*bu*(tu - t0);
601: td = x[k-1][j][i];
602: ad = 0.5*(t0 + td);
603: bd = PetscPowScalar(ad,bm1);
604: /* dd = bd * ad; */
605: dd = PetscPowScalar(ad,beta);
606: gd = coef*bd*(td - t0);
608: c[0].k = k-1; c[0].j = j; c[0].i = i;
609: v[0] = -hxhydhz*(dd - gd);
610: c[1].k = k; c[1].j = j; c[1].i = i;
611: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
612: c[2].k = k; c[2].j = j; c[2].i = i+1;
613: v[2] = -hyhzdhx*(de + ge);
614: c[3].k = k; c[3].j = j+1; c[3].i = i;
615: v[3] = -hzhxdhy*(dn + gn);
616: c[4].k = k+1; c[4].j = j; c[4].i = i;
617: v[4] = -hxhydhz*(du + gu);
618: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
620: /* left-hand bottom up corner */
621: } else {
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; c[0].j = j; c[0].i = i;
631: v[0] = -hxhydhz*(dd - gd);
632: c[1].k = k; c[1].j = j; c[1].i = i;
633: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
634: c[2].k = k; c[2].j = j; c[2].i = i+1;
635: v[2] = -hyhzdhx*(de + ge);
636: c[3].k = k; c[3].j = j+1; c[3].i = i;
637: v[3] = -hzhxdhy*(dn + gn);
638: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
639: }
641: /* left-hand top edge */
642: } else if (j == my-1) {
644: ts = x[k][j-1][i];
645: as = 0.5*(t0 + ts);
646: bs = PetscPowScalar(as,bm1);
647: /* ds = bs * as; */
648: ds = PetscPowScalar(as,beta);
649: gs = coef*bs*(ts - t0);
651: /* left-hand top down corner */
652: if (k == 0) {
654: tu = x[k+1][j][i];
655: au = 0.5*(t0 + tu);
656: bu = PetscPowScalar(au,bm1);
657: /* du = bu * au; */
658: du = PetscPowScalar(au,beta);
659: gu = coef*bu*(tu - t0);
661: c[0].k = k; c[0].j = j-1; c[0].i = i;
662: v[0] = -hzhxdhy*(ds - gs);
663: c[1].k = k; c[1].j = j; c[1].i = i;
664: v[1] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
665: c[2].k = k; c[2].j = j; c[2].i = i+1;
666: v[2] = -hyhzdhx*(de + ge);
667: c[3].k = k+1; c[3].j = j; c[3].i = i;
668: v[3] = -hxhydhz*(du + gu);
669: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
671: /* left-hand top interior edge */
672: } else if (k < mz-1) {
674: tu = x[k+1][j][i];
675: au = 0.5*(t0 + tu);
676: bu = PetscPowScalar(au,bm1);
677: /* du = bu * au; */
678: du = PetscPowScalar(au,beta);
679: gu = coef*bu*(tu - t0);
681: td = x[k-1][j][i];
682: ad = 0.5*(t0 + td);
683: bd = PetscPowScalar(ad,bm1);
684: /* dd = bd * ad; */
685: dd = PetscPowScalar(ad,beta);
686: gd = coef*bd*(td - t0);
688: c[0].k = k-1; c[0].j = j; c[0].i = i;
689: v[0] = -hxhydhz*(dd - gd);
690: c[1].k = k; c[1].j = j-1; c[1].i = i;
691: v[1] = -hzhxdhy*(ds - gs);
692: c[2].k = k; c[2].j = j; c[2].i = i;
693: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
694: c[3].k = k; c[3].j = j; c[3].i = i+1;
695: v[3] = -hyhzdhx*(de + ge);
696: c[4].k = k+1; c[4].j = j; c[4].i = i;
697: v[4] = -hxhydhz*(du + gu);
698: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
700: /* left-hand top up corner */
701: } else {
703: td = x[k-1][j][i];
704: ad = 0.5*(t0 + td);
705: bd = PetscPowScalar(ad,bm1);
706: /* dd = bd * ad; */
707: dd = PetscPowScalar(ad,beta);
708: gd = coef*bd*(td - t0);
710: c[0].k = k-1; c[0].j = j; c[0].i = i;
711: v[0] = -hxhydhz*(dd - gd);
712: c[1].k = k; c[1].j = j-1; c[1].i = i;
713: v[1] = -hzhxdhy*(ds - gs);
714: c[2].k = k; c[2].j = j; c[2].i = i;
715: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
716: c[3].k = k; c[3].j = j; c[3].i = i+1;
717: v[3] = -hyhzdhx*(de + ge);
718: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
719: }
721: } else {
723: ts = x[k][j-1][i];
724: as = 0.5*(t0 + ts);
725: bs = PetscPowScalar(as,bm1);
726: /* ds = bs * as; */
727: ds = PetscPowScalar(as,beta);
728: gs = coef*bs*(t0 - ts);
730: tn = x[k][j+1][i];
731: an = 0.5*(t0 + tn);
732: bn = PetscPowScalar(an,bm1);
733: /* dn = bn * an; */
734: dn = PetscPowScalar(an,beta);
735: gn = coef*bn*(tn - t0);
737: /* left-hand down interior edge */
738: if (k == 0) {
740: tu = x[k+1][j][i];
741: au = 0.5*(t0 + tu);
742: bu = PetscPowScalar(au,bm1);
743: /* du = bu * au; */
744: du = PetscPowScalar(au,beta);
745: gu = coef*bu*(tu - t0);
747: c[0].k = k; c[0].j = j-1; c[0].i = i;
748: v[0] = -hzhxdhy*(ds - gs);
749: c[1].k = k; c[1].j = j; c[1].i = i;
750: v[1] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
751: c[2].k = k; c[2].j = j; c[2].i = i+1;
752: v[2] = -hyhzdhx*(de + ge);
753: c[3].k = k; c[3].j = j+1; c[3].i = i;
754: v[3] = -hzhxdhy*(dn + gn);
755: c[4].k = k+1; c[4].j = j; c[4].i = i;
756: v[4] = -hxhydhz*(du + gu);
757: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
759: } else if (k == mz-1) { /* left-hand up interior edge */
761: td = x[k-1][j][i];
762: ad = 0.5*(t0 + td);
763: bd = PetscPowScalar(ad,bm1);
764: /* dd = bd * ad; */
765: dd = PetscPowScalar(ad,beta);
766: gd = coef*bd*(t0 - td);
768: c[0].k = k-1; c[0].j = j; c[0].i = i;
769: v[0] = -hxhydhz*(dd - gd);
770: c[1].k = k; c[1].j = j-1; c[1].i = i;
771: v[1] = -hzhxdhy*(ds - gs);
772: c[2].k = k; c[2].j = j; c[2].i = i;
773: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
774: c[3].k = k; c[3].j = j; c[3].i = i+1;
775: v[3] = -hyhzdhx*(de + ge);
776: c[4].k = k; c[4].j = j+1; c[4].i = i;
777: v[4] = -hzhxdhy*(dn + gn);
778: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
779: } else { /* left-hand interior plane */
781: td = x[k-1][j][i];
782: ad = 0.5*(t0 + td);
783: bd = PetscPowScalar(ad,bm1);
784: /* dd = bd * ad; */
785: dd = PetscPowScalar(ad,beta);
786: gd = coef*bd*(t0 - td);
788: tu = x[k+1][j][i];
789: au = 0.5*(t0 + tu);
790: bu = PetscPowScalar(au,bm1);
791: /* du = bu * au; */
792: du = PetscPowScalar(au,beta);
793: gu = coef*bu*(tu - t0);
795: c[0].k = k-1; c[0].j = j; c[0].i = i;
796: v[0] = -hxhydhz*(dd - gd);
797: c[1].k = k; c[1].j = j-1; c[1].i = i;
798: v[1] = -hzhxdhy*(ds - gs);
799: c[2].k = k; c[2].j = j; c[2].i = i;
800: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
801: c[3].k = k; c[3].j = j; c[3].i = i+1;
802: v[3] = -hyhzdhx*(de + ge);
803: c[4].k = k; c[4].j = j+1; c[4].i = i;
804: v[4] = -hzhxdhy*(dn + gn);
805: c[5].k = k+1; c[5].j = j; c[5].i = i;
806: v[5] = -hxhydhz*(du + gu);
807: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
808: }
809: }
811: } else if (i == mx-1) {
813: /* right-hand plane boundary */
814: tw = x[k][j][i-1];
815: aw = 0.5*(t0 + tw);
816: bw = PetscPowScalar(aw,bm1);
817: /* dw = bw * aw */
818: dw = PetscPowScalar(aw,beta);
819: gw = coef*bw*(t0 - tw);
821: te = tright;
822: ae = 0.5*(t0 + te);
823: be = PetscPowScalar(ae,bm1);
824: /* de = be * ae; */
825: de = PetscPowScalar(ae,beta);
826: ge = coef*be*(te - t0);
828: /* right-hand bottom edge */
829: if (j == 0) {
831: tn = x[k][j+1][i];
832: an = 0.5*(t0 + tn);
833: bn = PetscPowScalar(an,bm1);
834: /* dn = bn * an; */
835: dn = PetscPowScalar(an,beta);
836: gn = coef*bn*(tn - t0);
838: /* right-hand bottom down corner */
839: if (k == 0) {
841: tu = x[k+1][j][i];
842: au = 0.5*(t0 + tu);
843: bu = PetscPowScalar(au,bm1);
844: /* du = bu * au; */
845: du = PetscPowScalar(au,beta);
846: gu = coef*bu*(tu - t0);
848: c[0].k = k; c[0].j = j; c[0].i = i-1;
849: v[0] = -hyhzdhx*(dw - gw);
850: c[1].k = k; c[1].j = j; c[1].i = i;
851: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
852: c[2].k = k; c[2].j = j+1; c[2].i = i;
853: v[2] = -hzhxdhy*(dn + gn);
854: c[3].k = k+1; c[3].j = j; c[3].i = i;
855: v[3] = -hxhydhz*(du + gu);
856: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
858: /* right-hand bottom interior edge */
859: } else if (k < mz-1) {
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: td = x[k-1][j][i];
869: ad = 0.5*(t0 + td);
870: bd = PetscPowScalar(ad,bm1);
871: /* dd = bd * ad; */
872: dd = PetscPowScalar(ad,beta);
873: gd = coef*bd*(td - t0);
875: c[0].k = k-1; c[0].j = j; c[0].i = i;
876: v[0] = -hxhydhz*(dd - gd);
877: c[1].k = k; c[1].j = j; c[1].i = i-1;
878: v[1] = -hyhzdhx*(dw - gw);
879: c[2].k = k; c[2].j = j; c[2].i = i;
880: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
881: c[3].k = k; c[3].j = j+1; c[3].i = i;
882: v[3] = -hzhxdhy*(dn + gn);
883: c[4].k = k+1; c[4].j = j; c[4].i = i;
884: v[4] = -hxhydhz*(du + gu);
885: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
887: /* right-hand bottom up corner */
888: } else {
890: td = x[k-1][j][i];
891: ad = 0.5*(t0 + td);
892: bd = PetscPowScalar(ad,bm1);
893: /* dd = bd * ad; */
894: dd = PetscPowScalar(ad,beta);
895: gd = coef*bd*(td - t0);
897: c[0].k = k-1; c[0].j = j; c[0].i = i;
898: v[0] = -hxhydhz*(dd - gd);
899: c[1].k = k; c[1].j = j; c[1].i = i-1;
900: v[1] = -hyhzdhx*(dw - gw);
901: c[2].k = k; c[2].j = j; c[2].i = i;
902: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
903: c[3].k = k; c[3].j = j+1; c[3].i = i;
904: v[3] = -hzhxdhy*(dn + gn);
905: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
906: }
908: /* right-hand top edge */
909: } else if (j == my-1) {
911: ts = x[k][j-1][i];
912: as = 0.5*(t0 + ts);
913: bs = PetscPowScalar(as,bm1);
914: /* ds = bs * as; */
915: ds = PetscPowScalar(as,beta);
916: gs = coef*bs*(ts - t0);
918: /* right-hand top down corner */
919: if (k == 0) {
921: tu = x[k+1][j][i];
922: au = 0.5*(t0 + tu);
923: bu = PetscPowScalar(au,bm1);
924: /* du = bu * au; */
925: du = PetscPowScalar(au,beta);
926: gu = coef*bu*(tu - t0);
928: c[0].k = k; c[0].j = j-1; c[0].i = i;
929: v[0] = -hzhxdhy*(ds - gs);
930: c[1].k = k; c[1].j = j; c[1].i = i-1;
931: v[1] = -hyhzdhx*(dw - gw);
932: c[2].k = k; c[2].j = j; c[2].i = i;
933: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
934: c[3].k = k+1; c[3].j = j; c[3].i = i;
935: v[3] = -hxhydhz*(du + gu);
936: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
938: /* right-hand top interior edge */
939: } else if (k < mz-1) {
941: tu = x[k+1][j][i];
942: au = 0.5*(t0 + tu);
943: bu = PetscPowScalar(au,bm1);
944: /* du = bu * au; */
945: du = PetscPowScalar(au,beta);
946: gu = coef*bu*(tu - t0);
948: td = x[k-1][j][i];
949: ad = 0.5*(t0 + td);
950: bd = PetscPowScalar(ad,bm1);
951: /* dd = bd * ad; */
952: dd = PetscPowScalar(ad,beta);
953: gd = coef*bd*(td - t0);
955: c[0].k = k-1; c[0].j = j; c[0].i = i;
956: v[0] = -hxhydhz*(dd - gd);
957: c[1].k = k; c[1].j = j-1; c[1].i = i;
958: v[1] = -hzhxdhy*(ds - gs);
959: c[2].k = k; c[2].j = j; c[2].i = i-1;
960: v[2] = -hyhzdhx*(dw - gw);
961: c[3].k = k; c[3].j = j; c[3].i = i;
962: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
963: c[4].k = k+1; c[4].j = j; c[4].i = i;
964: v[4] = -hxhydhz*(du + gu);
965: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
967: /* right-hand top up corner */
968: } else {
970: td = x[k-1][j][i];
971: ad = 0.5*(t0 + td);
972: bd = PetscPowScalar(ad,bm1);
973: /* dd = bd * ad; */
974: dd = PetscPowScalar(ad,beta);
975: gd = coef*bd*(td - t0);
977: c[0].k = k-1; c[0].j = j; c[0].i = i;
978: v[0] = -hxhydhz*(dd - gd);
979: c[1].k = k; c[1].j = j-1; c[1].i = i;
980: v[1] = -hzhxdhy*(ds - gs);
981: c[2].k = k; c[2].j = j; c[2].i = i-1;
982: v[2] = -hyhzdhx*(dw - gw);
983: c[3].k = k; c[3].j = j; c[3].i = i;
984: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
985: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
986: }
988: } else {
990: ts = x[k][j-1][i];
991: as = 0.5*(t0 + ts);
992: bs = PetscPowScalar(as,bm1);
993: /* ds = bs * as; */
994: ds = PetscPowScalar(as,beta);
995: gs = coef*bs*(t0 - ts);
997: tn = x[k][j+1][i];
998: an = 0.5*(t0 + tn);
999: bn = PetscPowScalar(an,bm1);
1000: /* dn = bn * an; */
1001: dn = PetscPowScalar(an,beta);
1002: gn = coef*bn*(tn - t0);
1004: /* right-hand down interior edge */
1005: if (k == 0) {
1007: tu = x[k+1][j][i];
1008: au = 0.5*(t0 + tu);
1009: bu = PetscPowScalar(au,bm1);
1010: /* du = bu * au; */
1011: du = PetscPowScalar(au,beta);
1012: gu = coef*bu*(tu - t0);
1014: c[0].k = k; c[0].j = j-1; c[0].i = i;
1015: v[0] = -hzhxdhy*(ds - gs);
1016: c[1].k = k; c[1].j = j; c[1].i = i-1;
1017: v[1] = -hyhzdhx*(dw - gw);
1018: c[2].k = k; c[2].j = j; c[2].i = i;
1019: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1020: c[3].k = k; c[3].j = j+1; c[3].i = i;
1021: v[3] = -hzhxdhy*(dn + gn);
1022: c[4].k = k+1; c[4].j = j; c[4].i = i;
1023: v[4] = -hxhydhz*(du + gu);
1024: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1026: } else if (k == mz-1) { /* right-hand up interior edge */
1028: td = x[k-1][j][i];
1029: ad = 0.5*(t0 + td);
1030: bd = PetscPowScalar(ad,bm1);
1031: /* dd = bd * ad; */
1032: dd = PetscPowScalar(ad,beta);
1033: gd = coef*bd*(t0 - td);
1035: c[0].k = k-1; c[0].j = j; c[0].i = i;
1036: v[0] = -hxhydhz*(dd - gd);
1037: c[1].k = k; c[1].j = j-1; c[1].i = i;
1038: v[1] = -hzhxdhy*(ds - gs);
1039: c[2].k = k; c[2].j = j; c[2].i = i-1;
1040: v[2] = -hyhzdhx*(dw - gw);
1041: c[3].k = k; c[3].j = j; c[3].i = i;
1042: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1043: c[4].k = k; c[4].j = j+1; c[4].i = i;
1044: v[4] = -hzhxdhy*(dn + gn);
1045: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1047: } else { /* right-hand interior plane */
1049: td = x[k-1][j][i];
1050: ad = 0.5*(t0 + td);
1051: bd = PetscPowScalar(ad,bm1);
1052: /* dd = bd * ad; */
1053: dd = PetscPowScalar(ad,beta);
1054: gd = coef*bd*(t0 - td);
1056: tu = x[k+1][j][i];
1057: au = 0.5*(t0 + tu);
1058: bu = PetscPowScalar(au,bm1);
1059: /* du = bu * au; */
1060: du = PetscPowScalar(au,beta);
1061: gu = coef*bu*(tu - t0);
1063: c[0].k = k-1; c[0].j = j; c[0].i = i;
1064: v[0] = -hxhydhz*(dd - gd);
1065: c[1].k = k; c[1].j = j-1; c[1].i = i;
1066: v[1] = -hzhxdhy*(ds - gs);
1067: c[2].k = k; c[2].j = j; c[2].i = i-1;
1068: v[2] = -hyhzdhx*(dw - gw);
1069: c[3].k = k; c[3].j = j; c[3].i = i;
1070: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1071: c[4].k = k; c[4].j = j+1; c[4].i = i;
1072: v[4] = -hzhxdhy*(dn + gn);
1073: c[5].k = k+1; c[5].j = j; c[5].i = i;
1074: v[5] = -hxhydhz*(du + gu);
1075: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1076: }
1077: }
1079: } else if (j == 0) {
1081: tw = x[k][j][i-1];
1082: aw = 0.5*(t0 + tw);
1083: bw = PetscPowScalar(aw,bm1);
1084: /* dw = bw * aw */
1085: dw = PetscPowScalar(aw,beta);
1086: gw = coef*bw*(t0 - tw);
1088: te = x[k][j][i+1];
1089: ae = 0.5*(t0 + te);
1090: be = PetscPowScalar(ae,bm1);
1091: /* de = be * ae; */
1092: de = PetscPowScalar(ae,beta);
1093: ge = coef*be*(te - t0);
1095: tn = x[k][j+1][i];
1096: an = 0.5*(t0 + tn);
1097: bn = PetscPowScalar(an,bm1);
1098: /* dn = bn * an; */
1099: dn = PetscPowScalar(an,beta);
1100: gn = coef*bn*(tn - t0);
1102: /* bottom down interior edge */
1103: if (k == 0) {
1105: tu = x[k+1][j][i];
1106: au = 0.5*(t0 + tu);
1107: bu = PetscPowScalar(au,bm1);
1108: /* du = bu * au; */
1109: du = PetscPowScalar(au,beta);
1110: gu = coef*bu*(tu - t0);
1112: c[0].k = k; c[0].j = j; c[0].i = i-1;
1113: v[0] = -hyhzdhx*(dw - gw);
1114: c[1].k = k; c[1].j = j; c[1].i = i;
1115: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1116: c[2].k = k; c[2].j = j; c[2].i = i+1;
1117: v[2] = -hyhzdhx*(de + ge);
1118: c[3].k = k; c[3].j = j+1; c[3].i = i;
1119: v[3] = -hzhxdhy*(dn + gn);
1120: c[4].k = k+1; c[4].j = j; c[4].i = i;
1121: v[4] = -hxhydhz*(du + gu);
1122: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1124: } else if (k == mz-1) { /* bottom up interior edge */
1126: td = x[k-1][j][i];
1127: ad = 0.5*(t0 + td);
1128: bd = PetscPowScalar(ad,bm1);
1129: /* dd = bd * ad; */
1130: dd = PetscPowScalar(ad,beta);
1131: gd = coef*bd*(td - t0);
1133: c[0].k = k-1; c[0].j = j; c[0].i = i;
1134: v[0] = -hxhydhz*(dd - gd);
1135: c[1].k = k; c[1].j = j; c[1].i = i-1;
1136: v[1] = -hyhzdhx*(dw - gw);
1137: c[2].k = k; c[2].j = j; c[2].i = i;
1138: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1139: c[3].k = k; c[3].j = j; c[3].i = i+1;
1140: v[3] = -hyhzdhx*(de + ge);
1141: c[4].k = k; c[4].j = j+1; c[4].i = i;
1142: v[4] = -hzhxdhy*(dn + gn);
1143: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1145: } else { /* bottom interior plane */
1147: tu = x[k+1][j][i];
1148: au = 0.5*(t0 + tu);
1149: bu = PetscPowScalar(au,bm1);
1150: /* du = bu * au; */
1151: du = PetscPowScalar(au,beta);
1152: gu = coef*bu*(tu - t0);
1154: td = x[k-1][j][i];
1155: ad = 0.5*(t0 + td);
1156: bd = PetscPowScalar(ad,bm1);
1157: /* dd = bd * ad; */
1158: dd = PetscPowScalar(ad,beta);
1159: gd = coef*bd*(td - t0);
1161: c[0].k = k-1; c[0].j = j; c[0].i = i;
1162: v[0] = -hxhydhz*(dd - gd);
1163: c[1].k = k; c[1].j = j; c[1].i = i-1;
1164: v[1] = -hyhzdhx*(dw - gw);
1165: c[2].k = k; c[2].j = j; c[2].i = i;
1166: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1167: c[3].k = k; c[3].j = j; c[3].i = i+1;
1168: v[3] = -hyhzdhx*(de + ge);
1169: c[4].k = k; c[4].j = j+1; c[4].i = i;
1170: v[4] = -hzhxdhy*(dn + gn);
1171: c[5].k = k+1; c[5].j = j; c[5].i = i;
1172: v[5] = -hxhydhz*(du + gu);
1173: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1174: }
1176: } else if (j == my-1) {
1178: tw = x[k][j][i-1];
1179: aw = 0.5*(t0 + tw);
1180: bw = PetscPowScalar(aw,bm1);
1181: /* dw = bw * aw */
1182: dw = PetscPowScalar(aw,beta);
1183: gw = coef*bw*(t0 - tw);
1185: te = x[k][j][i+1];
1186: ae = 0.5*(t0 + te);
1187: be = PetscPowScalar(ae,bm1);
1188: /* de = be * ae; */
1189: de = PetscPowScalar(ae,beta);
1190: ge = coef*be*(te - t0);
1192: ts = x[k][j-1][i];
1193: as = 0.5*(t0 + ts);
1194: bs = PetscPowScalar(as,bm1);
1195: /* ds = bs * as; */
1196: ds = PetscPowScalar(as,beta);
1197: gs = coef*bs*(t0 - ts);
1199: /* top down interior edge */
1200: if (k == 0) {
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; c[0].j = j-1; c[0].i = i;
1210: v[0] = -hzhxdhy*(ds - gs);
1211: c[1].k = k; c[1].j = j; c[1].i = i-1;
1212: v[1] = -hyhzdhx*(dw - gw);
1213: c[2].k = k; c[2].j = j; c[2].i = i;
1214: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1215: c[3].k = k; c[3].j = j; c[3].i = i+1;
1216: v[3] = -hyhzdhx*(de + ge);
1217: c[4].k = k+1; c[4].j = j; c[4].i = i;
1218: v[4] = -hxhydhz*(du + gu);
1219: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1221: } else if (k == mz-1) { /* top up interior edge */
1223: td = x[k-1][j][i];
1224: ad = 0.5*(t0 + td);
1225: bd = PetscPowScalar(ad,bm1);
1226: /* dd = bd * ad; */
1227: dd = PetscPowScalar(ad,beta);
1228: gd = coef*bd*(td - t0);
1230: c[0].k = k-1; c[0].j = j; c[0].i = i;
1231: v[0] = -hxhydhz*(dd - gd);
1232: c[1].k = k; c[1].j = j-1; c[1].i = i;
1233: v[1] = -hzhxdhy*(ds - gs);
1234: c[2].k = k; c[2].j = j; c[2].i = i-1;
1235: v[2] = -hyhzdhx*(dw - gw);
1236: c[3].k = k; c[3].j = j; c[3].i = i;
1237: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1238: c[4].k = k; c[4].j = j; c[4].i = i+1;
1239: v[4] = -hyhzdhx*(de + ge);
1240: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1242: } else { /* top interior plane */
1244: tu = x[k+1][j][i];
1245: au = 0.5*(t0 + tu);
1246: bu = PetscPowScalar(au,bm1);
1247: /* du = bu * au; */
1248: du = PetscPowScalar(au,beta);
1249: gu = coef*bu*(tu - t0);
1251: td = x[k-1][j][i];
1252: ad = 0.5*(t0 + td);
1253: bd = PetscPowScalar(ad,bm1);
1254: /* dd = bd * ad; */
1255: dd = PetscPowScalar(ad,beta);
1256: gd = coef*bd*(td - t0);
1258: c[0].k = k-1; c[0].j = j; c[0].i = i;
1259: v[0] = -hxhydhz*(dd - gd);
1260: c[1].k = k; c[1].j = j-1; c[1].i = i;
1261: v[1] = -hzhxdhy*(ds - gs);
1262: c[2].k = k; c[2].j = j; c[2].i = i-1;
1263: v[2] = -hyhzdhx*(dw - gw);
1264: c[3].k = k; c[3].j = j; c[3].i = i;
1265: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1266: c[4].k = k; c[4].j = j; c[4].i = i+1;
1267: v[4] = -hyhzdhx*(de + ge);
1268: c[5].k = k+1; c[5].j = j; c[5].i = i;
1269: v[5] = -hxhydhz*(du + gu);
1270: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1271: }
1273: } else if (k == 0) {
1275: /* down interior plane */
1277: tw = x[k][j][i-1];
1278: aw = 0.5*(t0 + tw);
1279: bw = PetscPowScalar(aw,bm1);
1280: /* dw = bw * aw */
1281: dw = PetscPowScalar(aw,beta);
1282: gw = coef*bw*(t0 - tw);
1284: te = x[k][j][i+1];
1285: ae = 0.5*(t0 + te);
1286: be = PetscPowScalar(ae,bm1);
1287: /* de = be * ae; */
1288: de = PetscPowScalar(ae,beta);
1289: ge = coef*be*(te - t0);
1291: ts = x[k][j-1][i];
1292: as = 0.5*(t0 + ts);
1293: bs = PetscPowScalar(as,bm1);
1294: /* ds = bs * as; */
1295: ds = PetscPowScalar(as,beta);
1296: gs = coef*bs*(t0 - ts);
1298: tn = x[k][j+1][i];
1299: an = 0.5*(t0 + tn);
1300: bn = PetscPowScalar(an,bm1);
1301: /* dn = bn * an; */
1302: dn = PetscPowScalar(an,beta);
1303: gn = coef*bn*(tn - t0);
1305: tu = x[k+1][j][i];
1306: au = 0.5*(t0 + tu);
1307: bu = PetscPowScalar(au,bm1);
1308: /* du = bu * au; */
1309: du = PetscPowScalar(au,beta);
1310: gu = coef*bu*(tu - t0);
1312: c[0].k = k; c[0].j = j-1; c[0].i = i;
1313: v[0] = -hzhxdhy*(ds - gs);
1314: c[1].k = k; c[1].j = j; c[1].i = i-1;
1315: v[1] = -hyhzdhx*(dw - gw);
1316: c[2].k = k; c[2].j = j; c[2].i = i;
1317: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1318: c[3].k = k; c[3].j = j; c[3].i = i+1;
1319: v[3] = -hyhzdhx*(de + ge);
1320: c[4].k = k; c[4].j = j+1; c[4].i = i;
1321: v[4] = -hzhxdhy*(dn + gn);
1322: c[5].k = k+1; c[5].j = j; c[5].i = i;
1323: v[5] = -hxhydhz*(du + gu);
1324: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1326: } else if (k == mz-1) {
1328: /* up interior plane */
1330: tw = x[k][j][i-1];
1331: aw = 0.5*(t0 + tw);
1332: bw = PetscPowScalar(aw,bm1);
1333: /* dw = bw * aw */
1334: dw = PetscPowScalar(aw,beta);
1335: gw = coef*bw*(t0 - tw);
1337: te = x[k][j][i+1];
1338: ae = 0.5*(t0 + te);
1339: be = PetscPowScalar(ae,bm1);
1340: /* de = be * ae; */
1341: de = PetscPowScalar(ae,beta);
1342: ge = coef*be*(te - t0);
1344: ts = x[k][j-1][i];
1345: as = 0.5*(t0 + ts);
1346: bs = PetscPowScalar(as,bm1);
1347: /* ds = bs * as; */
1348: ds = PetscPowScalar(as,beta);
1349: gs = coef*bs*(t0 - ts);
1351: tn = x[k][j+1][i];
1352: an = 0.5*(t0 + tn);
1353: bn = PetscPowScalar(an,bm1);
1354: /* dn = bn * an; */
1355: dn = PetscPowScalar(an,beta);
1356: gn = coef*bn*(tn - t0);
1358: td = x[k-1][j][i];
1359: ad = 0.5*(t0 + td);
1360: bd = PetscPowScalar(ad,bm1);
1361: /* dd = bd * ad; */
1362: dd = PetscPowScalar(ad,beta);
1363: gd = coef*bd*(t0 - td);
1365: c[0].k = k-1; c[0].j = j; c[0].i = i;
1366: v[0] = -hxhydhz*(dd - gd);
1367: c[1].k = k; c[1].j = j-1; c[1].i = i;
1368: v[1] = -hzhxdhy*(ds - gs);
1369: c[2].k = k; c[2].j = j; c[2].i = i-1;
1370: v[2] = -hyhzdhx*(dw - gw);
1371: c[3].k = k; c[3].j = j; c[3].i = i;
1372: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1373: c[4].k = k; c[4].j = j; c[4].i = i+1;
1374: v[4] = -hyhzdhx*(de + ge);
1375: c[5].k = k; c[5].j = j+1; c[5].i = i;
1376: v[5] = -hzhxdhy*(dn + gn);
1377: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1378: }
1379: }
1380: }
1381: }
1382: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
1383: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
1384: if (jac != J) {
1385: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1386: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1387: }
1388: DMDAVecRestoreArray(da,localX,&x);
1389: DMRestoreLocalVector(da,&localX);
1391: PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
1392: return 0;
1393: }
1395: /*TEST
1397: test:
1398: nsize: 4
1399: args: -snes_monitor_short -pc_mg_type full -ksp_type fgmres -pc_type mg -snes_view -pc_mg_levels 2 -pc_mg_galerkin pmat
1400: requires: !single
1402: TEST*/