Actual source code: ex20.c
petsc-3.6.1 2015-08-06
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*);
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);
72: /* set problem parameters */
73: user.tleft = 1.0;
74: user.tright = 0.1;
75: user.beta = 2.5;
76: PetscOptionsGetReal(NULL,"-tleft",&user.tleft,NULL);
77: PetscOptionsGetReal(NULL,"-tright",&user.tright,NULL);
78: PetscOptionsGetReal(NULL,"-beta",&user.beta,NULL);
79: user.bm1 = user.beta - 1.0;
80: user.coef = user.beta/2.0;
82: /*
83: Set the DMDA (grid structure) for the grids.
84: */
85: 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);
86: DMSetApplicationContext(da,&user);
88: /*
89: Create the nonlinear solver
90: */
91: SNESCreate(PETSC_COMM_WORLD,&snes);
92: SNESSetDM(snes,da);
93: SNESSetFunction(snes,NULL,FormFunction,&user);
94: SNESSetJacobian(snes,NULL,NULL,FormJacobian,&user);
95: SNESSetFromOptions(snes);
96: SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);
98: SNESSolve(snes,NULL,NULL);
99: SNESGetIterationNumber(snes,&its);
100: SNESGetLinearSolveIterations(snes,&lits);
101: litspit = ((PetscReal)lits)/((PetscReal)its);
102: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
103: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
104: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",(double)litspit);
106: SNESDestroy(&snes);
107: DMDestroy(&da);
108: PetscFinalize();
110: return 0;
111: }
112: /* -------------------- Form initial approximation ----------------- */
115: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
116: {
117: AppCtx *user;
118: PetscInt i,j,k,xs,ys,xm,ym,zs,zm;
120: PetscScalar ***x;
121: DM da;
124: SNESGetDM(snes,&da);
125: DMGetApplicationContext(da,&user);
126: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
127: DMDAVecGetArray(da,X,&x);
129: /* Compute initial guess */
130: for (k=zs; k<zs+zm; k++) {
131: for (j=ys; j<ys+ym; j++) {
132: for (i=xs; i<xs+xm; i++) {
133: x[k][j][i] = user->tleft;
134: }
135: }
136: }
137: DMDAVecRestoreArray(da,X,&x);
138: return(0);
139: }
140: /* -------------------- Evaluate Function F(x) --------------------- */
143: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
144: {
145: AppCtx *user = (AppCtx*)ptr;
147: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
148: PetscScalar zero = 0.0,one = 1.0;
149: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
150: 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;
151: PetscScalar tleft,tright,beta,td,ad,dd,fd=0.0,tu,au,du=0.0,fu=0.0;
152: PetscScalar ***x,***f;
153: Vec localX;
154: DM da;
157: SNESGetDM(snes,&da);
158: DMGetLocalVector(da,&localX);
159: DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
160: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
161: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
162: tleft = user->tleft; tright = user->tright;
163: beta = user->beta;
165: /* Get ghost points */
166: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
167: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
168: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
169: DMDAVecGetArray(da,localX,&x);
170: DMDAVecGetArray(da,F,&f);
172: /* Evaluate function */
173: for (k=zs; k<zs+zm; k++) {
174: for (j=ys; j<ys+ym; j++) {
175: for (i=xs; i<xs+xm; i++) {
176: t0 = x[k][j][i];
178: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
180: /* general interior volume */
182: tw = x[k][j][i-1];
183: aw = 0.5*(t0 + tw);
184: dw = PetscPowScalar(aw,beta);
185: fw = dw*(t0 - tw);
187: te = x[k][j][i+1];
188: ae = 0.5*(t0 + te);
189: de = PetscPowScalar(ae,beta);
190: fe = de*(te - t0);
192: ts = x[k][j-1][i];
193: as = 0.5*(t0 + ts);
194: ds = PetscPowScalar(as,beta);
195: fs = ds*(t0 - ts);
197: tn = x[k][j+1][i];
198: an = 0.5*(t0 + tn);
199: dn = PetscPowScalar(an,beta);
200: fn = dn*(tn - t0);
202: td = x[k-1][j][i];
203: ad = 0.5*(t0 + td);
204: dd = PetscPowScalar(ad,beta);
205: fd = dd*(t0 - td);
207: tu = x[k+1][j][i];
208: au = 0.5*(t0 + tu);
209: du = PetscPowScalar(au,beta);
210: fu = du*(tu - t0);
212: } else if (i == 0) {
214: /* left-hand (west) boundary */
215: tw = tleft;
216: aw = 0.5*(t0 + tw);
217: dw = PetscPowScalar(aw,beta);
218: fw = dw*(t0 - tw);
220: te = x[k][j][i+1];
221: ae = 0.5*(t0 + te);
222: de = PetscPowScalar(ae,beta);
223: fe = de*(te - t0);
225: if (j > 0) {
226: ts = x[k][j-1][i];
227: as = 0.5*(t0 + ts);
228: ds = PetscPowScalar(as,beta);
229: fs = ds*(t0 - ts);
230: } else {
231: fs = zero;
232: }
234: if (j < my-1) {
235: tn = x[k][j+1][i];
236: an = 0.5*(t0 + tn);
237: dn = PetscPowScalar(an,beta);
238: fn = dn*(tn - t0);
239: } else {
240: fn = zero;
241: }
243: if (k > 0) {
244: td = x[k-1][j][i];
245: ad = 0.5*(t0 + td);
246: dd = PetscPowScalar(ad,beta);
247: fd = dd*(t0 - td);
248: } else {
249: fd = zero;
250: }
252: if (k < mz-1) {
253: tu = x[k+1][j][i];
254: au = 0.5*(t0 + tu);
255: du = PetscPowScalar(au,beta);
256: fu = du*(tu - t0);
257: } else {
258: fu = zero;
259: }
261: } else if (i == mx-1) {
263: /* right-hand (east) boundary */
264: tw = x[k][j][i-1];
265: aw = 0.5*(t0 + tw);
266: dw = PetscPowScalar(aw,beta);
267: fw = dw*(t0 - tw);
269: te = tright;
270: ae = 0.5*(t0 + te);
271: de = PetscPowScalar(ae,beta);
272: fe = de*(te - t0);
274: if (j > 0) {
275: ts = x[k][j-1][i];
276: as = 0.5*(t0 + ts);
277: ds = PetscPowScalar(as,beta);
278: fs = ds*(t0 - ts);
279: } else {
280: fs = zero;
281: }
283: if (j < my-1) {
284: tn = x[k][j+1][i];
285: an = 0.5*(t0 + tn);
286: dn = PetscPowScalar(an,beta);
287: fn = dn*(tn - t0);
288: } else {
289: fn = zero;
290: }
292: if (k > 0) {
293: td = x[k-1][j][i];
294: ad = 0.5*(t0 + td);
295: dd = PetscPowScalar(ad,beta);
296: fd = dd*(t0 - td);
297: } else {
298: fd = zero;
299: }
301: if (k < mz-1) {
302: tu = x[k+1][j][i];
303: au = 0.5*(t0 + tu);
304: du = PetscPowScalar(au,beta);
305: fu = du*(tu - t0);
306: } else {
307: fu = zero;
308: }
310: } else if (j == 0) {
312: /* bottom (south) boundary, and i <> 0 or mx-1 */
313: tw = x[k][j][i-1];
314: aw = 0.5*(t0 + tw);
315: dw = PetscPowScalar(aw,beta);
316: fw = dw*(t0 - tw);
318: te = x[k][j][i+1];
319: ae = 0.5*(t0 + te);
320: de = PetscPowScalar(ae,beta);
321: fe = de*(te - t0);
323: fs = zero;
325: tn = x[k][j+1][i];
326: an = 0.5*(t0 + tn);
327: dn = PetscPowScalar(an,beta);
328: fn = dn*(tn - t0);
330: if (k > 0) {
331: td = x[k-1][j][i];
332: ad = 0.5*(t0 + td);
333: dd = PetscPowScalar(ad,beta);
334: fd = dd*(t0 - td);
335: } else {
336: fd = zero;
337: }
339: if (k < mz-1) {
340: tu = x[k+1][j][i];
341: au = 0.5*(t0 + tu);
342: du = PetscPowScalar(au,beta);
343: fu = du*(tu - t0);
344: } else {
345: fu = zero;
346: }
348: } else if (j == my-1) {
350: /* top (north) boundary, and i <> 0 or mx-1 */
351: tw = x[k][j][i-1];
352: aw = 0.5*(t0 + tw);
353: dw = PetscPowScalar(aw,beta);
354: fw = dw*(t0 - tw);
356: te = x[k][j][i+1];
357: ae = 0.5*(t0 + te);
358: de = PetscPowScalar(ae,beta);
359: fe = de*(te - t0);
361: ts = x[k][j-1][i];
362: as = 0.5*(t0 + ts);
363: ds = PetscPowScalar(as,beta);
364: fs = ds*(t0 - ts);
366: fn = zero;
368: if (k > 0) {
369: td = x[k-1][j][i];
370: ad = 0.5*(t0 + td);
371: dd = PetscPowScalar(ad,beta);
372: fd = dd*(t0 - td);
373: } else {
374: fd = zero;
375: }
377: if (k < mz-1) {
378: tu = x[k+1][j][i];
379: au = 0.5*(t0 + tu);
380: du = PetscPowScalar(au,beta);
381: fu = du*(tu - t0);
382: } else {
383: fu = zero;
384: }
386: } else if (k == 0) {
388: /* down boundary (interior only) */
389: tw = x[k][j][i-1];
390: aw = 0.5*(t0 + tw);
391: dw = PetscPowScalar(aw,beta);
392: fw = dw*(t0 - tw);
394: te = x[k][j][i+1];
395: ae = 0.5*(t0 + te);
396: de = PetscPowScalar(ae,beta);
397: fe = de*(te - t0);
399: ts = x[k][j-1][i];
400: as = 0.5*(t0 + ts);
401: ds = PetscPowScalar(as,beta);
402: fs = ds*(t0 - ts);
404: tn = x[k][j+1][i];
405: an = 0.5*(t0 + tn);
406: dn = PetscPowScalar(an,beta);
407: fn = dn*(tn - t0);
409: fd = zero;
411: tu = x[k+1][j][i];
412: au = 0.5*(t0 + tu);
413: du = PetscPowScalar(au,beta);
414: fu = du*(tu - t0);
416: } else if (k == mz-1) {
418: /* up boundary (interior only) */
419: tw = x[k][j][i-1];
420: aw = 0.5*(t0 + tw);
421: dw = PetscPowScalar(aw,beta);
422: fw = dw*(t0 - tw);
424: te = x[k][j][i+1];
425: ae = 0.5*(t0 + te);
426: de = PetscPowScalar(ae,beta);
427: fe = de*(te - t0);
429: ts = x[k][j-1][i];
430: as = 0.5*(t0 + ts);
431: ds = PetscPowScalar(as,beta);
432: fs = ds*(t0 - ts);
434: tn = x[k][j+1][i];
435: an = 0.5*(t0 + tn);
436: dn = PetscPowScalar(an,beta);
437: fn = dn*(tn - t0);
439: td = x[k-1][j][i];
440: ad = 0.5*(t0 + td);
441: dd = PetscPowScalar(ad,beta);
442: fd = dd*(t0 - td);
444: fu = zero;
445: }
447: f[k][j][i] = -hyhzdhx*(fe-fw) - hzhxdhy*(fn-fs) - hxhydhz*(fu-fd);
448: }
449: }
450: }
451: DMDAVecRestoreArray(da,localX,&x);
452: DMDAVecRestoreArray(da,F,&f);
453: DMRestoreLocalVector(da,&localX);
454: PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
455: return(0);
456: }
457: /* -------------------- Evaluate Jacobian F(x) --------------------- */
460: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
461: {
462: AppCtx *user = (AppCtx*)ptr;
464: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
465: PetscScalar one = 1.0;
466: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
467: PetscScalar t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw;
468: PetscScalar tleft,tright,beta,td,ad,dd,tu,au,du,v[7],bm1,coef;
469: PetscScalar ***x,bn,bs,be,bw,bu,bd,gn,gs,ge,gw,gu,gd;
470: Vec localX;
471: MatStencil c[7],row;
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: }