Actual source code: ex20.c
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: DA^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 petscda.h
44: #include petscmg.h
45: #include petscdmmg.h
47: /* User-defined application context */
49: typedef struct {
50: PetscReal tleft,tright; /* Dirichlet boundary conditions */
51: PetscReal beta,bm1,coef; /* nonlinear diffusivity parameterizations */
52: } AppCtx;
54: #define POWFLOP 5 /* assume a pow() takes five flops */
62: int main(int argc,char **argv)
63: {
64: DMMG *dmmg;
65: SNES snes;
66: AppCtx user;
68: PetscInt its,lits;
69: PetscReal litspit;
70: DA da;
72: PetscInitialize(&argc,&argv,PETSC_NULL,help);
74: /* set problem parameters */
75: user.tleft = 1.0;
76: user.tright = 0.1;
77: user.beta = 2.5;
78: PetscOptionsGetReal(PETSC_NULL,"-tleft",&user.tleft,PETSC_NULL);
79: PetscOptionsGetReal(PETSC_NULL,"-tright",&user.tright,PETSC_NULL);
80: PetscOptionsGetReal(PETSC_NULL,"-beta",&user.beta,PETSC_NULL);
81: user.bm1 = user.beta - 1.0;
82: user.coef = user.beta/2.0;
85: /*
86: Create the multilevel DA data structure
87: */
88: DMMGCreate(PETSC_COMM_WORLD,2,&user,&dmmg);
90: /*
91: Set the DA (grid structure) for the grids.
92: */
93: DACreate3d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,-5,-5,-5,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);
94: DMMGSetDM(dmmg,(DM)da);
95: DADestroy(da);
97: /*
98: Create the nonlinear solver, and tell the DMMG structure to use it
99: */
100: DMMGSetSNES(dmmg,FormFunction,FormJacobian);
102: /*
103: PreLoadBegin() means that the following section of code is run twice. The first time
104: through the flag PreLoading is on this the nonlinear solver is only run for a single step.
105: The second time through (the actually timed code) the maximum iterations is set to 10
106: Preload of the executable is done to eliminate from the timing the time spent bring the
107: executable into memory from disk (paging in).
108: */
109: PreLoadBegin(PETSC_TRUE,"Solve");
110: DMMGSetInitialGuess(dmmg,FormInitialGuess);
111: DMMGSolve(dmmg);
112: PreLoadEnd();
113: snes = DMMGGetSNES(dmmg);
114: SNESGetIterationNumber(snes,&its);
115: SNESGetNumberLinearIterations(snes,&lits);
116: litspit = ((PetscReal)lits)/((PetscReal)its);
117: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);
118: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
119: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / Newton = %e\n",litspit);
121: DMMGDestroy(dmmg);
122: PetscFinalize();
124: return 0;
125: }
126: /* -------------------- Form initial approximation ----------------- */
129: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
130: {
131: AppCtx *user = (AppCtx*)dmmg->user;
132: PetscInt i,j,k,xs,ys,xm,ym,zs,zm;
134: PetscReal tleft = user->tleft;
135: PetscScalar ***x;
139: /* Get ghost points */
140: DAGetCorners((DA)dmmg->dm,&xs,&ys,&zs,&xm,&ym,&zm);
141: DAVecGetArray((DA)dmmg->dm,X,&x);
143: /* Compute initial guess */
144: for (k=zs; k<zs+zm; k++) {
145: for (j=ys; j<ys+ym; j++) {
146: for (i=xs; i<xs+xm; i++) {
147: x[k][j][i] = tleft;
148: }
149: }
150: }
151: DAVecRestoreArray((DA)dmmg->dm,X,&x);
152: return(0);
153: }
154: /* -------------------- Evaluate Function F(x) --------------------- */
157: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ptr)
158: {
159: DMMG dmmg = (DMMG)ptr;
160: AppCtx *user = (AppCtx*)dmmg->user;
162: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
163: PetscScalar zero = 0.0,one = 1.0;
164: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
165: 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;
166: PetscScalar tleft,tright,beta,td,ad,dd,fd,tu,au,du,fu;
167: PetscScalar ***x,***f;
168: Vec localX;
171: DAGetLocalVector((DA)dmmg->dm,&localX);
172: DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,&mz,0,0,0,0,0,0,0);
173: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
174: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
175: tleft = user->tleft; tright = user->tright;
176: beta = user->beta;
177:
178: /* Get ghost points */
179: DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
180: DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
181: DAGetCorners((DA)dmmg->dm,&xs,&ys,&zs,&xm,&ym,&zm);
182: DAVecGetArray((DA)dmmg->dm,localX,&x);
183: DAVecGetArray((DA)dmmg->dm,F,&f);
185: /* Evaluate function */
186: for (k=zs; k<zs+zm; k++) {
187: for (j=ys; j<ys+ym; j++) {
188: for (i=xs; i<xs+xm; i++) {
189: t0 = x[k][j][i];
191: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
193: /* general interior volume */
195: tw = x[k][j][i-1];
196: aw = 0.5*(t0 + tw);
197: dw = pow(aw,beta);
198: fw = dw*(t0 - tw);
200: te = x[k][j][i+1];
201: ae = 0.5*(t0 + te);
202: de = pow(ae,beta);
203: fe = de*(te - t0);
205: ts = x[k][j-1][i];
206: as = 0.5*(t0 + ts);
207: ds = pow(as,beta);
208: fs = ds*(t0 - ts);
209:
210: tn = x[k][j+1][i];
211: an = 0.5*(t0 + tn);
212: dn = pow(an,beta);
213: fn = dn*(tn - t0);
215: td = x[k-1][j][i];
216: ad = 0.5*(t0 + td);
217: dd = pow(ad,beta);
218: fd = dd*(t0 - td);
220: tu = x[k+1][j][i];
221: au = 0.5*(t0 + tu);
222: du = pow(au,beta);
223: fu = du*(tu - t0);
225: } else if (i == 0) {
227: /* left-hand (west) boundary */
228: tw = tleft;
229: aw = 0.5*(t0 + tw);
230: dw = pow(aw,beta);
231: fw = dw*(t0 - tw);
233: te = x[k][j][i+1];
234: ae = 0.5*(t0 + te);
235: de = pow(ae,beta);
236: fe = de*(te - t0);
238: if (j > 0) {
239: ts = x[k][j-1][i];
240: as = 0.5*(t0 + ts);
241: ds = pow(as,beta);
242: fs = ds*(t0 - ts);
243: } else {
244: fs = zero;
245: }
247: if (j < my-1) {
248: tn = x[k][j+1][i];
249: an = 0.5*(t0 + tn);
250: dn = pow(an,beta);
251: fn = dn*(tn - t0);
252: } else {
253: fn = zero;
254: }
256: if (k > 0) {
257: td = x[k-1][j][i];
258: ad = 0.5*(t0 + td);
259: dd = pow(ad,beta);
260: fd = dd*(t0 - td);
261: } else {
262: fd = zero;
263: }
265: if (k < mz-1) {
266: tu = x[k+1][j][i];
267: au = 0.5*(t0 + tu);
268: du = pow(au,beta);
269: fu = du*(tu - t0);
270: } else {
271: fu = zero;
272: }
274: } else if (i == mx-1) {
276: /* right-hand (east) boundary */
277: tw = x[k][j][i-1];
278: aw = 0.5*(t0 + tw);
279: dw = pow(aw,beta);
280: fw = dw*(t0 - tw);
281:
282: te = tright;
283: ae = 0.5*(t0 + te);
284: de = pow(ae,beta);
285: fe = de*(te - t0);
286:
287: if (j > 0) {
288: ts = x[k][j-1][i];
289: as = 0.5*(t0 + ts);
290: ds = pow(as,beta);
291: fs = ds*(t0 - ts);
292: } else {
293: fs = zero;
294: }
295:
296: if (j < my-1) {
297: tn = x[k][j+1][i];
298: an = 0.5*(t0 + tn);
299: dn = pow(an,beta);
300: fn = dn*(tn - t0);
301: } else {
302: fn = zero;
303: }
305: if (k > 0) {
306: td = x[k-1][j][i];
307: ad = 0.5*(t0 + td);
308: dd = pow(ad,beta);
309: fd = dd*(t0 - td);
310: } else {
311: fd = zero;
312: }
314: if (k < mz-1) {
315: tu = x[k+1][j][i];
316: au = 0.5*(t0 + tu);
317: du = pow(au,beta);
318: fu = du*(tu - t0);
319: } else {
320: fu = zero;
321: }
323: } else if (j == 0) {
325: /* bottom (south) boundary, and i <> 0 or mx-1 */
326: tw = x[k][j][i-1];
327: aw = 0.5*(t0 + tw);
328: dw = pow(aw,beta);
329: fw = dw*(t0 - tw);
331: te = x[k][j][i+1];
332: ae = 0.5*(t0 + te);
333: de = pow(ae,beta);
334: fe = de*(te - t0);
336: fs = zero;
338: tn = x[k][j+1][i];
339: an = 0.5*(t0 + tn);
340: dn = pow(an,beta);
341: fn = dn*(tn - t0);
343: if (k > 0) {
344: td = x[k-1][j][i];
345: ad = 0.5*(t0 + td);
346: dd = pow(ad,beta);
347: fd = dd*(t0 - td);
348: } else {
349: fd = zero;
350: }
352: if (k < mz-1) {
353: tu = x[k+1][j][i];
354: au = 0.5*(t0 + tu);
355: du = pow(au,beta);
356: fu = du*(tu - t0);
357: } else {
358: fu = zero;
359: }
361: } else if (j == my-1) {
363: /* top (north) boundary, and i <> 0 or mx-1 */
364: tw = x[k][j][i-1];
365: aw = 0.5*(t0 + tw);
366: dw = pow(aw,beta);
367: fw = dw*(t0 - tw);
369: te = x[k][j][i+1];
370: ae = 0.5*(t0 + te);
371: de = pow(ae,beta);
372: fe = de*(te - t0);
374: ts = x[k][j-1][i];
375: as = 0.5*(t0 + ts);
376: ds = pow(as,beta);
377: fs = ds*(t0 - ts);
379: fn = zero;
381: if (k > 0) {
382: td = x[k-1][j][i];
383: ad = 0.5*(t0 + td);
384: dd = pow(ad,beta);
385: fd = dd*(t0 - td);
386: } else {
387: fd = zero;
388: }
390: if (k < mz-1) {
391: tu = x[k+1][j][i];
392: au = 0.5*(t0 + tu);
393: du = pow(au,beta);
394: fu = du*(tu - t0);
395: } else {
396: fu = zero;
397: }
399: } else if (k == 0) {
401: /* down boundary (interior only) */
402: tw = x[k][j][i-1];
403: aw = 0.5*(t0 + tw);
404: dw = pow(aw,beta);
405: fw = dw*(t0 - tw);
407: te = x[k][j][i+1];
408: ae = 0.5*(t0 + te);
409: de = pow(ae,beta);
410: fe = de*(te - t0);
412: ts = x[k][j-1][i];
413: as = 0.5*(t0 + ts);
414: ds = pow(as,beta);
415: fs = ds*(t0 - ts);
417: tn = x[k][j+1][i];
418: an = 0.5*(t0 + tn);
419: dn = pow(an,beta);
420: fn = dn*(tn - t0);
422: fd = zero;
424: tu = x[k+1][j][i];
425: au = 0.5*(t0 + tu);
426: du = pow(au,beta);
427: fu = du*(tu - t0);
428:
429: } else if (k == mz-1) {
431: /* up boundary (interior only) */
432: tw = x[k][j][i-1];
433: aw = 0.5*(t0 + tw);
434: dw = pow(aw,beta);
435: fw = dw*(t0 - tw);
437: te = x[k][j][i+1];
438: ae = 0.5*(t0 + te);
439: de = pow(ae,beta);
440: fe = de*(te - t0);
442: ts = x[k][j-1][i];
443: as = 0.5*(t0 + ts);
444: ds = pow(as,beta);
445: fs = ds*(t0 - ts);
447: tn = x[k][j+1][i];
448: an = 0.5*(t0 + tn);
449: dn = pow(an,beta);
450: fn = dn*(tn - t0);
452: td = x[k-1][j][i];
453: ad = 0.5*(t0 + td);
454: dd = pow(ad,beta);
455: fd = dd*(t0 - td);
457: fu = zero;
458: }
460: f[k][j][i] = - hyhzdhx*(fe-fw) - hzhxdhy*(fn-fs) - hxhydhz*(fu-fd);
461: }
462: }
463: }
464: DAVecRestoreArray((DA)dmmg->dm,localX,&x);
465: DAVecRestoreArray((DA)dmmg->dm,F,&f);
466: DARestoreLocalVector((DA)dmmg->dm,&localX);
467: PetscLogFlops((22 + 4*POWFLOP)*ym*xm);
468: return(0);
469: }
470: /* -------------------- Evaluate Jacobian F(x) --------------------- */
473: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ptr)
474: {
475: DMMG dmmg = (DMMG)ptr;
476: AppCtx *user = (AppCtx*)dmmg->user;
478: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
479: PetscScalar one = 1.0;
480: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
481: PetscScalar t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw;
482: PetscScalar tleft,tright,beta,td,ad,dd,tu,au,du,v[7],bm1,coef;
483: PetscScalar ***x,bn,bs,be,bw,bu,bd,gn,gs,ge,gw,gu,gd;
484: Vec localX;
485: MatStencil c[7],row;
486: Mat jac = *B;
489: DAGetLocalVector((DA)dmmg->dm,&localX);
490: DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,&mz,0,0,0,0,0,0,0);
491: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
492: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
493: tleft = user->tleft; tright = user->tright;
494: beta = user->beta; bm1 = user->bm1; coef = user->coef;
496: /* Get ghost points */
497: DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
498: DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
499: DAGetCorners((DA)dmmg->dm,&xs,&ys,&zs,&xm,&ym,&zm);
500: DAVecGetArray((DA)dmmg->dm,localX,&x);
502: /* Evaluate Jacobian of function */
503: for (k=zs; k<zs+zm; k++) {
504: for (j=ys; j<ys+ym; j++) {
505: for (i=xs; i<xs+xm; i++) {
506: t0 = x[k][j][i];
507: row.k = k; row.j = j; row.i = i;
508: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
510: /* general interior volume */
512: tw = x[k][j][i-1];
513: aw = 0.5*(t0 + tw);
514: bw = pow(aw,bm1);
515: /* dw = bw * aw */
516: dw = pow(aw,beta);
517: gw = coef*bw*(t0 - tw);
519: te = x[k][j][i+1];
520: ae = 0.5*(t0 + te);
521: be = pow(ae,bm1);
522: /* de = be * ae; */
523: de = pow(ae,beta);
524: ge = coef*be*(te - t0);
526: ts = x[k][j-1][i];
527: as = 0.5*(t0 + ts);
528: bs = pow(as,bm1);
529: /* ds = bs * as; */
530: ds = pow(as,beta);
531: gs = coef*bs*(t0 - ts);
532:
533: tn = x[k][j+1][i];
534: an = 0.5*(t0 + tn);
535: bn = pow(an,bm1);
536: /* dn = bn * an; */
537: dn = pow(an,beta);
538: gn = coef*bn*(tn - t0);
540: td = x[k-1][j][i];
541: ad = 0.5*(t0 + td);
542: bd = pow(ad,bm1);
543: /* dd = bd * ad; */
544: dd = pow(ad,beta);
545: gd = coef*bd*(t0 - td);
546:
547: tu = x[k+1][j][i];
548: au = 0.5*(t0 + tu);
549: bu = pow(au,bm1);
550: /* du = bu * au; */
551: du = pow(au,beta);
552: gu = coef*bu*(tu - t0);
554: c[0].k = k-1; c[0].j = j; c[0].i = i; v[0] = - hxhydhz*(dd - gd);
555: c[1].k = k; c[1].j = j-1; c[1].i = i;
556: v[1] = - hzhxdhy*(ds - gs);
557: c[2].k = k; c[2].j = j; c[2].i = i-1;
558: v[2] = - hyhzdhx*(dw - gw);
559: c[3].k = k; c[3].j = j; c[3].i = i;
560: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
561: c[4].k = k; c[4].j = j; c[4].i = i+1;
562: v[4] = - hyhzdhx*(de + ge);
563: c[5].k = k; c[5].j = j+1; c[5].i = i;
564: v[5] = - hzhxdhy*(dn + gn);
565: c[6].k = k+1; c[6].j = j; c[6].i = i;
566: v[6] = - hxhydhz*(du + gu);
567: MatSetValuesStencil(jac,1,&row,7,c,v,INSERT_VALUES);
569: } else if (i == 0) {
571: /* left-hand plane boundary */
572: tw = tleft;
573: aw = 0.5*(t0 + tw);
574: bw = pow(aw,bm1);
575: /* dw = bw * aw */
576: dw = pow(aw,beta);
577: gw = coef*bw*(t0 - tw);
578:
579: te = x[k][j][i+1];
580: ae = 0.5*(t0 + te);
581: be = pow(ae,bm1);
582: /* de = be * ae; */
583: de = pow(ae,beta);
584: ge = coef*be*(te - t0);
585:
586: /* left-hand bottom edge */
587: if (j == 0) {
589: tn = x[k][j+1][i];
590: an = 0.5*(t0 + tn);
591: bn = pow(an,bm1);
592: /* dn = bn * an; */
593: dn = pow(an,beta);
594: gn = coef*bn*(tn - t0);
595:
596: /* left-hand bottom down corner */
597: if (k == 0) {
599: tu = x[k+1][j][i];
600: au = 0.5*(t0 + tu);
601: bu = pow(au,bm1);
602: /* du = bu * au; */
603: du = pow(au,beta);
604: gu = coef*bu*(tu - t0);
605:
606: c[0].k = k; c[0].j = j; c[0].i = i;
607: v[0] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
608: c[1].k = k; c[1].j = j; c[1].i = i+1;
609: v[1] = - hyhzdhx*(de + ge);
610: c[2].k = k; c[2].j = j+1; c[2].i = i;
611: v[2] = - hzhxdhy*(dn + gn);
612: c[3].k = k+1; c[3].j = j; c[3].i = i;
613: v[3] = - hxhydhz*(du + gu);
614: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
616: /* left-hand bottom interior edge */
617: } else if (k < mz-1) {
619: tu = x[k+1][j][i];
620: au = 0.5*(t0 + tu);
621: bu = pow(au,bm1);
622: /* du = bu * au; */
623: du = pow(au,beta);
624: gu = coef*bu*(tu - t0);
625:
626: td = x[k-1][j][i];
627: ad = 0.5*(t0 + td);
628: bd = pow(ad,bm1);
629: /* dd = bd * ad; */
630: dd = pow(ad,beta);
631: gd = coef*bd*(td - t0);
632:
633: c[0].k = k-1; c[0].j = j; c[0].i = i;
634: v[0] = - hxhydhz*(dd - gd);
635: c[1].k = k; c[1].j = j; c[1].i = i;
636: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
637: c[2].k = k; c[2].j = j; c[2].i = i+1;
638: v[2] = - hyhzdhx*(de + ge);
639: c[3].k = k; c[3].j = j+1; c[3].i = i;
640: v[3] = - hzhxdhy*(dn + gn);
641: c[4].k = k+1; c[4].j = j; c[4].i = i;
642: v[4] = - hxhydhz*(du + gu);
643: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
645: /* left-hand bottom up corner */
646: } else {
648: td = x[k-1][j][i];
649: ad = 0.5*(t0 + td);
650: bd = pow(ad,bm1);
651: /* dd = bd * ad; */
652: dd = pow(ad,beta);
653: gd = coef*bd*(td - t0);
654:
655: c[0].k = k-1; c[0].j = j; c[0].i = i;
656: v[0] = - hxhydhz*(dd - gd);
657: c[1].k = k; c[1].j = j; c[1].i = i;
658: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
659: c[2].k = k; c[2].j = j; c[2].i = i+1;
660: v[2] = - hyhzdhx*(de + ge);
661: c[3].k = k; c[3].j = j+1; c[3].i = i;
662: v[3] = - hzhxdhy*(dn + gn);
663: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
664: }
666: /* left-hand top edge */
667: } else if (j == my-1) {
669: ts = x[k][j-1][i];
670: as = 0.5*(t0 + ts);
671: bs = pow(as,bm1);
672: /* ds = bs * as; */
673: ds = pow(as,beta);
674: gs = coef*bs*(ts - t0);
675:
676: /* left-hand top down corner */
677: if (k == 0) {
679: tu = x[k+1][j][i];
680: au = 0.5*(t0 + tu);
681: bu = pow(au,bm1);
682: /* du = bu * au; */
683: du = pow(au,beta);
684: gu = coef*bu*(tu - t0);
685:
686: c[0].k = k; c[0].j = j-1; c[0].i = i;
687: v[0] = - hzhxdhy*(ds - gs);
688: c[1].k = k; c[1].j = j; c[1].i = i;
689: v[1] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
690: c[2].k = k; c[2].j = j; c[2].i = i+1;
691: v[2] = - hyhzdhx*(de + ge);
692: c[3].k = k+1; c[3].j = j; c[3].i = i;
693: v[3] = - hxhydhz*(du + gu);
694: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
696: /* left-hand top interior edge */
697: } else if (k < mz-1) {
699: tu = x[k+1][j][i];
700: au = 0.5*(t0 + tu);
701: bu = pow(au,bm1);
702: /* du = bu * au; */
703: du = pow(au,beta);
704: gu = coef*bu*(tu - t0);
705:
706: td = x[k-1][j][i];
707: ad = 0.5*(t0 + td);
708: bd = pow(ad,bm1);
709: /* dd = bd * ad; */
710: dd = pow(ad,beta);
711: gd = coef*bd*(td - t0);
712:
713: c[0].k = k-1; c[0].j = j; c[0].i = i;
714: v[0] = - hxhydhz*(dd - gd);
715: c[1].k = k; c[1].j = j-1; c[1].i = i;
716: v[1] = - hzhxdhy*(ds - gs);
717: c[2].k = k; c[2].j = j; c[2].i = i;
718: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
719: c[3].k = k; c[3].j = j; c[3].i = i+1;
720: v[3] = - hyhzdhx*(de + ge);
721: c[4].k = k+1; c[4].j = j; c[4].i = i;
722: v[4] = - hxhydhz*(du + gu);
723: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
725: /* left-hand top up corner */
726: } else {
728: td = x[k-1][j][i];
729: ad = 0.5*(t0 + td);
730: bd = pow(ad,bm1);
731: /* dd = bd * ad; */
732: dd = pow(ad,beta);
733: gd = coef*bd*(td - t0);
734:
735: c[0].k = k-1; c[0].j = j; c[0].i = i;
736: v[0] = - hxhydhz*(dd - gd);
737: c[1].k = k; c[1].j = j-1; c[1].i = i;
738: v[1] = - hzhxdhy*(ds - gs);
739: c[2].k = k; c[2].j = j; c[2].i = i;
740: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
741: c[3].k = k; c[3].j = j; c[3].i = i+1;
742: v[3] = - hyhzdhx*(de + ge);
743: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
744: }
746: } else {
748: ts = x[k][j-1][i];
749: as = 0.5*(t0 + ts);
750: bs = pow(as,bm1);
751: /* ds = bs * as; */
752: ds = pow(as,beta);
753: gs = coef*bs*(t0 - ts);
754:
755: tn = x[k][j+1][i];
756: an = 0.5*(t0 + tn);
757: bn = pow(an,bm1);
758: /* dn = bn * an; */
759: dn = pow(an,beta);
760: gn = coef*bn*(tn - t0);
762: /* left-hand down interior edge */
763: if (k == 0) {
765: tu = x[k+1][j][i];
766: au = 0.5*(t0 + tu);
767: bu = pow(au,bm1);
768: /* du = bu * au; */
769: du = pow(au,beta);
770: gu = coef*bu*(tu - t0);
772: c[0].k = k; c[0].j = j-1; c[0].i = i;
773: v[0] = - hzhxdhy*(ds - gs);
774: c[1].k = k; c[1].j = j; c[1].i = i;
775: v[1] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
776: c[2].k = k; c[2].j = j; c[2].i = i+1;
777: v[2] = - hyhzdhx*(de + ge);
778: c[3].k = k; c[3].j = j+1; c[3].i = i;
779: v[3] = - hzhxdhy*(dn + gn);
780: c[4].k = k+1; c[4].j = j; c[4].i = i;
781: v[4] = - hxhydhz*(du + gu);
782: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
783: }
785: /* left-hand up interior edge */
786: else if (k == mz-1) {
788: td = x[k-1][j][i];
789: ad = 0.5*(t0 + td);
790: bd = pow(ad,bm1);
791: /* dd = bd * ad; */
792: dd = pow(ad,beta);
793: gd = coef*bd*(t0 - td);
794:
795: c[0].k = k-1; c[0].j = j; c[0].i = i;
796: v[0] = - hxhydhz*(dd - gd);
797: c[1].k = k; c[1].j = j-1; c[1].i = i;
798: v[1] = - hzhxdhy*(ds - gs);
799: c[2].k = k; c[2].j = j; c[2].i = i;
800: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
801: c[3].k = k; c[3].j = j; c[3].i = i+1;
802: v[3] = - hyhzdhx*(de + ge);
803: c[4].k = k; c[4].j = j+1; c[4].i = i;
804: v[4] = - hzhxdhy*(dn + gn);
805: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
806: }
808: /* left-hand interior plane */
809: else {
811: td = x[k-1][j][i];
812: ad = 0.5*(t0 + td);
813: bd = pow(ad,bm1);
814: /* dd = bd * ad; */
815: dd = pow(ad,beta);
816: gd = coef*bd*(t0 - td);
817:
818: tu = x[k+1][j][i];
819: au = 0.5*(t0 + tu);
820: bu = pow(au,bm1);
821: /* du = bu * au; */
822: du = pow(au,beta);
823: gu = coef*bu*(tu - t0);
825: c[0].k = k-1; c[0].j = j; c[0].i = i;
826: v[0] = - hxhydhz*(dd - gd);
827: c[1].k = k; c[1].j = j-1; c[1].i = i;
828: v[1] = - hzhxdhy*(ds - gs);
829: c[2].k = k; c[2].j = j; c[2].i = i;
830: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
831: c[3].k = k; c[3].j = j; c[3].i = i+1;
832: v[3] = - hyhzdhx*(de + ge);
833: c[4].k = k; c[4].j = j+1; c[4].i = i;
834: v[4] = - hzhxdhy*(dn + gn);
835: c[5].k = k+1; c[5].j = j; c[5].i = i;
836: v[5] = - hxhydhz*(du + gu);
837: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
838: }
839: }
841: } else if (i == mx-1) {
843: /* right-hand plane boundary */
844: tw = x[k][j][i-1];
845: aw = 0.5*(t0 + tw);
846: bw = pow(aw,bm1);
847: /* dw = bw * aw */
848: dw = pow(aw,beta);
849: gw = coef*bw*(t0 - tw);
850:
851: te = tright;
852: ae = 0.5*(t0 + te);
853: be = pow(ae,bm1);
854: /* de = be * ae; */
855: de = pow(ae,beta);
856: ge = coef*be*(te - t0);
857:
858: /* right-hand bottom edge */
859: if (j == 0) {
861: tn = x[k][j+1][i];
862: an = 0.5*(t0 + tn);
863: bn = pow(an,bm1);
864: /* dn = bn * an; */
865: dn = pow(an,beta);
866: gn = coef*bn*(tn - t0);
867:
868: /* right-hand bottom down corner */
869: if (k == 0) {
871: tu = x[k+1][j][i];
872: au = 0.5*(t0 + tu);
873: bu = pow(au,bm1);
874: /* du = bu * au; */
875: du = pow(au,beta);
876: gu = coef*bu*(tu - t0);
877:
878: c[0].k = k; c[0].j = j; c[0].i = i-1;
879: v[0] = - hyhzdhx*(dw - gw);
880: c[1].k = k; c[1].j = j; c[1].i = i;
881: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
882: c[2].k = k; c[2].j = j+1; c[2].i = i;
883: v[2] = - hzhxdhy*(dn + gn);
884: c[3].k = k+1; c[3].j = j; c[3].i = i;
885: v[3] = - hxhydhz*(du + gu);
886: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
888: /* right-hand bottom interior edge */
889: } else if (k < mz-1) {
891: tu = x[k+1][j][i];
892: au = 0.5*(t0 + tu);
893: bu = pow(au,bm1);
894: /* du = bu * au; */
895: du = pow(au,beta);
896: gu = coef*bu*(tu - t0);
897:
898: td = x[k-1][j][i];
899: ad = 0.5*(t0 + td);
900: bd = pow(ad,bm1);
901: /* dd = bd * ad; */
902: dd = pow(ad,beta);
903: gd = coef*bd*(td - t0);
904:
905: c[0].k = k-1; c[0].j = j; c[0].i = i;
906: v[0] = - hxhydhz*(dd - gd);
907: c[1].k = k; c[1].j = j; c[1].i = i-1;
908: v[1] = - hyhzdhx*(dw - gw);
909: c[2].k = k; c[2].j = j; c[2].i = i;
910: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
911: c[3].k = k; c[3].j = j+1; c[3].i = i;
912: v[3] = - hzhxdhy*(dn + gn);
913: c[4].k = k+1; c[4].j = j; c[4].i = i;
914: v[4] = - hxhydhz*(du + gu);
915: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
917: /* right-hand bottom up corner */
918: } else {
920: td = x[k-1][j][i];
921: ad = 0.5*(t0 + td);
922: bd = pow(ad,bm1);
923: /* dd = bd * ad; */
924: dd = pow(ad,beta);
925: gd = coef*bd*(td - t0);
926:
927: c[0].k = k-1; c[0].j = j; c[0].i = i;
928: v[0] = - hxhydhz*(dd - gd);
929: c[1].k = k; c[1].j = j; c[1].i = i-1;
930: v[1] = - hyhzdhx*(dw - gw);
931: c[2].k = k; c[2].j = j; c[2].i = i;
932: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
933: c[3].k = k; c[3].j = j+1; c[3].i = i;
934: v[3] = - hzhxdhy*(dn + gn);
935: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
936: }
938: /* right-hand top edge */
939: } else if (j == my-1) {
941: ts = x[k][j-1][i];
942: as = 0.5*(t0 + ts);
943: bs = pow(as,bm1);
944: /* ds = bs * as; */
945: ds = pow(as,beta);
946: gs = coef*bs*(ts - t0);
947:
948: /* right-hand top down corner */
949: if (k == 0) {
951: tu = x[k+1][j][i];
952: au = 0.5*(t0 + tu);
953: bu = pow(au,bm1);
954: /* du = bu * au; */
955: du = pow(au,beta);
956: gu = coef*bu*(tu - t0);
957:
958: c[0].k = k; c[0].j = j-1; c[0].i = i;
959: v[0] = - hzhxdhy*(ds - gs);
960: c[1].k = k; c[1].j = j; c[1].i = i-1;
961: v[1] = - hyhzdhx*(dw - gw);
962: c[2].k = k; c[2].j = j; c[2].i = i;
963: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
964: c[3].k = k+1; c[3].j = j; c[3].i = i;
965: v[3] = - hxhydhz*(du + gu);
966: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
968: /* right-hand top interior edge */
969: } else if (k < mz-1) {
971: tu = x[k+1][j][i];
972: au = 0.5*(t0 + tu);
973: bu = pow(au,bm1);
974: /* du = bu * au; */
975: du = pow(au,beta);
976: gu = coef*bu*(tu - t0);
977:
978: td = x[k-1][j][i];
979: ad = 0.5*(t0 + td);
980: bd = pow(ad,bm1);
981: /* dd = bd * ad; */
982: dd = pow(ad,beta);
983: gd = coef*bd*(td - t0);
984:
985: c[0].k = k-1; c[0].j = j; c[0].i = i;
986: v[0] = - hxhydhz*(dd - gd);
987: c[1].k = k; c[1].j = j-1; c[1].i = i;
988: v[1] = - hzhxdhy*(ds - gs);
989: c[2].k = k; c[2].j = j; c[2].i = i-1;
990: v[2] = - hyhzdhx*(dw - gw);
991: c[3].k = k; c[3].j = j; c[3].i = i;
992: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
993: c[4].k = k+1; c[4].j = j; c[4].i = i;
994: v[4] = - hxhydhz*(du + gu);
995: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
997: /* right-hand top up corner */
998: } else {
1000: td = x[k-1][j][i];
1001: ad = 0.5*(t0 + td);
1002: bd = pow(ad,bm1);
1003: /* dd = bd * ad; */
1004: dd = pow(ad,beta);
1005: gd = coef*bd*(td - t0);
1006:
1007: c[0].k = k-1; c[0].j = j; c[0].i = i;
1008: v[0] = - hxhydhz*(dd - gd);
1009: c[1].k = k; c[1].j = j-1; c[1].i = i;
1010: v[1] = - hzhxdhy*(ds - gs);
1011: c[2].k = k; c[2].j = j; c[2].i = i-1;
1012: v[2] = - hyhzdhx*(dw - gw);
1013: c[3].k = k; c[3].j = j; c[3].i = i;
1014: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1015: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
1016: }
1018: } else {
1020: ts = x[k][j-1][i];
1021: as = 0.5*(t0 + ts);
1022: bs = pow(as,bm1);
1023: /* ds = bs * as; */
1024: ds = pow(as,beta);
1025: gs = coef*bs*(t0 - ts);
1026:
1027: tn = x[k][j+1][i];
1028: an = 0.5*(t0 + tn);
1029: bn = pow(an,bm1);
1030: /* dn = bn * an; */
1031: dn = pow(an,beta);
1032: gn = coef*bn*(tn - t0);
1034: /* right-hand down interior edge */
1035: if (k == 0) {
1037: tu = x[k+1][j][i];
1038: au = 0.5*(t0 + tu);
1039: bu = pow(au,bm1);
1040: /* du = bu * au; */
1041: du = pow(au,beta);
1042: gu = coef*bu*(tu - t0);
1044: c[0].k = k; c[0].j = j-1; c[0].i = i;
1045: v[0] = - hzhxdhy*(ds - gs);
1046: c[1].k = k; c[1].j = j; c[1].i = i-1;
1047: v[1] = - hyhzdhx*(dw - gw);
1048: c[2].k = k; c[2].j = j; c[2].i = i;
1049: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1050: c[3].k = k; c[3].j = j+1; c[3].i = i;
1051: v[3] = - hzhxdhy*(dn + gn);
1052: c[4].k = k+1; c[4].j = j; c[4].i = i;
1053: v[4] = - hxhydhz*(du + gu);
1054: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1055: }
1057: /* right-hand up interior edge */
1058: else if (k == mz-1) {
1060: td = x[k-1][j][i];
1061: ad = 0.5*(t0 + td);
1062: bd = pow(ad,bm1);
1063: /* dd = bd * ad; */
1064: dd = pow(ad,beta);
1065: gd = coef*bd*(t0 - td);
1066:
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 + gd);
1075: c[4].k = k; c[4].j = j+1; c[4].i = i;
1076: v[4] = - hzhxdhy*(dn + gn);
1077: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1078: }
1080: /* right-hand interior plane */
1081: else {
1083: td = x[k-1][j][i];
1084: ad = 0.5*(t0 + td);
1085: bd = pow(ad,bm1);
1086: /* dd = bd * ad; */
1087: dd = pow(ad,beta);
1088: gd = coef*bd*(t0 - td);
1089:
1090: tu = x[k+1][j][i];
1091: au = 0.5*(t0 + tu);
1092: bu = pow(au,bm1);
1093: /* du = bu * au; */
1094: du = pow(au,beta);
1095: gu = coef*bu*(tu - t0);
1097: c[0].k = k-1; c[0].j = j; c[0].i = i;
1098: v[0] = - hxhydhz*(dd - gd);
1099: c[1].k = k; c[1].j = j-1; c[1].i = i;
1100: v[1] = - hzhxdhy*(ds - gs);
1101: c[2].k = k; c[2].j = j; c[2].i = i-1;
1102: v[2] = - hyhzdhx*(dw - gw);
1103: c[3].k = k; c[3].j = j; c[3].i = i;
1104: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1105: c[4].k = k; c[4].j = j+1; c[4].i = i;
1106: v[4] = - hzhxdhy*(dn + gn);
1107: c[5].k = k+1; c[5].j = j; c[5].i = i;
1108: v[5] = - hxhydhz*(du + gu);
1109: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1110: }
1111: }
1113: } else if (j == 0) {
1115: tw = x[k][j][i-1];
1116: aw = 0.5*(t0 + tw);
1117: bw = pow(aw,bm1);
1118: /* dw = bw * aw */
1119: dw = pow(aw,beta);
1120: gw = coef*bw*(t0 - tw);
1122: te = x[k][j][i+1];
1123: ae = 0.5*(t0 + te);
1124: be = pow(ae,bm1);
1125: /* de = be * ae; */
1126: de = pow(ae,beta);
1127: ge = coef*be*(te - t0);
1129: tn = x[k][j+1][i];
1130: an = 0.5*(t0 + tn);
1131: bn = pow(an,bm1);
1132: /* dn = bn * an; */
1133: dn = pow(an,beta);
1134: gn = coef*bn*(tn - t0);
1137: /* bottom down interior edge */
1138: if (k == 0) {
1140: tu = x[k+1][j][i];
1141: au = 0.5*(t0 + tu);
1142: bu = pow(au,bm1);
1143: /* du = bu * au; */
1144: du = pow(au,beta);
1145: gu = coef*bu*(tu - t0);
1147: c[0].k = k; c[0].j = j; c[0].i = i-1;
1148: v[0] = - hyhzdhx*(dw - gw);
1149: c[1].k = k; c[1].j = j; c[1].i = i;
1150: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1151: c[2].k = k; c[2].j = j; c[2].i = i+1;
1152: v[2] = - hyhzdhx*(de + ge);
1153: c[3].k = k; c[3].j = j+1; c[3].i = i;
1154: v[3] = - hzhxdhy*(dn + gn);
1155: c[4].k = k+1; c[4].j = j; c[4].i = i;
1156: v[4] = - hxhydhz*(du + gu);
1157: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1158: }
1160: /* bottom up interior edge */
1161: else if (k == mz-1) {
1163: td = x[k-1][j][i];
1164: ad = 0.5*(t0 + td);
1165: bd = pow(ad,bm1);
1166: /* dd = bd * ad; */
1167: dd = pow(ad,beta);
1168: gd = coef*bd*(td - t0);
1170: c[0].k = k-1; c[0].j = j; c[0].i = i;
1171: v[0] = - hxhydhz*(dd - gd);
1172: c[1].k = k; c[1].j = j; c[1].i = i-1;
1173: v[1] = - hyhzdhx*(dw - gw);
1174: c[2].k = k; c[2].j = j; c[2].i = i;
1175: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1176: c[3].k = k; c[3].j = j; c[3].i = i+1;
1177: v[3] = - hyhzdhx*(de + ge);
1178: c[4].k = k; c[4].j = j+1; c[4].i = i;
1179: v[4] = - hzhxdhy*(dn + gn);
1180: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1181: }
1183: /* bottom interior plane */
1184: else {
1186: tu = x[k+1][j][i];
1187: au = 0.5*(t0 + tu);
1188: bu = pow(au,bm1);
1189: /* du = bu * au; */
1190: du = pow(au,beta);
1191: gu = coef*bu*(tu - t0);
1193: td = x[k-1][j][i];
1194: ad = 0.5*(t0 + td);
1195: bd = pow(ad,bm1);
1196: /* dd = bd * ad; */
1197: dd = pow(ad,beta);
1198: gd = coef*bd*(td - t0);
1200: c[0].k = k-1; c[0].j = j; c[0].i = i;
1201: v[0] = - hxhydhz*(dd - gd);
1202: c[1].k = k; c[1].j = j; c[1].i = i-1;
1203: v[1] = - hyhzdhx*(dw - gw);
1204: c[2].k = k; c[2].j = j; c[2].i = i;
1205: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1206: c[3].k = k; c[3].j = j; c[3].i = i+1;
1207: v[3] = - hyhzdhx*(de + ge);
1208: c[4].k = k; c[4].j = j+1; c[4].i = i;
1209: v[4] = - hzhxdhy*(dn + gn);
1210: c[5].k = k+1; c[5].j = j; c[5].i = i;
1211: v[5] = - hxhydhz*(du + gu);
1212: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1213: }
1215: } else if (j == my-1) {
1217: tw = x[k][j][i-1];
1218: aw = 0.5*(t0 + tw);
1219: bw = pow(aw,bm1);
1220: /* dw = bw * aw */
1221: dw = pow(aw,beta);
1222: gw = coef*bw*(t0 - tw);
1224: te = x[k][j][i+1];
1225: ae = 0.5*(t0 + te);
1226: be = pow(ae,bm1);
1227: /* de = be * ae; */
1228: de = pow(ae,beta);
1229: ge = coef*be*(te - t0);
1231: ts = x[k][j-1][i];
1232: as = 0.5*(t0 + ts);
1233: bs = pow(as,bm1);
1234: /* ds = bs * as; */
1235: ds = pow(as,beta);
1236: gs = coef*bs*(t0 - ts);
1237:
1238: /* top down interior edge */
1239: if (k == 0) {
1241: tu = x[k+1][j][i];
1242: au = 0.5*(t0 + tu);
1243: bu = pow(au,bm1);
1244: /* du = bu * au; */
1245: du = pow(au,beta);
1246: gu = coef*bu*(tu - t0);
1248: c[0].k = k; c[0].j = j-1; c[0].i = i;
1249: v[0] = - hzhxdhy*(ds - gs);
1250: c[1].k = k; c[1].j = j; c[1].i = i-1;
1251: v[1] = - hyhzdhx*(dw - gw);
1252: c[2].k = k; c[2].j = j; c[2].i = i;
1253: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1254: c[3].k = k; c[3].j = j; c[3].i = i+1;
1255: v[3] = - hyhzdhx*(de + ge);
1256: c[4].k = k+1; c[4].j = j; c[4].i = i;
1257: v[4] = - hxhydhz*(du + gu);
1258: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1259: }
1261: /* top up interior edge */
1262: else if (k == mz-1) {
1264: td = x[k-1][j][i];
1265: ad = 0.5*(t0 + td);
1266: bd = pow(ad,bm1);
1267: /* dd = bd * ad; */
1268: dd = pow(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 + gd);
1279: c[4].k = k; c[4].j = j; c[4].i = i+1;
1280: v[4] = - hyhzdhx*(de + ge);
1281: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1282: }
1284: /* top interior plane */
1285: else {
1287: tu = x[k+1][j][i];
1288: au = 0.5*(t0 + tu);
1289: bu = pow(au,bm1);
1290: /* du = bu * au; */
1291: du = pow(au,beta);
1292: gu = coef*bu*(tu - t0);
1294: td = x[k-1][j][i];
1295: ad = 0.5*(t0 + td);
1296: bd = pow(ad,bm1);
1297: /* dd = bd * ad; */
1298: dd = pow(ad,beta);
1299: gd = coef*bd*(td - t0);
1301: c[0].k = k-1; c[0].j = j; c[0].i = i;
1302: v[0] = - hxhydhz*(dd - gd);
1303: c[1].k = k; c[1].j = j-1; c[1].i = i;
1304: v[1] = - hzhxdhy*(ds - gs);
1305: c[2].k = k; c[2].j = j; c[2].i = i-1;
1306: v[2] = - hyhzdhx*(dw - gw);
1307: c[3].k = k; c[3].j = j; c[3].i = i;
1308: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1309: c[4].k = k; c[4].j = j; c[4].i = i+1;
1310: v[4] = - hyhzdhx*(de + ge);
1311: c[5].k = k+1; c[5].j = j; c[5].i = i;
1312: v[5] = - hxhydhz*(du + gu);
1313: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1314: }
1316: } else if (k == 0) {
1318: /* down interior plane */
1320: tw = x[k][j][i-1];
1321: aw = 0.5*(t0 + tw);
1322: bw = pow(aw,bm1);
1323: /* dw = bw * aw */
1324: dw = pow(aw,beta);
1325: gw = coef*bw*(t0 - tw);
1327: te = x[k][j][i+1];
1328: ae = 0.5*(t0 + te);
1329: be = pow(ae,bm1);
1330: /* de = be * ae; */
1331: de = pow(ae,beta);
1332: ge = coef*be*(te - t0);
1334: ts = x[k][j-1][i];
1335: as = 0.5*(t0 + ts);
1336: bs = pow(as,bm1);
1337: /* ds = bs * as; */
1338: ds = pow(as,beta);
1339: gs = coef*bs*(t0 - ts);
1340:
1341: tn = x[k][j+1][i];
1342: an = 0.5*(t0 + tn);
1343: bn = pow(an,bm1);
1344: /* dn = bn * an; */
1345: dn = pow(an,beta);
1346: gn = coef*bn*(tn - t0);
1347:
1348: tu = x[k+1][j][i];
1349: au = 0.5*(t0 + tu);
1350: bu = pow(au,bm1);
1351: /* du = bu * au; */
1352: du = pow(au,beta);
1353: gu = coef*bu*(tu - t0);
1355: c[0].k = k; c[0].j = j-1; c[0].i = i;
1356: v[0] = - hzhxdhy*(ds - gs);
1357: c[1].k = k; c[1].j = j; c[1].i = i-1;
1358: v[1] = - hyhzdhx*(dw - gw);
1359: c[2].k = k; c[2].j = j; c[2].i = i;
1360: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1361: c[3].k = k; c[3].j = j; c[3].i = i+1;
1362: v[3] = - hyhzdhx*(de + ge);
1363: c[4].k = k; c[4].j = j+1; c[4].i = i;
1364: v[4] = - hzhxdhy*(dn + gn);
1365: c[5].k = k+1; c[5].j = j; c[5].i = i;
1366: v[5] = - hxhydhz*(du + gu);
1367: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1368: }
1369:
1370: else if (k == mz-1) {
1372: /* up interior plane */
1374: tw = x[k][j][i-1];
1375: aw = 0.5*(t0 + tw);
1376: bw = pow(aw,bm1);
1377: /* dw = bw * aw */
1378: dw = pow(aw,beta);
1379: gw = coef*bw*(t0 - tw);
1381: te = x[k][j][i+1];
1382: ae = 0.5*(t0 + te);
1383: be = pow(ae,bm1);
1384: /* de = be * ae; */
1385: de = pow(ae,beta);
1386: ge = coef*be*(te - t0);
1388: ts = x[k][j-1][i];
1389: as = 0.5*(t0 + ts);
1390: bs = pow(as,bm1);
1391: /* ds = bs * as; */
1392: ds = pow(as,beta);
1393: gs = coef*bs*(t0 - ts);
1394:
1395: tn = x[k][j+1][i];
1396: an = 0.5*(t0 + tn);
1397: bn = pow(an,bm1);
1398: /* dn = bn * an; */
1399: dn = pow(an,beta);
1400: gn = coef*bn*(tn - t0);
1401:
1402: td = x[k-1][j][i];
1403: ad = 0.5*(t0 + td);
1404: bd = pow(ad,bm1);
1405: /* dd = bd * ad; */
1406: dd = pow(ad,beta);
1407: gd = coef*bd*(t0 - td);
1408:
1409: c[0].k = k-1; c[0].j = j; c[0].i = i;
1410: v[0] = - hxhydhz*(dd - gd);
1411: c[1].k = k; c[1].j = j-1; c[1].i = i;
1412: v[1] = - hzhxdhy*(ds - gs);
1413: c[2].k = k; c[2].j = j; c[2].i = i-1;
1414: v[2] = - hyhzdhx*(dw - gw);
1415: c[3].k = k; c[3].j = j; c[3].i = i;
1416: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1417: c[4].k = k; c[4].j = j; c[4].i = i+1;
1418: v[4] = - hyhzdhx*(de + ge);
1419: c[5].k = k; c[5].j = j+1; c[5].i = i;
1420: v[5] = - hzhxdhy*(dn + gn);
1421: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1422: }
1423: }
1424: }
1425: }
1426: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
1427: DAVecRestoreArray((DA)dmmg->dm,localX,&x);
1428: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
1429: DARestoreLocalVector((DA)dmmg->dm,&localX);
1431: PetscLogFlops((41 + 8*POWFLOP)*xm*ym);
1432: return(0);
1433: }