Actual source code: ex20.c
petsc-3.4.5 2014-06-29
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 <petscdmda.h>
45: /* User-defined application context */
47: typedef struct {
48: PetscReal tleft,tright; /* Dirichlet boundary conditions */
49: PetscReal beta,bm1,coef; /* nonlinear diffusivity parameterizations */
50: } AppCtx;
52: #define POWFLOP 5 /* assume a pow() takes five flops */
54: extern PetscErrorCode FormInitialGuess(SNES,Vec,void*);
55: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
56: extern PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
60: int main(int argc,char **argv)
61: {
62: SNES snes;
63: AppCtx user;
65: PetscInt its,lits;
66: PetscReal litspit;
67: DM da;
69: PetscInitialize(&argc,&argv,NULL,help);
71: /* set problem parameters */
72: user.tleft = 1.0;
73: user.tright = 0.1;
74: user.beta = 2.5;
75: PetscOptionsGetReal(NULL,"-tleft",&user.tleft,NULL);
76: PetscOptionsGetReal(NULL,"-tright",&user.tright,NULL);
77: PetscOptionsGetReal(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,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-5,-5,-5,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);
85: DMSetApplicationContext(da,&user);
87: /*
88: Create the nonlinear solver
89: */
90: SNESCreate(PETSC_COMM_WORLD,&snes);
91: SNESSetDM(snes,da);
92: SNESSetFunction(snes,NULL,FormFunction,&user);
93: SNESSetJacobian(snes,NULL,NULL,FormJacobian,&user);
94: SNESSetFromOptions(snes);
95: SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);
97: SNESSolve(snes,NULL,NULL);
98: SNESGetIterationNumber(snes,&its);
99: SNESGetLinearSolveIterations(snes,&lits);
100: litspit = ((PetscReal)lits)/((PetscReal)its);
101: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
102: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
103: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",litspit);
105: SNESDestroy(&snes);
106: DMDestroy(&da);
107: PetscFinalize();
109: return 0;
110: }
111: /* -------------------- Form initial approximation ----------------- */
114: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
115: {
116: AppCtx *user;
117: PetscInt i,j,k,xs,ys,xm,ym,zs,zm;
119: PetscScalar ***x;
120: DM da;
123: SNESGetDM(snes,&da);
124: DMGetApplicationContext(da,&user);
125: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
126: DMDAVecGetArray(da,X,&x);
128: /* Compute initial guess */
129: for (k=zs; k<zs+zm; k++) {
130: for (j=ys; j<ys+ym; j++) {
131: for (i=xs; i<xs+xm; i++) {
132: x[k][j][i] = user->tleft;
133: }
134: }
135: }
136: DMDAVecRestoreArray(da,X,&x);
137: return(0);
138: }
139: /* -------------------- Evaluate Function F(x) --------------------- */
142: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
143: {
144: AppCtx *user = (AppCtx*)ptr;
146: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
147: PetscScalar zero = 0.0,one = 1.0;
148: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
149: 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;
150: PetscScalar tleft,tright,beta,td,ad,dd,fd=0.0,tu,au,du=0.0,fu=0.0;
151: PetscScalar ***x,***f;
152: Vec localX;
153: DM da;
156: SNESGetDM(snes,&da);
157: DMGetLocalVector(da,&localX);
158: DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
159: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
160: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
161: tleft = user->tleft; tright = user->tright;
162: beta = user->beta;
164: /* Get ghost points */
165: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
166: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
167: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
168: DMDAVecGetArray(da,localX,&x);
169: DMDAVecGetArray(da,F,&f);
171: /* Evaluate function */
172: for (k=zs; k<zs+zm; k++) {
173: for (j=ys; j<ys+ym; j++) {
174: for (i=xs; i<xs+xm; i++) {
175: t0 = x[k][j][i];
177: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
179: /* general interior volume */
181: tw = x[k][j][i-1];
182: aw = 0.5*(t0 + tw);
183: dw = PetscPowScalar(aw,beta);
184: fw = dw*(t0 - tw);
186: te = x[k][j][i+1];
187: ae = 0.5*(t0 + te);
188: de = PetscPowScalar(ae,beta);
189: fe = de*(te - t0);
191: ts = x[k][j-1][i];
192: as = 0.5*(t0 + ts);
193: ds = PetscPowScalar(as,beta);
194: fs = ds*(t0 - ts);
196: tn = x[k][j+1][i];
197: an = 0.5*(t0 + tn);
198: dn = PetscPowScalar(an,beta);
199: fn = dn*(tn - t0);
201: td = x[k-1][j][i];
202: ad = 0.5*(t0 + td);
203: dd = PetscPowScalar(ad,beta);
204: fd = dd*(t0 - td);
206: tu = x[k+1][j][i];
207: au = 0.5*(t0 + tu);
208: du = PetscPowScalar(au,beta);
209: fu = du*(tu - t0);
211: } else if (i == 0) {
213: /* left-hand (west) boundary */
214: tw = tleft;
215: aw = 0.5*(t0 + tw);
216: dw = PetscPowScalar(aw,beta);
217: fw = dw*(t0 - tw);
219: te = x[k][j][i+1];
220: ae = 0.5*(t0 + te);
221: de = PetscPowScalar(ae,beta);
222: fe = de*(te - t0);
224: if (j > 0) {
225: ts = x[k][j-1][i];
226: as = 0.5*(t0 + ts);
227: ds = PetscPowScalar(as,beta);
228: fs = ds*(t0 - ts);
229: } else {
230: fs = zero;
231: }
233: if (j < my-1) {
234: tn = x[k][j+1][i];
235: an = 0.5*(t0 + tn);
236: dn = PetscPowScalar(an,beta);
237: fn = dn*(tn - t0);
238: } else {
239: fn = zero;
240: }
242: if (k > 0) {
243: td = x[k-1][j][i];
244: ad = 0.5*(t0 + td);
245: dd = PetscPowScalar(ad,beta);
246: fd = dd*(t0 - td);
247: } else {
248: fd = zero;
249: }
251: if (k < mz-1) {
252: tu = x[k+1][j][i];
253: au = 0.5*(t0 + tu);
254: du = PetscPowScalar(au,beta);
255: fu = du*(tu - t0);
256: } else {
257: fu = zero;
258: }
260: } else if (i == mx-1) {
262: /* right-hand (east) boundary */
263: tw = x[k][j][i-1];
264: aw = 0.5*(t0 + tw);
265: dw = PetscPowScalar(aw,beta);
266: fw = dw*(t0 - tw);
268: te = tright;
269: ae = 0.5*(t0 + te);
270: de = PetscPowScalar(ae,beta);
271: fe = de*(te - t0);
273: if (j > 0) {
274: ts = x[k][j-1][i];
275: as = 0.5*(t0 + ts);
276: ds = PetscPowScalar(as,beta);
277: fs = ds*(t0 - ts);
278: } else {
279: fs = zero;
280: }
282: if (j < my-1) {
283: tn = x[k][j+1][i];
284: an = 0.5*(t0 + tn);
285: dn = PetscPowScalar(an,beta);
286: fn = dn*(tn - t0);
287: } else {
288: fn = zero;
289: }
291: if (k > 0) {
292: td = x[k-1][j][i];
293: ad = 0.5*(t0 + td);
294: dd = PetscPowScalar(ad,beta);
295: fd = dd*(t0 - td);
296: } else {
297: fd = zero;
298: }
300: if (k < mz-1) {
301: tu = x[k+1][j][i];
302: au = 0.5*(t0 + tu);
303: du = PetscPowScalar(au,beta);
304: fu = du*(tu - t0);
305: } else {
306: fu = zero;
307: }
309: } else if (j == 0) {
311: /* bottom (south) boundary, and i <> 0 or mx-1 */
312: tw = x[k][j][i-1];
313: aw = 0.5*(t0 + tw);
314: dw = PetscPowScalar(aw,beta);
315: fw = dw*(t0 - tw);
317: te = x[k][j][i+1];
318: ae = 0.5*(t0 + te);
319: de = PetscPowScalar(ae,beta);
320: fe = de*(te - t0);
322: fs = zero;
324: tn = x[k][j+1][i];
325: an = 0.5*(t0 + tn);
326: dn = PetscPowScalar(an,beta);
327: fn = dn*(tn - t0);
329: if (k > 0) {
330: td = x[k-1][j][i];
331: ad = 0.5*(t0 + td);
332: dd = PetscPowScalar(ad,beta);
333: fd = dd*(t0 - td);
334: } else {
335: fd = zero;
336: }
338: if (k < mz-1) {
339: tu = x[k+1][j][i];
340: au = 0.5*(t0 + tu);
341: du = PetscPowScalar(au,beta);
342: fu = du*(tu - t0);
343: } else {
344: fu = zero;
345: }
347: } else if (j == my-1) {
349: /* top (north) boundary, and i <> 0 or mx-1 */
350: tw = x[k][j][i-1];
351: aw = 0.5*(t0 + tw);
352: dw = PetscPowScalar(aw,beta);
353: fw = dw*(t0 - tw);
355: te = x[k][j][i+1];
356: ae = 0.5*(t0 + te);
357: de = PetscPowScalar(ae,beta);
358: fe = de*(te - t0);
360: ts = x[k][j-1][i];
361: as = 0.5*(t0 + ts);
362: ds = PetscPowScalar(as,beta);
363: fs = ds*(t0 - ts);
365: fn = zero;
367: if (k > 0) {
368: td = x[k-1][j][i];
369: ad = 0.5*(t0 + td);
370: dd = PetscPowScalar(ad,beta);
371: fd = dd*(t0 - td);
372: } else {
373: fd = zero;
374: }
376: if (k < mz-1) {
377: tu = x[k+1][j][i];
378: au = 0.5*(t0 + tu);
379: du = PetscPowScalar(au,beta);
380: fu = du*(tu - t0);
381: } else {
382: fu = zero;
383: }
385: } else if (k == 0) {
387: /* down boundary (interior only) */
388: tw = x[k][j][i-1];
389: aw = 0.5*(t0 + tw);
390: dw = PetscPowScalar(aw,beta);
391: fw = dw*(t0 - tw);
393: te = x[k][j][i+1];
394: ae = 0.5*(t0 + te);
395: de = PetscPowScalar(ae,beta);
396: fe = de*(te - t0);
398: ts = x[k][j-1][i];
399: as = 0.5*(t0 + ts);
400: ds = PetscPowScalar(as,beta);
401: fs = ds*(t0 - ts);
403: tn = x[k][j+1][i];
404: an = 0.5*(t0 + tn);
405: dn = PetscPowScalar(an,beta);
406: fn = dn*(tn - t0);
408: fd = zero;
410: tu = x[k+1][j][i];
411: au = 0.5*(t0 + tu);
412: du = PetscPowScalar(au,beta);
413: fu = du*(tu - t0);
415: } else if (k == mz-1) {
417: /* up boundary (interior only) */
418: tw = x[k][j][i-1];
419: aw = 0.5*(t0 + tw);
420: dw = PetscPowScalar(aw,beta);
421: fw = dw*(t0 - tw);
423: te = x[k][j][i+1];
424: ae = 0.5*(t0 + te);
425: de = PetscPowScalar(ae,beta);
426: fe = de*(te - t0);
428: ts = x[k][j-1][i];
429: as = 0.5*(t0 + ts);
430: ds = PetscPowScalar(as,beta);
431: fs = ds*(t0 - ts);
433: tn = x[k][j+1][i];
434: an = 0.5*(t0 + tn);
435: dn = PetscPowScalar(an,beta);
436: fn = dn*(tn - t0);
438: td = x[k-1][j][i];
439: ad = 0.5*(t0 + td);
440: dd = PetscPowScalar(ad,beta);
441: fd = dd*(t0 - td);
443: fu = zero;
444: }
446: f[k][j][i] = -hyhzdhx*(fe-fw) - hzhxdhy*(fn-fs) - hxhydhz*(fu-fd);
447: }
448: }
449: }
450: DMDAVecRestoreArray(da,localX,&x);
451: DMDAVecRestoreArray(da,F,&f);
452: DMRestoreLocalVector(da,&localX);
453: PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
454: return(0);
455: }
456: /* -------------------- Evaluate Jacobian F(x) --------------------- */
459: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ptr)
460: {
461: AppCtx *user = (AppCtx*)ptr;
463: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
464: PetscScalar one = 1.0;
465: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
466: PetscScalar t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw;
467: PetscScalar tleft,tright,beta,td,ad,dd,tu,au,du,v[7],bm1,coef;
468: PetscScalar ***x,bn,bs,be,bw,bu,bd,gn,gs,ge,gw,gu,gd;
469: Vec localX;
470: MatStencil c[7],row;
471: Mat jac = *B;
472: DM da;
475: SNESGetDM(snes,&da);
476: DMGetLocalVector(da,&localX);
477: DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
478: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
479: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
480: tleft = user->tleft; tright = user->tright;
481: beta = user->beta; bm1 = user->bm1; coef = user->coef;
483: /* Get ghost points */
484: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
485: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
486: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
487: DMDAVecGetArray(da,localX,&x);
489: /* Evaluate Jacobian of function */
490: for (k=zs; k<zs+zm; k++) {
491: for (j=ys; j<ys+ym; j++) {
492: for (i=xs; i<xs+xm; i++) {
493: t0 = x[k][j][i];
494: row.k = k; row.j = j; row.i = i;
495: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
497: /* general interior volume */
499: tw = x[k][j][i-1];
500: aw = 0.5*(t0 + tw);
501: bw = PetscPowScalar(aw,bm1);
502: /* dw = bw * aw */
503: dw = PetscPowScalar(aw,beta);
504: gw = coef*bw*(t0 - tw);
506: te = x[k][j][i+1];
507: ae = 0.5*(t0 + te);
508: be = PetscPowScalar(ae,bm1);
509: /* de = be * ae; */
510: de = PetscPowScalar(ae,beta);
511: ge = coef*be*(te - t0);
513: ts = x[k][j-1][i];
514: as = 0.5*(t0 + ts);
515: bs = PetscPowScalar(as,bm1);
516: /* ds = bs * as; */
517: ds = PetscPowScalar(as,beta);
518: gs = coef*bs*(t0 - ts);
520: tn = x[k][j+1][i];
521: an = 0.5*(t0 + tn);
522: bn = PetscPowScalar(an,bm1);
523: /* dn = bn * an; */
524: dn = PetscPowScalar(an,beta);
525: gn = coef*bn*(tn - t0);
527: td = x[k-1][j][i];
528: ad = 0.5*(t0 + td);
529: bd = PetscPowScalar(ad,bm1);
530: /* dd = bd * ad; */
531: dd = PetscPowScalar(ad,beta);
532: gd = coef*bd*(t0 - td);
534: tu = x[k+1][j][i];
535: au = 0.5*(t0 + tu);
536: bu = PetscPowScalar(au,bm1);
537: /* du = bu * au; */
538: du = PetscPowScalar(au,beta);
539: gu = coef*bu*(tu - t0);
541: c[0].k = k-1; c[0].j = j; c[0].i = i; v[0] = -hxhydhz*(dd - gd);
542: c[1].k = k; c[1].j = j-1; c[1].i = i;
543: v[1] = -hzhxdhy*(ds - gs);
544: c[2].k = k; c[2].j = j; c[2].i = i-1;
545: v[2] = -hyhzdhx*(dw - gw);
546: c[3].k = k; c[3].j = j; c[3].i = i;
547: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
548: c[4].k = k; c[4].j = j; c[4].i = i+1;
549: v[4] = -hyhzdhx*(de + ge);
550: c[5].k = k; c[5].j = j+1; c[5].i = i;
551: v[5] = -hzhxdhy*(dn + gn);
552: c[6].k = k+1; c[6].j = j; c[6].i = i;
553: v[6] = -hxhydhz*(du + gu);
554: MatSetValuesStencil(jac,1,&row,7,c,v,INSERT_VALUES);
556: } 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) {
576: tn = x[k][j+1][i];
577: an = 0.5*(t0 + tn);
578: bn = PetscPowScalar(an,bm1);
579: /* dn = bn * an; */
580: dn = PetscPowScalar(an,beta);
581: gn = coef*bn*(tn - t0);
583: /* left-hand bottom down corner */
584: if (k == 0) {
586: tu = x[k+1][j][i];
587: au = 0.5*(t0 + tu);
588: bu = PetscPowScalar(au,bm1);
589: /* du = bu * au; */
590: du = PetscPowScalar(au,beta);
591: gu = coef*bu*(tu - t0);
593: c[0].k = k; c[0].j = j; c[0].i = i;
594: v[0] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
595: c[1].k = k; c[1].j = j; c[1].i = i+1;
596: v[1] = -hyhzdhx*(de + ge);
597: c[2].k = k; c[2].j = j+1; c[2].i = i;
598: v[2] = -hzhxdhy*(dn + gn);
599: c[3].k = k+1; c[3].j = j; c[3].i = i;
600: v[3] = -hxhydhz*(du + gu);
601: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
603: /* left-hand bottom interior edge */
604: } else if (k < mz-1) {
606: tu = x[k+1][j][i];
607: au = 0.5*(t0 + tu);
608: bu = PetscPowScalar(au,bm1);
609: /* du = bu * au; */
610: du = PetscPowScalar(au,beta);
611: gu = coef*bu*(tu - t0);
613: td = x[k-1][j][i];
614: ad = 0.5*(t0 + td);
615: bd = PetscPowScalar(ad,bm1);
616: /* dd = bd * ad; */
617: dd = PetscPowScalar(ad,beta);
618: gd = coef*bd*(td - t0);
620: c[0].k = k-1; c[0].j = j; c[0].i = i;
621: v[0] = -hxhydhz*(dd - gd);
622: c[1].k = k; c[1].j = j; c[1].i = i;
623: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
624: c[2].k = k; c[2].j = j; c[2].i = i+1;
625: v[2] = -hyhzdhx*(de + ge);
626: c[3].k = k; c[3].j = j+1; c[3].i = i;
627: v[3] = -hzhxdhy*(dn + gn);
628: c[4].k = k+1; c[4].j = j; c[4].i = i;
629: v[4] = -hxhydhz*(du + gu);
630: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
632: /* left-hand bottom up corner */
633: } else {
635: td = x[k-1][j][i];
636: ad = 0.5*(t0 + td);
637: bd = PetscPowScalar(ad,bm1);
638: /* dd = bd * ad; */
639: dd = PetscPowScalar(ad,beta);
640: gd = coef*bd*(td - t0);
642: c[0].k = k-1; c[0].j = j; c[0].i = i;
643: v[0] = -hxhydhz*(dd - gd);
644: c[1].k = k; c[1].j = j; c[1].i = i;
645: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
646: c[2].k = k; c[2].j = j; c[2].i = i+1;
647: v[2] = -hyhzdhx*(de + ge);
648: c[3].k = k; c[3].j = j+1; c[3].i = i;
649: v[3] = -hzhxdhy*(dn + gn);
650: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
651: }
653: /* left-hand top edge */
654: } else if (j == my-1) {
656: ts = x[k][j-1][i];
657: as = 0.5*(t0 + ts);
658: bs = PetscPowScalar(as,bm1);
659: /* ds = bs * as; */
660: ds = PetscPowScalar(as,beta);
661: gs = coef*bs*(ts - t0);
663: /* left-hand top down corner */
664: if (k == 0) {
666: tu = x[k+1][j][i];
667: au = 0.5*(t0 + tu);
668: bu = PetscPowScalar(au,bm1);
669: /* du = bu * au; */
670: du = PetscPowScalar(au,beta);
671: gu = coef*bu*(tu - t0);
673: c[0].k = k; c[0].j = j-1; c[0].i = i;
674: v[0] = -hzhxdhy*(ds - gs);
675: c[1].k = k; c[1].j = j; c[1].i = i;
676: v[1] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
677: c[2].k = k; c[2].j = j; c[2].i = i+1;
678: v[2] = -hyhzdhx*(de + ge);
679: c[3].k = k+1; c[3].j = j; c[3].i = i;
680: v[3] = -hxhydhz*(du + gu);
681: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
683: /* left-hand top interior edge */
684: } else if (k < mz-1) {
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: td = x[k-1][j][i];
694: ad = 0.5*(t0 + td);
695: bd = PetscPowScalar(ad,bm1);
696: /* dd = bd * ad; */
697: dd = PetscPowScalar(ad,beta);
698: gd = coef*bd*(td - t0);
700: c[0].k = k-1; c[0].j = j; c[0].i = i;
701: v[0] = -hxhydhz*(dd - gd);
702: c[1].k = k; c[1].j = j-1; c[1].i = i;
703: v[1] = -hzhxdhy*(ds - gs);
704: c[2].k = k; c[2].j = j; c[2].i = i;
705: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
706: c[3].k = k; c[3].j = j; c[3].i = i+1;
707: v[3] = -hyhzdhx*(de + ge);
708: c[4].k = k+1; c[4].j = j; c[4].i = i;
709: v[4] = -hxhydhz*(du + gu);
710: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
712: /* left-hand top up corner */
713: } else {
715: td = x[k-1][j][i];
716: ad = 0.5*(t0 + td);
717: bd = PetscPowScalar(ad,bm1);
718: /* dd = bd * ad; */
719: dd = PetscPowScalar(ad,beta);
720: gd = coef*bd*(td - t0);
722: c[0].k = k-1; c[0].j = j; c[0].i = i;
723: v[0] = -hxhydhz*(dd - gd);
724: c[1].k = k; c[1].j = j-1; c[1].i = i;
725: v[1] = -hzhxdhy*(ds - gs);
726: c[2].k = k; c[2].j = j; c[2].i = i;
727: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
728: c[3].k = k; c[3].j = j; c[3].i = i+1;
729: v[3] = -hyhzdhx*(de + ge);
730: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
731: }
733: } else {
735: ts = x[k][j-1][i];
736: as = 0.5*(t0 + ts);
737: bs = PetscPowScalar(as,bm1);
738: /* ds = bs * as; */
739: ds = PetscPowScalar(as,beta);
740: gs = coef*bs*(t0 - ts);
742: tn = x[k][j+1][i];
743: an = 0.5*(t0 + tn);
744: bn = PetscPowScalar(an,bm1);
745: /* dn = bn * an; */
746: dn = PetscPowScalar(an,beta);
747: gn = coef*bn*(tn - t0);
749: /* left-hand down interior edge */
750: if (k == 0) {
752: tu = x[k+1][j][i];
753: au = 0.5*(t0 + tu);
754: bu = PetscPowScalar(au,bm1);
755: /* du = bu * au; */
756: du = PetscPowScalar(au,beta);
757: gu = coef*bu*(tu - t0);
759: c[0].k = k; c[0].j = j-1; c[0].i = i;
760: v[0] = -hzhxdhy*(ds - gs);
761: c[1].k = k; c[1].j = j; c[1].i = i;
762: v[1] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
763: c[2].k = k; c[2].j = j; c[2].i = i+1;
764: v[2] = -hyhzdhx*(de + ge);
765: c[3].k = k; c[3].j = j+1; c[3].i = i;
766: v[3] = -hzhxdhy*(dn + gn);
767: c[4].k = k+1; c[4].j = j; c[4].i = i;
768: v[4] = -hxhydhz*(du + gu);
769: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
771: } else if (k == mz-1) { /* left-hand up interior edge */
773: td = x[k-1][j][i];
774: ad = 0.5*(t0 + td);
775: bd = PetscPowScalar(ad,bm1);
776: /* dd = bd * ad; */
777: dd = PetscPowScalar(ad,beta);
778: gd = coef*bd*(t0 - td);
780: c[0].k = k-1; c[0].j = j; c[0].i = i;
781: v[0] = -hxhydhz*(dd - gd);
782: c[1].k = k; c[1].j = j-1; c[1].i = i;
783: v[1] = -hzhxdhy*(ds - gs);
784: c[2].k = k; c[2].j = j; c[2].i = i;
785: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
786: c[3].k = k; c[3].j = j; c[3].i = i+1;
787: v[3] = -hyhzdhx*(de + ge);
788: c[4].k = k; c[4].j = j+1; c[4].i = i;
789: v[4] = -hzhxdhy*(dn + gn);
790: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
791: } else { /* left-hand interior plane */
793: td = x[k-1][j][i];
794: ad = 0.5*(t0 + td);
795: bd = PetscPowScalar(ad,bm1);
796: /* dd = bd * ad; */
797: dd = PetscPowScalar(ad,beta);
798: gd = coef*bd*(t0 - td);
800: tu = x[k+1][j][i];
801: au = 0.5*(t0 + tu);
802: bu = PetscPowScalar(au,bm1);
803: /* du = bu * au; */
804: du = PetscPowScalar(au,beta);
805: gu = coef*bu*(tu - t0);
807: c[0].k = k-1; c[0].j = j; c[0].i = i;
808: v[0] = -hxhydhz*(dd - gd);
809: c[1].k = k; c[1].j = j-1; c[1].i = i;
810: v[1] = -hzhxdhy*(ds - gs);
811: c[2].k = k; c[2].j = j; c[2].i = i;
812: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
813: c[3].k = k; c[3].j = j; c[3].i = i+1;
814: v[3] = -hyhzdhx*(de + ge);
815: c[4].k = k; c[4].j = j+1; c[4].i = i;
816: v[4] = -hzhxdhy*(dn + gn);
817: c[5].k = k+1; c[5].j = j; c[5].i = i;
818: v[5] = -hxhydhz*(du + gu);
819: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
820: }
821: }
823: } else if (i == mx-1) {
825: /* right-hand plane boundary */
826: tw = x[k][j][i-1];
827: aw = 0.5*(t0 + tw);
828: bw = PetscPowScalar(aw,bm1);
829: /* dw = bw * aw */
830: dw = PetscPowScalar(aw,beta);
831: gw = coef*bw*(t0 - tw);
833: te = tright;
834: ae = 0.5*(t0 + te);
835: be = PetscPowScalar(ae,bm1);
836: /* de = be * ae; */
837: de = PetscPowScalar(ae,beta);
838: ge = coef*be*(te - t0);
840: /* right-hand bottom edge */
841: if (j == 0) {
843: tn = x[k][j+1][i];
844: an = 0.5*(t0 + tn);
845: bn = PetscPowScalar(an,bm1);
846: /* dn = bn * an; */
847: dn = PetscPowScalar(an,beta);
848: gn = coef*bn*(tn - t0);
850: /* right-hand bottom down corner */
851: if (k == 0) {
853: tu = x[k+1][j][i];
854: au = 0.5*(t0 + tu);
855: bu = PetscPowScalar(au,bm1);
856: /* du = bu * au; */
857: du = PetscPowScalar(au,beta);
858: gu = coef*bu*(tu - t0);
860: c[0].k = k; c[0].j = j; c[0].i = i-1;
861: v[0] = -hyhzdhx*(dw - gw);
862: c[1].k = k; c[1].j = j; c[1].i = i;
863: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
864: c[2].k = k; c[2].j = j+1; c[2].i = i;
865: v[2] = -hzhxdhy*(dn + gn);
866: c[3].k = k+1; c[3].j = j; c[3].i = i;
867: v[3] = -hxhydhz*(du + gu);
868: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
870: /* right-hand bottom interior edge */
871: } else if (k < mz-1) {
873: tu = x[k+1][j][i];
874: au = 0.5*(t0 + tu);
875: bu = PetscPowScalar(au,bm1);
876: /* du = bu * au; */
877: du = PetscPowScalar(au,beta);
878: gu = coef*bu*(tu - t0);
880: td = x[k-1][j][i];
881: ad = 0.5*(t0 + td);
882: bd = PetscPowScalar(ad,bm1);
883: /* dd = bd * ad; */
884: dd = PetscPowScalar(ad,beta);
885: gd = coef*bd*(td - t0);
887: c[0].k = k-1; c[0].j = j; c[0].i = i;
888: v[0] = -hxhydhz*(dd - gd);
889: c[1].k = k; c[1].j = j; c[1].i = i-1;
890: v[1] = -hyhzdhx*(dw - gw);
891: c[2].k = k; c[2].j = j; c[2].i = i;
892: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
893: c[3].k = k; c[3].j = j+1; c[3].i = i;
894: v[3] = -hzhxdhy*(dn + gn);
895: c[4].k = k+1; c[4].j = j; c[4].i = i;
896: v[4] = -hxhydhz*(du + gu);
897: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
899: /* right-hand bottom up corner */
900: } else {
902: td = x[k-1][j][i];
903: ad = 0.5*(t0 + td);
904: bd = PetscPowScalar(ad,bm1);
905: /* dd = bd * ad; */
906: dd = PetscPowScalar(ad,beta);
907: gd = coef*bd*(td - t0);
909: c[0].k = k-1; c[0].j = j; c[0].i = i;
910: v[0] = -hxhydhz*(dd - gd);
911: c[1].k = k; c[1].j = j; c[1].i = i-1;
912: v[1] = -hyhzdhx*(dw - gw);
913: c[2].k = k; c[2].j = j; c[2].i = i;
914: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
915: c[3].k = k; c[3].j = j+1; c[3].i = i;
916: v[3] = -hzhxdhy*(dn + gn);
917: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
918: }
920: /* right-hand top edge */
921: } else if (j == my-1) {
923: ts = x[k][j-1][i];
924: as = 0.5*(t0 + ts);
925: bs = PetscPowScalar(as,bm1);
926: /* ds = bs * as; */
927: ds = PetscPowScalar(as,beta);
928: gs = coef*bs*(ts - t0);
930: /* right-hand top down corner */
931: if (k == 0) {
933: tu = x[k+1][j][i];
934: au = 0.5*(t0 + tu);
935: bu = PetscPowScalar(au,bm1);
936: /* du = bu * au; */
937: du = PetscPowScalar(au,beta);
938: gu = coef*bu*(tu - t0);
940: c[0].k = k; c[0].j = j-1; c[0].i = i;
941: v[0] = -hzhxdhy*(ds - gs);
942: c[1].k = k; c[1].j = j; c[1].i = i-1;
943: v[1] = -hyhzdhx*(dw - gw);
944: c[2].k = k; c[2].j = j; c[2].i = i;
945: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
946: c[3].k = k+1; c[3].j = j; c[3].i = i;
947: v[3] = -hxhydhz*(du + gu);
948: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
950: /* right-hand top interior edge */
951: } else if (k < mz-1) {
953: tu = x[k+1][j][i];
954: au = 0.5*(t0 + tu);
955: bu = PetscPowScalar(au,bm1);
956: /* du = bu * au; */
957: du = PetscPowScalar(au,beta);
958: gu = coef*bu*(tu - t0);
960: td = x[k-1][j][i];
961: ad = 0.5*(t0 + td);
962: bd = PetscPowScalar(ad,bm1);
963: /* dd = bd * ad; */
964: dd = PetscPowScalar(ad,beta);
965: gd = coef*bd*(td - t0);
967: c[0].k = k-1; c[0].j = j; c[0].i = i;
968: v[0] = -hxhydhz*(dd - gd);
969: c[1].k = k; c[1].j = j-1; c[1].i = i;
970: v[1] = -hzhxdhy*(ds - gs);
971: c[2].k = k; c[2].j = j; c[2].i = i-1;
972: v[2] = -hyhzdhx*(dw - gw);
973: c[3].k = k; c[3].j = j; c[3].i = i;
974: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
975: c[4].k = k+1; c[4].j = j; c[4].i = i;
976: v[4] = -hxhydhz*(du + gu);
977: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
979: /* right-hand top up corner */
980: } else {
982: td = x[k-1][j][i];
983: ad = 0.5*(t0 + td);
984: bd = PetscPowScalar(ad,bm1);
985: /* dd = bd * ad; */
986: dd = PetscPowScalar(ad,beta);
987: gd = coef*bd*(td - t0);
989: c[0].k = k-1; c[0].j = j; c[0].i = i;
990: v[0] = -hxhydhz*(dd - gd);
991: c[1].k = k; c[1].j = j-1; c[1].i = i;
992: v[1] = -hzhxdhy*(ds - gs);
993: c[2].k = k; c[2].j = j; c[2].i = i-1;
994: v[2] = -hyhzdhx*(dw - gw);
995: c[3].k = k; c[3].j = j; c[3].i = i;
996: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
997: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
998: }
1000: } else {
1002: ts = x[k][j-1][i];
1003: as = 0.5*(t0 + ts);
1004: bs = PetscPowScalar(as,bm1);
1005: /* ds = bs * as; */
1006: ds = PetscPowScalar(as,beta);
1007: gs = coef*bs*(t0 - ts);
1009: tn = x[k][j+1][i];
1010: an = 0.5*(t0 + tn);
1011: bn = PetscPowScalar(an,bm1);
1012: /* dn = bn * an; */
1013: dn = PetscPowScalar(an,beta);
1014: gn = coef*bn*(tn - t0);
1016: /* right-hand down interior edge */
1017: if (k == 0) {
1019: tu = x[k+1][j][i];
1020: au = 0.5*(t0 + tu);
1021: bu = PetscPowScalar(au,bm1);
1022: /* du = bu * au; */
1023: du = PetscPowScalar(au,beta);
1024: gu = coef*bu*(tu - t0);
1026: c[0].k = k; c[0].j = j-1; c[0].i = i;
1027: v[0] = -hzhxdhy*(ds - gs);
1028: c[1].k = k; c[1].j = j; c[1].i = i-1;
1029: v[1] = -hyhzdhx*(dw - gw);
1030: c[2].k = k; c[2].j = j; c[2].i = i;
1031: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1032: c[3].k = k; c[3].j = j+1; c[3].i = i;
1033: v[3] = -hzhxdhy*(dn + gn);
1034: c[4].k = k+1; c[4].j = j; c[4].i = i;
1035: v[4] = -hxhydhz*(du + gu);
1036: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1038: } else if (k == mz-1) { /* right-hand up interior edge */
1040: td = x[k-1][j][i];
1041: ad = 0.5*(t0 + td);
1042: bd = PetscPowScalar(ad,bm1);
1043: /* dd = bd * ad; */
1044: dd = PetscPowScalar(ad,beta);
1045: gd = coef*bd*(t0 - td);
1047: c[0].k = k-1; c[0].j = j; c[0].i = i;
1048: v[0] = -hxhydhz*(dd - gd);
1049: c[1].k = k; c[1].j = j-1; c[1].i = i;
1050: v[1] = -hzhxdhy*(ds - gs);
1051: c[2].k = k; c[2].j = j; c[2].i = i-1;
1052: v[2] = -hyhzdhx*(dw - gw);
1053: c[3].k = k; c[3].j = j; c[3].i = i;
1054: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1055: c[4].k = k; c[4].j = j+1; c[4].i = i;
1056: v[4] = -hzhxdhy*(dn + gn);
1057: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1059: } else { /* right-hand interior plane */
1061: td = x[k-1][j][i];
1062: ad = 0.5*(t0 + td);
1063: bd = PetscPowScalar(ad,bm1);
1064: /* dd = bd * ad; */
1065: dd = PetscPowScalar(ad,beta);
1066: gd = coef*bd*(t0 - td);
1068: tu = x[k+1][j][i];
1069: au = 0.5*(t0 + tu);
1070: bu = PetscPowScalar(au,bm1);
1071: /* du = bu * au; */
1072: du = PetscPowScalar(au,beta);
1073: gu = coef*bu*(tu - t0);
1075: c[0].k = k-1; c[0].j = j; c[0].i = i;
1076: v[0] = -hxhydhz*(dd - gd);
1077: c[1].k = k; c[1].j = j-1; c[1].i = i;
1078: v[1] = -hzhxdhy*(ds - gs);
1079: c[2].k = k; c[2].j = j; c[2].i = i-1;
1080: v[2] = -hyhzdhx*(dw - gw);
1081: c[3].k = k; c[3].j = j; c[3].i = i;
1082: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1083: c[4].k = k; c[4].j = j+1; c[4].i = i;
1084: v[4] = -hzhxdhy*(dn + gn);
1085: c[5].k = k+1; c[5].j = j; c[5].i = i;
1086: v[5] = -hxhydhz*(du + gu);
1087: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1088: }
1089: }
1091: } else if (j == 0) {
1093: tw = x[k][j][i-1];
1094: aw = 0.5*(t0 + tw);
1095: bw = PetscPowScalar(aw,bm1);
1096: /* dw = bw * aw */
1097: dw = PetscPowScalar(aw,beta);
1098: gw = coef*bw*(t0 - tw);
1100: te = x[k][j][i+1];
1101: ae = 0.5*(t0 + te);
1102: be = PetscPowScalar(ae,bm1);
1103: /* de = be * ae; */
1104: de = PetscPowScalar(ae,beta);
1105: ge = coef*be*(te - t0);
1107: tn = x[k][j+1][i];
1108: an = 0.5*(t0 + tn);
1109: bn = PetscPowScalar(an,bm1);
1110: /* dn = bn * an; */
1111: dn = PetscPowScalar(an,beta);
1112: gn = coef*bn*(tn - t0);
1115: /* bottom down interior edge */
1116: if (k == 0) {
1118: tu = x[k+1][j][i];
1119: au = 0.5*(t0 + tu);
1120: bu = PetscPowScalar(au,bm1);
1121: /* du = bu * au; */
1122: du = PetscPowScalar(au,beta);
1123: gu = coef*bu*(tu - t0);
1125: c[0].k = k; c[0].j = j; c[0].i = i-1;
1126: v[0] = -hyhzdhx*(dw - gw);
1127: c[1].k = k; c[1].j = j; c[1].i = i;
1128: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1129: c[2].k = k; c[2].j = j; c[2].i = i+1;
1130: v[2] = -hyhzdhx*(de + ge);
1131: c[3].k = k; c[3].j = j+1; c[3].i = i;
1132: v[3] = -hzhxdhy*(dn + gn);
1133: c[4].k = k+1; c[4].j = j; c[4].i = i;
1134: v[4] = -hxhydhz*(du + gu);
1135: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1137: } else if (k == mz-1) { /* bottom up interior edge */
1139: td = x[k-1][j][i];
1140: ad = 0.5*(t0 + td);
1141: bd = PetscPowScalar(ad,bm1);
1142: /* dd = bd * ad; */
1143: dd = PetscPowScalar(ad,beta);
1144: gd = coef*bd*(td - t0);
1146: c[0].k = k-1; c[0].j = j; c[0].i = i;
1147: v[0] = -hxhydhz*(dd - gd);
1148: c[1].k = k; c[1].j = j; c[1].i = i-1;
1149: v[1] = -hyhzdhx*(dw - gw);
1150: c[2].k = k; c[2].j = j; c[2].i = i;
1151: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1152: c[3].k = k; c[3].j = j; c[3].i = i+1;
1153: v[3] = -hyhzdhx*(de + ge);
1154: c[4].k = k; c[4].j = j+1; c[4].i = i;
1155: v[4] = -hzhxdhy*(dn + gn);
1156: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1158: } else { /* bottom interior plane */
1160: tu = x[k+1][j][i];
1161: au = 0.5*(t0 + tu);
1162: bu = PetscPowScalar(au,bm1);
1163: /* du = bu * au; */
1164: du = PetscPowScalar(au,beta);
1165: gu = coef*bu*(tu - t0);
1167: td = x[k-1][j][i];
1168: ad = 0.5*(t0 + td);
1169: bd = PetscPowScalar(ad,bm1);
1170: /* dd = bd * ad; */
1171: dd = PetscPowScalar(ad,beta);
1172: gd = coef*bd*(td - t0);
1174: c[0].k = k-1; c[0].j = j; c[0].i = i;
1175: v[0] = -hxhydhz*(dd - gd);
1176: c[1].k = k; c[1].j = j; c[1].i = i-1;
1177: v[1] = -hyhzdhx*(dw - gw);
1178: c[2].k = k; c[2].j = j; c[2].i = i;
1179: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1180: c[3].k = k; c[3].j = j; c[3].i = i+1;
1181: v[3] = -hyhzdhx*(de + ge);
1182: c[4].k = k; c[4].j = j+1; c[4].i = i;
1183: v[4] = -hzhxdhy*(dn + gn);
1184: c[5].k = k+1; c[5].j = j; c[5].i = i;
1185: v[5] = -hxhydhz*(du + gu);
1186: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1187: }
1189: } else if (j == my-1) {
1191: tw = x[k][j][i-1];
1192: aw = 0.5*(t0 + tw);
1193: bw = PetscPowScalar(aw,bm1);
1194: /* dw = bw * aw */
1195: dw = PetscPowScalar(aw,beta);
1196: gw = coef*bw*(t0 - tw);
1198: te = x[k][j][i+1];
1199: ae = 0.5*(t0 + te);
1200: be = PetscPowScalar(ae,bm1);
1201: /* de = be * ae; */
1202: de = PetscPowScalar(ae,beta);
1203: ge = coef*be*(te - t0);
1205: ts = x[k][j-1][i];
1206: as = 0.5*(t0 + ts);
1207: bs = PetscPowScalar(as,bm1);
1208: /* ds = bs * as; */
1209: ds = PetscPowScalar(as,beta);
1210: gs = coef*bs*(t0 - ts);
1212: /* top down interior edge */
1213: if (k == 0) {
1215: tu = x[k+1][j][i];
1216: au = 0.5*(t0 + tu);
1217: bu = PetscPowScalar(au,bm1);
1218: /* du = bu * au; */
1219: du = PetscPowScalar(au,beta);
1220: gu = coef*bu*(tu - t0);
1222: c[0].k = k; c[0].j = j-1; c[0].i = i;
1223: v[0] = -hzhxdhy*(ds - gs);
1224: c[1].k = k; c[1].j = j; c[1].i = i-1;
1225: v[1] = -hyhzdhx*(dw - gw);
1226: c[2].k = k; c[2].j = j; c[2].i = i;
1227: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1228: c[3].k = k; c[3].j = j; c[3].i = i+1;
1229: v[3] = -hyhzdhx*(de + ge);
1230: c[4].k = k+1; c[4].j = j; c[4].i = i;
1231: v[4] = -hxhydhz*(du + gu);
1232: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1234: } else if (k == mz-1) { /* top up interior edge */
1236: td = x[k-1][j][i];
1237: ad = 0.5*(t0 + td);
1238: bd = PetscPowScalar(ad,bm1);
1239: /* dd = bd * ad; */
1240: dd = PetscPowScalar(ad,beta);
1241: gd = coef*bd*(td - t0);
1243: c[0].k = k-1; c[0].j = j; c[0].i = i;
1244: v[0] = -hxhydhz*(dd - gd);
1245: c[1].k = k; c[1].j = j-1; c[1].i = i;
1246: v[1] = -hzhxdhy*(ds - gs);
1247: c[2].k = k; c[2].j = j; c[2].i = i-1;
1248: v[2] = -hyhzdhx*(dw - gw);
1249: c[3].k = k; c[3].j = j; c[3].i = i;
1250: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1251: c[4].k = k; c[4].j = j; c[4].i = i+1;
1252: v[4] = -hyhzdhx*(de + ge);
1253: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1255: } else { /* top interior plane */
1257: tu = x[k+1][j][i];
1258: au = 0.5*(t0 + tu);
1259: bu = PetscPowScalar(au,bm1);
1260: /* du = bu * au; */
1261: du = PetscPowScalar(au,beta);
1262: gu = coef*bu*(tu - t0);
1264: td = x[k-1][j][i];
1265: ad = 0.5*(t0 + td);
1266: bd = PetscPowScalar(ad,bm1);
1267: /* dd = bd * ad; */
1268: dd = PetscPowScalar(ad,beta);
1269: gd = coef*bd*(td - t0);
1271: c[0].k = k-1; c[0].j = j; c[0].i = i;
1272: v[0] = -hxhydhz*(dd - gd);
1273: c[1].k = k; c[1].j = j-1; c[1].i = i;
1274: v[1] = -hzhxdhy*(ds - gs);
1275: c[2].k = k; c[2].j = j; c[2].i = i-1;
1276: v[2] = -hyhzdhx*(dw - gw);
1277: c[3].k = k; c[3].j = j; c[3].i = i;
1278: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1279: c[4].k = k; c[4].j = j; c[4].i = i+1;
1280: v[4] = -hyhzdhx*(de + ge);
1281: c[5].k = k+1; c[5].j = j; c[5].i = i;
1282: v[5] = -hxhydhz*(du + gu);
1283: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1284: }
1286: } else if (k == 0) {
1288: /* down interior plane */
1290: tw = x[k][j][i-1];
1291: aw = 0.5*(t0 + tw);
1292: bw = PetscPowScalar(aw,bm1);
1293: /* dw = bw * aw */
1294: dw = PetscPowScalar(aw,beta);
1295: gw = coef*bw*(t0 - tw);
1297: te = x[k][j][i+1];
1298: ae = 0.5*(t0 + te);
1299: be = PetscPowScalar(ae,bm1);
1300: /* de = be * ae; */
1301: de = PetscPowScalar(ae,beta);
1302: ge = coef*be*(te - t0);
1304: ts = x[k][j-1][i];
1305: as = 0.5*(t0 + ts);
1306: bs = PetscPowScalar(as,bm1);
1307: /* ds = bs * as; */
1308: ds = PetscPowScalar(as,beta);
1309: gs = coef*bs*(t0 - ts);
1311: tn = x[k][j+1][i];
1312: an = 0.5*(t0 + tn);
1313: bn = PetscPowScalar(an,bm1);
1314: /* dn = bn * an; */
1315: dn = PetscPowScalar(an,beta);
1316: gn = coef*bn*(tn - t0);
1318: tu = x[k+1][j][i];
1319: au = 0.5*(t0 + tu);
1320: bu = PetscPowScalar(au,bm1);
1321: /* du = bu * au; */
1322: du = PetscPowScalar(au,beta);
1323: gu = coef*bu*(tu - t0);
1325: c[0].k = k; c[0].j = j-1; c[0].i = i;
1326: v[0] = -hzhxdhy*(ds - gs);
1327: c[1].k = k; c[1].j = j; c[1].i = i-1;
1328: v[1] = -hyhzdhx*(dw - gw);
1329: c[2].k = k; c[2].j = j; c[2].i = i;
1330: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1331: c[3].k = k; c[3].j = j; c[3].i = i+1;
1332: v[3] = -hyhzdhx*(de + ge);
1333: c[4].k = k; c[4].j = j+1; c[4].i = i;
1334: v[4] = -hzhxdhy*(dn + gn);
1335: c[5].k = k+1; c[5].j = j; c[5].i = i;
1336: v[5] = -hxhydhz*(du + gu);
1337: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1339: } else if (k == mz-1) {
1341: /* up interior plane */
1343: tw = x[k][j][i-1];
1344: aw = 0.5*(t0 + tw);
1345: bw = PetscPowScalar(aw,bm1);
1346: /* dw = bw * aw */
1347: dw = PetscPowScalar(aw,beta);
1348: gw = coef*bw*(t0 - tw);
1350: te = x[k][j][i+1];
1351: ae = 0.5*(t0 + te);
1352: be = PetscPowScalar(ae,bm1);
1353: /* de = be * ae; */
1354: de = PetscPowScalar(ae,beta);
1355: ge = coef*be*(te - t0);
1357: ts = x[k][j-1][i];
1358: as = 0.5*(t0 + ts);
1359: bs = PetscPowScalar(as,bm1);
1360: /* ds = bs * as; */
1361: ds = PetscPowScalar(as,beta);
1362: gs = coef*bs*(t0 - ts);
1364: tn = x[k][j+1][i];
1365: an = 0.5*(t0 + tn);
1366: bn = PetscPowScalar(an,bm1);
1367: /* dn = bn * an; */
1368: dn = PetscPowScalar(an,beta);
1369: gn = coef*bn*(tn - t0);
1371: td = x[k-1][j][i];
1372: ad = 0.5*(t0 + td);
1373: bd = PetscPowScalar(ad,bm1);
1374: /* dd = bd * ad; */
1375: dd = PetscPowScalar(ad,beta);
1376: gd = coef*bd*(t0 - td);
1378: c[0].k = k-1; c[0].j = j; c[0].i = i;
1379: v[0] = -hxhydhz*(dd - gd);
1380: c[1].k = k; c[1].j = j-1; c[1].i = i;
1381: v[1] = -hzhxdhy*(ds - gs);
1382: c[2].k = k; c[2].j = j; c[2].i = i-1;
1383: v[2] = -hyhzdhx*(dw - gw);
1384: c[3].k = k; c[3].j = j; c[3].i = i;
1385: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1386: c[4].k = k; c[4].j = j; c[4].i = i+1;
1387: v[4] = -hyhzdhx*(de + ge);
1388: c[5].k = k; c[5].j = j+1; c[5].i = i;
1389: v[5] = -hzhxdhy*(dn + gn);
1390: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1391: }
1392: }
1393: }
1394: }
1395: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
1396: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
1397: if (jac != *J) {
1398: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1399: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1400: }
1401: DMDAVecRestoreArray(da,localX,&x);
1402: DMRestoreLocalVector(da,&localX);
1404: PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
1405: return(0);
1406: }