Actual source code: ex20.c
petsc-3.11.4 2019-09-28
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*/
23: /*
25: This example models the partial differential equation
27: - Div(alpha* T^beta (GRAD T)) = 0.
29: where beta = 2.5 and alpha = 1.0
31: BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = dT/dn_up = dT/dn_down = 0.
33: in the unit square, which is uniformly discretized in each of x and
34: y in this simple encoding. The degrees of freedom are cell centered.
36: A finite volume approximation with the usual 7-point stencil
37: is used to discretize the boundary value problem to obtain a
38: nonlinear system of equations.
40: This code was contributed by Nickolas Jovanovic based on ex18.c
42: */
44: #include <petscsnes.h>
45: #include <petscdm.h>
46: #include <petscdmda.h>
48: /* User-defined application context */
50: typedef struct {
51: PetscReal tleft,tright; /* Dirichlet boundary conditions */
52: PetscReal beta,bm1,coef; /* nonlinear diffusivity parameterizations */
53: } AppCtx;
55: #define POWFLOP 5 /* assume a pow() takes five flops */
57: extern PetscErrorCode FormInitialGuess(SNES,Vec,void*);
58: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
59: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
61: int main(int argc,char **argv)
62: {
63: SNES snes;
64: AppCtx user;
66: PetscInt its,lits;
67: PetscReal litspit;
68: DM da;
70: PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
71: /* set problem parameters */
72: user.tleft = 1.0;
73: user.tright = 0.1;
74: user.beta = 2.5;
75: PetscOptionsGetReal(NULL,NULL,"-tleft",&user.tleft,NULL);
76: PetscOptionsGetReal(NULL,NULL,"-tright",&user.tright,NULL);
77: PetscOptionsGetReal(NULL,NULL,"-beta",&user.beta,NULL);
78: user.bm1 = user.beta - 1.0;
79: user.coef = user.beta/2.0;
81: /*
82: Set the DMDA (grid structure) for the grids.
83: */
84: 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);
85: DMSetFromOptions(da);
86: DMSetUp(da);
87: DMSetApplicationContext(da,&user);
89: /*
90: Create the nonlinear solver
91: */
92: SNESCreate(PETSC_COMM_WORLD,&snes);
93: SNESSetDM(snes,da);
94: SNESSetFunction(snes,NULL,FormFunction,&user);
95: SNESSetJacobian(snes,NULL,NULL,FormJacobian,&user);
96: SNESSetFromOptions(snes);
97: SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);
99: SNESSolve(snes,NULL,NULL);
100: SNESGetIterationNumber(snes,&its);
101: SNESGetLinearSolveIterations(snes,&lits);
102: litspit = ((PetscReal)lits)/((PetscReal)its);
103: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
104: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
105: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",(double)litspit);
107: SNESDestroy(&snes);
108: DMDestroy(&da);
109: PetscFinalize();
110: return ierr;
111: }
112: /* -------------------- Form initial approximation ----------------- */
113: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
114: {
115: AppCtx *user;
116: PetscInt i,j,k,xs,ys,xm,ym,zs,zm;
118: PetscScalar ***x;
119: DM da;
122: SNESGetDM(snes,&da);
123: DMGetApplicationContext(da,&user);
124: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
125: DMDAVecGetArray(da,X,&x);
127: /* Compute initial guess */
128: for (k=zs; k<zs+zm; k++) {
129: for (j=ys; j<ys+ym; j++) {
130: for (i=xs; i<xs+xm; i++) {
131: x[k][j][i] = user->tleft;
132: }
133: }
134: }
135: DMDAVecRestoreArray(da,X,&x);
136: return(0);
137: }
138: /* -------------------- Evaluate Function F(x) --------------------- */
139: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
140: {
141: AppCtx *user = (AppCtx*)ptr;
143: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
144: PetscScalar zero = 0.0,one = 1.0;
145: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
146: 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;
147: PetscScalar tleft,tright,beta,td,ad,dd,fd=0.0,tu,au,du=0.0,fu=0.0;
148: PetscScalar ***x,***f;
149: Vec localX;
150: DM da;
153: SNESGetDM(snes,&da);
154: DMGetLocalVector(da,&localX);
155: DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
156: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
157: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
158: tleft = user->tleft; tright = user->tright;
159: beta = user->beta;
161: /* Get ghost points */
162: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
163: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
164: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
165: DMDAVecGetArray(da,localX,&x);
166: DMDAVecGetArray(da,F,&f);
168: /* Evaluate function */
169: for (k=zs; k<zs+zm; k++) {
170: for (j=ys; j<ys+ym; j++) {
171: for (i=xs; i<xs+xm; i++) {
172: t0 = x[k][j][i];
174: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
176: /* general interior volume */
178: tw = x[k][j][i-1];
179: aw = 0.5*(t0 + tw);
180: dw = PetscPowScalar(aw,beta);
181: fw = dw*(t0 - tw);
183: te = x[k][j][i+1];
184: ae = 0.5*(t0 + te);
185: de = PetscPowScalar(ae,beta);
186: fe = de*(te - t0);
188: ts = x[k][j-1][i];
189: as = 0.5*(t0 + ts);
190: ds = PetscPowScalar(as,beta);
191: fs = ds*(t0 - ts);
193: tn = x[k][j+1][i];
194: an = 0.5*(t0 + tn);
195: dn = PetscPowScalar(an,beta);
196: fn = dn*(tn - t0);
198: td = x[k-1][j][i];
199: ad = 0.5*(t0 + td);
200: dd = PetscPowScalar(ad,beta);
201: fd = dd*(t0 - td);
203: tu = x[k+1][j][i];
204: au = 0.5*(t0 + tu);
205: du = PetscPowScalar(au,beta);
206: fu = du*(tu - t0);
208: } else if (i == 0) {
210: /* left-hand (west) boundary */
211: tw = tleft;
212: aw = 0.5*(t0 + tw);
213: dw = PetscPowScalar(aw,beta);
214: fw = dw*(t0 - tw);
216: te = x[k][j][i+1];
217: ae = 0.5*(t0 + te);
218: de = PetscPowScalar(ae,beta);
219: fe = de*(te - t0);
221: if (j > 0) {
222: ts = x[k][j-1][i];
223: as = 0.5*(t0 + ts);
224: ds = PetscPowScalar(as,beta);
225: fs = ds*(t0 - ts);
226: } else {
227: fs = zero;
228: }
230: if (j < my-1) {
231: tn = x[k][j+1][i];
232: an = 0.5*(t0 + tn);
233: dn = PetscPowScalar(an,beta);
234: fn = dn*(tn - t0);
235: } else {
236: fn = zero;
237: }
239: if (k > 0) {
240: td = x[k-1][j][i];
241: ad = 0.5*(t0 + td);
242: dd = PetscPowScalar(ad,beta);
243: fd = dd*(t0 - td);
244: } else {
245: fd = zero;
246: }
248: if (k < mz-1) {
249: tu = x[k+1][j][i];
250: au = 0.5*(t0 + tu);
251: du = PetscPowScalar(au,beta);
252: fu = du*(tu - t0);
253: } else {
254: fu = zero;
255: }
257: } else if (i == mx-1) {
259: /* right-hand (east) boundary */
260: tw = x[k][j][i-1];
261: aw = 0.5*(t0 + tw);
262: dw = PetscPowScalar(aw,beta);
263: fw = dw*(t0 - tw);
265: te = tright;
266: ae = 0.5*(t0 + te);
267: de = PetscPowScalar(ae,beta);
268: fe = de*(te - t0);
270: if (j > 0) {
271: ts = x[k][j-1][i];
272: as = 0.5*(t0 + ts);
273: ds = PetscPowScalar(as,beta);
274: fs = ds*(t0 - ts);
275: } else {
276: fs = zero;
277: }
279: if (j < my-1) {
280: tn = x[k][j+1][i];
281: an = 0.5*(t0 + tn);
282: dn = PetscPowScalar(an,beta);
283: fn = dn*(tn - t0);
284: } else {
285: fn = zero;
286: }
288: if (k > 0) {
289: td = x[k-1][j][i];
290: ad = 0.5*(t0 + td);
291: dd = PetscPowScalar(ad,beta);
292: fd = dd*(t0 - td);
293: } else {
294: fd = zero;
295: }
297: if (k < mz-1) {
298: tu = x[k+1][j][i];
299: au = 0.5*(t0 + tu);
300: du = PetscPowScalar(au,beta);
301: fu = du*(tu - t0);
302: } else {
303: fu = zero;
304: }
306: } else if (j == 0) {
308: /* bottom (south) boundary, and i <> 0 or mx-1 */
309: tw = x[k][j][i-1];
310: aw = 0.5*(t0 + tw);
311: dw = PetscPowScalar(aw,beta);
312: fw = dw*(t0 - tw);
314: te = x[k][j][i+1];
315: ae = 0.5*(t0 + te);
316: de = PetscPowScalar(ae,beta);
317: fe = de*(te - t0);
319: fs = zero;
321: tn = x[k][j+1][i];
322: an = 0.5*(t0 + tn);
323: dn = PetscPowScalar(an,beta);
324: fn = dn*(tn - t0);
326: if (k > 0) {
327: td = x[k-1][j][i];
328: ad = 0.5*(t0 + td);
329: dd = PetscPowScalar(ad,beta);
330: fd = dd*(t0 - td);
331: } else {
332: fd = zero;
333: }
335: if (k < mz-1) {
336: tu = x[k+1][j][i];
337: au = 0.5*(t0 + tu);
338: du = PetscPowScalar(au,beta);
339: fu = du*(tu - t0);
340: } else {
341: fu = zero;
342: }
344: } else if (j == my-1) {
346: /* top (north) boundary, and i <> 0 or mx-1 */
347: tw = x[k][j][i-1];
348: aw = 0.5*(t0 + tw);
349: dw = PetscPowScalar(aw,beta);
350: fw = dw*(t0 - tw);
352: te = x[k][j][i+1];
353: ae = 0.5*(t0 + te);
354: de = PetscPowScalar(ae,beta);
355: fe = de*(te - t0);
357: ts = x[k][j-1][i];
358: as = 0.5*(t0 + ts);
359: ds = PetscPowScalar(as,beta);
360: fs = ds*(t0 - ts);
362: fn = zero;
364: if (k > 0) {
365: td = x[k-1][j][i];
366: ad = 0.5*(t0 + td);
367: dd = PetscPowScalar(ad,beta);
368: fd = dd*(t0 - td);
369: } else {
370: fd = zero;
371: }
373: if (k < mz-1) {
374: tu = x[k+1][j][i];
375: au = 0.5*(t0 + tu);
376: du = PetscPowScalar(au,beta);
377: fu = du*(tu - t0);
378: } else {
379: fu = zero;
380: }
382: } else if (k == 0) {
384: /* down boundary (interior only) */
385: tw = x[k][j][i-1];
386: aw = 0.5*(t0 + tw);
387: dw = PetscPowScalar(aw,beta);
388: fw = dw*(t0 - tw);
390: te = x[k][j][i+1];
391: ae = 0.5*(t0 + te);
392: de = PetscPowScalar(ae,beta);
393: fe = de*(te - t0);
395: ts = x[k][j-1][i];
396: as = 0.5*(t0 + ts);
397: ds = PetscPowScalar(as,beta);
398: fs = ds*(t0 - ts);
400: tn = x[k][j+1][i];
401: an = 0.5*(t0 + tn);
402: dn = PetscPowScalar(an,beta);
403: fn = dn*(tn - t0);
405: fd = zero;
407: tu = x[k+1][j][i];
408: au = 0.5*(t0 + tu);
409: du = PetscPowScalar(au,beta);
410: fu = du*(tu - t0);
412: } else if (k == mz-1) {
414: /* up boundary (interior only) */
415: tw = x[k][j][i-1];
416: aw = 0.5*(t0 + tw);
417: dw = PetscPowScalar(aw,beta);
418: fw = dw*(t0 - tw);
420: te = x[k][j][i+1];
421: ae = 0.5*(t0 + te);
422: de = PetscPowScalar(ae,beta);
423: fe = de*(te - t0);
425: ts = x[k][j-1][i];
426: as = 0.5*(t0 + ts);
427: ds = PetscPowScalar(as,beta);
428: fs = ds*(t0 - ts);
430: tn = x[k][j+1][i];
431: an = 0.5*(t0 + tn);
432: dn = PetscPowScalar(an,beta);
433: fn = dn*(tn - t0);
435: td = x[k-1][j][i];
436: ad = 0.5*(t0 + td);
437: dd = PetscPowScalar(ad,beta);
438: fd = dd*(t0 - td);
440: fu = zero;
441: }
443: f[k][j][i] = -hyhzdhx*(fe-fw) - hzhxdhy*(fn-fs) - hxhydhz*(fu-fd);
444: }
445: }
446: }
447: DMDAVecRestoreArray(da,localX,&x);
448: DMDAVecRestoreArray(da,F,&f);
449: DMRestoreLocalVector(da,&localX);
450: PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
451: return(0);
452: }
453: /* -------------------- Evaluate Jacobian F(x) --------------------- */
454: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
455: {
456: AppCtx *user = (AppCtx*)ptr;
458: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
459: PetscScalar one = 1.0;
460: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
461: PetscScalar t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw;
462: PetscScalar tleft,tright,beta,td,ad,dd,tu,au,du,v[7],bm1,coef;
463: PetscScalar ***x,bn,bs,be,bw,bu,bd,gn,gs,ge,gw,gu,gd;
464: Vec localX;
465: MatStencil c[7],row;
466: DM da;
469: SNESGetDM(snes,&da);
470: DMGetLocalVector(da,&localX);
471: DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
472: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
473: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
474: tleft = user->tleft; tright = user->tright;
475: beta = user->beta; bm1 = user->bm1; coef = user->coef;
477: /* Get ghost points */
478: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
479: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
480: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
481: DMDAVecGetArray(da,localX,&x);
483: /* Evaluate Jacobian of function */
484: for (k=zs; k<zs+zm; k++) {
485: for (j=ys; j<ys+ym; j++) {
486: for (i=xs; i<xs+xm; i++) {
487: t0 = x[k][j][i];
488: row.k = k; row.j = j; row.i = i;
489: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
491: /* general interior volume */
493: tw = x[k][j][i-1];
494: aw = 0.5*(t0 + tw);
495: bw = PetscPowScalar(aw,bm1);
496: /* dw = bw * aw */
497: dw = PetscPowScalar(aw,beta);
498: gw = coef*bw*(t0 - tw);
500: te = x[k][j][i+1];
501: ae = 0.5*(t0 + te);
502: be = PetscPowScalar(ae,bm1);
503: /* de = be * ae; */
504: de = PetscPowScalar(ae,beta);
505: ge = coef*be*(te - t0);
507: ts = x[k][j-1][i];
508: as = 0.5*(t0 + ts);
509: bs = PetscPowScalar(as,bm1);
510: /* ds = bs * as; */
511: ds = PetscPowScalar(as,beta);
512: gs = coef*bs*(t0 - ts);
514: tn = x[k][j+1][i];
515: an = 0.5*(t0 + tn);
516: bn = PetscPowScalar(an,bm1);
517: /* dn = bn * an; */
518: dn = PetscPowScalar(an,beta);
519: gn = coef*bn*(tn - t0);
521: td = x[k-1][j][i];
522: ad = 0.5*(t0 + td);
523: bd = PetscPowScalar(ad,bm1);
524: /* dd = bd * ad; */
525: dd = PetscPowScalar(ad,beta);
526: gd = coef*bd*(t0 - td);
528: tu = x[k+1][j][i];
529: au = 0.5*(t0 + tu);
530: bu = PetscPowScalar(au,bm1);
531: /* du = bu * au; */
532: du = PetscPowScalar(au,beta);
533: gu = coef*bu*(tu - t0);
535: c[0].k = k-1; c[0].j = j; c[0].i = i; v[0] = -hxhydhz*(dd - gd);
536: c[1].k = k; c[1].j = j-1; c[1].i = i;
537: v[1] = -hzhxdhy*(ds - gs);
538: c[2].k = k; c[2].j = j; c[2].i = i-1;
539: v[2] = -hyhzdhx*(dw - gw);
540: c[3].k = k; c[3].j = j; c[3].i = i;
541: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
542: c[4].k = k; c[4].j = j; c[4].i = i+1;
543: v[4] = -hyhzdhx*(de + ge);
544: c[5].k = k; c[5].j = j+1; c[5].i = i;
545: v[5] = -hzhxdhy*(dn + gn);
546: c[6].k = k+1; c[6].j = j; c[6].i = i;
547: v[6] = -hxhydhz*(du + gu);
548: MatSetValuesStencil(jac,1,&row,7,c,v,INSERT_VALUES);
550: } else if (i == 0) {
552: /* left-hand plane boundary */
553: tw = tleft;
554: aw = 0.5*(t0 + tw);
555: bw = PetscPowScalar(aw,bm1);
556: /* dw = bw * aw */
557: dw = PetscPowScalar(aw,beta);
558: gw = coef*bw*(t0 - tw);
560: te = x[k][j][i+1];
561: ae = 0.5*(t0 + te);
562: be = PetscPowScalar(ae,bm1);
563: /* de = be * ae; */
564: de = PetscPowScalar(ae,beta);
565: ge = coef*be*(te - t0);
567: /* left-hand bottom edge */
568: if (j == 0) {
570: tn = x[k][j+1][i];
571: an = 0.5*(t0 + tn);
572: bn = PetscPowScalar(an,bm1);
573: /* dn = bn * an; */
574: dn = PetscPowScalar(an,beta);
575: gn = coef*bn*(tn - t0);
577: /* left-hand bottom down corner */
578: if (k == 0) {
580: tu = x[k+1][j][i];
581: au = 0.5*(t0 + tu);
582: bu = PetscPowScalar(au,bm1);
583: /* du = bu * au; */
584: du = PetscPowScalar(au,beta);
585: gu = coef*bu*(tu - t0);
587: c[0].k = k; c[0].j = j; c[0].i = i;
588: v[0] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
589: c[1].k = k; c[1].j = j; c[1].i = i+1;
590: v[1] = -hyhzdhx*(de + ge);
591: c[2].k = k; c[2].j = j+1; c[2].i = i;
592: v[2] = -hzhxdhy*(dn + gn);
593: c[3].k = k+1; c[3].j = j; c[3].i = i;
594: v[3] = -hxhydhz*(du + gu);
595: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
597: /* left-hand bottom interior edge */
598: } else if (k < mz-1) {
600: tu = x[k+1][j][i];
601: au = 0.5*(t0 + tu);
602: bu = PetscPowScalar(au,bm1);
603: /* du = bu * au; */
604: du = PetscPowScalar(au,beta);
605: gu = coef*bu*(tu - t0);
607: td = x[k-1][j][i];
608: ad = 0.5*(t0 + td);
609: bd = PetscPowScalar(ad,bm1);
610: /* dd = bd * ad; */
611: dd = PetscPowScalar(ad,beta);
612: gd = coef*bd*(td - t0);
614: c[0].k = k-1; c[0].j = j; c[0].i = i;
615: v[0] = -hxhydhz*(dd - gd);
616: c[1].k = k; c[1].j = j; c[1].i = i;
617: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
618: c[2].k = k; c[2].j = j; c[2].i = i+1;
619: v[2] = -hyhzdhx*(de + ge);
620: c[3].k = k; c[3].j = j+1; c[3].i = i;
621: v[3] = -hzhxdhy*(dn + gn);
622: c[4].k = k+1; c[4].j = j; c[4].i = i;
623: v[4] = -hxhydhz*(du + gu);
624: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
626: /* left-hand bottom up corner */
627: } else {
629: td = x[k-1][j][i];
630: ad = 0.5*(t0 + td);
631: bd = PetscPowScalar(ad,bm1);
632: /* dd = bd * ad; */
633: dd = PetscPowScalar(ad,beta);
634: gd = coef*bd*(td - t0);
636: c[0].k = k-1; c[0].j = j; c[0].i = i;
637: v[0] = -hxhydhz*(dd - gd);
638: c[1].k = k; c[1].j = j; c[1].i = i;
639: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
640: c[2].k = k; c[2].j = j; c[2].i = i+1;
641: v[2] = -hyhzdhx*(de + ge);
642: c[3].k = k; c[3].j = j+1; c[3].i = i;
643: v[3] = -hzhxdhy*(dn + gn);
644: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
645: }
647: /* left-hand top edge */
648: } else if (j == my-1) {
650: ts = x[k][j-1][i];
651: as = 0.5*(t0 + ts);
652: bs = PetscPowScalar(as,bm1);
653: /* ds = bs * as; */
654: ds = PetscPowScalar(as,beta);
655: gs = coef*bs*(ts - t0);
657: /* left-hand top down corner */
658: if (k == 0) {
660: tu = x[k+1][j][i];
661: au = 0.5*(t0 + tu);
662: bu = PetscPowScalar(au,bm1);
663: /* du = bu * au; */
664: du = PetscPowScalar(au,beta);
665: gu = coef*bu*(tu - t0);
667: c[0].k = k; c[0].j = j-1; c[0].i = i;
668: v[0] = -hzhxdhy*(ds - gs);
669: c[1].k = k; c[1].j = j; c[1].i = i;
670: v[1] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
671: c[2].k = k; c[2].j = j; c[2].i = i+1;
672: v[2] = -hyhzdhx*(de + ge);
673: c[3].k = k+1; c[3].j = j; c[3].i = i;
674: v[3] = -hxhydhz*(du + gu);
675: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
677: /* left-hand top interior edge */
678: } else if (k < mz-1) {
680: tu = x[k+1][j][i];
681: au = 0.5*(t0 + tu);
682: bu = PetscPowScalar(au,bm1);
683: /* du = bu * au; */
684: du = PetscPowScalar(au,beta);
685: gu = coef*bu*(tu - t0);
687: td = x[k-1][j][i];
688: ad = 0.5*(t0 + td);
689: bd = PetscPowScalar(ad,bm1);
690: /* dd = bd * ad; */
691: dd = PetscPowScalar(ad,beta);
692: gd = coef*bd*(td - t0);
694: c[0].k = k-1; c[0].j = j; c[0].i = i;
695: v[0] = -hxhydhz*(dd - gd);
696: c[1].k = k; c[1].j = j-1; c[1].i = i;
697: v[1] = -hzhxdhy*(ds - gs);
698: c[2].k = k; c[2].j = j; c[2].i = i;
699: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
700: c[3].k = k; c[3].j = j; c[3].i = i+1;
701: v[3] = -hyhzdhx*(de + ge);
702: c[4].k = k+1; c[4].j = j; c[4].i = i;
703: v[4] = -hxhydhz*(du + gu);
704: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
706: /* left-hand top up corner */
707: } else {
709: td = x[k-1][j][i];
710: ad = 0.5*(t0 + td);
711: bd = PetscPowScalar(ad,bm1);
712: /* dd = bd * ad; */
713: dd = PetscPowScalar(ad,beta);
714: gd = coef*bd*(td - t0);
716: c[0].k = k-1; c[0].j = j; c[0].i = i;
717: v[0] = -hxhydhz*(dd - gd);
718: c[1].k = k; c[1].j = j-1; c[1].i = i;
719: v[1] = -hzhxdhy*(ds - gs);
720: c[2].k = k; c[2].j = j; c[2].i = i;
721: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
722: c[3].k = k; c[3].j = j; c[3].i = i+1;
723: v[3] = -hyhzdhx*(de + ge);
724: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
725: }
727: } else {
729: ts = x[k][j-1][i];
730: as = 0.5*(t0 + ts);
731: bs = PetscPowScalar(as,bm1);
732: /* ds = bs * as; */
733: ds = PetscPowScalar(as,beta);
734: gs = coef*bs*(t0 - ts);
736: tn = x[k][j+1][i];
737: an = 0.5*(t0 + tn);
738: bn = PetscPowScalar(an,bm1);
739: /* dn = bn * an; */
740: dn = PetscPowScalar(an,beta);
741: gn = coef*bn*(tn - t0);
743: /* left-hand down interior edge */
744: if (k == 0) {
746: tu = x[k+1][j][i];
747: au = 0.5*(t0 + tu);
748: bu = PetscPowScalar(au,bm1);
749: /* du = bu * au; */
750: du = PetscPowScalar(au,beta);
751: gu = coef*bu*(tu - t0);
753: c[0].k = k; c[0].j = j-1; c[0].i = i;
754: v[0] = -hzhxdhy*(ds - gs);
755: c[1].k = k; c[1].j = j; c[1].i = i;
756: v[1] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
757: c[2].k = k; c[2].j = j; c[2].i = i+1;
758: v[2] = -hyhzdhx*(de + ge);
759: c[3].k = k; c[3].j = j+1; c[3].i = i;
760: v[3] = -hzhxdhy*(dn + gn);
761: c[4].k = k+1; c[4].j = j; c[4].i = i;
762: v[4] = -hxhydhz*(du + gu);
763: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
765: } else if (k == mz-1) { /* left-hand up interior edge */
767: td = x[k-1][j][i];
768: ad = 0.5*(t0 + td);
769: bd = PetscPowScalar(ad,bm1);
770: /* dd = bd * ad; */
771: dd = PetscPowScalar(ad,beta);
772: gd = coef*bd*(t0 - td);
774: c[0].k = k-1; c[0].j = j; c[0].i = i;
775: v[0] = -hxhydhz*(dd - gd);
776: c[1].k = k; c[1].j = j-1; c[1].i = i;
777: v[1] = -hzhxdhy*(ds - gs);
778: c[2].k = k; c[2].j = j; c[2].i = i;
779: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
780: c[3].k = k; c[3].j = j; c[3].i = i+1;
781: v[3] = -hyhzdhx*(de + ge);
782: c[4].k = k; c[4].j = j+1; c[4].i = i;
783: v[4] = -hzhxdhy*(dn + gn);
784: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
785: } else { /* left-hand interior plane */
787: td = x[k-1][j][i];
788: ad = 0.5*(t0 + td);
789: bd = PetscPowScalar(ad,bm1);
790: /* dd = bd * ad; */
791: dd = PetscPowScalar(ad,beta);
792: gd = coef*bd*(t0 - td);
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-1; c[0].j = j; c[0].i = i;
802: v[0] = -hxhydhz*(dd - gd);
803: c[1].k = k; c[1].j = j-1; c[1].i = i;
804: v[1] = -hzhxdhy*(ds - gs);
805: c[2].k = k; c[2].j = j; c[2].i = i;
806: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
807: c[3].k = k; c[3].j = j; c[3].i = i+1;
808: v[3] = -hyhzdhx*(de + ge);
809: c[4].k = k; c[4].j = j+1; c[4].i = i;
810: v[4] = -hzhxdhy*(dn + gn);
811: c[5].k = k+1; c[5].j = j; c[5].i = i;
812: v[5] = -hxhydhz*(du + gu);
813: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
814: }
815: }
817: } else if (i == mx-1) {
819: /* right-hand plane boundary */
820: tw = x[k][j][i-1];
821: aw = 0.5*(t0 + tw);
822: bw = PetscPowScalar(aw,bm1);
823: /* dw = bw * aw */
824: dw = PetscPowScalar(aw,beta);
825: gw = coef*bw*(t0 - tw);
827: te = tright;
828: ae = 0.5*(t0 + te);
829: be = PetscPowScalar(ae,bm1);
830: /* de = be * ae; */
831: de = PetscPowScalar(ae,beta);
832: ge = coef*be*(te - t0);
834: /* right-hand bottom edge */
835: if (j == 0) {
837: tn = x[k][j+1][i];
838: an = 0.5*(t0 + tn);
839: bn = PetscPowScalar(an,bm1);
840: /* dn = bn * an; */
841: dn = PetscPowScalar(an,beta);
842: gn = coef*bn*(tn - t0);
844: /* right-hand bottom down corner */
845: if (k == 0) {
847: tu = x[k+1][j][i];
848: au = 0.5*(t0 + tu);
849: bu = PetscPowScalar(au,bm1);
850: /* du = bu * au; */
851: du = PetscPowScalar(au,beta);
852: gu = coef*bu*(tu - t0);
854: c[0].k = k; c[0].j = j; c[0].i = i-1;
855: v[0] = -hyhzdhx*(dw - gw);
856: c[1].k = k; c[1].j = j; c[1].i = i;
857: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
858: c[2].k = k; c[2].j = j+1; c[2].i = i;
859: v[2] = -hzhxdhy*(dn + gn);
860: c[3].k = k+1; c[3].j = j; c[3].i = i;
861: v[3] = -hxhydhz*(du + gu);
862: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
864: /* right-hand bottom interior edge */
865: } else if (k < mz-1) {
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: td = x[k-1][j][i];
875: ad = 0.5*(t0 + td);
876: bd = PetscPowScalar(ad,bm1);
877: /* dd = bd * ad; */
878: dd = PetscPowScalar(ad,beta);
879: gd = coef*bd*(td - t0);
881: c[0].k = k-1; c[0].j = j; c[0].i = i;
882: v[0] = -hxhydhz*(dd - gd);
883: c[1].k = k; c[1].j = j; c[1].i = i-1;
884: v[1] = -hyhzdhx*(dw - gw);
885: c[2].k = k; c[2].j = j; c[2].i = i;
886: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
887: c[3].k = k; c[3].j = j+1; c[3].i = i;
888: v[3] = -hzhxdhy*(dn + gn);
889: c[4].k = k+1; c[4].j = j; c[4].i = i;
890: v[4] = -hxhydhz*(du + gu);
891: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
893: /* right-hand bottom up corner */
894: } else {
896: td = x[k-1][j][i];
897: ad = 0.5*(t0 + td);
898: bd = PetscPowScalar(ad,bm1);
899: /* dd = bd * ad; */
900: dd = PetscPowScalar(ad,beta);
901: gd = coef*bd*(td - t0);
903: c[0].k = k-1; c[0].j = j; c[0].i = i;
904: v[0] = -hxhydhz*(dd - gd);
905: c[1].k = k; c[1].j = j; c[1].i = i-1;
906: v[1] = -hyhzdhx*(dw - gw);
907: c[2].k = k; c[2].j = j; c[2].i = i;
908: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
909: c[3].k = k; c[3].j = j+1; c[3].i = i;
910: v[3] = -hzhxdhy*(dn + gn);
911: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
912: }
914: /* right-hand top edge */
915: } else if (j == my-1) {
917: ts = x[k][j-1][i];
918: as = 0.5*(t0 + ts);
919: bs = PetscPowScalar(as,bm1);
920: /* ds = bs * as; */
921: ds = PetscPowScalar(as,beta);
922: gs = coef*bs*(ts - t0);
924: /* right-hand top down corner */
925: if (k == 0) {
927: tu = x[k+1][j][i];
928: au = 0.5*(t0 + tu);
929: bu = PetscPowScalar(au,bm1);
930: /* du = bu * au; */
931: du = PetscPowScalar(au,beta);
932: gu = coef*bu*(tu - t0);
934: c[0].k = k; c[0].j = j-1; c[0].i = i;
935: v[0] = -hzhxdhy*(ds - gs);
936: c[1].k = k; c[1].j = j; c[1].i = i-1;
937: v[1] = -hyhzdhx*(dw - gw);
938: c[2].k = k; c[2].j = j; c[2].i = i;
939: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
940: c[3].k = k+1; c[3].j = j; c[3].i = i;
941: v[3] = -hxhydhz*(du + gu);
942: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
944: /* right-hand top interior edge */
945: } else if (k < mz-1) {
947: tu = x[k+1][j][i];
948: au = 0.5*(t0 + tu);
949: bu = PetscPowScalar(au,bm1);
950: /* du = bu * au; */
951: du = PetscPowScalar(au,beta);
952: gu = coef*bu*(tu - t0);
954: td = x[k-1][j][i];
955: ad = 0.5*(t0 + td);
956: bd = PetscPowScalar(ad,bm1);
957: /* dd = bd * ad; */
958: dd = PetscPowScalar(ad,beta);
959: gd = coef*bd*(td - t0);
961: c[0].k = k-1; c[0].j = j; c[0].i = i;
962: v[0] = -hxhydhz*(dd - gd);
963: c[1].k = k; c[1].j = j-1; c[1].i = i;
964: v[1] = -hzhxdhy*(ds - gs);
965: c[2].k = k; c[2].j = j; c[2].i = i-1;
966: v[2] = -hyhzdhx*(dw - gw);
967: c[3].k = k; c[3].j = j; c[3].i = i;
968: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
969: c[4].k = k+1; c[4].j = j; c[4].i = i;
970: v[4] = -hxhydhz*(du + gu);
971: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
973: /* right-hand top up corner */
974: } else {
976: td = x[k-1][j][i];
977: ad = 0.5*(t0 + td);
978: bd = PetscPowScalar(ad,bm1);
979: /* dd = bd * ad; */
980: dd = PetscPowScalar(ad,beta);
981: gd = coef*bd*(td - t0);
983: c[0].k = k-1; c[0].j = j; c[0].i = i;
984: v[0] = -hxhydhz*(dd - gd);
985: c[1].k = k; c[1].j = j-1; c[1].i = i;
986: v[1] = -hzhxdhy*(ds - gs);
987: c[2].k = k; c[2].j = j; c[2].i = i-1;
988: v[2] = -hyhzdhx*(dw - gw);
989: c[3].k = k; c[3].j = j; c[3].i = i;
990: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
991: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
992: }
994: } else {
996: ts = x[k][j-1][i];
997: as = 0.5*(t0 + ts);
998: bs = PetscPowScalar(as,bm1);
999: /* ds = bs * as; */
1000: ds = PetscPowScalar(as,beta);
1001: gs = coef*bs*(t0 - ts);
1003: tn = x[k][j+1][i];
1004: an = 0.5*(t0 + tn);
1005: bn = PetscPowScalar(an,bm1);
1006: /* dn = bn * an; */
1007: dn = PetscPowScalar(an,beta);
1008: gn = coef*bn*(tn - t0);
1010: /* right-hand down interior edge */
1011: if (k == 0) {
1013: tu = x[k+1][j][i];
1014: au = 0.5*(t0 + tu);
1015: bu = PetscPowScalar(au,bm1);
1016: /* du = bu * au; */
1017: du = PetscPowScalar(au,beta);
1018: gu = coef*bu*(tu - t0);
1020: c[0].k = k; c[0].j = j-1; c[0].i = i;
1021: v[0] = -hzhxdhy*(ds - gs);
1022: c[1].k = k; c[1].j = j; c[1].i = i-1;
1023: v[1] = -hyhzdhx*(dw - gw);
1024: c[2].k = k; c[2].j = j; c[2].i = i;
1025: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1026: c[3].k = k; c[3].j = j+1; c[3].i = i;
1027: v[3] = -hzhxdhy*(dn + gn);
1028: c[4].k = k+1; c[4].j = j; c[4].i = i;
1029: v[4] = -hxhydhz*(du + gu);
1030: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1032: } else if (k == mz-1) { /* right-hand up interior edge */
1034: td = x[k-1][j][i];
1035: ad = 0.5*(t0 + td);
1036: bd = PetscPowScalar(ad,bm1);
1037: /* dd = bd * ad; */
1038: dd = PetscPowScalar(ad,beta);
1039: gd = coef*bd*(t0 - td);
1041: c[0].k = k-1; c[0].j = j; c[0].i = i;
1042: v[0] = -hxhydhz*(dd - gd);
1043: c[1].k = k; c[1].j = j-1; c[1].i = i;
1044: v[1] = -hzhxdhy*(ds - gs);
1045: c[2].k = k; c[2].j = j; c[2].i = i-1;
1046: v[2] = -hyhzdhx*(dw - gw);
1047: c[3].k = k; c[3].j = j; c[3].i = i;
1048: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1049: c[4].k = k; c[4].j = j+1; c[4].i = i;
1050: v[4] = -hzhxdhy*(dn + gn);
1051: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1053: } else { /* right-hand interior plane */
1055: td = x[k-1][j][i];
1056: ad = 0.5*(t0 + td);
1057: bd = PetscPowScalar(ad,bm1);
1058: /* dd = bd * ad; */
1059: dd = PetscPowScalar(ad,beta);
1060: gd = coef*bd*(t0 - td);
1062: tu = x[k+1][j][i];
1063: au = 0.5*(t0 + tu);
1064: bu = PetscPowScalar(au,bm1);
1065: /* du = bu * au; */
1066: du = PetscPowScalar(au,beta);
1067: gu = coef*bu*(tu - t0);
1069: c[0].k = k-1; c[0].j = j; c[0].i = i;
1070: v[0] = -hxhydhz*(dd - gd);
1071: c[1].k = k; c[1].j = j-1; c[1].i = i;
1072: v[1] = -hzhxdhy*(ds - gs);
1073: c[2].k = k; c[2].j = j; c[2].i = i-1;
1074: v[2] = -hyhzdhx*(dw - gw);
1075: c[3].k = k; c[3].j = j; c[3].i = i;
1076: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1077: c[4].k = k; c[4].j = j+1; c[4].i = i;
1078: v[4] = -hzhxdhy*(dn + gn);
1079: c[5].k = k+1; c[5].j = j; c[5].i = i;
1080: v[5] = -hxhydhz*(du + gu);
1081: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1082: }
1083: }
1085: } else if (j == 0) {
1087: tw = x[k][j][i-1];
1088: aw = 0.5*(t0 + tw);
1089: bw = PetscPowScalar(aw,bm1);
1090: /* dw = bw * aw */
1091: dw = PetscPowScalar(aw,beta);
1092: gw = coef*bw*(t0 - tw);
1094: te = x[k][j][i+1];
1095: ae = 0.5*(t0 + te);
1096: be = PetscPowScalar(ae,bm1);
1097: /* de = be * ae; */
1098: de = PetscPowScalar(ae,beta);
1099: ge = coef*be*(te - t0);
1101: tn = x[k][j+1][i];
1102: an = 0.5*(t0 + tn);
1103: bn = PetscPowScalar(an,bm1);
1104: /* dn = bn * an; */
1105: dn = PetscPowScalar(an,beta);
1106: gn = coef*bn*(tn - t0);
1109: /* bottom down interior edge */
1110: if (k == 0) {
1112: tu = x[k+1][j][i];
1113: au = 0.5*(t0 + tu);
1114: bu = PetscPowScalar(au,bm1);
1115: /* du = bu * au; */
1116: du = PetscPowScalar(au,beta);
1117: gu = coef*bu*(tu - t0);
1119: c[0].k = k; c[0].j = j; c[0].i = i-1;
1120: v[0] = -hyhzdhx*(dw - gw);
1121: c[1].k = k; c[1].j = j; c[1].i = i;
1122: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1123: c[2].k = k; c[2].j = j; c[2].i = i+1;
1124: v[2] = -hyhzdhx*(de + ge);
1125: c[3].k = k; c[3].j = j+1; c[3].i = i;
1126: v[3] = -hzhxdhy*(dn + gn);
1127: c[4].k = k+1; c[4].j = j; c[4].i = i;
1128: v[4] = -hxhydhz*(du + gu);
1129: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1131: } else if (k == mz-1) { /* bottom up interior edge */
1133: td = x[k-1][j][i];
1134: ad = 0.5*(t0 + td);
1135: bd = PetscPowScalar(ad,bm1);
1136: /* dd = bd * ad; */
1137: dd = PetscPowScalar(ad,beta);
1138: gd = coef*bd*(td - t0);
1140: c[0].k = k-1; c[0].j = j; c[0].i = i;
1141: v[0] = -hxhydhz*(dd - gd);
1142: c[1].k = k; c[1].j = j; c[1].i = i-1;
1143: v[1] = -hyhzdhx*(dw - gw);
1144: c[2].k = k; c[2].j = j; c[2].i = i;
1145: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1146: c[3].k = k; c[3].j = j; c[3].i = i+1;
1147: v[3] = -hyhzdhx*(de + ge);
1148: c[4].k = k; c[4].j = j+1; c[4].i = i;
1149: v[4] = -hzhxdhy*(dn + gn);
1150: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1152: } else { /* bottom interior plane */
1154: tu = x[k+1][j][i];
1155: au = 0.5*(t0 + tu);
1156: bu = PetscPowScalar(au,bm1);
1157: /* du = bu * au; */
1158: du = PetscPowScalar(au,beta);
1159: gu = coef*bu*(tu - t0);
1161: td = x[k-1][j][i];
1162: ad = 0.5*(t0 + td);
1163: bd = PetscPowScalar(ad,bm1);
1164: /* dd = bd * ad; */
1165: dd = PetscPowScalar(ad,beta);
1166: gd = coef*bd*(td - t0);
1168: c[0].k = k-1; c[0].j = j; c[0].i = i;
1169: v[0] = -hxhydhz*(dd - gd);
1170: c[1].k = k; c[1].j = j; c[1].i = i-1;
1171: v[1] = -hyhzdhx*(dw - gw);
1172: c[2].k = k; c[2].j = j; c[2].i = i;
1173: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1174: c[3].k = k; c[3].j = j; c[3].i = i+1;
1175: v[3] = -hyhzdhx*(de + ge);
1176: c[4].k = k; c[4].j = j+1; c[4].i = i;
1177: v[4] = -hzhxdhy*(dn + gn);
1178: c[5].k = k+1; c[5].j = j; c[5].i = i;
1179: v[5] = -hxhydhz*(du + gu);
1180: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1181: }
1183: } else if (j == my-1) {
1185: tw = x[k][j][i-1];
1186: aw = 0.5*(t0 + tw);
1187: bw = PetscPowScalar(aw,bm1);
1188: /* dw = bw * aw */
1189: dw = PetscPowScalar(aw,beta);
1190: gw = coef*bw*(t0 - tw);
1192: te = x[k][j][i+1];
1193: ae = 0.5*(t0 + te);
1194: be = PetscPowScalar(ae,bm1);
1195: /* de = be * ae; */
1196: de = PetscPowScalar(ae,beta);
1197: ge = coef*be*(te - t0);
1199: ts = x[k][j-1][i];
1200: as = 0.5*(t0 + ts);
1201: bs = PetscPowScalar(as,bm1);
1202: /* ds = bs * as; */
1203: ds = PetscPowScalar(as,beta);
1204: gs = coef*bs*(t0 - ts);
1206: /* top down interior edge */
1207: if (k == 0) {
1209: tu = x[k+1][j][i];
1210: au = 0.5*(t0 + tu);
1211: bu = PetscPowScalar(au,bm1);
1212: /* du = bu * au; */
1213: du = PetscPowScalar(au,beta);
1214: gu = coef*bu*(tu - t0);
1216: c[0].k = k; c[0].j = j-1; c[0].i = i;
1217: v[0] = -hzhxdhy*(ds - gs);
1218: c[1].k = k; c[1].j = j; c[1].i = i-1;
1219: v[1] = -hyhzdhx*(dw - gw);
1220: c[2].k = k; c[2].j = j; c[2].i = i;
1221: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1222: c[3].k = k; c[3].j = j; c[3].i = i+1;
1223: v[3] = -hyhzdhx*(de + ge);
1224: c[4].k = k+1; c[4].j = j; c[4].i = i;
1225: v[4] = -hxhydhz*(du + gu);
1226: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1228: } else if (k == mz-1) { /* top up interior edge */
1230: td = x[k-1][j][i];
1231: ad = 0.5*(t0 + td);
1232: bd = PetscPowScalar(ad,bm1);
1233: /* dd = bd * ad; */
1234: dd = PetscPowScalar(ad,beta);
1235: gd = coef*bd*(td - t0);
1237: c[0].k = k-1; c[0].j = j; c[0].i = i;
1238: v[0] = -hxhydhz*(dd - gd);
1239: c[1].k = k; c[1].j = j-1; c[1].i = i;
1240: v[1] = -hzhxdhy*(ds - gs);
1241: c[2].k = k; c[2].j = j; c[2].i = i-1;
1242: v[2] = -hyhzdhx*(dw - gw);
1243: c[3].k = k; c[3].j = j; c[3].i = i;
1244: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1245: c[4].k = k; c[4].j = j; c[4].i = i+1;
1246: v[4] = -hyhzdhx*(de + ge);
1247: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1249: } else { /* top interior plane */
1251: tu = x[k+1][j][i];
1252: au = 0.5*(t0 + tu);
1253: bu = PetscPowScalar(au,bm1);
1254: /* du = bu * au; */
1255: du = PetscPowScalar(au,beta);
1256: gu = coef*bu*(tu - t0);
1258: td = x[k-1][j][i];
1259: ad = 0.5*(t0 + td);
1260: bd = PetscPowScalar(ad,bm1);
1261: /* dd = bd * ad; */
1262: dd = PetscPowScalar(ad,beta);
1263: gd = coef*bd*(td - t0);
1265: c[0].k = k-1; c[0].j = j; c[0].i = i;
1266: v[0] = -hxhydhz*(dd - gd);
1267: c[1].k = k; c[1].j = j-1; c[1].i = i;
1268: v[1] = -hzhxdhy*(ds - gs);
1269: c[2].k = k; c[2].j = j; c[2].i = i-1;
1270: v[2] = -hyhzdhx*(dw - gw);
1271: c[3].k = k; c[3].j = j; c[3].i = i;
1272: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1273: c[4].k = k; c[4].j = j; c[4].i = i+1;
1274: v[4] = -hyhzdhx*(de + ge);
1275: c[5].k = k+1; c[5].j = j; c[5].i = i;
1276: v[5] = -hxhydhz*(du + gu);
1277: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1278: }
1280: } else if (k == 0) {
1282: /* down interior plane */
1284: tw = x[k][j][i-1];
1285: aw = 0.5*(t0 + tw);
1286: bw = PetscPowScalar(aw,bm1);
1287: /* dw = bw * aw */
1288: dw = PetscPowScalar(aw,beta);
1289: gw = coef*bw*(t0 - tw);
1291: te = x[k][j][i+1];
1292: ae = 0.5*(t0 + te);
1293: be = PetscPowScalar(ae,bm1);
1294: /* de = be * ae; */
1295: de = PetscPowScalar(ae,beta);
1296: ge = coef*be*(te - t0);
1298: ts = x[k][j-1][i];
1299: as = 0.5*(t0 + ts);
1300: bs = PetscPowScalar(as,bm1);
1301: /* ds = bs * as; */
1302: ds = PetscPowScalar(as,beta);
1303: gs = coef*bs*(t0 - ts);
1305: tn = x[k][j+1][i];
1306: an = 0.5*(t0 + tn);
1307: bn = PetscPowScalar(an,bm1);
1308: /* dn = bn * an; */
1309: dn = PetscPowScalar(an,beta);
1310: gn = coef*bn*(tn - t0);
1312: tu = x[k+1][j][i];
1313: au = 0.5*(t0 + tu);
1314: bu = PetscPowScalar(au,bm1);
1315: /* du = bu * au; */
1316: du = PetscPowScalar(au,beta);
1317: gu = coef*bu*(tu - t0);
1319: c[0].k = k; c[0].j = j-1; c[0].i = i;
1320: v[0] = -hzhxdhy*(ds - gs);
1321: c[1].k = k; c[1].j = j; c[1].i = i-1;
1322: v[1] = -hyhzdhx*(dw - gw);
1323: c[2].k = k; c[2].j = j; c[2].i = i;
1324: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1325: c[3].k = k; c[3].j = j; c[3].i = i+1;
1326: v[3] = -hyhzdhx*(de + ge);
1327: c[4].k = k; c[4].j = j+1; c[4].i = i;
1328: v[4] = -hzhxdhy*(dn + gn);
1329: c[5].k = k+1; c[5].j = j; c[5].i = i;
1330: v[5] = -hxhydhz*(du + gu);
1331: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1333: } else if (k == mz-1) {
1335: /* up interior plane */
1337: tw = x[k][j][i-1];
1338: aw = 0.5*(t0 + tw);
1339: bw = PetscPowScalar(aw,bm1);
1340: /* dw = bw * aw */
1341: dw = PetscPowScalar(aw,beta);
1342: gw = coef*bw*(t0 - tw);
1344: te = x[k][j][i+1];
1345: ae = 0.5*(t0 + te);
1346: be = PetscPowScalar(ae,bm1);
1347: /* de = be * ae; */
1348: de = PetscPowScalar(ae,beta);
1349: ge = coef*be*(te - t0);
1351: ts = x[k][j-1][i];
1352: as = 0.5*(t0 + ts);
1353: bs = PetscPowScalar(as,bm1);
1354: /* ds = bs * as; */
1355: ds = PetscPowScalar(as,beta);
1356: gs = coef*bs*(t0 - ts);
1358: tn = x[k][j+1][i];
1359: an = 0.5*(t0 + tn);
1360: bn = PetscPowScalar(an,bm1);
1361: /* dn = bn * an; */
1362: dn = PetscPowScalar(an,beta);
1363: gn = coef*bn*(tn - t0);
1365: td = x[k-1][j][i];
1366: ad = 0.5*(t0 + td);
1367: bd = PetscPowScalar(ad,bm1);
1368: /* dd = bd * ad; */
1369: dd = PetscPowScalar(ad,beta);
1370: gd = coef*bd*(t0 - td);
1372: c[0].k = k-1; c[0].j = j; c[0].i = i;
1373: v[0] = -hxhydhz*(dd - gd);
1374: c[1].k = k; c[1].j = j-1; c[1].i = i;
1375: v[1] = -hzhxdhy*(ds - gs);
1376: c[2].k = k; c[2].j = j; c[2].i = i-1;
1377: v[2] = -hyhzdhx*(dw - gw);
1378: c[3].k = k; c[3].j = j; c[3].i = i;
1379: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1380: c[4].k = k; c[4].j = j; c[4].i = i+1;
1381: v[4] = -hyhzdhx*(de + ge);
1382: c[5].k = k; c[5].j = j+1; c[5].i = i;
1383: v[5] = -hzhxdhy*(dn + gn);
1384: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1385: }
1386: }
1387: }
1388: }
1389: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
1390: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
1391: if (jac != J) {
1392: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1393: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1394: }
1395: DMDAVecRestoreArray(da,localX,&x);
1396: DMRestoreLocalVector(da,&localX);
1398: PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
1399: return(0);
1400: }
1403: /*TEST
1405: test:
1406: nsize: 4
1407: args: -snes_monitor_short -pc_mg_type full -ksp_type fgmres -pc_type mg -snes_view -pc_mg_levels 2 -pc_mg_galerkin pmat
1408: requires: !single
1410: TEST*/