Actual source code: ex20.c
petsc-3.8.4 2018-03-24
2: static char help[] ="Nonlinear Radiative Transport PDE with multigrid in 3d.\n\
3: Uses 3-dimensional distributed arrays.\n\
4: A 3-dim simplified Radiative Transport test problem is used, with analytic Jacobian. \n\
5: \n\
6: Solves the linear systems via multilevel methods \n\
7: \n\
8: The command line\n\
9: options are:\n\
10: -tleft <tl>, where <tl> indicates the left Diriclet BC \n\
11: -tright <tr>, where <tr> indicates the right Diriclet BC \n\
12: -beta <beta>, where <beta> indicates the exponent in T \n\n";
14: /*T
15: Concepts: SNES^solving a system of nonlinear equations
16: Concepts: DM^using distributed arrays
17: Concepts: multigrid;
18: Processors: n
19: T*/
21: /*
23: This example models the partial differential equation
25: - Div(alpha* T^beta (GRAD T)) = 0.
27: where beta = 2.5 and alpha = 1.0
29: BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = dT/dn_up = dT/dn_down = 0.
31: in the unit square, which is uniformly discretized in each of x and
32: y in this simple encoding. The degrees of freedom are cell centered.
34: A finite volume approximation with the usual 7-point stencil
35: is used to discretize the boundary value problem to obtain a
36: nonlinear system of equations.
38: This code was contributed by Nickolas Jovanovic based on ex18.c
40: */
42: #include <petscsnes.h>
43: #include <petscdm.h>
44: #include <petscdmda.h>
46: /* User-defined application context */
48: typedef struct {
49: PetscReal tleft,tright; /* Dirichlet boundary conditions */
50: PetscReal beta,bm1,coef; /* nonlinear diffusivity parameterizations */
51: } AppCtx;
53: #define POWFLOP 5 /* assume a pow() takes five flops */
55: extern PetscErrorCode FormInitialGuess(SNES,Vec,void*);
56: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
57: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
59: int main(int argc,char **argv)
60: {
61: SNES snes;
62: AppCtx user;
64: PetscInt its,lits;
65: PetscReal litspit;
66: DM da;
68: PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
69: /* set problem parameters */
70: user.tleft = 1.0;
71: user.tright = 0.1;
72: user.beta = 2.5;
73: PetscOptionsGetReal(NULL,NULL,"-tleft",&user.tleft,NULL);
74: PetscOptionsGetReal(NULL,NULL,"-tright",&user.tright,NULL);
75: PetscOptionsGetReal(NULL,NULL,"-beta",&user.beta,NULL);
76: user.bm1 = user.beta - 1.0;
77: user.coef = user.beta/2.0;
79: /*
80: Set the DMDA (grid structure) for the grids.
81: */
82: 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);
83: DMSetFromOptions(da);
84: DMSetUp(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",(double)litspit);
105: SNESDestroy(&snes);
106: DMDestroy(&da);
107: PetscFinalize();
108: return ierr;
109: }
110: /* -------------------- Form initial approximation ----------------- */
111: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
112: {
113: AppCtx *user;
114: PetscInt i,j,k,xs,ys,xm,ym,zs,zm;
116: PetscScalar ***x;
117: DM da;
120: SNESGetDM(snes,&da);
121: DMGetApplicationContext(da,&user);
122: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
123: DMDAVecGetArray(da,X,&x);
125: /* Compute initial guess */
126: for (k=zs; k<zs+zm; k++) {
127: for (j=ys; j<ys+ym; j++) {
128: for (i=xs; i<xs+xm; i++) {
129: x[k][j][i] = user->tleft;
130: }
131: }
132: }
133: DMDAVecRestoreArray(da,X,&x);
134: return(0);
135: }
136: /* -------------------- Evaluate Function F(x) --------------------- */
137: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
138: {
139: AppCtx *user = (AppCtx*)ptr;
141: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
142: PetscScalar zero = 0.0,one = 1.0;
143: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
144: 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;
145: PetscScalar tleft,tright,beta,td,ad,dd,fd=0.0,tu,au,du=0.0,fu=0.0;
146: PetscScalar ***x,***f;
147: Vec localX;
148: DM da;
151: SNESGetDM(snes,&da);
152: DMGetLocalVector(da,&localX);
153: DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
154: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
155: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
156: tleft = user->tleft; tright = user->tright;
157: beta = user->beta;
159: /* Get ghost points */
160: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
161: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
162: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
163: DMDAVecGetArray(da,localX,&x);
164: DMDAVecGetArray(da,F,&f);
166: /* Evaluate function */
167: for (k=zs; k<zs+zm; k++) {
168: for (j=ys; j<ys+ym; j++) {
169: for (i=xs; i<xs+xm; i++) {
170: t0 = x[k][j][i];
172: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
174: /* general interior volume */
176: tw = x[k][j][i-1];
177: aw = 0.5*(t0 + tw);
178: dw = PetscPowScalar(aw,beta);
179: fw = dw*(t0 - tw);
181: te = x[k][j][i+1];
182: ae = 0.5*(t0 + te);
183: de = PetscPowScalar(ae,beta);
184: fe = de*(te - t0);
186: ts = x[k][j-1][i];
187: as = 0.5*(t0 + ts);
188: ds = PetscPowScalar(as,beta);
189: fs = ds*(t0 - ts);
191: tn = x[k][j+1][i];
192: an = 0.5*(t0 + tn);
193: dn = PetscPowScalar(an,beta);
194: fn = dn*(tn - t0);
196: td = x[k-1][j][i];
197: ad = 0.5*(t0 + td);
198: dd = PetscPowScalar(ad,beta);
199: fd = dd*(t0 - td);
201: tu = x[k+1][j][i];
202: au = 0.5*(t0 + tu);
203: du = PetscPowScalar(au,beta);
204: fu = du*(tu - t0);
206: } else if (i == 0) {
208: /* left-hand (west) boundary */
209: tw = tleft;
210: aw = 0.5*(t0 + tw);
211: dw = PetscPowScalar(aw,beta);
212: fw = dw*(t0 - tw);
214: te = x[k][j][i+1];
215: ae = 0.5*(t0 + te);
216: de = PetscPowScalar(ae,beta);
217: fe = de*(te - t0);
219: if (j > 0) {
220: ts = x[k][j-1][i];
221: as = 0.5*(t0 + ts);
222: ds = PetscPowScalar(as,beta);
223: fs = ds*(t0 - ts);
224: } else {
225: fs = zero;
226: }
228: if (j < my-1) {
229: tn = x[k][j+1][i];
230: an = 0.5*(t0 + tn);
231: dn = PetscPowScalar(an,beta);
232: fn = dn*(tn - t0);
233: } else {
234: fn = zero;
235: }
237: if (k > 0) {
238: td = x[k-1][j][i];
239: ad = 0.5*(t0 + td);
240: dd = PetscPowScalar(ad,beta);
241: fd = dd*(t0 - td);
242: } else {
243: fd = zero;
244: }
246: if (k < mz-1) {
247: tu = x[k+1][j][i];
248: au = 0.5*(t0 + tu);
249: du = PetscPowScalar(au,beta);
250: fu = du*(tu - t0);
251: } else {
252: fu = zero;
253: }
255: } else if (i == mx-1) {
257: /* right-hand (east) boundary */
258: tw = x[k][j][i-1];
259: aw = 0.5*(t0 + tw);
260: dw = PetscPowScalar(aw,beta);
261: fw = dw*(t0 - tw);
263: te = tright;
264: ae = 0.5*(t0 + te);
265: de = PetscPowScalar(ae,beta);
266: fe = de*(te - t0);
268: if (j > 0) {
269: ts = x[k][j-1][i];
270: as = 0.5*(t0 + ts);
271: ds = PetscPowScalar(as,beta);
272: fs = ds*(t0 - ts);
273: } else {
274: fs = zero;
275: }
277: if (j < my-1) {
278: tn = x[k][j+1][i];
279: an = 0.5*(t0 + tn);
280: dn = PetscPowScalar(an,beta);
281: fn = dn*(tn - t0);
282: } else {
283: fn = zero;
284: }
286: if (k > 0) {
287: td = x[k-1][j][i];
288: ad = 0.5*(t0 + td);
289: dd = PetscPowScalar(ad,beta);
290: fd = dd*(t0 - td);
291: } else {
292: fd = zero;
293: }
295: if (k < mz-1) {
296: tu = x[k+1][j][i];
297: au = 0.5*(t0 + tu);
298: du = PetscPowScalar(au,beta);
299: fu = du*(tu - t0);
300: } else {
301: fu = zero;
302: }
304: } else if (j == 0) {
306: /* bottom (south) boundary, and i <> 0 or mx-1 */
307: tw = x[k][j][i-1];
308: aw = 0.5*(t0 + tw);
309: dw = PetscPowScalar(aw,beta);
310: fw = dw*(t0 - tw);
312: te = x[k][j][i+1];
313: ae = 0.5*(t0 + te);
314: de = PetscPowScalar(ae,beta);
315: fe = de*(te - t0);
317: fs = zero;
319: tn = x[k][j+1][i];
320: an = 0.5*(t0 + tn);
321: dn = PetscPowScalar(an,beta);
322: fn = dn*(tn - t0);
324: if (k > 0) {
325: td = x[k-1][j][i];
326: ad = 0.5*(t0 + td);
327: dd = PetscPowScalar(ad,beta);
328: fd = dd*(t0 - td);
329: } else {
330: fd = zero;
331: }
333: if (k < mz-1) {
334: tu = x[k+1][j][i];
335: au = 0.5*(t0 + tu);
336: du = PetscPowScalar(au,beta);
337: fu = du*(tu - t0);
338: } else {
339: fu = zero;
340: }
342: } else if (j == my-1) {
344: /* top (north) boundary, and i <> 0 or mx-1 */
345: tw = x[k][j][i-1];
346: aw = 0.5*(t0 + tw);
347: dw = PetscPowScalar(aw,beta);
348: fw = dw*(t0 - tw);
350: te = x[k][j][i+1];
351: ae = 0.5*(t0 + te);
352: de = PetscPowScalar(ae,beta);
353: fe = de*(te - t0);
355: ts = x[k][j-1][i];
356: as = 0.5*(t0 + ts);
357: ds = PetscPowScalar(as,beta);
358: fs = ds*(t0 - ts);
360: fn = zero;
362: if (k > 0) {
363: td = x[k-1][j][i];
364: ad = 0.5*(t0 + td);
365: dd = PetscPowScalar(ad,beta);
366: fd = dd*(t0 - td);
367: } else {
368: fd = zero;
369: }
371: if (k < mz-1) {
372: tu = x[k+1][j][i];
373: au = 0.5*(t0 + tu);
374: du = PetscPowScalar(au,beta);
375: fu = du*(tu - t0);
376: } else {
377: fu = zero;
378: }
380: } else if (k == 0) {
382: /* down boundary (interior only) */
383: tw = x[k][j][i-1];
384: aw = 0.5*(t0 + tw);
385: dw = PetscPowScalar(aw,beta);
386: fw = dw*(t0 - tw);
388: te = x[k][j][i+1];
389: ae = 0.5*(t0 + te);
390: de = PetscPowScalar(ae,beta);
391: fe = de*(te - t0);
393: ts = x[k][j-1][i];
394: as = 0.5*(t0 + ts);
395: ds = PetscPowScalar(as,beta);
396: fs = ds*(t0 - ts);
398: tn = x[k][j+1][i];
399: an = 0.5*(t0 + tn);
400: dn = PetscPowScalar(an,beta);
401: fn = dn*(tn - t0);
403: fd = zero;
405: tu = x[k+1][j][i];
406: au = 0.5*(t0 + tu);
407: du = PetscPowScalar(au,beta);
408: fu = du*(tu - t0);
410: } else if (k == mz-1) {
412: /* up boundary (interior only) */
413: tw = x[k][j][i-1];
414: aw = 0.5*(t0 + tw);
415: dw = PetscPowScalar(aw,beta);
416: fw = dw*(t0 - tw);
418: te = x[k][j][i+1];
419: ae = 0.5*(t0 + te);
420: de = PetscPowScalar(ae,beta);
421: fe = de*(te - t0);
423: ts = x[k][j-1][i];
424: as = 0.5*(t0 + ts);
425: ds = PetscPowScalar(as,beta);
426: fs = ds*(t0 - ts);
428: tn = x[k][j+1][i];
429: an = 0.5*(t0 + tn);
430: dn = PetscPowScalar(an,beta);
431: fn = dn*(tn - t0);
433: td = x[k-1][j][i];
434: ad = 0.5*(t0 + td);
435: dd = PetscPowScalar(ad,beta);
436: fd = dd*(t0 - td);
438: fu = zero;
439: }
441: f[k][j][i] = -hyhzdhx*(fe-fw) - hzhxdhy*(fn-fs) - hxhydhz*(fu-fd);
442: }
443: }
444: }
445: DMDAVecRestoreArray(da,localX,&x);
446: DMDAVecRestoreArray(da,F,&f);
447: DMRestoreLocalVector(da,&localX);
448: PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
449: return(0);
450: }
451: /* -------------------- Evaluate Jacobian F(x) --------------------- */
452: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
453: {
454: AppCtx *user = (AppCtx*)ptr;
456: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
457: PetscScalar one = 1.0;
458: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
459: PetscScalar t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw;
460: PetscScalar tleft,tright,beta,td,ad,dd,tu,au,du,v[7],bm1,coef;
461: PetscScalar ***x,bn,bs,be,bw,bu,bd,gn,gs,ge,gw,gu,gd;
462: Vec localX;
463: MatStencil c[7],row;
464: DM da;
467: SNESGetDM(snes,&da);
468: DMGetLocalVector(da,&localX);
469: DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
470: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
471: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
472: tleft = user->tleft; tright = user->tright;
473: beta = user->beta; bm1 = user->bm1; coef = user->coef;
475: /* Get ghost points */
476: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
477: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
478: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
479: DMDAVecGetArray(da,localX,&x);
481: /* Evaluate Jacobian of function */
482: for (k=zs; k<zs+zm; k++) {
483: for (j=ys; j<ys+ym; j++) {
484: for (i=xs; i<xs+xm; i++) {
485: t0 = x[k][j][i];
486: row.k = k; row.j = j; row.i = i;
487: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
489: /* general interior volume */
491: tw = x[k][j][i-1];
492: aw = 0.5*(t0 + tw);
493: bw = PetscPowScalar(aw,bm1);
494: /* dw = bw * aw */
495: dw = PetscPowScalar(aw,beta);
496: gw = coef*bw*(t0 - tw);
498: te = x[k][j][i+1];
499: ae = 0.5*(t0 + te);
500: be = PetscPowScalar(ae,bm1);
501: /* de = be * ae; */
502: de = PetscPowScalar(ae,beta);
503: ge = coef*be*(te - t0);
505: ts = x[k][j-1][i];
506: as = 0.5*(t0 + ts);
507: bs = PetscPowScalar(as,bm1);
508: /* ds = bs * as; */
509: ds = PetscPowScalar(as,beta);
510: gs = coef*bs*(t0 - ts);
512: tn = x[k][j+1][i];
513: an = 0.5*(t0 + tn);
514: bn = PetscPowScalar(an,bm1);
515: /* dn = bn * an; */
516: dn = PetscPowScalar(an,beta);
517: gn = coef*bn*(tn - t0);
519: td = x[k-1][j][i];
520: ad = 0.5*(t0 + td);
521: bd = PetscPowScalar(ad,bm1);
522: /* dd = bd * ad; */
523: dd = PetscPowScalar(ad,beta);
524: gd = coef*bd*(t0 - td);
526: tu = x[k+1][j][i];
527: au = 0.5*(t0 + tu);
528: bu = PetscPowScalar(au,bm1);
529: /* du = bu * au; */
530: du = PetscPowScalar(au,beta);
531: gu = coef*bu*(tu - t0);
533: c[0].k = k-1; c[0].j = j; c[0].i = i; v[0] = -hxhydhz*(dd - gd);
534: c[1].k = k; c[1].j = j-1; c[1].i = i;
535: v[1] = -hzhxdhy*(ds - gs);
536: c[2].k = k; c[2].j = j; c[2].i = i-1;
537: v[2] = -hyhzdhx*(dw - gw);
538: c[3].k = k; c[3].j = j; c[3].i = i;
539: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
540: c[4].k = k; c[4].j = j; c[4].i = i+1;
541: v[4] = -hyhzdhx*(de + ge);
542: c[5].k = k; c[5].j = j+1; c[5].i = i;
543: v[5] = -hzhxdhy*(dn + gn);
544: c[6].k = k+1; c[6].j = j; c[6].i = i;
545: v[6] = -hxhydhz*(du + gu);
546: MatSetValuesStencil(jac,1,&row,7,c,v,INSERT_VALUES);
548: } else if (i == 0) {
550: /* left-hand plane boundary */
551: tw = tleft;
552: aw = 0.5*(t0 + tw);
553: bw = PetscPowScalar(aw,bm1);
554: /* dw = bw * aw */
555: dw = PetscPowScalar(aw,beta);
556: gw = coef*bw*(t0 - tw);
558: te = x[k][j][i+1];
559: ae = 0.5*(t0 + te);
560: be = PetscPowScalar(ae,bm1);
561: /* de = be * ae; */
562: de = PetscPowScalar(ae,beta);
563: ge = coef*be*(te - t0);
565: /* left-hand bottom edge */
566: if (j == 0) {
568: tn = x[k][j+1][i];
569: an = 0.5*(t0 + tn);
570: bn = PetscPowScalar(an,bm1);
571: /* dn = bn * an; */
572: dn = PetscPowScalar(an,beta);
573: gn = coef*bn*(tn - t0);
575: /* left-hand bottom down corner */
576: if (k == 0) {
578: tu = x[k+1][j][i];
579: au = 0.5*(t0 + tu);
580: bu = PetscPowScalar(au,bm1);
581: /* du = bu * au; */
582: du = PetscPowScalar(au,beta);
583: gu = coef*bu*(tu - t0);
585: c[0].k = k; c[0].j = j; c[0].i = i;
586: v[0] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
587: c[1].k = k; c[1].j = j; c[1].i = i+1;
588: v[1] = -hyhzdhx*(de + ge);
589: c[2].k = k; c[2].j = j+1; c[2].i = i;
590: v[2] = -hzhxdhy*(dn + gn);
591: c[3].k = k+1; c[3].j = j; c[3].i = i;
592: v[3] = -hxhydhz*(du + gu);
593: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
595: /* left-hand bottom interior edge */
596: } else if (k < mz-1) {
598: tu = x[k+1][j][i];
599: au = 0.5*(t0 + tu);
600: bu = PetscPowScalar(au,bm1);
601: /* du = bu * au; */
602: du = PetscPowScalar(au,beta);
603: gu = coef*bu*(tu - t0);
605: td = x[k-1][j][i];
606: ad = 0.5*(t0 + td);
607: bd = PetscPowScalar(ad,bm1);
608: /* dd = bd * ad; */
609: dd = PetscPowScalar(ad,beta);
610: gd = coef*bd*(td - t0);
612: c[0].k = k-1; c[0].j = j; c[0].i = i;
613: v[0] = -hxhydhz*(dd - gd);
614: c[1].k = k; c[1].j = j; c[1].i = i;
615: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
616: c[2].k = k; c[2].j = j; c[2].i = i+1;
617: v[2] = -hyhzdhx*(de + ge);
618: c[3].k = k; c[3].j = j+1; c[3].i = i;
619: v[3] = -hzhxdhy*(dn + gn);
620: c[4].k = k+1; c[4].j = j; c[4].i = i;
621: v[4] = -hxhydhz*(du + gu);
622: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
624: /* left-hand bottom up corner */
625: } else {
627: td = x[k-1][j][i];
628: ad = 0.5*(t0 + td);
629: bd = PetscPowScalar(ad,bm1);
630: /* dd = bd * ad; */
631: dd = PetscPowScalar(ad,beta);
632: gd = coef*bd*(td - t0);
634: c[0].k = k-1; c[0].j = j; c[0].i = i;
635: v[0] = -hxhydhz*(dd - gd);
636: c[1].k = k; c[1].j = j; c[1].i = i;
637: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
638: c[2].k = k; c[2].j = j; c[2].i = i+1;
639: v[2] = -hyhzdhx*(de + ge);
640: c[3].k = k; c[3].j = j+1; c[3].i = i;
641: v[3] = -hzhxdhy*(dn + gn);
642: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
643: }
645: /* left-hand top edge */
646: } else if (j == my-1) {
648: ts = x[k][j-1][i];
649: as = 0.5*(t0 + ts);
650: bs = PetscPowScalar(as,bm1);
651: /* ds = bs * as; */
652: ds = PetscPowScalar(as,beta);
653: gs = coef*bs*(ts - t0);
655: /* left-hand top down corner */
656: if (k == 0) {
658: tu = x[k+1][j][i];
659: au = 0.5*(t0 + tu);
660: bu = PetscPowScalar(au,bm1);
661: /* du = bu * au; */
662: du = PetscPowScalar(au,beta);
663: gu = coef*bu*(tu - t0);
665: c[0].k = k; c[0].j = j-1; c[0].i = i;
666: v[0] = -hzhxdhy*(ds - gs);
667: c[1].k = k; c[1].j = j; c[1].i = i;
668: v[1] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
669: c[2].k = k; c[2].j = j; c[2].i = i+1;
670: v[2] = -hyhzdhx*(de + ge);
671: c[3].k = k+1; c[3].j = j; c[3].i = i;
672: v[3] = -hxhydhz*(du + gu);
673: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
675: /* left-hand top interior edge */
676: } else if (k < mz-1) {
678: tu = x[k+1][j][i];
679: au = 0.5*(t0 + tu);
680: bu = PetscPowScalar(au,bm1);
681: /* du = bu * au; */
682: du = PetscPowScalar(au,beta);
683: gu = coef*bu*(tu - t0);
685: td = x[k-1][j][i];
686: ad = 0.5*(t0 + td);
687: bd = PetscPowScalar(ad,bm1);
688: /* dd = bd * ad; */
689: dd = PetscPowScalar(ad,beta);
690: gd = coef*bd*(td - t0);
692: c[0].k = k-1; c[0].j = j; c[0].i = i;
693: v[0] = -hxhydhz*(dd - gd);
694: c[1].k = k; c[1].j = j-1; c[1].i = i;
695: v[1] = -hzhxdhy*(ds - gs);
696: c[2].k = k; c[2].j = j; c[2].i = i;
697: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
698: c[3].k = k; c[3].j = j; c[3].i = i+1;
699: v[3] = -hyhzdhx*(de + ge);
700: c[4].k = k+1; c[4].j = j; c[4].i = i;
701: v[4] = -hxhydhz*(du + gu);
702: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
704: /* left-hand top up corner */
705: } else {
707: td = x[k-1][j][i];
708: ad = 0.5*(t0 + td);
709: bd = PetscPowScalar(ad,bm1);
710: /* dd = bd * ad; */
711: dd = PetscPowScalar(ad,beta);
712: gd = coef*bd*(td - t0);
714: c[0].k = k-1; c[0].j = j; c[0].i = i;
715: v[0] = -hxhydhz*(dd - gd);
716: c[1].k = k; c[1].j = j-1; c[1].i = i;
717: v[1] = -hzhxdhy*(ds - gs);
718: c[2].k = k; c[2].j = j; c[2].i = i;
719: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
720: c[3].k = k; c[3].j = j; c[3].i = i+1;
721: v[3] = -hyhzdhx*(de + ge);
722: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
723: }
725: } else {
727: ts = x[k][j-1][i];
728: as = 0.5*(t0 + ts);
729: bs = PetscPowScalar(as,bm1);
730: /* ds = bs * as; */
731: ds = PetscPowScalar(as,beta);
732: gs = coef*bs*(t0 - ts);
734: tn = x[k][j+1][i];
735: an = 0.5*(t0 + tn);
736: bn = PetscPowScalar(an,bm1);
737: /* dn = bn * an; */
738: dn = PetscPowScalar(an,beta);
739: gn = coef*bn*(tn - t0);
741: /* left-hand down interior edge */
742: if (k == 0) {
744: tu = x[k+1][j][i];
745: au = 0.5*(t0 + tu);
746: bu = PetscPowScalar(au,bm1);
747: /* du = bu * au; */
748: du = PetscPowScalar(au,beta);
749: gu = coef*bu*(tu - t0);
751: c[0].k = k; c[0].j = j-1; c[0].i = i;
752: v[0] = -hzhxdhy*(ds - gs);
753: c[1].k = k; c[1].j = j; c[1].i = i;
754: v[1] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
755: c[2].k = k; c[2].j = j; c[2].i = i+1;
756: v[2] = -hyhzdhx*(de + ge);
757: c[3].k = k; c[3].j = j+1; c[3].i = i;
758: v[3] = -hzhxdhy*(dn + gn);
759: c[4].k = k+1; c[4].j = j; c[4].i = i;
760: v[4] = -hxhydhz*(du + gu);
761: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
763: } else if (k == mz-1) { /* left-hand up interior edge */
765: td = x[k-1][j][i];
766: ad = 0.5*(t0 + td);
767: bd = PetscPowScalar(ad,bm1);
768: /* dd = bd * ad; */
769: dd = PetscPowScalar(ad,beta);
770: gd = coef*bd*(t0 - td);
772: c[0].k = k-1; c[0].j = j; c[0].i = i;
773: v[0] = -hxhydhz*(dd - gd);
774: c[1].k = k; c[1].j = j-1; c[1].i = i;
775: v[1] = -hzhxdhy*(ds - gs);
776: c[2].k = k; c[2].j = j; c[2].i = i;
777: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
778: c[3].k = k; c[3].j = j; c[3].i = i+1;
779: v[3] = -hyhzdhx*(de + ge);
780: c[4].k = k; c[4].j = j+1; c[4].i = i;
781: v[4] = -hzhxdhy*(dn + gn);
782: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
783: } else { /* left-hand interior plane */
785: td = x[k-1][j][i];
786: ad = 0.5*(t0 + td);
787: bd = PetscPowScalar(ad,bm1);
788: /* dd = bd * ad; */
789: dd = PetscPowScalar(ad,beta);
790: gd = coef*bd*(t0 - td);
792: tu = x[k+1][j][i];
793: au = 0.5*(t0 + tu);
794: bu = PetscPowScalar(au,bm1);
795: /* du = bu * au; */
796: du = PetscPowScalar(au,beta);
797: gu = coef*bu*(tu - t0);
799: c[0].k = k-1; c[0].j = j; c[0].i = i;
800: v[0] = -hxhydhz*(dd - gd);
801: c[1].k = k; c[1].j = j-1; c[1].i = i;
802: v[1] = -hzhxdhy*(ds - gs);
803: c[2].k = k; c[2].j = j; c[2].i = i;
804: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
805: c[3].k = k; c[3].j = j; c[3].i = i+1;
806: v[3] = -hyhzdhx*(de + ge);
807: c[4].k = k; c[4].j = j+1; c[4].i = i;
808: v[4] = -hzhxdhy*(dn + gn);
809: c[5].k = k+1; c[5].j = j; c[5].i = i;
810: v[5] = -hxhydhz*(du + gu);
811: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
812: }
813: }
815: } else if (i == mx-1) {
817: /* right-hand plane boundary */
818: tw = x[k][j][i-1];
819: aw = 0.5*(t0 + tw);
820: bw = PetscPowScalar(aw,bm1);
821: /* dw = bw * aw */
822: dw = PetscPowScalar(aw,beta);
823: gw = coef*bw*(t0 - tw);
825: te = tright;
826: ae = 0.5*(t0 + te);
827: be = PetscPowScalar(ae,bm1);
828: /* de = be * ae; */
829: de = PetscPowScalar(ae,beta);
830: ge = coef*be*(te - t0);
832: /* right-hand bottom edge */
833: if (j == 0) {
835: tn = x[k][j+1][i];
836: an = 0.5*(t0 + tn);
837: bn = PetscPowScalar(an,bm1);
838: /* dn = bn * an; */
839: dn = PetscPowScalar(an,beta);
840: gn = coef*bn*(tn - t0);
842: /* right-hand bottom down corner */
843: if (k == 0) {
845: tu = x[k+1][j][i];
846: au = 0.5*(t0 + tu);
847: bu = PetscPowScalar(au,bm1);
848: /* du = bu * au; */
849: du = PetscPowScalar(au,beta);
850: gu = coef*bu*(tu - t0);
852: c[0].k = k; c[0].j = j; c[0].i = i-1;
853: v[0] = -hyhzdhx*(dw - gw);
854: c[1].k = k; c[1].j = j; c[1].i = i;
855: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
856: c[2].k = k; c[2].j = j+1; c[2].i = i;
857: v[2] = -hzhxdhy*(dn + gn);
858: c[3].k = k+1; c[3].j = j; c[3].i = i;
859: v[3] = -hxhydhz*(du + gu);
860: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
862: /* right-hand bottom interior edge */
863: } else if (k < mz-1) {
865: tu = x[k+1][j][i];
866: au = 0.5*(t0 + tu);
867: bu = PetscPowScalar(au,bm1);
868: /* du = bu * au; */
869: du = PetscPowScalar(au,beta);
870: gu = coef*bu*(tu - t0);
872: td = x[k-1][j][i];
873: ad = 0.5*(t0 + td);
874: bd = PetscPowScalar(ad,bm1);
875: /* dd = bd * ad; */
876: dd = PetscPowScalar(ad,beta);
877: gd = coef*bd*(td - t0);
879: c[0].k = k-1; c[0].j = j; c[0].i = i;
880: v[0] = -hxhydhz*(dd - gd);
881: c[1].k = k; c[1].j = j; c[1].i = i-1;
882: v[1] = -hyhzdhx*(dw - gw);
883: c[2].k = k; c[2].j = j; c[2].i = i;
884: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
885: c[3].k = k; c[3].j = j+1; c[3].i = i;
886: v[3] = -hzhxdhy*(dn + gn);
887: c[4].k = k+1; c[4].j = j; c[4].i = i;
888: v[4] = -hxhydhz*(du + gu);
889: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
891: /* right-hand bottom up corner */
892: } else {
894: td = x[k-1][j][i];
895: ad = 0.5*(t0 + td);
896: bd = PetscPowScalar(ad,bm1);
897: /* dd = bd * ad; */
898: dd = PetscPowScalar(ad,beta);
899: gd = coef*bd*(td - t0);
901: c[0].k = k-1; c[0].j = j; c[0].i = i;
902: v[0] = -hxhydhz*(dd - gd);
903: c[1].k = k; c[1].j = j; c[1].i = i-1;
904: v[1] = -hyhzdhx*(dw - gw);
905: c[2].k = k; c[2].j = j; c[2].i = i;
906: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
907: c[3].k = k; c[3].j = j+1; c[3].i = i;
908: v[3] = -hzhxdhy*(dn + gn);
909: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
910: }
912: /* right-hand top edge */
913: } else if (j == my-1) {
915: ts = x[k][j-1][i];
916: as = 0.5*(t0 + ts);
917: bs = PetscPowScalar(as,bm1);
918: /* ds = bs * as; */
919: ds = PetscPowScalar(as,beta);
920: gs = coef*bs*(ts - t0);
922: /* right-hand top down corner */
923: if (k == 0) {
925: tu = x[k+1][j][i];
926: au = 0.5*(t0 + tu);
927: bu = PetscPowScalar(au,bm1);
928: /* du = bu * au; */
929: du = PetscPowScalar(au,beta);
930: gu = coef*bu*(tu - t0);
932: c[0].k = k; c[0].j = j-1; c[0].i = i;
933: v[0] = -hzhxdhy*(ds - gs);
934: c[1].k = k; c[1].j = j; c[1].i = i-1;
935: v[1] = -hyhzdhx*(dw - gw);
936: c[2].k = k; c[2].j = j; c[2].i = i;
937: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
938: c[3].k = k+1; c[3].j = j; c[3].i = i;
939: v[3] = -hxhydhz*(du + gu);
940: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
942: /* right-hand top interior edge */
943: } else if (k < mz-1) {
945: tu = x[k+1][j][i];
946: au = 0.5*(t0 + tu);
947: bu = PetscPowScalar(au,bm1);
948: /* du = bu * au; */
949: du = PetscPowScalar(au,beta);
950: gu = coef*bu*(tu - t0);
952: td = x[k-1][j][i];
953: ad = 0.5*(t0 + td);
954: bd = PetscPowScalar(ad,bm1);
955: /* dd = bd * ad; */
956: dd = PetscPowScalar(ad,beta);
957: gd = coef*bd*(td - t0);
959: c[0].k = k-1; c[0].j = j; c[0].i = i;
960: v[0] = -hxhydhz*(dd - gd);
961: c[1].k = k; c[1].j = j-1; c[1].i = i;
962: v[1] = -hzhxdhy*(ds - gs);
963: c[2].k = k; c[2].j = j; c[2].i = i-1;
964: v[2] = -hyhzdhx*(dw - gw);
965: c[3].k = k; c[3].j = j; c[3].i = i;
966: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
967: c[4].k = k+1; c[4].j = j; c[4].i = i;
968: v[4] = -hxhydhz*(du + gu);
969: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
971: /* right-hand top up corner */
972: } else {
974: td = x[k-1][j][i];
975: ad = 0.5*(t0 + td);
976: bd = PetscPowScalar(ad,bm1);
977: /* dd = bd * ad; */
978: dd = PetscPowScalar(ad,beta);
979: gd = coef*bd*(td - t0);
981: c[0].k = k-1; c[0].j = j; c[0].i = i;
982: v[0] = -hxhydhz*(dd - gd);
983: c[1].k = k; c[1].j = j-1; c[1].i = i;
984: v[1] = -hzhxdhy*(ds - gs);
985: c[2].k = k; c[2].j = j; c[2].i = i-1;
986: v[2] = -hyhzdhx*(dw - gw);
987: c[3].k = k; c[3].j = j; c[3].i = i;
988: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
989: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
990: }
992: } else {
994: ts = x[k][j-1][i];
995: as = 0.5*(t0 + ts);
996: bs = PetscPowScalar(as,bm1);
997: /* ds = bs * as; */
998: ds = PetscPowScalar(as,beta);
999: gs = coef*bs*(t0 - ts);
1001: tn = x[k][j+1][i];
1002: an = 0.5*(t0 + tn);
1003: bn = PetscPowScalar(an,bm1);
1004: /* dn = bn * an; */
1005: dn = PetscPowScalar(an,beta);
1006: gn = coef*bn*(tn - t0);
1008: /* right-hand down interior edge */
1009: if (k == 0) {
1011: tu = x[k+1][j][i];
1012: au = 0.5*(t0 + tu);
1013: bu = PetscPowScalar(au,bm1);
1014: /* du = bu * au; */
1015: du = PetscPowScalar(au,beta);
1016: gu = coef*bu*(tu - t0);
1018: c[0].k = k; c[0].j = j-1; c[0].i = i;
1019: v[0] = -hzhxdhy*(ds - gs);
1020: c[1].k = k; c[1].j = j; c[1].i = i-1;
1021: v[1] = -hyhzdhx*(dw - gw);
1022: c[2].k = k; c[2].j = j; c[2].i = i;
1023: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1024: c[3].k = k; c[3].j = j+1; c[3].i = i;
1025: v[3] = -hzhxdhy*(dn + gn);
1026: c[4].k = k+1; c[4].j = j; c[4].i = i;
1027: v[4] = -hxhydhz*(du + gu);
1028: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1030: } else if (k == mz-1) { /* right-hand up interior edge */
1032: td = x[k-1][j][i];
1033: ad = 0.5*(t0 + td);
1034: bd = PetscPowScalar(ad,bm1);
1035: /* dd = bd * ad; */
1036: dd = PetscPowScalar(ad,beta);
1037: gd = coef*bd*(t0 - td);
1039: c[0].k = k-1; c[0].j = j; c[0].i = i;
1040: v[0] = -hxhydhz*(dd - gd);
1041: c[1].k = k; c[1].j = j-1; c[1].i = i;
1042: v[1] = -hzhxdhy*(ds - gs);
1043: c[2].k = k; c[2].j = j; c[2].i = i-1;
1044: v[2] = -hyhzdhx*(dw - gw);
1045: c[3].k = k; c[3].j = j; c[3].i = i;
1046: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1047: c[4].k = k; c[4].j = j+1; c[4].i = i;
1048: v[4] = -hzhxdhy*(dn + gn);
1049: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1051: } else { /* right-hand interior plane */
1053: td = x[k-1][j][i];
1054: ad = 0.5*(t0 + td);
1055: bd = PetscPowScalar(ad,bm1);
1056: /* dd = bd * ad; */
1057: dd = PetscPowScalar(ad,beta);
1058: gd = coef*bd*(t0 - td);
1060: tu = x[k+1][j][i];
1061: au = 0.5*(t0 + tu);
1062: bu = PetscPowScalar(au,bm1);
1063: /* du = bu * au; */
1064: du = PetscPowScalar(au,beta);
1065: gu = coef*bu*(tu - t0);
1067: c[0].k = k-1; c[0].j = j; c[0].i = i;
1068: v[0] = -hxhydhz*(dd - gd);
1069: c[1].k = k; c[1].j = j-1; c[1].i = i;
1070: v[1] = -hzhxdhy*(ds - gs);
1071: c[2].k = k; c[2].j = j; c[2].i = i-1;
1072: v[2] = -hyhzdhx*(dw - gw);
1073: c[3].k = k; c[3].j = j; c[3].i = i;
1074: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1075: c[4].k = k; c[4].j = j+1; c[4].i = i;
1076: v[4] = -hzhxdhy*(dn + gn);
1077: c[5].k = k+1; c[5].j = j; c[5].i = i;
1078: v[5] = -hxhydhz*(du + gu);
1079: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1080: }
1081: }
1083: } else if (j == 0) {
1085: tw = x[k][j][i-1];
1086: aw = 0.5*(t0 + tw);
1087: bw = PetscPowScalar(aw,bm1);
1088: /* dw = bw * aw */
1089: dw = PetscPowScalar(aw,beta);
1090: gw = coef*bw*(t0 - tw);
1092: te = x[k][j][i+1];
1093: ae = 0.5*(t0 + te);
1094: be = PetscPowScalar(ae,bm1);
1095: /* de = be * ae; */
1096: de = PetscPowScalar(ae,beta);
1097: ge = coef*be*(te - t0);
1099: tn = x[k][j+1][i];
1100: an = 0.5*(t0 + tn);
1101: bn = PetscPowScalar(an,bm1);
1102: /* dn = bn * an; */
1103: dn = PetscPowScalar(an,beta);
1104: gn = coef*bn*(tn - t0);
1107: /* bottom down interior edge */
1108: if (k == 0) {
1110: tu = x[k+1][j][i];
1111: au = 0.5*(t0 + tu);
1112: bu = PetscPowScalar(au,bm1);
1113: /* du = bu * au; */
1114: du = PetscPowScalar(au,beta);
1115: gu = coef*bu*(tu - t0);
1117: c[0].k = k; c[0].j = j; c[0].i = i-1;
1118: v[0] = -hyhzdhx*(dw - gw);
1119: c[1].k = k; c[1].j = j; c[1].i = i;
1120: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1121: c[2].k = k; c[2].j = j; c[2].i = i+1;
1122: v[2] = -hyhzdhx*(de + ge);
1123: c[3].k = k; c[3].j = j+1; c[3].i = i;
1124: v[3] = -hzhxdhy*(dn + gn);
1125: c[4].k = k+1; c[4].j = j; c[4].i = i;
1126: v[4] = -hxhydhz*(du + gu);
1127: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1129: } else if (k == mz-1) { /* bottom up interior edge */
1131: td = x[k-1][j][i];
1132: ad = 0.5*(t0 + td);
1133: bd = PetscPowScalar(ad,bm1);
1134: /* dd = bd * ad; */
1135: dd = PetscPowScalar(ad,beta);
1136: gd = coef*bd*(td - t0);
1138: c[0].k = k-1; c[0].j = j; c[0].i = i;
1139: v[0] = -hxhydhz*(dd - gd);
1140: c[1].k = k; c[1].j = j; c[1].i = i-1;
1141: v[1] = -hyhzdhx*(dw - gw);
1142: c[2].k = k; c[2].j = j; c[2].i = i;
1143: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1144: c[3].k = k; c[3].j = j; c[3].i = i+1;
1145: v[3] = -hyhzdhx*(de + ge);
1146: c[4].k = k; c[4].j = j+1; c[4].i = i;
1147: v[4] = -hzhxdhy*(dn + gn);
1148: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1150: } else { /* bottom interior plane */
1152: tu = x[k+1][j][i];
1153: au = 0.5*(t0 + tu);
1154: bu = PetscPowScalar(au,bm1);
1155: /* du = bu * au; */
1156: du = PetscPowScalar(au,beta);
1157: gu = coef*bu*(tu - t0);
1159: td = x[k-1][j][i];
1160: ad = 0.5*(t0 + td);
1161: bd = PetscPowScalar(ad,bm1);
1162: /* dd = bd * ad; */
1163: dd = PetscPowScalar(ad,beta);
1164: gd = coef*bd*(td - t0);
1166: c[0].k = k-1; c[0].j = j; c[0].i = i;
1167: v[0] = -hxhydhz*(dd - gd);
1168: c[1].k = k; c[1].j = j; c[1].i = i-1;
1169: v[1] = -hyhzdhx*(dw - gw);
1170: c[2].k = k; c[2].j = j; c[2].i = i;
1171: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1172: c[3].k = k; c[3].j = j; c[3].i = i+1;
1173: v[3] = -hyhzdhx*(de + ge);
1174: c[4].k = k; c[4].j = j+1; c[4].i = i;
1175: v[4] = -hzhxdhy*(dn + gn);
1176: c[5].k = k+1; c[5].j = j; c[5].i = i;
1177: v[5] = -hxhydhz*(du + gu);
1178: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1179: }
1181: } else if (j == my-1) {
1183: tw = x[k][j][i-1];
1184: aw = 0.5*(t0 + tw);
1185: bw = PetscPowScalar(aw,bm1);
1186: /* dw = bw * aw */
1187: dw = PetscPowScalar(aw,beta);
1188: gw = coef*bw*(t0 - tw);
1190: te = x[k][j][i+1];
1191: ae = 0.5*(t0 + te);
1192: be = PetscPowScalar(ae,bm1);
1193: /* de = be * ae; */
1194: de = PetscPowScalar(ae,beta);
1195: ge = coef*be*(te - t0);
1197: ts = x[k][j-1][i];
1198: as = 0.5*(t0 + ts);
1199: bs = PetscPowScalar(as,bm1);
1200: /* ds = bs * as; */
1201: ds = PetscPowScalar(as,beta);
1202: gs = coef*bs*(t0 - ts);
1204: /* top down interior edge */
1205: if (k == 0) {
1207: tu = x[k+1][j][i];
1208: au = 0.5*(t0 + tu);
1209: bu = PetscPowScalar(au,bm1);
1210: /* du = bu * au; */
1211: du = PetscPowScalar(au,beta);
1212: gu = coef*bu*(tu - t0);
1214: c[0].k = k; c[0].j = j-1; c[0].i = i;
1215: v[0] = -hzhxdhy*(ds - gs);
1216: c[1].k = k; c[1].j = j; c[1].i = i-1;
1217: v[1] = -hyhzdhx*(dw - gw);
1218: c[2].k = k; c[2].j = j; c[2].i = i;
1219: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1220: c[3].k = k; c[3].j = j; c[3].i = i+1;
1221: v[3] = -hyhzdhx*(de + ge);
1222: c[4].k = k+1; c[4].j = j; c[4].i = i;
1223: v[4] = -hxhydhz*(du + gu);
1224: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1226: } else if (k == mz-1) { /* top up interior edge */
1228: td = x[k-1][j][i];
1229: ad = 0.5*(t0 + td);
1230: bd = PetscPowScalar(ad,bm1);
1231: /* dd = bd * ad; */
1232: dd = PetscPowScalar(ad,beta);
1233: gd = coef*bd*(td - t0);
1235: c[0].k = k-1; c[0].j = j; c[0].i = i;
1236: v[0] = -hxhydhz*(dd - gd);
1237: c[1].k = k; c[1].j = j-1; c[1].i = i;
1238: v[1] = -hzhxdhy*(ds - gs);
1239: c[2].k = k; c[2].j = j; c[2].i = i-1;
1240: v[2] = -hyhzdhx*(dw - gw);
1241: c[3].k = k; c[3].j = j; c[3].i = i;
1242: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1243: c[4].k = k; c[4].j = j; c[4].i = i+1;
1244: v[4] = -hyhzdhx*(de + ge);
1245: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1247: } else { /* top interior plane */
1249: tu = x[k+1][j][i];
1250: au = 0.5*(t0 + tu);
1251: bu = PetscPowScalar(au,bm1);
1252: /* du = bu * au; */
1253: du = PetscPowScalar(au,beta);
1254: gu = coef*bu*(tu - t0);
1256: td = x[k-1][j][i];
1257: ad = 0.5*(t0 + td);
1258: bd = PetscPowScalar(ad,bm1);
1259: /* dd = bd * ad; */
1260: dd = PetscPowScalar(ad,beta);
1261: gd = coef*bd*(td - t0);
1263: c[0].k = k-1; c[0].j = j; c[0].i = i;
1264: v[0] = -hxhydhz*(dd - gd);
1265: c[1].k = k; c[1].j = j-1; c[1].i = i;
1266: v[1] = -hzhxdhy*(ds - gs);
1267: c[2].k = k; c[2].j = j; c[2].i = i-1;
1268: v[2] = -hyhzdhx*(dw - gw);
1269: c[3].k = k; c[3].j = j; c[3].i = i;
1270: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1271: c[4].k = k; c[4].j = j; c[4].i = i+1;
1272: v[4] = -hyhzdhx*(de + ge);
1273: c[5].k = k+1; c[5].j = j; c[5].i = i;
1274: v[5] = -hxhydhz*(du + gu);
1275: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1276: }
1278: } else if (k == 0) {
1280: /* down interior plane */
1282: tw = x[k][j][i-1];
1283: aw = 0.5*(t0 + tw);
1284: bw = PetscPowScalar(aw,bm1);
1285: /* dw = bw * aw */
1286: dw = PetscPowScalar(aw,beta);
1287: gw = coef*bw*(t0 - tw);
1289: te = x[k][j][i+1];
1290: ae = 0.5*(t0 + te);
1291: be = PetscPowScalar(ae,bm1);
1292: /* de = be * ae; */
1293: de = PetscPowScalar(ae,beta);
1294: ge = coef*be*(te - t0);
1296: ts = x[k][j-1][i];
1297: as = 0.5*(t0 + ts);
1298: bs = PetscPowScalar(as,bm1);
1299: /* ds = bs * as; */
1300: ds = PetscPowScalar(as,beta);
1301: gs = coef*bs*(t0 - ts);
1303: tn = x[k][j+1][i];
1304: an = 0.5*(t0 + tn);
1305: bn = PetscPowScalar(an,bm1);
1306: /* dn = bn * an; */
1307: dn = PetscPowScalar(an,beta);
1308: gn = coef*bn*(tn - t0);
1310: tu = x[k+1][j][i];
1311: au = 0.5*(t0 + tu);
1312: bu = PetscPowScalar(au,bm1);
1313: /* du = bu * au; */
1314: du = PetscPowScalar(au,beta);
1315: gu = coef*bu*(tu - t0);
1317: c[0].k = k; c[0].j = j-1; c[0].i = i;
1318: v[0] = -hzhxdhy*(ds - gs);
1319: c[1].k = k; c[1].j = j; c[1].i = i-1;
1320: v[1] = -hyhzdhx*(dw - gw);
1321: c[2].k = k; c[2].j = j; c[2].i = i;
1322: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1323: c[3].k = k; c[3].j = j; c[3].i = i+1;
1324: v[3] = -hyhzdhx*(de + ge);
1325: c[4].k = k; c[4].j = j+1; c[4].i = i;
1326: v[4] = -hzhxdhy*(dn + gn);
1327: c[5].k = k+1; c[5].j = j; c[5].i = i;
1328: v[5] = -hxhydhz*(du + gu);
1329: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1331: } else if (k == mz-1) {
1333: /* up interior plane */
1335: tw = x[k][j][i-1];
1336: aw = 0.5*(t0 + tw);
1337: bw = PetscPowScalar(aw,bm1);
1338: /* dw = bw * aw */
1339: dw = PetscPowScalar(aw,beta);
1340: gw = coef*bw*(t0 - tw);
1342: te = x[k][j][i+1];
1343: ae = 0.5*(t0 + te);
1344: be = PetscPowScalar(ae,bm1);
1345: /* de = be * ae; */
1346: de = PetscPowScalar(ae,beta);
1347: ge = coef*be*(te - t0);
1349: ts = x[k][j-1][i];
1350: as = 0.5*(t0 + ts);
1351: bs = PetscPowScalar(as,bm1);
1352: /* ds = bs * as; */
1353: ds = PetscPowScalar(as,beta);
1354: gs = coef*bs*(t0 - ts);
1356: tn = x[k][j+1][i];
1357: an = 0.5*(t0 + tn);
1358: bn = PetscPowScalar(an,bm1);
1359: /* dn = bn * an; */
1360: dn = PetscPowScalar(an,beta);
1361: gn = coef*bn*(tn - t0);
1363: td = x[k-1][j][i];
1364: ad = 0.5*(t0 + td);
1365: bd = PetscPowScalar(ad,bm1);
1366: /* dd = bd * ad; */
1367: dd = PetscPowScalar(ad,beta);
1368: gd = coef*bd*(t0 - td);
1370: c[0].k = k-1; c[0].j = j; c[0].i = i;
1371: v[0] = -hxhydhz*(dd - gd);
1372: c[1].k = k; c[1].j = j-1; c[1].i = i;
1373: v[1] = -hzhxdhy*(ds - gs);
1374: c[2].k = k; c[2].j = j; c[2].i = i-1;
1375: v[2] = -hyhzdhx*(dw - gw);
1376: c[3].k = k; c[3].j = j; c[3].i = i;
1377: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1378: c[4].k = k; c[4].j = j; c[4].i = i+1;
1379: v[4] = -hyhzdhx*(de + ge);
1380: c[5].k = k; c[5].j = j+1; c[5].i = i;
1381: v[5] = -hzhxdhy*(dn + gn);
1382: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1383: }
1384: }
1385: }
1386: }
1387: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
1388: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
1389: if (jac != J) {
1390: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1391: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1392: }
1393: DMDAVecRestoreArray(da,localX,&x);
1394: DMRestoreLocalVector(da,&localX);
1396: PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
1397: return(0);
1398: }