Actual source code: ex20.c
petsc-3.3-p7 2013-05-11
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: /*
22:
23: This example models the partial differential equation
24:
25: - Div(alpha* T^beta (GRAD T)) = 0.
26:
27: where beta = 2.5 and alpha = 1.0
28:
29: BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = dT/dn_up = dT/dn_down = 0.
30:
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.
33:
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
39:
40: */
42: #include <petscsnes.h>
43: #include <petscdmda.h>
44: #include <petscpcmg.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(DM,Vec);
56: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
57: extern PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,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,PETSC_NULL,help);
72: /* set problem parameters */
73: user.tleft = 1.0;
74: user.tright = 0.1;
75: user.beta = 2.5;
76: PetscOptionsGetReal(PETSC_NULL,"-tleft",&user.tleft,PETSC_NULL);
77: PetscOptionsGetReal(PETSC_NULL,"-tright",&user.tright,PETSC_NULL);
78: PetscOptionsGetReal(PETSC_NULL,"-beta",&user.beta,PETSC_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,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);
86: DMSetApplicationContext(da,&user);
88: /*
89: Create the nonlinear solver
90: */
91: SNESCreate(PETSC_COMM_WORLD,&snes);
92: SNESSetDM(snes,da);
93: SNESSetFunction(snes,PETSC_NULL,FormFunction,&user);
94: SNESSetJacobian(snes,PETSC_NULL,PETSC_NULL,FormJacobian,&user);
95: SNESSetFromOptions(snes);
96: DMSetInitialGuess(da,FormInitialGuess);
98: SNESSolve(snes,PETSC_NULL,PETSC_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",litspit);
106: SNESDestroy(&snes);
107: DMDestroy(&da);
108: PetscFinalize();
110: return 0;
111: }
112: /* -------------------- Form initial approximation ----------------- */
115: PetscErrorCode FormInitialGuess(DM da,Vec X)
116: {
117: AppCtx *user;
118: PetscInt i,j,k,xs,ys,xm,ym,zs,zm;
120: PetscScalar ***x;
123: DMGetApplicationContext(da,&user);
124: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
125: DMDAVecGetArray(da,X,&x);
127: /* Compute initial guess */
128: for (k=zs; k<zs+zm; k++) {
129: for (j=ys; j<ys+ym; j++) {
130: for (i=xs; i<xs+xm; i++) {
131: x[k][j][i] = user->tleft;
132: }
133: }
134: }
135: DMDAVecRestoreArray(da,X,&x);
136: return(0);
137: }
138: /* -------------------- Evaluate Function F(x) --------------------- */
141: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ptr)
142: {
143: AppCtx *user = (AppCtx*)ptr;
145: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
146: PetscScalar zero = 0.0,one = 1.0;
147: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
148: 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;
149: PetscScalar tleft,tright,beta,td,ad,dd,fd=0.0,tu,au,du=0.0,fu=0.0;
150: PetscScalar ***x,***f;
151: Vec localX;
152: DM da;
155: SNESGetDM(snes,&da);
156: DMGetLocalVector(da,&localX);
157: DMDAGetInfo(da,PETSC_NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
158: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
159: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
160: tleft = user->tleft; tright = user->tright;
161: beta = user->beta;
162:
163: /* Get ghost points */
164: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
165: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
166: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
167: DMDAVecGetArray(da,localX,&x);
168: DMDAVecGetArray(da,F,&f);
170: /* Evaluate function */
171: for (k=zs; k<zs+zm; k++) {
172: for (j=ys; j<ys+ym; j++) {
173: for (i=xs; i<xs+xm; i++) {
174: t0 = x[k][j][i];
176: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
178: /* general interior volume */
180: tw = x[k][j][i-1];
181: aw = 0.5*(t0 + tw);
182: dw = pow(aw,beta);
183: fw = dw*(t0 - tw);
185: te = x[k][j][i+1];
186: ae = 0.5*(t0 + te);
187: de = pow(ae,beta);
188: fe = de*(te - t0);
190: ts = x[k][j-1][i];
191: as = 0.5*(t0 + ts);
192: ds = pow(as,beta);
193: fs = ds*(t0 - ts);
194:
195: tn = x[k][j+1][i];
196: an = 0.5*(t0 + tn);
197: dn = pow(an,beta);
198: fn = dn*(tn - t0);
200: td = x[k-1][j][i];
201: ad = 0.5*(t0 + td);
202: dd = pow(ad,beta);
203: fd = dd*(t0 - td);
205: tu = x[k+1][j][i];
206: au = 0.5*(t0 + tu);
207: du = pow(au,beta);
208: fu = du*(tu - t0);
210: } else if (i == 0) {
212: /* left-hand (west) boundary */
213: tw = tleft;
214: aw = 0.5*(t0 + tw);
215: dw = pow(aw,beta);
216: fw = dw*(t0 - tw);
218: te = x[k][j][i+1];
219: ae = 0.5*(t0 + te);
220: de = pow(ae,beta);
221: fe = de*(te - t0);
223: if (j > 0) {
224: ts = x[k][j-1][i];
225: as = 0.5*(t0 + ts);
226: ds = pow(as,beta);
227: fs = ds*(t0 - ts);
228: } else {
229: fs = zero;
230: }
232: if (j < my-1) {
233: tn = x[k][j+1][i];
234: an = 0.5*(t0 + tn);
235: dn = pow(an,beta);
236: fn = dn*(tn - t0);
237: } else {
238: fn = zero;
239: }
241: if (k > 0) {
242: td = x[k-1][j][i];
243: ad = 0.5*(t0 + td);
244: dd = pow(ad,beta);
245: fd = dd*(t0 - td);
246: } else {
247: fd = zero;
248: }
250: if (k < mz-1) {
251: tu = x[k+1][j][i];
252: au = 0.5*(t0 + tu);
253: du = pow(au,beta);
254: fu = du*(tu - t0);
255: } else {
256: fu = zero;
257: }
259: } else if (i == mx-1) {
261: /* right-hand (east) boundary */
262: tw = x[k][j][i-1];
263: aw = 0.5*(t0 + tw);
264: dw = pow(aw,beta);
265: fw = dw*(t0 - tw);
266:
267: te = tright;
268: ae = 0.5*(t0 + te);
269: de = pow(ae,beta);
270: fe = de*(te - t0);
271:
272: if (j > 0) {
273: ts = x[k][j-1][i];
274: as = 0.5*(t0 + ts);
275: ds = pow(as,beta);
276: fs = ds*(t0 - ts);
277: } else {
278: fs = zero;
279: }
280:
281: if (j < my-1) {
282: tn = x[k][j+1][i];
283: an = 0.5*(t0 + tn);
284: dn = pow(an,beta);
285: fn = dn*(tn - t0);
286: } else {
287: fn = zero;
288: }
290: if (k > 0) {
291: td = x[k-1][j][i];
292: ad = 0.5*(t0 + td);
293: dd = pow(ad,beta);
294: fd = dd*(t0 - td);
295: } else {
296: fd = zero;
297: }
299: if (k < mz-1) {
300: tu = x[k+1][j][i];
301: au = 0.5*(t0 + tu);
302: du = pow(au,beta);
303: fu = du*(tu - t0);
304: } else {
305: fu = zero;
306: }
308: } else if (j == 0) {
310: /* bottom (south) boundary, and i <> 0 or mx-1 */
311: tw = x[k][j][i-1];
312: aw = 0.5*(t0 + tw);
313: dw = pow(aw,beta);
314: fw = dw*(t0 - tw);
316: te = x[k][j][i+1];
317: ae = 0.5*(t0 + te);
318: de = pow(ae,beta);
319: fe = de*(te - t0);
321: fs = zero;
323: tn = x[k][j+1][i];
324: an = 0.5*(t0 + tn);
325: dn = pow(an,beta);
326: fn = dn*(tn - t0);
328: if (k > 0) {
329: td = x[k-1][j][i];
330: ad = 0.5*(t0 + td);
331: dd = pow(ad,beta);
332: fd = dd*(t0 - td);
333: } else {
334: fd = zero;
335: }
337: if (k < mz-1) {
338: tu = x[k+1][j][i];
339: au = 0.5*(t0 + tu);
340: du = pow(au,beta);
341: fu = du*(tu - t0);
342: } else {
343: fu = zero;
344: }
346: } else if (j == my-1) {
348: /* top (north) boundary, and i <> 0 or mx-1 */
349: tw = x[k][j][i-1];
350: aw = 0.5*(t0 + tw);
351: dw = pow(aw,beta);
352: fw = dw*(t0 - tw);
354: te = x[k][j][i+1];
355: ae = 0.5*(t0 + te);
356: de = pow(ae,beta);
357: fe = de*(te - t0);
359: ts = x[k][j-1][i];
360: as = 0.5*(t0 + ts);
361: ds = pow(as,beta);
362: fs = ds*(t0 - ts);
364: fn = zero;
366: if (k > 0) {
367: td = x[k-1][j][i];
368: ad = 0.5*(t0 + td);
369: dd = pow(ad,beta);
370: fd = dd*(t0 - td);
371: } else {
372: fd = zero;
373: }
375: if (k < mz-1) {
376: tu = x[k+1][j][i];
377: au = 0.5*(t0 + tu);
378: du = pow(au,beta);
379: fu = du*(tu - t0);
380: } else {
381: fu = zero;
382: }
384: } else if (k == 0) {
386: /* down boundary (interior only) */
387: tw = x[k][j][i-1];
388: aw = 0.5*(t0 + tw);
389: dw = pow(aw,beta);
390: fw = dw*(t0 - tw);
392: te = x[k][j][i+1];
393: ae = 0.5*(t0 + te);
394: de = pow(ae,beta);
395: fe = de*(te - t0);
397: ts = x[k][j-1][i];
398: as = 0.5*(t0 + ts);
399: ds = pow(as,beta);
400: fs = ds*(t0 - ts);
402: tn = x[k][j+1][i];
403: an = 0.5*(t0 + tn);
404: dn = pow(an,beta);
405: fn = dn*(tn - t0);
407: fd = zero;
409: tu = x[k+1][j][i];
410: au = 0.5*(t0 + tu);
411: du = pow(au,beta);
412: fu = du*(tu - t0);
413:
414: } else if (k == mz-1) {
416: /* up boundary (interior only) */
417: tw = x[k][j][i-1];
418: aw = 0.5*(t0 + tw);
419: dw = pow(aw,beta);
420: fw = dw*(t0 - tw);
422: te = x[k][j][i+1];
423: ae = 0.5*(t0 + te);
424: de = pow(ae,beta);
425: fe = de*(te - t0);
427: ts = x[k][j-1][i];
428: as = 0.5*(t0 + ts);
429: ds = pow(as,beta);
430: fs = ds*(t0 - ts);
432: tn = x[k][j+1][i];
433: an = 0.5*(t0 + tn);
434: dn = pow(an,beta);
435: fn = dn*(tn - t0);
437: td = x[k-1][j][i];
438: ad = 0.5*(t0 + td);
439: dd = pow(ad,beta);
440: fd = dd*(t0 - td);
442: fu = zero;
443: }
445: f[k][j][i] = - hyhzdhx*(fe-fw) - hzhxdhy*(fn-fs) - hxhydhz*(fu-fd);
446: }
447: }
448: }
449: DMDAVecRestoreArray(da,localX,&x);
450: DMDAVecRestoreArray(da,F,&f);
451: DMRestoreLocalVector(da,&localX);
452: PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
453: return(0);
454: }
455: /* -------------------- Evaluate Jacobian F(x) --------------------- */
458: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ptr)
459: {
460: AppCtx *user = (AppCtx*)ptr;
462: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
463: PetscScalar one = 1.0;
464: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
465: PetscScalar t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw;
466: PetscScalar tleft,tright,beta,td,ad,dd,tu,au,du,v[7],bm1,coef;
467: PetscScalar ***x,bn,bs,be,bw,bu,bd,gn,gs,ge,gw,gu,gd;
468: Vec localX;
469: MatStencil c[7],row;
470: Mat jac = *B;
471: DM da;
474: SNESGetDM(snes,&da);
475: DMGetLocalVector(da,&localX);
476: DMDAGetInfo(da,PETSC_NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
477: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
478: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
479: tleft = user->tleft; tright = user->tright;
480: beta = user->beta; bm1 = user->bm1; coef = user->coef;
482: /* Get ghost points */
483: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
484: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
485: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
486: DMDAVecGetArray(da,localX,&x);
488: /* Evaluate Jacobian of function */
489: for (k=zs; k<zs+zm; k++) {
490: for (j=ys; j<ys+ym; j++) {
491: for (i=xs; i<xs+xm; i++) {
492: t0 = x[k][j][i];
493: row.k = k; row.j = j; row.i = i;
494: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
496: /* general interior volume */
498: tw = x[k][j][i-1];
499: aw = 0.5*(t0 + tw);
500: bw = pow(aw,bm1);
501: /* dw = bw * aw */
502: dw = pow(aw,beta);
503: gw = coef*bw*(t0 - tw);
505: te = x[k][j][i+1];
506: ae = 0.5*(t0 + te);
507: be = pow(ae,bm1);
508: /* de = be * ae; */
509: de = pow(ae,beta);
510: ge = coef*be*(te - t0);
512: ts = x[k][j-1][i];
513: as = 0.5*(t0 + ts);
514: bs = pow(as,bm1);
515: /* ds = bs * as; */
516: ds = pow(as,beta);
517: gs = coef*bs*(t0 - ts);
518:
519: tn = x[k][j+1][i];
520: an = 0.5*(t0 + tn);
521: bn = pow(an,bm1);
522: /* dn = bn * an; */
523: dn = pow(an,beta);
524: gn = coef*bn*(tn - t0);
526: td = x[k-1][j][i];
527: ad = 0.5*(t0 + td);
528: bd = pow(ad,bm1);
529: /* dd = bd * ad; */
530: dd = pow(ad,beta);
531: gd = coef*bd*(t0 - td);
532:
533: tu = x[k+1][j][i];
534: au = 0.5*(t0 + tu);
535: bu = pow(au,bm1);
536: /* du = bu * au; */
537: du = pow(au,beta);
538: gu = coef*bu*(tu - t0);
540: c[0].k = k-1; c[0].j = j; c[0].i = i; v[0] = - hxhydhz*(dd - gd);
541: c[1].k = k; c[1].j = j-1; c[1].i = i;
542: v[1] = - hzhxdhy*(ds - gs);
543: c[2].k = k; c[2].j = j; c[2].i = i-1;
544: v[2] = - hyhzdhx*(dw - gw);
545: c[3].k = k; c[3].j = j; c[3].i = i;
546: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
547: c[4].k = k; c[4].j = j; c[4].i = i+1;
548: v[4] = - hyhzdhx*(de + ge);
549: c[5].k = k; c[5].j = j+1; c[5].i = i;
550: v[5] = - hzhxdhy*(dn + gn);
551: c[6].k = k+1; c[6].j = j; c[6].i = i;
552: v[6] = - hxhydhz*(du + gu);
553: MatSetValuesStencil(jac,1,&row,7,c,v,INSERT_VALUES);
555: } else if (i == 0) {
557: /* left-hand plane boundary */
558: tw = tleft;
559: aw = 0.5*(t0 + tw);
560: bw = pow(aw,bm1);
561: /* dw = bw * aw */
562: dw = pow(aw,beta);
563: gw = coef*bw*(t0 - tw);
564:
565: te = x[k][j][i+1];
566: ae = 0.5*(t0 + te);
567: be = pow(ae,bm1);
568: /* de = be * ae; */
569: de = pow(ae,beta);
570: ge = coef*be*(te - t0);
571:
572: /* left-hand bottom edge */
573: if (j == 0) {
575: tn = x[k][j+1][i];
576: an = 0.5*(t0 + tn);
577: bn = pow(an,bm1);
578: /* dn = bn * an; */
579: dn = pow(an,beta);
580: gn = coef*bn*(tn - t0);
581:
582: /* left-hand bottom down corner */
583: if (k == 0) {
585: tu = x[k+1][j][i];
586: au = 0.5*(t0 + tu);
587: bu = pow(au,bm1);
588: /* du = bu * au; */
589: du = pow(au,beta);
590: gu = coef*bu*(tu - t0);
591:
592: c[0].k = k; c[0].j = j; c[0].i = i;
593: v[0] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
594: c[1].k = k; c[1].j = j; c[1].i = i+1;
595: v[1] = - hyhzdhx*(de + ge);
596: c[2].k = k; c[2].j = j+1; c[2].i = i;
597: v[2] = - hzhxdhy*(dn + gn);
598: c[3].k = k+1; c[3].j = j; c[3].i = i;
599: v[3] = - hxhydhz*(du + gu);
600: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
602: /* left-hand bottom interior edge */
603: } else if (k < mz-1) {
605: tu = x[k+1][j][i];
606: au = 0.5*(t0 + tu);
607: bu = pow(au,bm1);
608: /* du = bu * au; */
609: du = pow(au,beta);
610: gu = coef*bu*(tu - t0);
611:
612: td = x[k-1][j][i];
613: ad = 0.5*(t0 + td);
614: bd = pow(ad,bm1);
615: /* dd = bd * ad; */
616: dd = pow(ad,beta);
617: gd = coef*bd*(td - t0);
618:
619: c[0].k = k-1; c[0].j = j; c[0].i = i;
620: v[0] = - hxhydhz*(dd - gd);
621: c[1].k = k; c[1].j = j; c[1].i = i;
622: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
623: c[2].k = k; c[2].j = j; c[2].i = i+1;
624: v[2] = - hyhzdhx*(de + ge);
625: c[3].k = k; c[3].j = j+1; c[3].i = i;
626: v[3] = - hzhxdhy*(dn + gn);
627: c[4].k = k+1; c[4].j = j; c[4].i = i;
628: v[4] = - hxhydhz*(du + gu);
629: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
631: /* left-hand bottom up corner */
632: } else {
634: td = x[k-1][j][i];
635: ad = 0.5*(t0 + td);
636: bd = pow(ad,bm1);
637: /* dd = bd * ad; */
638: dd = pow(ad,beta);
639: gd = coef*bd*(td - t0);
640:
641: c[0].k = k-1; c[0].j = j; c[0].i = i;
642: v[0] = - hxhydhz*(dd - gd);
643: c[1].k = k; c[1].j = j; c[1].i = i;
644: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
645: c[2].k = k; c[2].j = j; c[2].i = i+1;
646: v[2] = - hyhzdhx*(de + ge);
647: c[3].k = k; c[3].j = j+1; c[3].i = i;
648: v[3] = - hzhxdhy*(dn + gn);
649: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
650: }
652: /* left-hand top edge */
653: } else if (j == my-1) {
655: ts = x[k][j-1][i];
656: as = 0.5*(t0 + ts);
657: bs = pow(as,bm1);
658: /* ds = bs * as; */
659: ds = pow(as,beta);
660: gs = coef*bs*(ts - t0);
661:
662: /* left-hand top down corner */
663: if (k == 0) {
665: tu = x[k+1][j][i];
666: au = 0.5*(t0 + tu);
667: bu = pow(au,bm1);
668: /* du = bu * au; */
669: du = pow(au,beta);
670: gu = coef*bu*(tu - t0);
671:
672: c[0].k = k; c[0].j = j-1; c[0].i = i;
673: v[0] = - hzhxdhy*(ds - gs);
674: c[1].k = k; c[1].j = j; c[1].i = i;
675: v[1] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
676: c[2].k = k; c[2].j = j; c[2].i = i+1;
677: v[2] = - hyhzdhx*(de + ge);
678: c[3].k = k+1; c[3].j = j; c[3].i = i;
679: v[3] = - hxhydhz*(du + gu);
680: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
682: /* left-hand top interior edge */
683: } else if (k < mz-1) {
685: tu = x[k+1][j][i];
686: au = 0.5*(t0 + tu);
687: bu = pow(au,bm1);
688: /* du = bu * au; */
689: du = pow(au,beta);
690: gu = coef*bu*(tu - t0);
691:
692: td = x[k-1][j][i];
693: ad = 0.5*(t0 + td);
694: bd = pow(ad,bm1);
695: /* dd = bd * ad; */
696: dd = pow(ad,beta);
697: gd = coef*bd*(td - t0);
698:
699: c[0].k = k-1; c[0].j = j; c[0].i = i;
700: v[0] = - hxhydhz*(dd - gd);
701: c[1].k = k; c[1].j = j-1; c[1].i = i;
702: v[1] = - hzhxdhy*(ds - gs);
703: c[2].k = k; c[2].j = j; c[2].i = i;
704: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
705: c[3].k = k; c[3].j = j; c[3].i = i+1;
706: v[3] = - hyhzdhx*(de + ge);
707: c[4].k = k+1; c[4].j = j; c[4].i = i;
708: v[4] = - hxhydhz*(du + gu);
709: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
711: /* left-hand top up corner */
712: } else {
714: td = x[k-1][j][i];
715: ad = 0.5*(t0 + td);
716: bd = pow(ad,bm1);
717: /* dd = bd * ad; */
718: dd = pow(ad,beta);
719: gd = coef*bd*(td - t0);
720:
721: c[0].k = k-1; c[0].j = j; c[0].i = i;
722: v[0] = - hxhydhz*(dd - gd);
723: c[1].k = k; c[1].j = j-1; c[1].i = i;
724: v[1] = - hzhxdhy*(ds - gs);
725: c[2].k = k; c[2].j = j; c[2].i = i;
726: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
727: c[3].k = k; c[3].j = j; c[3].i = i+1;
728: v[3] = - hyhzdhx*(de + ge);
729: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
730: }
732: } else {
734: ts = x[k][j-1][i];
735: as = 0.5*(t0 + ts);
736: bs = pow(as,bm1);
737: /* ds = bs * as; */
738: ds = pow(as,beta);
739: gs = coef*bs*(t0 - ts);
740:
741: tn = x[k][j+1][i];
742: an = 0.5*(t0 + tn);
743: bn = pow(an,bm1);
744: /* dn = bn * an; */
745: dn = pow(an,beta);
746: gn = coef*bn*(tn - t0);
748: /* left-hand down interior edge */
749: if (k == 0) {
751: tu = x[k+1][j][i];
752: au = 0.5*(t0 + tu);
753: bu = pow(au,bm1);
754: /* du = bu * au; */
755: du = pow(au,beta);
756: gu = coef*bu*(tu - t0);
758: c[0].k = k; c[0].j = j-1; c[0].i = i;
759: v[0] = - hzhxdhy*(ds - gs);
760: c[1].k = k; c[1].j = j; c[1].i = i;
761: v[1] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
762: c[2].k = k; c[2].j = j; c[2].i = i+1;
763: v[2] = - hyhzdhx*(de + ge);
764: c[3].k = k; c[3].j = j+1; c[3].i = i;
765: v[3] = - hzhxdhy*(dn + gn);
766: c[4].k = k+1; c[4].j = j; c[4].i = i;
767: v[4] = - hxhydhz*(du + gu);
768: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
769: }
771: /* left-hand up interior edge */
772: else if (k == mz-1) {
774: td = x[k-1][j][i];
775: ad = 0.5*(t0 + td);
776: bd = pow(ad,bm1);
777: /* dd = bd * ad; */
778: dd = pow(ad,beta);
779: gd = coef*bd*(t0 - td);
780:
781: c[0].k = k-1; c[0].j = j; c[0].i = i;
782: v[0] = - hxhydhz*(dd - gd);
783: c[1].k = k; c[1].j = j-1; c[1].i = i;
784: v[1] = - hzhxdhy*(ds - gs);
785: c[2].k = k; c[2].j = j; c[2].i = i;
786: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
787: c[3].k = k; c[3].j = j; c[3].i = i+1;
788: v[3] = - hyhzdhx*(de + ge);
789: c[4].k = k; c[4].j = j+1; c[4].i = i;
790: v[4] = - hzhxdhy*(dn + gn);
791: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
792: }
794: /* left-hand interior plane */
795: else {
797: td = x[k-1][j][i];
798: ad = 0.5*(t0 + td);
799: bd = pow(ad,bm1);
800: /* dd = bd * ad; */
801: dd = pow(ad,beta);
802: gd = coef*bd*(t0 - td);
803:
804: tu = x[k+1][j][i];
805: au = 0.5*(t0 + tu);
806: bu = pow(au,bm1);
807: /* du = bu * au; */
808: du = pow(au,beta);
809: gu = coef*bu*(tu - t0);
811: c[0].k = k-1; c[0].j = j; c[0].i = i;
812: v[0] = - hxhydhz*(dd - gd);
813: c[1].k = k; c[1].j = j-1; c[1].i = i;
814: v[1] = - hzhxdhy*(ds - gs);
815: c[2].k = k; c[2].j = j; c[2].i = i;
816: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
817: c[3].k = k; c[3].j = j; c[3].i = i+1;
818: v[3] = - hyhzdhx*(de + ge);
819: c[4].k = k; c[4].j = j+1; c[4].i = i;
820: v[4] = - hzhxdhy*(dn + gn);
821: c[5].k = k+1; c[5].j = j; c[5].i = i;
822: v[5] = - hxhydhz*(du + gu);
823: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
824: }
825: }
827: } else if (i == mx-1) {
829: /* right-hand plane boundary */
830: tw = x[k][j][i-1];
831: aw = 0.5*(t0 + tw);
832: bw = pow(aw,bm1);
833: /* dw = bw * aw */
834: dw = pow(aw,beta);
835: gw = coef*bw*(t0 - tw);
836:
837: te = tright;
838: ae = 0.5*(t0 + te);
839: be = pow(ae,bm1);
840: /* de = be * ae; */
841: de = pow(ae,beta);
842: ge = coef*be*(te - t0);
843:
844: /* right-hand bottom edge */
845: if (j == 0) {
847: tn = x[k][j+1][i];
848: an = 0.5*(t0 + tn);
849: bn = pow(an,bm1);
850: /* dn = bn * an; */
851: dn = pow(an,beta);
852: gn = coef*bn*(tn - t0);
853:
854: /* right-hand bottom down corner */
855: if (k == 0) {
857: tu = x[k+1][j][i];
858: au = 0.5*(t0 + tu);
859: bu = pow(au,bm1);
860: /* du = bu * au; */
861: du = pow(au,beta);
862: gu = coef*bu*(tu - t0);
863:
864: c[0].k = k; c[0].j = j; c[0].i = i-1;
865: v[0] = - hyhzdhx*(dw - gw);
866: c[1].k = k; c[1].j = j; c[1].i = i;
867: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
868: c[2].k = k; c[2].j = j+1; c[2].i = i;
869: v[2] = - hzhxdhy*(dn + gn);
870: c[3].k = k+1; c[3].j = j; c[3].i = i;
871: v[3] = - hxhydhz*(du + gu);
872: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
874: /* right-hand bottom interior edge */
875: } else if (k < mz-1) {
877: tu = x[k+1][j][i];
878: au = 0.5*(t0 + tu);
879: bu = pow(au,bm1);
880: /* du = bu * au; */
881: du = pow(au,beta);
882: gu = coef*bu*(tu - t0);
883:
884: td = x[k-1][j][i];
885: ad = 0.5*(t0 + td);
886: bd = pow(ad,bm1);
887: /* dd = bd * ad; */
888: dd = pow(ad,beta);
889: gd = coef*bd*(td - t0);
890:
891: c[0].k = k-1; c[0].j = j; c[0].i = i;
892: v[0] = - hxhydhz*(dd - gd);
893: c[1].k = k; c[1].j = j; c[1].i = i-1;
894: v[1] = - hyhzdhx*(dw - gw);
895: c[2].k = k; c[2].j = j; c[2].i = i;
896: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
897: c[3].k = k; c[3].j = j+1; c[3].i = i;
898: v[3] = - hzhxdhy*(dn + gn);
899: c[4].k = k+1; c[4].j = j; c[4].i = i;
900: v[4] = - hxhydhz*(du + gu);
901: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
903: /* right-hand bottom up corner */
904: } else {
906: td = x[k-1][j][i];
907: ad = 0.5*(t0 + td);
908: bd = pow(ad,bm1);
909: /* dd = bd * ad; */
910: dd = pow(ad,beta);
911: gd = coef*bd*(td - t0);
912:
913: c[0].k = k-1; c[0].j = j; c[0].i = i;
914: v[0] = - hxhydhz*(dd - gd);
915: c[1].k = k; c[1].j = j; c[1].i = i-1;
916: v[1] = - hyhzdhx*(dw - gw);
917: c[2].k = k; c[2].j = j; c[2].i = i;
918: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
919: c[3].k = k; c[3].j = j+1; c[3].i = i;
920: v[3] = - hzhxdhy*(dn + gn);
921: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
922: }
924: /* right-hand top edge */
925: } else if (j == my-1) {
927: ts = x[k][j-1][i];
928: as = 0.5*(t0 + ts);
929: bs = pow(as,bm1);
930: /* ds = bs * as; */
931: ds = pow(as,beta);
932: gs = coef*bs*(ts - t0);
933:
934: /* right-hand top down corner */
935: if (k == 0) {
937: tu = x[k+1][j][i];
938: au = 0.5*(t0 + tu);
939: bu = pow(au,bm1);
940: /* du = bu * au; */
941: du = pow(au,beta);
942: gu = coef*bu*(tu - t0);
943:
944: c[0].k = k; c[0].j = j-1; c[0].i = i;
945: v[0] = - hzhxdhy*(ds - gs);
946: c[1].k = k; c[1].j = j; c[1].i = i-1;
947: v[1] = - hyhzdhx*(dw - gw);
948: c[2].k = k; c[2].j = j; c[2].i = i;
949: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
950: c[3].k = k+1; c[3].j = j; c[3].i = i;
951: v[3] = - hxhydhz*(du + gu);
952: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
954: /* right-hand top interior edge */
955: } else if (k < mz-1) {
957: tu = x[k+1][j][i];
958: au = 0.5*(t0 + tu);
959: bu = pow(au,bm1);
960: /* du = bu * au; */
961: du = pow(au,beta);
962: gu = coef*bu*(tu - t0);
963:
964: td = x[k-1][j][i];
965: ad = 0.5*(t0 + td);
966: bd = pow(ad,bm1);
967: /* dd = bd * ad; */
968: dd = pow(ad,beta);
969: gd = coef*bd*(td - t0);
970:
971: c[0].k = k-1; c[0].j = j; c[0].i = i;
972: v[0] = - hxhydhz*(dd - gd);
973: c[1].k = k; c[1].j = j-1; c[1].i = i;
974: v[1] = - hzhxdhy*(ds - gs);
975: c[2].k = k; c[2].j = j; c[2].i = i-1;
976: v[2] = - hyhzdhx*(dw - gw);
977: c[3].k = k; c[3].j = j; c[3].i = i;
978: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
979: c[4].k = k+1; c[4].j = j; c[4].i = i;
980: v[4] = - hxhydhz*(du + gu);
981: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
983: /* right-hand top up corner */
984: } else {
986: td = x[k-1][j][i];
987: ad = 0.5*(t0 + td);
988: bd = pow(ad,bm1);
989: /* dd = bd * ad; */
990: dd = pow(ad,beta);
991: gd = coef*bd*(td - t0);
992:
993: c[0].k = k-1; c[0].j = j; c[0].i = i;
994: v[0] = - hxhydhz*(dd - gd);
995: c[1].k = k; c[1].j = j-1; c[1].i = i;
996: v[1] = - hzhxdhy*(ds - gs);
997: c[2].k = k; c[2].j = j; c[2].i = i-1;
998: v[2] = - hyhzdhx*(dw - gw);
999: c[3].k = k; c[3].j = j; c[3].i = i;
1000: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1001: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
1002: }
1004: } else {
1006: ts = x[k][j-1][i];
1007: as = 0.5*(t0 + ts);
1008: bs = pow(as,bm1);
1009: /* ds = bs * as; */
1010: ds = pow(as,beta);
1011: gs = coef*bs*(t0 - ts);
1012:
1013: tn = x[k][j+1][i];
1014: an = 0.5*(t0 + tn);
1015: bn = pow(an,bm1);
1016: /* dn = bn * an; */
1017: dn = pow(an,beta);
1018: gn = coef*bn*(tn - t0);
1020: /* right-hand down interior edge */
1021: if (k == 0) {
1023: tu = x[k+1][j][i];
1024: au = 0.5*(t0 + tu);
1025: bu = pow(au,bm1);
1026: /* du = bu * au; */
1027: du = pow(au,beta);
1028: gu = coef*bu*(tu - t0);
1030: c[0].k = k; c[0].j = j-1; c[0].i = i;
1031: v[0] = - hzhxdhy*(ds - gs);
1032: c[1].k = k; c[1].j = j; c[1].i = i-1;
1033: v[1] = - hyhzdhx*(dw - gw);
1034: c[2].k = k; c[2].j = j; c[2].i = i;
1035: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1036: c[3].k = k; c[3].j = j+1; c[3].i = i;
1037: v[3] = - hzhxdhy*(dn + gn);
1038: c[4].k = k+1; c[4].j = j; c[4].i = i;
1039: v[4] = - hxhydhz*(du + gu);
1040: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1041: }
1043: /* right-hand up interior edge */
1044: else if (k == mz-1) {
1046: td = x[k-1][j][i];
1047: ad = 0.5*(t0 + td);
1048: bd = pow(ad,bm1);
1049: /* dd = bd * ad; */
1050: dd = pow(ad,beta);
1051: gd = coef*bd*(t0 - td);
1052:
1053: c[0].k = k-1; c[0].j = j; c[0].i = i;
1054: v[0] = - hxhydhz*(dd - gd);
1055: c[1].k = k; c[1].j = j-1; c[1].i = i;
1056: v[1] = - hzhxdhy*(ds - gs);
1057: c[2].k = k; c[2].j = j; c[2].i = i-1;
1058: v[2] = - hyhzdhx*(dw - gw);
1059: c[3].k = k; c[3].j = j; c[3].i = i;
1060: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1061: c[4].k = k; c[4].j = j+1; c[4].i = i;
1062: v[4] = - hzhxdhy*(dn + gn);
1063: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1064: }
1066: /* right-hand interior plane */
1067: else {
1069: td = x[k-1][j][i];
1070: ad = 0.5*(t0 + td);
1071: bd = pow(ad,bm1);
1072: /* dd = bd * ad; */
1073: dd = pow(ad,beta);
1074: gd = coef*bd*(t0 - td);
1075:
1076: tu = x[k+1][j][i];
1077: au = 0.5*(t0 + tu);
1078: bu = pow(au,bm1);
1079: /* du = bu * au; */
1080: du = pow(au,beta);
1081: gu = coef*bu*(tu - t0);
1083: c[0].k = k-1; c[0].j = j; c[0].i = i;
1084: v[0] = - hxhydhz*(dd - gd);
1085: c[1].k = k; c[1].j = j-1; c[1].i = i;
1086: v[1] = - hzhxdhy*(ds - gs);
1087: c[2].k = k; c[2].j = j; c[2].i = i-1;
1088: v[2] = - hyhzdhx*(dw - gw);
1089: c[3].k = k; c[3].j = j; c[3].i = i;
1090: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1091: c[4].k = k; c[4].j = j+1; c[4].i = i;
1092: v[4] = - hzhxdhy*(dn + gn);
1093: c[5].k = k+1; c[5].j = j; c[5].i = i;
1094: v[5] = - hxhydhz*(du + gu);
1095: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1096: }
1097: }
1099: } else if (j == 0) {
1101: tw = x[k][j][i-1];
1102: aw = 0.5*(t0 + tw);
1103: bw = pow(aw,bm1);
1104: /* dw = bw * aw */
1105: dw = pow(aw,beta);
1106: gw = coef*bw*(t0 - tw);
1108: te = x[k][j][i+1];
1109: ae = 0.5*(t0 + te);
1110: be = pow(ae,bm1);
1111: /* de = be * ae; */
1112: de = pow(ae,beta);
1113: ge = coef*be*(te - t0);
1115: tn = x[k][j+1][i];
1116: an = 0.5*(t0 + tn);
1117: bn = pow(an,bm1);
1118: /* dn = bn * an; */
1119: dn = pow(an,beta);
1120: gn = coef*bn*(tn - t0);
1123: /* bottom down interior edge */
1124: if (k == 0) {
1126: tu = x[k+1][j][i];
1127: au = 0.5*(t0 + tu);
1128: bu = pow(au,bm1);
1129: /* du = bu * au; */
1130: du = pow(au,beta);
1131: gu = coef*bu*(tu - t0);
1133: c[0].k = k; c[0].j = j; c[0].i = i-1;
1134: v[0] = - hyhzdhx*(dw - gw);
1135: c[1].k = k; c[1].j = j; c[1].i = i;
1136: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1137: c[2].k = k; c[2].j = j; c[2].i = i+1;
1138: v[2] = - hyhzdhx*(de + ge);
1139: c[3].k = k; c[3].j = j+1; c[3].i = i;
1140: v[3] = - hzhxdhy*(dn + gn);
1141: c[4].k = k+1; c[4].j = j; c[4].i = i;
1142: v[4] = - hxhydhz*(du + gu);
1143: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1144: }
1146: /* bottom up interior edge */
1147: else if (k == mz-1) {
1149: td = x[k-1][j][i];
1150: ad = 0.5*(t0 + td);
1151: bd = pow(ad,bm1);
1152: /* dd = bd * ad; */
1153: dd = pow(ad,beta);
1154: gd = coef*bd*(td - t0);
1156: c[0].k = k-1; c[0].j = j; c[0].i = i;
1157: v[0] = - hxhydhz*(dd - gd);
1158: c[1].k = k; c[1].j = j; c[1].i = i-1;
1159: v[1] = - hyhzdhx*(dw - gw);
1160: c[2].k = k; c[2].j = j; c[2].i = i;
1161: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1162: c[3].k = k; c[3].j = j; c[3].i = i+1;
1163: v[3] = - hyhzdhx*(de + ge);
1164: c[4].k = k; c[4].j = j+1; c[4].i = i;
1165: v[4] = - hzhxdhy*(dn + gn);
1166: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1167: }
1169: /* bottom interior plane */
1170: else {
1172: tu = x[k+1][j][i];
1173: au = 0.5*(t0 + tu);
1174: bu = pow(au,bm1);
1175: /* du = bu * au; */
1176: du = pow(au,beta);
1177: gu = coef*bu*(tu - t0);
1179: td = x[k-1][j][i];
1180: ad = 0.5*(t0 + td);
1181: bd = pow(ad,bm1);
1182: /* dd = bd * ad; */
1183: dd = pow(ad,beta);
1184: gd = coef*bd*(td - t0);
1186: c[0].k = k-1; c[0].j = j; c[0].i = i;
1187: v[0] = - hxhydhz*(dd - gd);
1188: c[1].k = k; c[1].j = j; c[1].i = i-1;
1189: v[1] = - hyhzdhx*(dw - gw);
1190: c[2].k = k; c[2].j = j; c[2].i = i;
1191: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1192: c[3].k = k; c[3].j = j; c[3].i = i+1;
1193: v[3] = - hyhzdhx*(de + ge);
1194: c[4].k = k; c[4].j = j+1; c[4].i = i;
1195: v[4] = - hzhxdhy*(dn + gn);
1196: c[5].k = k+1; c[5].j = j; c[5].i = i;
1197: v[5] = - hxhydhz*(du + gu);
1198: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1199: }
1201: } else if (j == my-1) {
1203: tw = x[k][j][i-1];
1204: aw = 0.5*(t0 + tw);
1205: bw = pow(aw,bm1);
1206: /* dw = bw * aw */
1207: dw = pow(aw,beta);
1208: gw = coef*bw*(t0 - tw);
1210: te = x[k][j][i+1];
1211: ae = 0.5*(t0 + te);
1212: be = pow(ae,bm1);
1213: /* de = be * ae; */
1214: de = pow(ae,beta);
1215: ge = coef*be*(te - t0);
1217: ts = x[k][j-1][i];
1218: as = 0.5*(t0 + ts);
1219: bs = pow(as,bm1);
1220: /* ds = bs * as; */
1221: ds = pow(as,beta);
1222: gs = coef*bs*(t0 - ts);
1223:
1224: /* top down interior edge */
1225: if (k == 0) {
1227: tu = x[k+1][j][i];
1228: au = 0.5*(t0 + tu);
1229: bu = pow(au,bm1);
1230: /* du = bu * au; */
1231: du = pow(au,beta);
1232: gu = coef*bu*(tu - t0);
1234: c[0].k = k; c[0].j = j-1; c[0].i = i;
1235: v[0] = - hzhxdhy*(ds - gs);
1236: c[1].k = k; c[1].j = j; c[1].i = i-1;
1237: v[1] = - hyhzdhx*(dw - gw);
1238: c[2].k = k; c[2].j = j; c[2].i = i;
1239: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1240: c[3].k = k; c[3].j = j; c[3].i = i+1;
1241: v[3] = - hyhzdhx*(de + ge);
1242: c[4].k = k+1; c[4].j = j; c[4].i = i;
1243: v[4] = - hxhydhz*(du + gu);
1244: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1245: }
1247: /* top up interior edge */
1248: else if (k == mz-1) {
1250: td = x[k-1][j][i];
1251: ad = 0.5*(t0 + td);
1252: bd = pow(ad,bm1);
1253: /* dd = bd * ad; */
1254: dd = pow(ad,beta);
1255: gd = coef*bd*(td - t0);
1257: c[0].k = k-1; c[0].j = j; c[0].i = i;
1258: v[0] = - hxhydhz*(dd - gd);
1259: c[1].k = k; c[1].j = j-1; c[1].i = i;
1260: v[1] = - hzhxdhy*(ds - gs);
1261: c[2].k = k; c[2].j = j; c[2].i = i-1;
1262: v[2] = - hyhzdhx*(dw - gw);
1263: c[3].k = k; c[3].j = j; c[3].i = i;
1264: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1265: c[4].k = k; c[4].j = j; c[4].i = i+1;
1266: v[4] = - hyhzdhx*(de + ge);
1267: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1268: }
1270: /* top interior plane */
1271: else {
1273: tu = x[k+1][j][i];
1274: au = 0.5*(t0 + tu);
1275: bu = pow(au,bm1);
1276: /* du = bu * au; */
1277: du = pow(au,beta);
1278: gu = coef*bu*(tu - t0);
1280: td = x[k-1][j][i];
1281: ad = 0.5*(t0 + td);
1282: bd = pow(ad,bm1);
1283: /* dd = bd * ad; */
1284: dd = pow(ad,beta);
1285: gd = coef*bd*(td - t0);
1287: c[0].k = k-1; c[0].j = j; c[0].i = i;
1288: v[0] = - hxhydhz*(dd - gd);
1289: c[1].k = k; c[1].j = j-1; c[1].i = i;
1290: v[1] = - hzhxdhy*(ds - gs);
1291: c[2].k = k; c[2].j = j; c[2].i = i-1;
1292: v[2] = - hyhzdhx*(dw - gw);
1293: c[3].k = k; c[3].j = j; c[3].i = i;
1294: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1295: c[4].k = k; c[4].j = j; c[4].i = i+1;
1296: v[4] = - hyhzdhx*(de + ge);
1297: c[5].k = k+1; c[5].j = j; c[5].i = i;
1298: v[5] = - hxhydhz*(du + gu);
1299: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1300: }
1302: } else if (k == 0) {
1304: /* down interior plane */
1306: tw = x[k][j][i-1];
1307: aw = 0.5*(t0 + tw);
1308: bw = pow(aw,bm1);
1309: /* dw = bw * aw */
1310: dw = pow(aw,beta);
1311: gw = coef*bw*(t0 - tw);
1313: te = x[k][j][i+1];
1314: ae = 0.5*(t0 + te);
1315: be = pow(ae,bm1);
1316: /* de = be * ae; */
1317: de = pow(ae,beta);
1318: ge = coef*be*(te - t0);
1320: ts = x[k][j-1][i];
1321: as = 0.5*(t0 + ts);
1322: bs = pow(as,bm1);
1323: /* ds = bs * as; */
1324: ds = pow(as,beta);
1325: gs = coef*bs*(t0 - ts);
1326:
1327: tn = x[k][j+1][i];
1328: an = 0.5*(t0 + tn);
1329: bn = pow(an,bm1);
1330: /* dn = bn * an; */
1331: dn = pow(an,beta);
1332: gn = coef*bn*(tn - t0);
1333:
1334: tu = x[k+1][j][i];
1335: au = 0.5*(t0 + tu);
1336: bu = pow(au,bm1);
1337: /* du = bu * au; */
1338: du = pow(au,beta);
1339: gu = coef*bu*(tu - t0);
1341: c[0].k = k; c[0].j = j-1; c[0].i = i;
1342: v[0] = - hzhxdhy*(ds - gs);
1343: c[1].k = k; c[1].j = j; c[1].i = i-1;
1344: v[1] = - hyhzdhx*(dw - gw);
1345: c[2].k = k; c[2].j = j; c[2].i = i;
1346: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1347: c[3].k = k; c[3].j = j; c[3].i = i+1;
1348: v[3] = - hyhzdhx*(de + ge);
1349: c[4].k = k; c[4].j = j+1; c[4].i = i;
1350: v[4] = - hzhxdhy*(dn + gn);
1351: c[5].k = k+1; c[5].j = j; c[5].i = i;
1352: v[5] = - hxhydhz*(du + gu);
1353: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1354: }
1355:
1356: else if (k == mz-1) {
1358: /* up interior plane */
1360: tw = x[k][j][i-1];
1361: aw = 0.5*(t0 + tw);
1362: bw = pow(aw,bm1);
1363: /* dw = bw * aw */
1364: dw = pow(aw,beta);
1365: gw = coef*bw*(t0 - tw);
1367: te = x[k][j][i+1];
1368: ae = 0.5*(t0 + te);
1369: be = pow(ae,bm1);
1370: /* de = be * ae; */
1371: de = pow(ae,beta);
1372: ge = coef*be*(te - t0);
1374: ts = x[k][j-1][i];
1375: as = 0.5*(t0 + ts);
1376: bs = pow(as,bm1);
1377: /* ds = bs * as; */
1378: ds = pow(as,beta);
1379: gs = coef*bs*(t0 - ts);
1380:
1381: tn = x[k][j+1][i];
1382: an = 0.5*(t0 + tn);
1383: bn = pow(an,bm1);
1384: /* dn = bn * an; */
1385: dn = pow(an,beta);
1386: gn = coef*bn*(tn - t0);
1387:
1388: td = x[k-1][j][i];
1389: ad = 0.5*(t0 + td);
1390: bd = pow(ad,bm1);
1391: /* dd = bd * ad; */
1392: dd = pow(ad,beta);
1393: gd = coef*bd*(t0 - td);
1394:
1395: c[0].k = k-1; c[0].j = j; c[0].i = i;
1396: v[0] = - hxhydhz*(dd - gd);
1397: c[1].k = k; c[1].j = j-1; c[1].i = i;
1398: v[1] = - hzhxdhy*(ds - gs);
1399: c[2].k = k; c[2].j = j; c[2].i = i-1;
1400: v[2] = - hyhzdhx*(dw - gw);
1401: c[3].k = k; c[3].j = j; c[3].i = i;
1402: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1403: c[4].k = k; c[4].j = j; c[4].i = i+1;
1404: v[4] = - hyhzdhx*(de + ge);
1405: c[5].k = k; c[5].j = j+1; c[5].i = i;
1406: v[5] = - hzhxdhy*(dn + gn);
1407: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1408: }
1409: }
1410: }
1411: }
1412: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
1413: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
1414: if (jac != *J) {
1415: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1416: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1417: }
1418: DMDAVecRestoreArray(da,localX,&x);
1419: DMRestoreLocalVector(da,&localX);
1421: PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
1422: return(0);
1423: }