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);
101: DMMGSetFromOptions(dmmg);
103: /*
104: PreLoadBegin() means that the following section of code is run twice. The first time
105: through the flag PreLoading is on this the nonlinear solver is only run for a single step.
106: The second time through (the actually timed code) the maximum iterations is set to 10
107: Preload of the executable is done to eliminate from the timing the time spent bring the
108: executable into memory from disk (paging in).
109: */
110: PreLoadBegin(PETSC_TRUE,"Solve");
111: DMMGSetInitialGuess(dmmg,FormInitialGuess);
112: DMMGSolve(dmmg);
113: PreLoadEnd();
114: snes = DMMGGetSNES(dmmg);
115: SNESGetIterationNumber(snes,&its);
116: SNESGetLinearSolveIterations(snes,&lits);
117: litspit = ((PetscReal)lits)/((PetscReal)its);
118: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);
119: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
120: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / Newton = %e\n",litspit);
122: DMMGDestroy(dmmg);
123: PetscFinalize();
125: return 0;
126: }
127: /* -------------------- Form initial approximation ----------------- */
130: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
131: {
132: AppCtx *user = (AppCtx*)dmmg->user;
133: PetscInt i,j,k,xs,ys,xm,ym,zs,zm;
135: PetscReal tleft = user->tleft;
136: PetscScalar ***x;
140: /* Get ghost points */
141: DAGetCorners((DA)dmmg->dm,&xs,&ys,&zs,&xm,&ym,&zm);
142: DAVecGetArray((DA)dmmg->dm,X,&x);
144: /* Compute initial guess */
145: for (k=zs; k<zs+zm; k++) {
146: for (j=ys; j<ys+ym; j++) {
147: for (i=xs; i<xs+xm; i++) {
148: x[k][j][i] = tleft;
149: }
150: }
151: }
152: DAVecRestoreArray((DA)dmmg->dm,X,&x);
153: return(0);
154: }
155: /* -------------------- Evaluate Function F(x) --------------------- */
158: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ptr)
159: {
160: DMMG dmmg = (DMMG)ptr;
161: AppCtx *user = (AppCtx*)dmmg->user;
163: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
164: PetscScalar zero = 0.0,one = 1.0;
165: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
166: 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;
167: PetscScalar tleft,tright,beta,td,ad,dd,fd,tu,au,du,fu;
168: PetscScalar ***x,***f;
169: Vec localX;
172: DAGetLocalVector((DA)dmmg->dm,&localX);
173: DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,&mz,0,0,0,0,0,0,0);
174: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
175: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
176: tleft = user->tleft; tright = user->tright;
177: beta = user->beta;
178:
179: /* Get ghost points */
180: DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
181: DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
182: DAGetCorners((DA)dmmg->dm,&xs,&ys,&zs,&xm,&ym,&zm);
183: DAVecGetArray((DA)dmmg->dm,localX,&x);
184: DAVecGetArray((DA)dmmg->dm,F,&f);
186: /* Evaluate function */
187: for (k=zs; k<zs+zm; k++) {
188: for (j=ys; j<ys+ym; j++) {
189: for (i=xs; i<xs+xm; i++) {
190: t0 = x[k][j][i];
192: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
194: /* general interior volume */
196: tw = x[k][j][i-1];
197: aw = 0.5*(t0 + tw);
198: dw = pow(aw,beta);
199: fw = dw*(t0 - tw);
201: te = x[k][j][i+1];
202: ae = 0.5*(t0 + te);
203: de = pow(ae,beta);
204: fe = de*(te - t0);
206: ts = x[k][j-1][i];
207: as = 0.5*(t0 + ts);
208: ds = pow(as,beta);
209: fs = ds*(t0 - ts);
210:
211: tn = x[k][j+1][i];
212: an = 0.5*(t0 + tn);
213: dn = pow(an,beta);
214: fn = dn*(tn - t0);
216: td = x[k-1][j][i];
217: ad = 0.5*(t0 + td);
218: dd = pow(ad,beta);
219: fd = dd*(t0 - td);
221: tu = x[k+1][j][i];
222: au = 0.5*(t0 + tu);
223: du = pow(au,beta);
224: fu = du*(tu - t0);
226: } else if (i == 0) {
228: /* left-hand (west) boundary */
229: tw = tleft;
230: aw = 0.5*(t0 + tw);
231: dw = pow(aw,beta);
232: fw = dw*(t0 - tw);
234: te = x[k][j][i+1];
235: ae = 0.5*(t0 + te);
236: de = pow(ae,beta);
237: fe = de*(te - t0);
239: if (j > 0) {
240: ts = x[k][j-1][i];
241: as = 0.5*(t0 + ts);
242: ds = pow(as,beta);
243: fs = ds*(t0 - ts);
244: } else {
245: fs = zero;
246: }
248: if (j < my-1) {
249: tn = x[k][j+1][i];
250: an = 0.5*(t0 + tn);
251: dn = pow(an,beta);
252: fn = dn*(tn - t0);
253: } else {
254: fn = zero;
255: }
257: if (k > 0) {
258: td = x[k-1][j][i];
259: ad = 0.5*(t0 + td);
260: dd = pow(ad,beta);
261: fd = dd*(t0 - td);
262: } else {
263: fd = zero;
264: }
266: if (k < mz-1) {
267: tu = x[k+1][j][i];
268: au = 0.5*(t0 + tu);
269: du = pow(au,beta);
270: fu = du*(tu - t0);
271: } else {
272: fu = zero;
273: }
275: } else if (i == mx-1) {
277: /* right-hand (east) boundary */
278: tw = x[k][j][i-1];
279: aw = 0.5*(t0 + tw);
280: dw = pow(aw,beta);
281: fw = dw*(t0 - tw);
282:
283: te = tright;
284: ae = 0.5*(t0 + te);
285: de = pow(ae,beta);
286: fe = de*(te - t0);
287:
288: if (j > 0) {
289: ts = x[k][j-1][i];
290: as = 0.5*(t0 + ts);
291: ds = pow(as,beta);
292: fs = ds*(t0 - ts);
293: } else {
294: fs = zero;
295: }
296:
297: if (j < my-1) {
298: tn = x[k][j+1][i];
299: an = 0.5*(t0 + tn);
300: dn = pow(an,beta);
301: fn = dn*(tn - t0);
302: } else {
303: fn = zero;
304: }
306: if (k > 0) {
307: td = x[k-1][j][i];
308: ad = 0.5*(t0 + td);
309: dd = pow(ad,beta);
310: fd = dd*(t0 - td);
311: } else {
312: fd = zero;
313: }
315: if (k < mz-1) {
316: tu = x[k+1][j][i];
317: au = 0.5*(t0 + tu);
318: du = pow(au,beta);
319: fu = du*(tu - t0);
320: } else {
321: fu = zero;
322: }
324: } else if (j == 0) {
326: /* bottom (south) boundary, and i <> 0 or mx-1 */
327: tw = x[k][j][i-1];
328: aw = 0.5*(t0 + tw);
329: dw = pow(aw,beta);
330: fw = dw*(t0 - tw);
332: te = x[k][j][i+1];
333: ae = 0.5*(t0 + te);
334: de = pow(ae,beta);
335: fe = de*(te - t0);
337: fs = zero;
339: tn = x[k][j+1][i];
340: an = 0.5*(t0 + tn);
341: dn = pow(an,beta);
342: fn = dn*(tn - t0);
344: if (k > 0) {
345: td = x[k-1][j][i];
346: ad = 0.5*(t0 + td);
347: dd = pow(ad,beta);
348: fd = dd*(t0 - td);
349: } else {
350: fd = zero;
351: }
353: if (k < mz-1) {
354: tu = x[k+1][j][i];
355: au = 0.5*(t0 + tu);
356: du = pow(au,beta);
357: fu = du*(tu - t0);
358: } else {
359: fu = zero;
360: }
362: } else if (j == my-1) {
364: /* top (north) boundary, and i <> 0 or mx-1 */
365: tw = x[k][j][i-1];
366: aw = 0.5*(t0 + tw);
367: dw = pow(aw,beta);
368: fw = dw*(t0 - tw);
370: te = x[k][j][i+1];
371: ae = 0.5*(t0 + te);
372: de = pow(ae,beta);
373: fe = de*(te - t0);
375: ts = x[k][j-1][i];
376: as = 0.5*(t0 + ts);
377: ds = pow(as,beta);
378: fs = ds*(t0 - ts);
380: fn = zero;
382: if (k > 0) {
383: td = x[k-1][j][i];
384: ad = 0.5*(t0 + td);
385: dd = pow(ad,beta);
386: fd = dd*(t0 - td);
387: } else {
388: fd = zero;
389: }
391: if (k < mz-1) {
392: tu = x[k+1][j][i];
393: au = 0.5*(t0 + tu);
394: du = pow(au,beta);
395: fu = du*(tu - t0);
396: } else {
397: fu = zero;
398: }
400: } else if (k == 0) {
402: /* down boundary (interior only) */
403: tw = x[k][j][i-1];
404: aw = 0.5*(t0 + tw);
405: dw = pow(aw,beta);
406: fw = dw*(t0 - tw);
408: te = x[k][j][i+1];
409: ae = 0.5*(t0 + te);
410: de = pow(ae,beta);
411: fe = de*(te - t0);
413: ts = x[k][j-1][i];
414: as = 0.5*(t0 + ts);
415: ds = pow(as,beta);
416: fs = ds*(t0 - ts);
418: tn = x[k][j+1][i];
419: an = 0.5*(t0 + tn);
420: dn = pow(an,beta);
421: fn = dn*(tn - t0);
423: fd = zero;
425: tu = x[k+1][j][i];
426: au = 0.5*(t0 + tu);
427: du = pow(au,beta);
428: fu = du*(tu - t0);
429:
430: } else if (k == mz-1) {
432: /* up boundary (interior only) */
433: tw = x[k][j][i-1];
434: aw = 0.5*(t0 + tw);
435: dw = pow(aw,beta);
436: fw = dw*(t0 - tw);
438: te = x[k][j][i+1];
439: ae = 0.5*(t0 + te);
440: de = pow(ae,beta);
441: fe = de*(te - t0);
443: ts = x[k][j-1][i];
444: as = 0.5*(t0 + ts);
445: ds = pow(as,beta);
446: fs = ds*(t0 - ts);
448: tn = x[k][j+1][i];
449: an = 0.5*(t0 + tn);
450: dn = pow(an,beta);
451: fn = dn*(tn - t0);
453: td = x[k-1][j][i];
454: ad = 0.5*(t0 + td);
455: dd = pow(ad,beta);
456: fd = dd*(t0 - td);
458: fu = zero;
459: }
461: f[k][j][i] = - hyhzdhx*(fe-fw) - hzhxdhy*(fn-fs) - hxhydhz*(fu-fd);
462: }
463: }
464: }
465: DAVecRestoreArray((DA)dmmg->dm,localX,&x);
466: DAVecRestoreArray((DA)dmmg->dm,F,&f);
467: DARestoreLocalVector((DA)dmmg->dm,&localX);
468: PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
469: return(0);
470: }
471: /* -------------------- Evaluate Jacobian F(x) --------------------- */
474: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ptr)
475: {
476: DMMG dmmg = (DMMG)ptr;
477: AppCtx *user = (AppCtx*)dmmg->user;
479: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
480: PetscScalar one = 1.0;
481: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
482: PetscScalar t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw;
483: PetscScalar tleft,tright,beta,td,ad,dd,tu,au,du,v[7],bm1,coef;
484: PetscScalar ***x,bn,bs,be,bw,bu,bd,gn,gs,ge,gw,gu,gd;
485: Vec localX;
486: MatStencil c[7],row;
487: Mat jac = *B;
490: DAGetLocalVector((DA)dmmg->dm,&localX);
491: DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,&mz,0,0,0,0,0,0,0);
492: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
493: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
494: tleft = user->tleft; tright = user->tright;
495: beta = user->beta; bm1 = user->bm1; coef = user->coef;
497: /* Get ghost points */
498: DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
499: DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
500: DAGetCorners((DA)dmmg->dm,&xs,&ys,&zs,&xm,&ym,&zm);
501: DAVecGetArray((DA)dmmg->dm,localX,&x);
503: /* Evaluate Jacobian of function */
504: for (k=zs; k<zs+zm; k++) {
505: for (j=ys; j<ys+ym; j++) {
506: for (i=xs; i<xs+xm; i++) {
507: t0 = x[k][j][i];
508: row.k = k; row.j = j; row.i = i;
509: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
511: /* general interior volume */
513: tw = x[k][j][i-1];
514: aw = 0.5*(t0 + tw);
515: bw = pow(aw,bm1);
516: /* dw = bw * aw */
517: dw = pow(aw,beta);
518: gw = coef*bw*(t0 - tw);
520: te = x[k][j][i+1];
521: ae = 0.5*(t0 + te);
522: be = pow(ae,bm1);
523: /* de = be * ae; */
524: de = pow(ae,beta);
525: ge = coef*be*(te - t0);
527: ts = x[k][j-1][i];
528: as = 0.5*(t0 + ts);
529: bs = pow(as,bm1);
530: /* ds = bs * as; */
531: ds = pow(as,beta);
532: gs = coef*bs*(t0 - ts);
533:
534: tn = x[k][j+1][i];
535: an = 0.5*(t0 + tn);
536: bn = pow(an,bm1);
537: /* dn = bn * an; */
538: dn = pow(an,beta);
539: gn = coef*bn*(tn - t0);
541: td = x[k-1][j][i];
542: ad = 0.5*(t0 + td);
543: bd = pow(ad,bm1);
544: /* dd = bd * ad; */
545: dd = pow(ad,beta);
546: gd = coef*bd*(t0 - td);
547:
548: tu = x[k+1][j][i];
549: au = 0.5*(t0 + tu);
550: bu = pow(au,bm1);
551: /* du = bu * au; */
552: du = pow(au,beta);
553: gu = coef*bu*(tu - t0);
555: c[0].k = k-1; c[0].j = j; c[0].i = i; v[0] = - hxhydhz*(dd - gd);
556: c[1].k = k; c[1].j = j-1; c[1].i = i;
557: v[1] = - hzhxdhy*(ds - gs);
558: c[2].k = k; c[2].j = j; c[2].i = i-1;
559: v[2] = - hyhzdhx*(dw - gw);
560: c[3].k = k; c[3].j = j; c[3].i = i;
561: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
562: c[4].k = k; c[4].j = j; c[4].i = i+1;
563: v[4] = - hyhzdhx*(de + ge);
564: c[5].k = k; c[5].j = j+1; c[5].i = i;
565: v[5] = - hzhxdhy*(dn + gn);
566: c[6].k = k+1; c[6].j = j; c[6].i = i;
567: v[6] = - hxhydhz*(du + gu);
568: MatSetValuesStencil(jac,1,&row,7,c,v,INSERT_VALUES);
570: } else if (i == 0) {
572: /* left-hand plane boundary */
573: tw = tleft;
574: aw = 0.5*(t0 + tw);
575: bw = pow(aw,bm1);
576: /* dw = bw * aw */
577: dw = pow(aw,beta);
578: gw = coef*bw*(t0 - tw);
579:
580: te = x[k][j][i+1];
581: ae = 0.5*(t0 + te);
582: be = pow(ae,bm1);
583: /* de = be * ae; */
584: de = pow(ae,beta);
585: ge = coef*be*(te - t0);
586:
587: /* left-hand bottom edge */
588: if (j == 0) {
590: tn = x[k][j+1][i];
591: an = 0.5*(t0 + tn);
592: bn = pow(an,bm1);
593: /* dn = bn * an; */
594: dn = pow(an,beta);
595: gn = coef*bn*(tn - t0);
596:
597: /* left-hand bottom down corner */
598: if (k == 0) {
600: tu = x[k+1][j][i];
601: au = 0.5*(t0 + tu);
602: bu = pow(au,bm1);
603: /* du = bu * au; */
604: du = pow(au,beta);
605: gu = coef*bu*(tu - t0);
606:
607: c[0].k = k; c[0].j = j; c[0].i = i;
608: v[0] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
609: c[1].k = k; c[1].j = j; c[1].i = i+1;
610: v[1] = - hyhzdhx*(de + ge);
611: c[2].k = k; c[2].j = j+1; c[2].i = i;
612: v[2] = - hzhxdhy*(dn + gn);
613: c[3].k = k+1; c[3].j = j; c[3].i = i;
614: v[3] = - hxhydhz*(du + gu);
615: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
617: /* left-hand bottom interior edge */
618: } else if (k < mz-1) {
620: tu = x[k+1][j][i];
621: au = 0.5*(t0 + tu);
622: bu = pow(au,bm1);
623: /* du = bu * au; */
624: du = pow(au,beta);
625: gu = coef*bu*(tu - t0);
626:
627: td = x[k-1][j][i];
628: ad = 0.5*(t0 + td);
629: bd = pow(ad,bm1);
630: /* dd = bd * ad; */
631: dd = pow(ad,beta);
632: gd = coef*bd*(td - t0);
633:
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 + du + gd - gu);
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: c[4].k = k+1; c[4].j = j; c[4].i = i;
643: v[4] = - hxhydhz*(du + gu);
644: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
646: /* left-hand bottom up corner */
647: } else {
649: td = x[k-1][j][i];
650: ad = 0.5*(t0 + td);
651: bd = pow(ad,bm1);
652: /* dd = bd * ad; */
653: dd = pow(ad,beta);
654: gd = coef*bd*(td - t0);
655:
656: c[0].k = k-1; c[0].j = j; c[0].i = i;
657: v[0] = - hxhydhz*(dd - gd);
658: c[1].k = k; c[1].j = j; c[1].i = i;
659: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
660: c[2].k = k; c[2].j = j; c[2].i = i+1;
661: v[2] = - hyhzdhx*(de + ge);
662: c[3].k = k; c[3].j = j+1; c[3].i = i;
663: v[3] = - hzhxdhy*(dn + gn);
664: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
665: }
667: /* left-hand top edge */
668: } else if (j == my-1) {
670: ts = x[k][j-1][i];
671: as = 0.5*(t0 + ts);
672: bs = pow(as,bm1);
673: /* ds = bs * as; */
674: ds = pow(as,beta);
675: gs = coef*bs*(ts - t0);
676:
677: /* left-hand top down corner */
678: if (k == 0) {
680: tu = x[k+1][j][i];
681: au = 0.5*(t0 + tu);
682: bu = pow(au,bm1);
683: /* du = bu * au; */
684: du = pow(au,beta);
685: gu = coef*bu*(tu - t0);
686:
687: c[0].k = k; c[0].j = j-1; c[0].i = i;
688: v[0] = - hzhxdhy*(ds - gs);
689: c[1].k = k; c[1].j = j; c[1].i = i;
690: v[1] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
691: c[2].k = k; c[2].j = j; c[2].i = i+1;
692: v[2] = - hyhzdhx*(de + ge);
693: c[3].k = k+1; c[3].j = j; c[3].i = i;
694: v[3] = - hxhydhz*(du + gu);
695: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
697: /* left-hand top interior edge */
698: } else if (k < mz-1) {
700: tu = x[k+1][j][i];
701: au = 0.5*(t0 + tu);
702: bu = pow(au,bm1);
703: /* du = bu * au; */
704: du = pow(au,beta);
705: gu = coef*bu*(tu - t0);
706:
707: td = x[k-1][j][i];
708: ad = 0.5*(t0 + td);
709: bd = pow(ad,bm1);
710: /* dd = bd * ad; */
711: dd = pow(ad,beta);
712: gd = coef*bd*(td - t0);
713:
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 + du + gd - gu);
720: c[3].k = k; c[3].j = j; c[3].i = i+1;
721: v[3] = - hyhzdhx*(de + ge);
722: c[4].k = k+1; c[4].j = j; c[4].i = i;
723: v[4] = - hxhydhz*(du + gu);
724: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
726: /* left-hand top up corner */
727: } else {
729: td = x[k-1][j][i];
730: ad = 0.5*(t0 + td);
731: bd = pow(ad,bm1);
732: /* dd = bd * ad; */
733: dd = pow(ad,beta);
734: gd = coef*bd*(td - t0);
735:
736: c[0].k = k-1; c[0].j = j; c[0].i = i;
737: v[0] = - hxhydhz*(dd - gd);
738: c[1].k = k; c[1].j = j-1; c[1].i = i;
739: v[1] = - hzhxdhy*(ds - gs);
740: c[2].k = k; c[2].j = j; c[2].i = i;
741: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
742: c[3].k = k; c[3].j = j; c[3].i = i+1;
743: v[3] = - hyhzdhx*(de + ge);
744: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
745: }
747: } else {
749: ts = x[k][j-1][i];
750: as = 0.5*(t0 + ts);
751: bs = pow(as,bm1);
752: /* ds = bs * as; */
753: ds = pow(as,beta);
754: gs = coef*bs*(t0 - ts);
755:
756: tn = x[k][j+1][i];
757: an = 0.5*(t0 + tn);
758: bn = pow(an,bm1);
759: /* dn = bn * an; */
760: dn = pow(an,beta);
761: gn = coef*bn*(tn - t0);
763: /* left-hand down interior edge */
764: if (k == 0) {
766: tu = x[k+1][j][i];
767: au = 0.5*(t0 + tu);
768: bu = pow(au,bm1);
769: /* du = bu * au; */
770: du = pow(au,beta);
771: gu = coef*bu*(tu - t0);
773: c[0].k = k; c[0].j = j-1; c[0].i = i;
774: v[0] = - hzhxdhy*(ds - gs);
775: c[1].k = k; c[1].j = j; c[1].i = i;
776: v[1] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
777: c[2].k = k; c[2].j = j; c[2].i = i+1;
778: v[2] = - hyhzdhx*(de + ge);
779: c[3].k = k; c[3].j = j+1; c[3].i = i;
780: v[3] = - hzhxdhy*(dn + gn);
781: c[4].k = k+1; c[4].j = j; c[4].i = i;
782: v[4] = - hxhydhz*(du + gu);
783: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
784: }
786: /* left-hand up interior edge */
787: else if (k == mz-1) {
789: td = x[k-1][j][i];
790: ad = 0.5*(t0 + td);
791: bd = pow(ad,bm1);
792: /* dd = bd * ad; */
793: dd = pow(ad,beta);
794: gd = coef*bd*(t0 - td);
795:
796: c[0].k = k-1; c[0].j = j; c[0].i = i;
797: v[0] = - hxhydhz*(dd - gd);
798: c[1].k = k; c[1].j = j-1; c[1].i = i;
799: v[1] = - hzhxdhy*(ds - gs);
800: c[2].k = k; c[2].j = j; c[2].i = i;
801: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
802: c[3].k = k; c[3].j = j; c[3].i = i+1;
803: v[3] = - hyhzdhx*(de + ge);
804: c[4].k = k; c[4].j = j+1; c[4].i = i;
805: v[4] = - hzhxdhy*(dn + gn);
806: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
807: }
809: /* left-hand interior plane */
810: else {
812: td = x[k-1][j][i];
813: ad = 0.5*(t0 + td);
814: bd = pow(ad,bm1);
815: /* dd = bd * ad; */
816: dd = pow(ad,beta);
817: gd = coef*bd*(t0 - td);
818:
819: tu = x[k+1][j][i];
820: au = 0.5*(t0 + tu);
821: bu = pow(au,bm1);
822: /* du = bu * au; */
823: du = pow(au,beta);
824: gu = coef*bu*(tu - t0);
826: c[0].k = k-1; c[0].j = j; c[0].i = i;
827: v[0] = - hxhydhz*(dd - gd);
828: c[1].k = k; c[1].j = j-1; c[1].i = i;
829: v[1] = - hzhxdhy*(ds - gs);
830: c[2].k = k; c[2].j = j; c[2].i = i;
831: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
832: c[3].k = k; c[3].j = j; c[3].i = i+1;
833: v[3] = - hyhzdhx*(de + ge);
834: c[4].k = k; c[4].j = j+1; c[4].i = i;
835: v[4] = - hzhxdhy*(dn + gn);
836: c[5].k = k+1; c[5].j = j; c[5].i = i;
837: v[5] = - hxhydhz*(du + gu);
838: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
839: }
840: }
842: } else if (i == mx-1) {
844: /* right-hand plane boundary */
845: tw = x[k][j][i-1];
846: aw = 0.5*(t0 + tw);
847: bw = pow(aw,bm1);
848: /* dw = bw * aw */
849: dw = pow(aw,beta);
850: gw = coef*bw*(t0 - tw);
851:
852: te = tright;
853: ae = 0.5*(t0 + te);
854: be = pow(ae,bm1);
855: /* de = be * ae; */
856: de = pow(ae,beta);
857: ge = coef*be*(te - t0);
858:
859: /* right-hand bottom edge */
860: if (j == 0) {
862: tn = x[k][j+1][i];
863: an = 0.5*(t0 + tn);
864: bn = pow(an,bm1);
865: /* dn = bn * an; */
866: dn = pow(an,beta);
867: gn = coef*bn*(tn - t0);
868:
869: /* right-hand bottom down corner */
870: if (k == 0) {
872: tu = x[k+1][j][i];
873: au = 0.5*(t0 + tu);
874: bu = pow(au,bm1);
875: /* du = bu * au; */
876: du = pow(au,beta);
877: gu = coef*bu*(tu - t0);
878:
879: c[0].k = k; c[0].j = j; c[0].i = i-1;
880: v[0] = - hyhzdhx*(dw - gw);
881: c[1].k = k; c[1].j = j; c[1].i = i;
882: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
883: c[2].k = k; c[2].j = j+1; c[2].i = i;
884: v[2] = - hzhxdhy*(dn + gn);
885: c[3].k = k+1; c[3].j = j; c[3].i = i;
886: v[3] = - hxhydhz*(du + gu);
887: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
889: /* right-hand bottom interior edge */
890: } else if (k < mz-1) {
892: tu = x[k+1][j][i];
893: au = 0.5*(t0 + tu);
894: bu = pow(au,bm1);
895: /* du = bu * au; */
896: du = pow(au,beta);
897: gu = coef*bu*(tu - t0);
898:
899: td = x[k-1][j][i];
900: ad = 0.5*(t0 + td);
901: bd = pow(ad,bm1);
902: /* dd = bd * ad; */
903: dd = pow(ad,beta);
904: gd = coef*bd*(td - t0);
905:
906: c[0].k = k-1; c[0].j = j; c[0].i = i;
907: v[0] = - hxhydhz*(dd - gd);
908: c[1].k = k; c[1].j = j; c[1].i = i-1;
909: v[1] = - hyhzdhx*(dw - gw);
910: c[2].k = k; c[2].j = j; c[2].i = i;
911: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
912: c[3].k = k; c[3].j = j+1; c[3].i = i;
913: v[3] = - hzhxdhy*(dn + gn);
914: c[4].k = k+1; c[4].j = j; c[4].i = i;
915: v[4] = - hxhydhz*(du + gu);
916: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
918: /* right-hand bottom up corner */
919: } else {
921: td = x[k-1][j][i];
922: ad = 0.5*(t0 + td);
923: bd = pow(ad,bm1);
924: /* dd = bd * ad; */
925: dd = pow(ad,beta);
926: gd = coef*bd*(td - t0);
927:
928: c[0].k = k-1; c[0].j = j; c[0].i = i;
929: v[0] = - hxhydhz*(dd - gd);
930: c[1].k = k; c[1].j = j; c[1].i = i-1;
931: v[1] = - hyhzdhx*(dw - gw);
932: c[2].k = k; c[2].j = j; c[2].i = i;
933: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
934: c[3].k = k; c[3].j = j+1; c[3].i = i;
935: v[3] = - hzhxdhy*(dn + gn);
936: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
937: }
939: /* right-hand top edge */
940: } else if (j == my-1) {
942: ts = x[k][j-1][i];
943: as = 0.5*(t0 + ts);
944: bs = pow(as,bm1);
945: /* ds = bs * as; */
946: ds = pow(as,beta);
947: gs = coef*bs*(ts - t0);
948:
949: /* right-hand top down corner */
950: if (k == 0) {
952: tu = x[k+1][j][i];
953: au = 0.5*(t0 + tu);
954: bu = pow(au,bm1);
955: /* du = bu * au; */
956: du = pow(au,beta);
957: gu = coef*bu*(tu - t0);
958:
959: c[0].k = k; c[0].j = j-1; c[0].i = i;
960: v[0] = - hzhxdhy*(ds - gs);
961: c[1].k = k; c[1].j = j; c[1].i = i-1;
962: v[1] = - hyhzdhx*(dw - gw);
963: c[2].k = k; c[2].j = j; c[2].i = i;
964: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
965: c[3].k = k+1; c[3].j = j; c[3].i = i;
966: v[3] = - hxhydhz*(du + gu);
967: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
969: /* right-hand top interior edge */
970: } else if (k < mz-1) {
972: tu = x[k+1][j][i];
973: au = 0.5*(t0 + tu);
974: bu = pow(au,bm1);
975: /* du = bu * au; */
976: du = pow(au,beta);
977: gu = coef*bu*(tu - t0);
978:
979: td = x[k-1][j][i];
980: ad = 0.5*(t0 + td);
981: bd = pow(ad,bm1);
982: /* dd = bd * ad; */
983: dd = pow(ad,beta);
984: gd = coef*bd*(td - t0);
985:
986: c[0].k = k-1; c[0].j = j; c[0].i = i;
987: v[0] = - hxhydhz*(dd - gd);
988: c[1].k = k; c[1].j = j-1; c[1].i = i;
989: v[1] = - hzhxdhy*(ds - gs);
990: c[2].k = k; c[2].j = j; c[2].i = i-1;
991: v[2] = - hyhzdhx*(dw - gw);
992: c[3].k = k; c[3].j = j; c[3].i = i;
993: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
994: c[4].k = k+1; c[4].j = j; c[4].i = i;
995: v[4] = - hxhydhz*(du + gu);
996: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
998: /* right-hand top up corner */
999: } else {
1001: td = x[k-1][j][i];
1002: ad = 0.5*(t0 + td);
1003: bd = pow(ad,bm1);
1004: /* dd = bd * ad; */
1005: dd = pow(ad,beta);
1006: gd = coef*bd*(td - t0);
1007:
1008: c[0].k = k-1; c[0].j = j; c[0].i = i;
1009: v[0] = - hxhydhz*(dd - gd);
1010: c[1].k = k; c[1].j = j-1; c[1].i = i;
1011: v[1] = - hzhxdhy*(ds - gs);
1012: c[2].k = k; c[2].j = j; c[2].i = i-1;
1013: v[2] = - hyhzdhx*(dw - gw);
1014: c[3].k = k; c[3].j = j; c[3].i = i;
1015: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1016: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
1017: }
1019: } else {
1021: ts = x[k][j-1][i];
1022: as = 0.5*(t0 + ts);
1023: bs = pow(as,bm1);
1024: /* ds = bs * as; */
1025: ds = pow(as,beta);
1026: gs = coef*bs*(t0 - ts);
1027:
1028: tn = x[k][j+1][i];
1029: an = 0.5*(t0 + tn);
1030: bn = pow(an,bm1);
1031: /* dn = bn * an; */
1032: dn = pow(an,beta);
1033: gn = coef*bn*(tn - t0);
1035: /* right-hand down interior edge */
1036: if (k == 0) {
1038: tu = x[k+1][j][i];
1039: au = 0.5*(t0 + tu);
1040: bu = pow(au,bm1);
1041: /* du = bu * au; */
1042: du = pow(au,beta);
1043: gu = coef*bu*(tu - t0);
1045: c[0].k = k; c[0].j = j-1; c[0].i = i;
1046: v[0] = - hzhxdhy*(ds - gs);
1047: c[1].k = k; c[1].j = j; c[1].i = i-1;
1048: v[1] = - hyhzdhx*(dw - gw);
1049: c[2].k = k; c[2].j = j; c[2].i = i;
1050: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1051: c[3].k = k; c[3].j = j+1; c[3].i = i;
1052: v[3] = - hzhxdhy*(dn + gn);
1053: c[4].k = k+1; c[4].j = j; c[4].i = i;
1054: v[4] = - hxhydhz*(du + gu);
1055: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1056: }
1058: /* right-hand up interior edge */
1059: else if (k == mz-1) {
1061: td = x[k-1][j][i];
1062: ad = 0.5*(t0 + td);
1063: bd = pow(ad,bm1);
1064: /* dd = bd * ad; */
1065: dd = pow(ad,beta);
1066: gd = coef*bd*(t0 - td);
1067:
1068: c[0].k = k-1; c[0].j = j; c[0].i = i;
1069: v[0] = - hxhydhz*(dd - gd);
1070: c[1].k = k; c[1].j = j-1; c[1].i = i;
1071: v[1] = - hzhxdhy*(ds - gs);
1072: c[2].k = k; c[2].j = j; c[2].i = i-1;
1073: v[2] = - hyhzdhx*(dw - gw);
1074: c[3].k = k; c[3].j = j; c[3].i = i;
1075: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1076: c[4].k = k; c[4].j = j+1; c[4].i = i;
1077: v[4] = - hzhxdhy*(dn + gn);
1078: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1079: }
1081: /* right-hand interior plane */
1082: else {
1084: td = x[k-1][j][i];
1085: ad = 0.5*(t0 + td);
1086: bd = pow(ad,bm1);
1087: /* dd = bd * ad; */
1088: dd = pow(ad,beta);
1089: gd = coef*bd*(t0 - td);
1090:
1091: tu = x[k+1][j][i];
1092: au = 0.5*(t0 + tu);
1093: bu = pow(au,bm1);
1094: /* du = bu * au; */
1095: du = pow(au,beta);
1096: gu = coef*bu*(tu - t0);
1098: c[0].k = k-1; c[0].j = j; c[0].i = i;
1099: v[0] = - hxhydhz*(dd - gd);
1100: c[1].k = k; c[1].j = j-1; c[1].i = i;
1101: v[1] = - hzhxdhy*(ds - gs);
1102: c[2].k = k; c[2].j = j; c[2].i = i-1;
1103: v[2] = - hyhzdhx*(dw - gw);
1104: c[3].k = k; c[3].j = j; c[3].i = i;
1105: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1106: c[4].k = k; c[4].j = j+1; c[4].i = i;
1107: v[4] = - hzhxdhy*(dn + gn);
1108: c[5].k = k+1; c[5].j = j; c[5].i = i;
1109: v[5] = - hxhydhz*(du + gu);
1110: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1111: }
1112: }
1114: } else if (j == 0) {
1116: tw = x[k][j][i-1];
1117: aw = 0.5*(t0 + tw);
1118: bw = pow(aw,bm1);
1119: /* dw = bw * aw */
1120: dw = pow(aw,beta);
1121: gw = coef*bw*(t0 - tw);
1123: te = x[k][j][i+1];
1124: ae = 0.5*(t0 + te);
1125: be = pow(ae,bm1);
1126: /* de = be * ae; */
1127: de = pow(ae,beta);
1128: ge = coef*be*(te - t0);
1130: tn = x[k][j+1][i];
1131: an = 0.5*(t0 + tn);
1132: bn = pow(an,bm1);
1133: /* dn = bn * an; */
1134: dn = pow(an,beta);
1135: gn = coef*bn*(tn - t0);
1138: /* bottom down interior edge */
1139: if (k == 0) {
1141: tu = x[k+1][j][i];
1142: au = 0.5*(t0 + tu);
1143: bu = pow(au,bm1);
1144: /* du = bu * au; */
1145: du = pow(au,beta);
1146: gu = coef*bu*(tu - t0);
1148: c[0].k = k; c[0].j = j; c[0].i = i-1;
1149: v[0] = - hyhzdhx*(dw - gw);
1150: c[1].k = k; c[1].j = j; c[1].i = i;
1151: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1152: c[2].k = k; c[2].j = j; c[2].i = i+1;
1153: v[2] = - hyhzdhx*(de + ge);
1154: c[3].k = k; c[3].j = j+1; c[3].i = i;
1155: v[3] = - hzhxdhy*(dn + gn);
1156: c[4].k = k+1; c[4].j = j; c[4].i = i;
1157: v[4] = - hxhydhz*(du + gu);
1158: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1159: }
1161: /* bottom up interior edge */
1162: else if (k == mz-1) {
1164: td = x[k-1][j][i];
1165: ad = 0.5*(t0 + td);
1166: bd = pow(ad,bm1);
1167: /* dd = bd * ad; */
1168: dd = pow(ad,beta);
1169: gd = coef*bd*(td - t0);
1171: c[0].k = k-1; c[0].j = j; c[0].i = i;
1172: v[0] = - hxhydhz*(dd - gd);
1173: c[1].k = k; c[1].j = j; c[1].i = i-1;
1174: v[1] = - hyhzdhx*(dw - gw);
1175: c[2].k = k; c[2].j = j; c[2].i = i;
1176: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1177: c[3].k = k; c[3].j = j; c[3].i = i+1;
1178: v[3] = - hyhzdhx*(de + ge);
1179: c[4].k = k; c[4].j = j+1; c[4].i = i;
1180: v[4] = - hzhxdhy*(dn + gn);
1181: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1182: }
1184: /* bottom interior plane */
1185: else {
1187: tu = x[k+1][j][i];
1188: au = 0.5*(t0 + tu);
1189: bu = pow(au,bm1);
1190: /* du = bu * au; */
1191: du = pow(au,beta);
1192: gu = coef*bu*(tu - t0);
1194: td = x[k-1][j][i];
1195: ad = 0.5*(t0 + td);
1196: bd = pow(ad,bm1);
1197: /* dd = bd * ad; */
1198: dd = pow(ad,beta);
1199: gd = coef*bd*(td - t0);
1201: c[0].k = k-1; c[0].j = j; c[0].i = i;
1202: v[0] = - hxhydhz*(dd - gd);
1203: c[1].k = k; c[1].j = j; c[1].i = i-1;
1204: v[1] = - hyhzdhx*(dw - gw);
1205: c[2].k = k; c[2].j = j; c[2].i = i;
1206: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1207: c[3].k = k; c[3].j = j; c[3].i = i+1;
1208: v[3] = - hyhzdhx*(de + ge);
1209: c[4].k = k; c[4].j = j+1; c[4].i = i;
1210: v[4] = - hzhxdhy*(dn + gn);
1211: c[5].k = k+1; c[5].j = j; c[5].i = i;
1212: v[5] = - hxhydhz*(du + gu);
1213: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1214: }
1216: } else if (j == my-1) {
1218: tw = x[k][j][i-1];
1219: aw = 0.5*(t0 + tw);
1220: bw = pow(aw,bm1);
1221: /* dw = bw * aw */
1222: dw = pow(aw,beta);
1223: gw = coef*bw*(t0 - tw);
1225: te = x[k][j][i+1];
1226: ae = 0.5*(t0 + te);
1227: be = pow(ae,bm1);
1228: /* de = be * ae; */
1229: de = pow(ae,beta);
1230: ge = coef*be*(te - t0);
1232: ts = x[k][j-1][i];
1233: as = 0.5*(t0 + ts);
1234: bs = pow(as,bm1);
1235: /* ds = bs * as; */
1236: ds = pow(as,beta);
1237: gs = coef*bs*(t0 - ts);
1238:
1239: /* top down interior edge */
1240: if (k == 0) {
1242: tu = x[k+1][j][i];
1243: au = 0.5*(t0 + tu);
1244: bu = pow(au,bm1);
1245: /* du = bu * au; */
1246: du = pow(au,beta);
1247: gu = coef*bu*(tu - t0);
1249: c[0].k = k; c[0].j = j-1; c[0].i = i;
1250: v[0] = - hzhxdhy*(ds - gs);
1251: c[1].k = k; c[1].j = j; c[1].i = i-1;
1252: v[1] = - hyhzdhx*(dw - gw);
1253: c[2].k = k; c[2].j = j; c[2].i = i;
1254: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1255: c[3].k = k; c[3].j = j; c[3].i = i+1;
1256: v[3] = - hyhzdhx*(de + ge);
1257: c[4].k = k+1; c[4].j = j; c[4].i = i;
1258: v[4] = - hxhydhz*(du + gu);
1259: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1260: }
1262: /* top up interior edge */
1263: else if (k == mz-1) {
1265: td = x[k-1][j][i];
1266: ad = 0.5*(t0 + td);
1267: bd = pow(ad,bm1);
1268: /* dd = bd * ad; */
1269: dd = pow(ad,beta);
1270: gd = coef*bd*(td - t0);
1272: c[0].k = k-1; c[0].j = j; c[0].i = i;
1273: v[0] = - hxhydhz*(dd - gd);
1274: c[1].k = k; c[1].j = j-1; c[1].i = i;
1275: v[1] = - hzhxdhy*(ds - gs);
1276: c[2].k = k; c[2].j = j; c[2].i = i-1;
1277: v[2] = - hyhzdhx*(dw - gw);
1278: c[3].k = k; c[3].j = j; c[3].i = i;
1279: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1280: c[4].k = k; c[4].j = j; c[4].i = i+1;
1281: v[4] = - hyhzdhx*(de + ge);
1282: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1283: }
1285: /* top interior plane */
1286: else {
1288: tu = x[k+1][j][i];
1289: au = 0.5*(t0 + tu);
1290: bu = pow(au,bm1);
1291: /* du = bu * au; */
1292: du = pow(au,beta);
1293: gu = coef*bu*(tu - t0);
1295: td = x[k-1][j][i];
1296: ad = 0.5*(t0 + td);
1297: bd = pow(ad,bm1);
1298: /* dd = bd * ad; */
1299: dd = pow(ad,beta);
1300: gd = coef*bd*(td - t0);
1302: c[0].k = k-1; c[0].j = j; c[0].i = i;
1303: v[0] = - hxhydhz*(dd - gd);
1304: c[1].k = k; c[1].j = j-1; c[1].i = i;
1305: v[1] = - hzhxdhy*(ds - gs);
1306: c[2].k = k; c[2].j = j; c[2].i = i-1;
1307: v[2] = - hyhzdhx*(dw - gw);
1308: c[3].k = k; c[3].j = j; c[3].i = i;
1309: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1310: c[4].k = k; c[4].j = j; c[4].i = i+1;
1311: v[4] = - hyhzdhx*(de + ge);
1312: c[5].k = k+1; c[5].j = j; c[5].i = i;
1313: v[5] = - hxhydhz*(du + gu);
1314: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1315: }
1317: } else if (k == 0) {
1319: /* down interior plane */
1321: tw = x[k][j][i-1];
1322: aw = 0.5*(t0 + tw);
1323: bw = pow(aw,bm1);
1324: /* dw = bw * aw */
1325: dw = pow(aw,beta);
1326: gw = coef*bw*(t0 - tw);
1328: te = x[k][j][i+1];
1329: ae = 0.5*(t0 + te);
1330: be = pow(ae,bm1);
1331: /* de = be * ae; */
1332: de = pow(ae,beta);
1333: ge = coef*be*(te - t0);
1335: ts = x[k][j-1][i];
1336: as = 0.5*(t0 + ts);
1337: bs = pow(as,bm1);
1338: /* ds = bs * as; */
1339: ds = pow(as,beta);
1340: gs = coef*bs*(t0 - ts);
1341:
1342: tn = x[k][j+1][i];
1343: an = 0.5*(t0 + tn);
1344: bn = pow(an,bm1);
1345: /* dn = bn * an; */
1346: dn = pow(an,beta);
1347: gn = coef*bn*(tn - t0);
1348:
1349: tu = x[k+1][j][i];
1350: au = 0.5*(t0 + tu);
1351: bu = pow(au,bm1);
1352: /* du = bu * au; */
1353: du = pow(au,beta);
1354: gu = coef*bu*(tu - t0);
1356: c[0].k = k; c[0].j = j-1; c[0].i = i;
1357: v[0] = - hzhxdhy*(ds - gs);
1358: c[1].k = k; c[1].j = j; c[1].i = i-1;
1359: v[1] = - hyhzdhx*(dw - gw);
1360: c[2].k = k; c[2].j = j; c[2].i = i;
1361: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1362: c[3].k = k; c[3].j = j; c[3].i = i+1;
1363: v[3] = - hyhzdhx*(de + ge);
1364: c[4].k = k; c[4].j = j+1; c[4].i = i;
1365: v[4] = - hzhxdhy*(dn + gn);
1366: c[5].k = k+1; c[5].j = j; c[5].i = i;
1367: v[5] = - hxhydhz*(du + gu);
1368: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1369: }
1370:
1371: else if (k == mz-1) {
1373: /* up interior plane */
1375: tw = x[k][j][i-1];
1376: aw = 0.5*(t0 + tw);
1377: bw = pow(aw,bm1);
1378: /* dw = bw * aw */
1379: dw = pow(aw,beta);
1380: gw = coef*bw*(t0 - tw);
1382: te = x[k][j][i+1];
1383: ae = 0.5*(t0 + te);
1384: be = pow(ae,bm1);
1385: /* de = be * ae; */
1386: de = pow(ae,beta);
1387: ge = coef*be*(te - t0);
1389: ts = x[k][j-1][i];
1390: as = 0.5*(t0 + ts);
1391: bs = pow(as,bm1);
1392: /* ds = bs * as; */
1393: ds = pow(as,beta);
1394: gs = coef*bs*(t0 - ts);
1395:
1396: tn = x[k][j+1][i];
1397: an = 0.5*(t0 + tn);
1398: bn = pow(an,bm1);
1399: /* dn = bn * an; */
1400: dn = pow(an,beta);
1401: gn = coef*bn*(tn - t0);
1402:
1403: td = x[k-1][j][i];
1404: ad = 0.5*(t0 + td);
1405: bd = pow(ad,bm1);
1406: /* dd = bd * ad; */
1407: dd = pow(ad,beta);
1408: gd = coef*bd*(t0 - td);
1409:
1410: c[0].k = k-1; c[0].j = j; c[0].i = i;
1411: v[0] = - hxhydhz*(dd - gd);
1412: c[1].k = k; c[1].j = j-1; c[1].i = i;
1413: v[1] = - hzhxdhy*(ds - gs);
1414: c[2].k = k; c[2].j = j; c[2].i = i-1;
1415: v[2] = - hyhzdhx*(dw - gw);
1416: c[3].k = k; c[3].j = j; c[3].i = i;
1417: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1418: c[4].k = k; c[4].j = j; c[4].i = i+1;
1419: v[4] = - hyhzdhx*(de + ge);
1420: c[5].k = k; c[5].j = j+1; c[5].i = i;
1421: v[5] = - hzhxdhy*(dn + gn);
1422: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1423: }
1424: }
1425: }
1426: }
1427: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
1428: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
1429: if (jac != *J) {
1430: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1431: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1432: }
1433: DAVecRestoreArray((DA)dmmg->dm,localX,&x);
1434: DARestoreLocalVector((DA)dmmg->dm,&localX);
1436: PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
1437: return(0);
1438: }