Actual source code: ex20.c
1: /* $Id: ex20.c,v 1.20 2001/08/07 21:31:17 bsmith Exp $ */
4: static char help[] ="Nonlinear Radiative Transport PDE with multigrid in 3d.n
5: Uses 3-dimensional distributed arrays.n
6: A 3-dim simplified Radiative Transport test problem is used, with analytic Jacobian. n
7: n
8: Solves the linear systems via multilevel methods n
9: n
10: The command linen
11: options are:n
12: -tleft <tl>, where <tl> indicates the left Diriclet BC n
13: -tright <tr>, where <tr> indicates the right Diriclet BC n
14: -beta <beta>, where <beta> indicates the exponent in T nn";
16: /*T
17: Concepts: SNES^solving a system of nonlinear equations
18: Concepts: DA^using distributed arrays
19: Concepts: multigrid;
20: Processors: n
21: T*/
23: /*
24:
25: This example models the partial differential equation
26:
27: - Div(alpha* T^beta (GRAD T)) = 0.
28:
29: where beta = 2.5 and alpha = 1.0
30:
31: BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = dT/dn_up = dT/dn_down = 0.
32:
33: in the unit square, which is uniformly discretized in each of x and
34: y in this simple encoding. The degrees of freedom are cell centered.
35:
36: A finite volume approximation with the usual 7-point stencil
37: is used to discretize the boundary value problem to obtain a
38: nonlinear system of equations.
40: This code was contributed by Nickolas Jovanovic based on ex18.c
41:
42: */
44: #include petscsnes.h
45: #include petscda.h
46: #include petscmg.h
48: /* User-defined application context */
50: typedef struct {
51: PetscReal tleft,tright; /* Dirichlet boundary conditions */
52: PetscReal beta,bm1,coef; /* nonlinear diffusivity parameterizations */
53: } AppCtx;
55: #define POWFLOP 5 /* assume a pow() takes five flops */
57: extern int FormInitialGuess(SNES,Vec,void*);
58: extern int FormFunction(SNES,Vec,Vec,void*);
59: extern int FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
61: #undef __FUNCT__
63: int main(int argc,char **argv)
64: {
65: DMMG *dmmg;
66: SNES snes;
67: AppCtx user;
68: int ierr,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 = %dn",its);
118: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %dn",lits);
119: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / Newton = %en",litspit);
121: DMMGDestroy(dmmg);
122: PetscFinalize();
124: return 0;
125: }
126: /* -------------------- Form initial approximation ----------------- */
127: #undef __FUNCT__
129: int FormInitialGuess(SNES snes,Vec X,void *ptr)
130: {
131: DMMG dmmg = (DMMG)ptr;
132: AppCtx *user = (AppCtx*)dmmg->user;
133: int i,j,k,ierr,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,(void**)&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,(void**)&x);
152: return(0);
153: }
154: /* -------------------- Evaluate Function F(x) --------------------- */
155: #undef __FUNCT__
157: int FormFunction(SNES snes,Vec X,Vec F,void* ptr)
158: {
159: DMMG dmmg = (DMMG)ptr;
160: AppCtx *user = (AppCtx*)dmmg->user;
161: int ierr,i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
162: PetscScalar zero = 0.0,one = 1.0;
163: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
164: 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;
165: PetscScalar tleft,tright,beta,td,ad,dd,fd,tu,au,du,fu;
166: PetscScalar ***x,***f;
167: Vec localX;
170: DAGetLocalVector((DA)dmmg->dm,&localX);
171: DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,&mz,0,0,0,0,0,0,0);
172: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
173: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
174: tleft = user->tleft; tright = user->tright;
175: beta = user->beta;
176:
177: /* Get ghost points */
178: DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
179: DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
180: DAGetCorners((DA)dmmg->dm,&xs,&ys,&zs,&xm,&ym,&zm);
181: DAVecGetArray((DA)dmmg->dm,localX,(void**)&x);
182: DAVecGetArray((DA)dmmg->dm,F,(void**)&f);
184: /* Evaluate function */
185: for (k=zs; k<zs+zm; k++) {
186: for (j=ys; j<ys+ym; j++) {
187: for (i=xs; i<xs+xm; i++) {
188: t0 = x[k][j][i];
190: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
192: /* general interior volume */
194: tw = x[k][j][i-1];
195: aw = 0.5*(t0 + tw);
196: dw = pow(aw,beta);
197: fw = dw*(t0 - tw);
199: te = x[k][j][i+1];
200: ae = 0.5*(t0 + te);
201: de = pow(ae,beta);
202: fe = de*(te - t0);
204: ts = x[k][j-1][i];
205: as = 0.5*(t0 + ts);
206: ds = pow(as,beta);
207: fs = ds*(t0 - ts);
208:
209: tn = x[k][j+1][i];
210: an = 0.5*(t0 + tn);
211: dn = pow(an,beta);
212: fn = dn*(tn - t0);
214: td = x[k-1][j][i];
215: ad = 0.5*(t0 + td);
216: dd = pow(ad,beta);
217: fd = dd*(t0 - td);
219: tu = x[k+1][j][i];
220: au = 0.5*(t0 + tu);
221: du = pow(au,beta);
222: fu = du*(tu - t0);
224: } else if (i == 0) {
226: /* left-hand (west) boundary */
227: tw = tleft;
228: aw = 0.5*(t0 + tw);
229: dw = pow(aw,beta);
230: fw = dw*(t0 - tw);
232: te = x[k][j][i+1];
233: ae = 0.5*(t0 + te);
234: de = pow(ae,beta);
235: fe = de*(te - t0);
237: if (j > 0) {
238: ts = x[k][j-1][i];
239: as = 0.5*(t0 + ts);
240: ds = pow(as,beta);
241: fs = ds*(t0 - ts);
242: } else {
243: fs = zero;
244: }
246: if (j < my-1) {
247: tn = x[k][j+1][i];
248: an = 0.5*(t0 + tn);
249: dn = pow(an,beta);
250: fn = dn*(tn - t0);
251: } else {
252: fn = zero;
253: }
255: if (k > 0) {
256: td = x[k-1][j][i];
257: ad = 0.5*(t0 + td);
258: dd = pow(ad,beta);
259: fd = dd*(t0 - td);
260: } else {
261: fd = zero;
262: }
264: if (k < mz-1) {
265: tu = x[k+1][j][i];
266: au = 0.5*(t0 + tu);
267: du = pow(au,beta);
268: fu = du*(tu - t0);
269: } else {
270: fu = zero;
271: }
273: } else if (i == mx-1) {
275: /* right-hand (east) boundary */
276: tw = x[k][j][i-1];
277: aw = 0.5*(t0 + tw);
278: dw = pow(aw,beta);
279: fw = dw*(t0 - tw);
280:
281: te = tright;
282: ae = 0.5*(t0 + te);
283: de = pow(ae,beta);
284: fe = de*(te - t0);
285:
286: if (j > 0) {
287: ts = x[k][j-1][i];
288: as = 0.5*(t0 + ts);
289: ds = pow(as,beta);
290: fs = ds*(t0 - ts);
291: } else {
292: fs = zero;
293: }
294:
295: if (j < my-1) {
296: tn = x[k][j+1][i];
297: an = 0.5*(t0 + tn);
298: dn = pow(an,beta);
299: fn = dn*(tn - t0);
300: } else {
301: fn = zero;
302: }
304: if (k > 0) {
305: td = x[k-1][j][i];
306: ad = 0.5*(t0 + td);
307: dd = pow(ad,beta);
308: fd = dd*(t0 - td);
309: } else {
310: fd = zero;
311: }
313: if (k < mz-1) {
314: tu = x[k+1][j][i];
315: au = 0.5*(t0 + tu);
316: du = pow(au,beta);
317: fu = du*(tu - t0);
318: } else {
319: fu = zero;
320: }
322: } else if (j == 0) {
324: /* bottom (south) boundary, and i <> 0 or mx-1 */
325: tw = x[k][j][i-1];
326: aw = 0.5*(t0 + tw);
327: dw = pow(aw,beta);
328: fw = dw*(t0 - tw);
330: te = x[k][j][i+1];
331: ae = 0.5*(t0 + te);
332: de = pow(ae,beta);
333: fe = de*(te - t0);
335: fs = zero;
337: tn = x[k][j+1][i];
338: an = 0.5*(t0 + tn);
339: dn = pow(an,beta);
340: fn = dn*(tn - t0);
342: if (k > 0) {
343: td = x[k-1][j][i];
344: ad = 0.5*(t0 + td);
345: dd = pow(ad,beta);
346: fd = dd*(t0 - td);
347: } else {
348: fd = zero;
349: }
351: if (k < mz-1) {
352: tu = x[k+1][j][i];
353: au = 0.5*(t0 + tu);
354: du = pow(au,beta);
355: fu = du*(tu - t0);
356: } else {
357: fu = zero;
358: }
360: } else if (j == my-1) {
362: /* top (north) boundary, and i <> 0 or mx-1 */
363: tw = x[k][j][i-1];
364: aw = 0.5*(t0 + tw);
365: dw = pow(aw,beta);
366: fw = dw*(t0 - tw);
368: te = x[k][j][i+1];
369: ae = 0.5*(t0 + te);
370: de = pow(ae,beta);
371: fe = de*(te - t0);
373: ts = x[k][j-1][i];
374: as = 0.5*(t0 + ts);
375: ds = pow(as,beta);
376: fs = ds*(t0 - ts);
378: fn = zero;
380: if (k > 0) {
381: td = x[k-1][j][i];
382: ad = 0.5*(t0 + td);
383: dd = pow(ad,beta);
384: fd = dd*(t0 - td);
385: } else {
386: fd = zero;
387: }
389: if (k < mz-1) {
390: tu = x[k+1][j][i];
391: au = 0.5*(t0 + tu);
392: du = pow(au,beta);
393: fu = du*(tu - t0);
394: } else {
395: fu = zero;
396: }
398: } else if (k == 0) {
400: /* down boundary (interior only) */
401: tw = x[k][j][i-1];
402: aw = 0.5*(t0 + tw);
403: dw = pow(aw,beta);
404: fw = dw*(t0 - tw);
406: te = x[k][j][i+1];
407: ae = 0.5*(t0 + te);
408: de = pow(ae,beta);
409: fe = de*(te - t0);
411: ts = x[k][j-1][i];
412: as = 0.5*(t0 + ts);
413: ds = pow(as,beta);
414: fs = ds*(t0 - ts);
416: tn = x[k][j+1][i];
417: an = 0.5*(t0 + tn);
418: dn = pow(an,beta);
419: fn = dn*(tn - t0);
421: fd = zero;
423: tu = x[k+1][j][i];
424: au = 0.5*(t0 + tu);
425: du = pow(au,beta);
426: fu = du*(tu - t0);
427:
428: } else if (k == mz-1) {
430: /* up boundary (interior only) */
431: tw = x[k][j][i-1];
432: aw = 0.5*(t0 + tw);
433: dw = pow(aw,beta);
434: fw = dw*(t0 - tw);
436: te = x[k][j][i+1];
437: ae = 0.5*(t0 + te);
438: de = pow(ae,beta);
439: fe = de*(te - t0);
441: ts = x[k][j-1][i];
442: as = 0.5*(t0 + ts);
443: ds = pow(as,beta);
444: fs = ds*(t0 - ts);
446: tn = x[k][j+1][i];
447: an = 0.5*(t0 + tn);
448: dn = pow(an,beta);
449: fn = dn*(tn - t0);
451: td = x[k-1][j][i];
452: ad = 0.5*(t0 + td);
453: dd = pow(ad,beta);
454: fd = dd*(t0 - td);
456: fu = zero;
457: }
459: f[k][j][i] = - hyhzdhx*(fe-fw) - hzhxdhy*(fn-fs) - hxhydhz*(fu-fd);
460: }
461: }
462: }
463: DAVecRestoreArray((DA)dmmg->dm,localX,(void**)&x);
464: DAVecRestoreArray((DA)dmmg->dm,F,(void**)&f);
465: DARestoreLocalVector((DA)dmmg->dm,&localX);
466: PetscLogFlops((22 + 4*POWFLOP)*ym*xm);
467: return(0);
468: }
469: /* -------------------- Evaluate Jacobian F(x) --------------------- */
470: #undef __FUNCT__
472: int FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ptr)
473: {
474: DMMG dmmg = (DMMG)ptr;
475: AppCtx *user = (AppCtx*)dmmg->user;
476: int ierr,i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
477: PetscScalar zero = 0.0,one = 1.0;
478: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
479: 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;
480: PetscScalar tleft,tright,beta,td,ad,dd,fd,tu,au,du,fu,v[7],bm1,coef;
481: PetscScalar ***x,bn,bs,be,bw,bu,bd,gn,gs,ge,gw,gu,gd;
482: Vec localX;
483: MatStencil c[7],row;
484: Mat jac = *B;
487: DAGetLocalVector((DA)dmmg->dm,&localX);
488: DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,&mz,0,0,0,0,0,0,0);
489: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
490: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
491: tleft = user->tleft; tright = user->tright;
492: beta = user->beta; bm1 = user->bm1; coef = user->coef;
494: /* Get ghost points */
495: DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
496: DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
497: DAGetCorners((DA)dmmg->dm,&xs,&ys,&zs,&xm,&ym,&zm);
498: DAVecGetArray((DA)dmmg->dm,localX,(void**)&x);
500: /* Evaluate Jacobian of function */
501: for (k=zs; k<zs+zm; k++) {
502: for (j=ys; j<ys+ym; j++) {
503: for (i=xs; i<xs+xm; i++) {
504: t0 = x[k][j][i];
505: row.k = k; row.j = j; row.i = i;
506: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
508: /* general interior volume */
510: tw = x[k][j][i-1];
511: aw = 0.5*(t0 + tw);
512: bw = pow(aw,bm1);
513: /* dw = bw * aw */
514: dw = pow(aw,beta);
515: gw = coef*bw*(t0 - tw);
517: te = x[k][j][i+1];
518: ae = 0.5*(t0 + te);
519: be = pow(ae,bm1);
520: /* de = be * ae; */
521: de = pow(ae,beta);
522: ge = coef*be*(te - t0);
524: ts = x[k][j-1][i];
525: as = 0.5*(t0 + ts);
526: bs = pow(as,bm1);
527: /* ds = bs * as; */
528: ds = pow(as,beta);
529: gs = coef*bs*(t0 - ts);
530:
531: tn = x[k][j+1][i];
532: an = 0.5*(t0 + tn);
533: bn = pow(an,bm1);
534: /* dn = bn * an; */
535: dn = pow(an,beta);
536: gn = coef*bn*(tn - t0);
538: td = x[k-1][j][i];
539: ad = 0.5*(t0 + td);
540: bd = pow(ad,bm1);
541: /* dd = bd * ad; */
542: dd = pow(ad,beta);
543: gd = coef*bd*(t0 - td);
544:
545: tu = x[k+1][j][i];
546: au = 0.5*(t0 + tu);
547: bu = pow(au,bm1);
548: /* du = bu * au; */
549: du = pow(au,beta);
550: gu = coef*bu*(tu - t0);
552: c[0].k = k-1; c[0].j = j; c[0].i = i; v[0] = - hxhydhz*(dd - gd);
553: c[1].k = k; c[1].j = j-1; c[1].i = i;
554: v[1] = - hzhxdhy*(ds - gs);
555: c[2].k = k; c[2].j = j; c[2].i = i-1;
556: v[2] = - hyhzdhx*(dw - gw);
557: c[3].k = k; c[3].j = j; c[3].i = i;
558: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
559: c[4].k = k; c[4].j = j; c[4].i = i+1;
560: v[4] = - hyhzdhx*(de + ge);
561: c[5].k = k; c[5].j = j+1; c[5].i = i;
562: v[5] = - hzhxdhy*(dn + gn);
563: c[6].k = k+1; c[6].j = j; c[6].i = i;
564: v[6] = - hxhydhz*(du + gu);
565: ierr = MatSetValuesStencil(jac,1,&row,7,c,v,INSERT_VALUES);
567: } else if (i == 0) {
569: /* left-hand plane boundary */
570: tw = tleft;
571: aw = 0.5*(t0 + tw);
572: bw = pow(aw,bm1);
573: /* dw = bw * aw */
574: dw = pow(aw,beta);
575: gw = coef*bw*(t0 - tw);
576:
577: te = x[k][j][i+1];
578: ae = 0.5*(t0 + te);
579: be = pow(ae,bm1);
580: /* de = be * ae; */
581: de = pow(ae,beta);
582: ge = coef*be*(te - t0);
583:
584: /* left-hand bottom edge */
585: if (j == 0) {
587: tn = x[k][j+1][i];
588: an = 0.5*(t0 + tn);
589: bn = pow(an,bm1);
590: /* dn = bn * an; */
591: dn = pow(an,beta);
592: gn = coef*bn*(tn - t0);
593:
594: /* left-hand bottom down corner */
595: if (k == 0) {
597: tu = x[k+1][j][i];
598: au = 0.5*(t0 + tu);
599: bu = pow(au,bm1);
600: /* du = bu * au; */
601: du = pow(au,beta);
602: gu = coef*bu*(tu - t0);
603:
604: c[0].k = k; c[0].j = j; c[0].i = i;
605: v[0] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
606: c[1].k = k; c[1].j = j; c[1].i = i+1;
607: v[1] = - hyhzdhx*(de + ge);
608: c[2].k = k; c[2].j = j+1; c[2].i = i;
609: v[2] = - hzhxdhy*(dn + gn);
610: c[3].k = k+1; c[3].j = j; c[3].i = i;
611: v[3] = - hxhydhz*(du + gu);
612: ierr = MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
614: /* left-hand bottom interior edge */
615: } else if (k < mz-1) {
617: tu = x[k+1][j][i];
618: au = 0.5*(t0 + tu);
619: bu = pow(au,bm1);
620: /* du = bu * au; */
621: du = pow(au,beta);
622: gu = coef*bu*(tu - t0);
623:
624: td = x[k-1][j][i];
625: ad = 0.5*(t0 + td);
626: bd = pow(ad,bm1);
627: /* dd = bd * ad; */
628: dd = pow(ad,beta);
629: gd = coef*bd*(td - t0);
630:
631: c[0].k = k-1; c[0].j = j; c[0].i = i;
632: v[0] = - hxhydhz*(dd - gd);
633: c[1].k = k; c[1].j = j; c[1].i = i;
634: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
635: c[2].k = k; c[2].j = j; c[2].i = i+1;
636: v[2] = - hyhzdhx*(de + ge);
637: c[3].k = k; c[3].j = j+1; c[3].i = i;
638: v[3] = - hzhxdhy*(dn + gn);
639: c[4].k = k+1; c[4].j = j; c[4].i = i;
640: v[4] = - hxhydhz*(du + gu);
641: ierr = MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
643: /* left-hand bottom up corner */
644: } else {
646: td = x[k-1][j][i];
647: ad = 0.5*(t0 + td);
648: bd = pow(ad,bm1);
649: /* dd = bd * ad; */
650: dd = pow(ad,beta);
651: gd = coef*bd*(td - t0);
652:
653: c[0].k = k-1; c[0].j = j; c[0].i = i;
654: v[0] = - hxhydhz*(dd - gd);
655: c[1].k = k; c[1].j = j; c[1].i = i;
656: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
657: c[2].k = k; c[2].j = j; c[2].i = i+1;
658: v[2] = - hyhzdhx*(de + ge);
659: c[3].k = k; c[3].j = j+1; c[3].i = i;
660: v[3] = - hzhxdhy*(dn + gn);
661: ierr = MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
662: }
664: /* left-hand top edge */
665: } else if (j == my-1) {
667: ts = x[k][j-1][i];
668: as = 0.5*(t0 + ts);
669: bs = pow(as,bm1);
670: /* ds = bs * as; */
671: ds = pow(as,beta);
672: gs = coef*bs*(ts - t0);
673:
674: /* left-hand top down corner */
675: if (k == 0) {
677: tu = x[k+1][j][i];
678: au = 0.5*(t0 + tu);
679: bu = pow(au,bm1);
680: /* du = bu * au; */
681: du = pow(au,beta);
682: gu = coef*bu*(tu - t0);
683:
684: c[0].k = k; c[0].j = j-1; c[0].i = i;
685: v[0] = - hzhxdhy*(ds - gs);
686: c[1].k = k; c[1].j = j; c[1].i = i;
687: v[1] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
688: c[2].k = k; c[2].j = j; c[2].i = i+1;
689: v[2] = - hyhzdhx*(de + ge);
690: c[3].k = k+1; c[3].j = j; c[3].i = i;
691: v[3] = - hxhydhz*(du + gu);
692: ierr = MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
694: /* left-hand top interior edge */
695: } else if (k < mz-1) {
697: tu = x[k+1][j][i];
698: au = 0.5*(t0 + tu);
699: bu = pow(au,bm1);
700: /* du = bu * au; */
701: du = pow(au,beta);
702: gu = coef*bu*(tu - t0);
703:
704: td = x[k-1][j][i];
705: ad = 0.5*(t0 + td);
706: bd = pow(ad,bm1);
707: /* dd = bd * ad; */
708: dd = pow(ad,beta);
709: gd = coef*bd*(td - t0);
710:
711: c[0].k = k-1; c[0].j = j; c[0].i = i;
712: v[0] = - hxhydhz*(dd - gd);
713: c[1].k = k; c[1].j = j-1; c[1].i = i;
714: v[1] = - hzhxdhy*(ds - gs);
715: c[2].k = k; c[2].j = j; c[2].i = i;
716: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
717: c[3].k = k; c[3].j = j; c[3].i = i+1;
718: v[3] = - hyhzdhx*(de + ge);
719: c[4].k = k+1; c[4].j = j; c[4].i = i;
720: v[4] = - hxhydhz*(du + gu);
721: ierr = MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
723: /* left-hand top up corner */
724: } else {
726: td = x[k-1][j][i];
727: ad = 0.5*(t0 + td);
728: bd = pow(ad,bm1);
729: /* dd = bd * ad; */
730: dd = pow(ad,beta);
731: gd = coef*bd*(td - t0);
732:
733: c[0].k = k-1; c[0].j = j; c[0].i = i;
734: v[0] = - hxhydhz*(dd - gd);
735: c[1].k = k; c[1].j = j-1; c[1].i = i;
736: v[1] = - hzhxdhy*(ds - gs);
737: c[2].k = k; c[2].j = j; c[2].i = i;
738: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
739: c[3].k = k; c[3].j = j; c[3].i = i+1;
740: v[3] = - hyhzdhx*(de + ge);
741: ierr = MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
742: }
744: } else {
746: ts = x[k][j-1][i];
747: as = 0.5*(t0 + ts);
748: bs = pow(as,bm1);
749: /* ds = bs * as; */
750: ds = pow(as,beta);
751: gs = coef*bs*(t0 - ts);
752:
753: tn = x[k][j+1][i];
754: an = 0.5*(t0 + tn);
755: bn = pow(an,bm1);
756: /* dn = bn * an; */
757: dn = pow(an,beta);
758: gn = coef*bn*(tn - t0);
760: /* left-hand down interior edge */
761: if (k == 0) {
763: tu = x[k+1][j][i];
764: au = 0.5*(t0 + tu);
765: bu = pow(au,bm1);
766: /* du = bu * au; */
767: du = pow(au,beta);
768: gu = coef*bu*(tu - t0);
770: c[0].k = k; c[0].j = j-1; c[0].i = i;
771: v[0] = - hzhxdhy*(ds - gs);
772: c[1].k = k; c[1].j = j; c[1].i = i;
773: v[1] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
774: c[2].k = k; c[2].j = j; c[2].i = i+1;
775: v[2] = - hyhzdhx*(de + ge);
776: c[3].k = k; c[3].j = j+1; c[3].i = i;
777: v[3] = - hzhxdhy*(dn + gn);
778: c[4].k = k+1; c[4].j = j; c[4].i = i;
779: v[4] = - hxhydhz*(du + gu);
780: ierr = MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
781: }
783: /* left-hand up interior edge */
784: else if (k == mz-1) {
786: td = x[k-1][j][i];
787: ad = 0.5*(t0 + td);
788: bd = pow(ad,bm1);
789: /* dd = bd * ad; */
790: dd = pow(ad,beta);
791: gd = coef*bd*(t0 - td);
792:
793: c[0].k = k-1; c[0].j = j; c[0].i = i;
794: v[0] = - hxhydhz*(dd - gd);
795: c[1].k = k; c[1].j = j-1; c[1].i = i;
796: v[1] = - hzhxdhy*(ds - gs);
797: c[2].k = k; c[2].j = j; c[2].i = i;
798: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
799: c[3].k = k; c[3].j = j; c[3].i = i+1;
800: v[3] = - hyhzdhx*(de + ge);
801: c[4].k = k; c[4].j = j+1; c[4].i = i;
802: v[4] = - hzhxdhy*(dn + gn);
803: ierr = MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
804: }
806: /* left-hand interior plane */
807: else {
809: td = x[k-1][j][i];
810: ad = 0.5*(t0 + td);
811: bd = pow(ad,bm1);
812: /* dd = bd * ad; */
813: dd = pow(ad,beta);
814: gd = coef*bd*(t0 - td);
815:
816: tu = x[k+1][j][i];
817: au = 0.5*(t0 + tu);
818: bu = pow(au,bm1);
819: /* du = bu * au; */
820: du = pow(au,beta);
821: gu = coef*bu*(tu - t0);
823: c[0].k = k-1; c[0].j = j; c[0].i = i;
824: v[0] = - hxhydhz*(dd - gd);
825: c[1].k = k; c[1].j = j-1; c[1].i = i;
826: v[1] = - hzhxdhy*(ds - gs);
827: c[2].k = k; c[2].j = j; c[2].i = i;
828: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
829: c[3].k = k; c[3].j = j; c[3].i = i+1;
830: v[3] = - hyhzdhx*(de + ge);
831: c[4].k = k; c[4].j = j+1; c[4].i = i;
832: v[4] = - hzhxdhy*(dn + gn);
833: c[5].k = k+1; c[5].j = j; c[5].i = i;
834: v[5] = - hxhydhz*(du + gu);
835: ierr = MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
836: }
837: }
839: } else if (i == mx-1) {
841: /* right-hand plane boundary */
842: tw = x[k][j][i-1];
843: aw = 0.5*(t0 + tw);
844: bw = pow(aw,bm1);
845: /* dw = bw * aw */
846: dw = pow(aw,beta);
847: gw = coef*bw*(t0 - tw);
848:
849: te = tright;
850: ae = 0.5*(t0 + te);
851: be = pow(ae,bm1);
852: /* de = be * ae; */
853: de = pow(ae,beta);
854: ge = coef*be*(te - t0);
855:
856: /* right-hand bottom edge */
857: if (j == 0) {
859: tn = x[k][j+1][i];
860: an = 0.5*(t0 + tn);
861: bn = pow(an,bm1);
862: /* dn = bn * an; */
863: dn = pow(an,beta);
864: gn = coef*bn*(tn - t0);
865:
866: /* right-hand bottom down corner */
867: if (k == 0) {
869: tu = x[k+1][j][i];
870: au = 0.5*(t0 + tu);
871: bu = pow(au,bm1);
872: /* du = bu * au; */
873: du = pow(au,beta);
874: gu = coef*bu*(tu - t0);
875:
876: c[0].k = k; c[0].j = j; c[0].i = i-1;
877: v[0] = - hyhzdhx*(dw - gw);
878: c[1].k = k; c[1].j = j; c[1].i = i;
879: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
880: c[2].k = k; c[2].j = j+1; c[2].i = i;
881: v[2] = - hzhxdhy*(dn + gn);
882: c[3].k = k+1; c[3].j = j; c[3].i = i;
883: v[3] = - hxhydhz*(du + gu);
884: ierr = MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
886: /* right-hand bottom interior edge */
887: } else if (k < mz-1) {
889: tu = x[k+1][j][i];
890: au = 0.5*(t0 + tu);
891: bu = pow(au,bm1);
892: /* du = bu * au; */
893: du = pow(au,beta);
894: gu = coef*bu*(tu - t0);
895:
896: td = x[k-1][j][i];
897: ad = 0.5*(t0 + td);
898: bd = pow(ad,bm1);
899: /* dd = bd * ad; */
900: dd = pow(ad,beta);
901: gd = coef*bd*(td - t0);
902:
903: c[0].k = k-1; c[0].j = j; c[0].i = i;
904: v[0] = - hxhydhz*(dd - gd);
905: c[1].k = k; c[1].j = j; c[1].i = i-1;
906: v[1] = - hyhzdhx*(dw - gw);
907: c[2].k = k; c[2].j = j; c[2].i = i;
908: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
909: c[3].k = k; c[3].j = j+1; c[3].i = i;
910: v[3] = - hzhxdhy*(dn + gn);
911: c[4].k = k+1; c[4].j = j; c[4].i = i;
912: v[4] = - hxhydhz*(du + gu);
913: ierr = MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
915: /* right-hand bottom up corner */
916: } else {
918: td = x[k-1][j][i];
919: ad = 0.5*(t0 + td);
920: bd = pow(ad,bm1);
921: /* dd = bd * ad; */
922: dd = pow(ad,beta);
923: gd = coef*bd*(td - t0);
924:
925: c[0].k = k-1; c[0].j = j; c[0].i = i;
926: v[0] = - hxhydhz*(dd - gd);
927: c[1].k = k; c[1].j = j; c[1].i = i-1;
928: v[1] = - hyhzdhx*(dw - gw);
929: c[2].k = k; c[2].j = j; c[2].i = i;
930: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
931: c[3].k = k; c[3].j = j+1; c[3].i = i;
932: v[3] = - hzhxdhy*(dn + gn);
933: ierr = MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
934: }
936: /* right-hand top edge */
937: } else if (j == my-1) {
939: ts = x[k][j-1][i];
940: as = 0.5*(t0 + ts);
941: bs = pow(as,bm1);
942: /* ds = bs * as; */
943: ds = pow(as,beta);
944: gs = coef*bs*(ts - t0);
945:
946: /* right-hand top down corner */
947: if (k == 0) {
949: tu = x[k+1][j][i];
950: au = 0.5*(t0 + tu);
951: bu = pow(au,bm1);
952: /* du = bu * au; */
953: du = pow(au,beta);
954: gu = coef*bu*(tu - t0);
955:
956: c[0].k = k; c[0].j = j-1; c[0].i = i;
957: v[0] = - hzhxdhy*(ds - gs);
958: c[1].k = k; c[1].j = j; c[1].i = i-1;
959: v[1] = - hyhzdhx*(dw - gw);
960: c[2].k = k; c[2].j = j; c[2].i = i;
961: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
962: c[3].k = k+1; c[3].j = j; c[3].i = i;
963: v[3] = - hxhydhz*(du + gu);
964: ierr = MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
966: /* right-hand top interior edge */
967: } else if (k < mz-1) {
969: tu = x[k+1][j][i];
970: au = 0.5*(t0 + tu);
971: bu = pow(au,bm1);
972: /* du = bu * au; */
973: du = pow(au,beta);
974: gu = coef*bu*(tu - t0);
975:
976: td = x[k-1][j][i];
977: ad = 0.5*(t0 + td);
978: bd = pow(ad,bm1);
979: /* dd = bd * ad; */
980: dd = pow(ad,beta);
981: gd = coef*bd*(td - t0);
982:
983: c[0].k = k-1; c[0].j = j; c[0].i = i;
984: v[0] = - hxhydhz*(dd - gd);
985: c[1].k = k; c[1].j = j-1; c[1].i = i;
986: v[1] = - hzhxdhy*(ds - gs);
987: c[2].k = k; c[2].j = j; c[2].i = i-1;
988: v[2] = - hyhzdhx*(dw - gw);
989: c[3].k = k; c[3].j = j; c[3].i = i;
990: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
991: c[4].k = k+1; c[4].j = j; c[4].i = i;
992: v[4] = - hxhydhz*(du + gu);
993: ierr = MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
995: /* right-hand top up corner */
996: } else {
998: td = x[k-1][j][i];
999: ad = 0.5*(t0 + td);
1000: bd = pow(ad,bm1);
1001: /* dd = bd * ad; */
1002: dd = pow(ad,beta);
1003: gd = coef*bd*(td - t0);
1004:
1005: c[0].k = k-1; c[0].j = j; c[0].i = i;
1006: v[0] = - hxhydhz*(dd - gd);
1007: c[1].k = k; c[1].j = j-1; c[1].i = i;
1008: v[1] = - hzhxdhy*(ds - gs);
1009: c[2].k = k; c[2].j = j; c[2].i = i-1;
1010: v[2] = - hyhzdhx*(dw - gw);
1011: c[3].k = k; c[3].j = j; c[3].i = i;
1012: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1013: ierr = MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
1014: }
1016: } else {
1018: ts = x[k][j-1][i];
1019: as = 0.5*(t0 + ts);
1020: bs = pow(as,bm1);
1021: /* ds = bs * as; */
1022: ds = pow(as,beta);
1023: gs = coef*bs*(t0 - ts);
1024:
1025: tn = x[k][j+1][i];
1026: an = 0.5*(t0 + tn);
1027: bn = pow(an,bm1);
1028: /* dn = bn * an; */
1029: dn = pow(an,beta);
1030: gn = coef*bn*(tn - t0);
1032: /* right-hand down interior edge */
1033: if (k == 0) {
1035: tu = x[k+1][j][i];
1036: au = 0.5*(t0 + tu);
1037: bu = pow(au,bm1);
1038: /* du = bu * au; */
1039: du = pow(au,beta);
1040: gu = coef*bu*(tu - t0);
1042: c[0].k = k; c[0].j = j-1; c[0].i = i;
1043: v[0] = - hzhxdhy*(ds - gs);
1044: c[1].k = k; c[1].j = j; c[1].i = i-1;
1045: v[1] = - hyhzdhx*(dw - gw);
1046: c[2].k = k; c[2].j = j; c[2].i = i;
1047: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1048: c[3].k = k; c[3].j = j+1; c[3].i = i;
1049: v[3] = - hzhxdhy*(dn + gn);
1050: c[4].k = k+1; c[4].j = j; c[4].i = i;
1051: v[4] = - hxhydhz*(du + gu);
1052: ierr = MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1053: }
1055: /* right-hand up interior edge */
1056: else if (k == mz-1) {
1058: td = x[k-1][j][i];
1059: ad = 0.5*(t0 + td);
1060: bd = pow(ad,bm1);
1061: /* dd = bd * ad; */
1062: dd = pow(ad,beta);
1063: gd = coef*bd*(t0 - td);
1064:
1065: c[0].k = k-1; c[0].j = j; c[0].i = i;
1066: v[0] = - hxhydhz*(dd - gd);
1067: c[1].k = k; c[1].j = j-1; c[1].i = i;
1068: v[1] = - hzhxdhy*(ds - gs);
1069: c[2].k = k; c[2].j = j; c[2].i = i-1;
1070: v[2] = - hyhzdhx*(dw - gw);
1071: c[3].k = k; c[3].j = j; c[3].i = i;
1072: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1073: c[4].k = k; c[4].j = j+1; c[4].i = i;
1074: v[4] = - hzhxdhy*(dn + gn);
1075: ierr = MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1076: }
1078: /* right-hand interior plane */
1079: else {
1081: td = x[k-1][j][i];
1082: ad = 0.5*(t0 + td);
1083: bd = pow(ad,bm1);
1084: /* dd = bd * ad; */
1085: dd = pow(ad,beta);
1086: gd = coef*bd*(t0 - td);
1087:
1088: tu = x[k+1][j][i];
1089: au = 0.5*(t0 + tu);
1090: bu = pow(au,bm1);
1091: /* du = bu * au; */
1092: du = pow(au,beta);
1093: gu = coef*bu*(tu - t0);
1095: c[0].k = k-1; c[0].j = j; c[0].i = i;
1096: v[0] = - hxhydhz*(dd - gd);
1097: c[1].k = k; c[1].j = j-1; c[1].i = i;
1098: v[1] = - hzhxdhy*(ds - gs);
1099: c[2].k = k; c[2].j = j; c[2].i = i-1;
1100: v[2] = - hyhzdhx*(dw - gw);
1101: c[3].k = k; c[3].j = j; c[3].i = i;
1102: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1103: c[4].k = k; c[4].j = j+1; c[4].i = i;
1104: v[4] = - hzhxdhy*(dn + gn);
1105: c[5].k = k+1; c[5].j = j; c[5].i = i;
1106: v[5] = - hxhydhz*(du + gu);
1107: ierr = MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1108: }
1109: }
1111: } else if (j == 0) {
1113: tw = x[k][j][i-1];
1114: aw = 0.5*(t0 + tw);
1115: bw = pow(aw,bm1);
1116: /* dw = bw * aw */
1117: dw = pow(aw,beta);
1118: gw = coef*bw*(t0 - tw);
1120: te = x[k][j][i+1];
1121: ae = 0.5*(t0 + te);
1122: be = pow(ae,bm1);
1123: /* de = be * ae; */
1124: de = pow(ae,beta);
1125: ge = coef*be*(te - t0);
1127: tn = x[k][j+1][i];
1128: an = 0.5*(t0 + tn);
1129: bn = pow(an,bm1);
1130: /* dn = bn * an; */
1131: dn = pow(an,beta);
1132: gn = coef*bn*(tn - t0);
1135: /* bottom down interior edge */
1136: if (k == 0) {
1138: tu = x[k+1][j][i];
1139: au = 0.5*(t0 + tu);
1140: bu = pow(au,bm1);
1141: /* du = bu * au; */
1142: du = pow(au,beta);
1143: gu = coef*bu*(tu - t0);
1145: c[0].k = k; c[0].j = j; c[0].i = i-1;
1146: v[0] = - hyhzdhx*(dw - gw);
1147: c[1].k = k; c[1].j = j; c[1].i = i;
1148: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1149: c[2].k = k; c[2].j = j; c[2].i = i+1;
1150: v[2] = - hyhzdhx*(de + ge);
1151: c[3].k = k; c[3].j = j+1; c[3].i = i;
1152: v[3] = - hzhxdhy*(dn + gn);
1153: c[4].k = k+1; c[4].j = j; c[4].i = i;
1154: v[4] = - hxhydhz*(du + gu);
1155: ierr = MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1156: }
1158: /* bottom up interior edge */
1159: else if (k == mz-1) {
1161: td = x[k-1][j][i];
1162: ad = 0.5*(t0 + td);
1163: bd = pow(ad,bm1);
1164: /* dd = bd * ad; */
1165: dd = pow(ad,beta);
1166: gd = coef*bd*(td - t0);
1168: c[0].k = k-1; c[0].j = j; c[0].i = i;
1169: v[0] = - hxhydhz*(dd - gd);
1170: c[1].k = k; c[1].j = j; c[1].i = i-1;
1171: v[1] = - hyhzdhx*(dw - gw);
1172: c[2].k = k; c[2].j = j; c[2].i = i;
1173: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1174: c[3].k = k; c[3].j = j; c[3].i = i+1;
1175: v[3] = - hyhzdhx*(de + ge);
1176: c[4].k = k; c[4].j = j+1; c[4].i = i;
1177: v[4] = - hzhxdhy*(dn + gn);
1178: ierr = MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1179: }
1181: /* bottom interior plane */
1182: else {
1184: tu = x[k+1][j][i];
1185: au = 0.5*(t0 + tu);
1186: bu = pow(au,bm1);
1187: /* du = bu * au; */
1188: du = pow(au,beta);
1189: gu = coef*bu*(tu - t0);
1191: td = x[k-1][j][i];
1192: ad = 0.5*(t0 + td);
1193: bd = pow(ad,bm1);
1194: /* dd = bd * ad; */
1195: dd = pow(ad,beta);
1196: gd = coef*bd*(td - t0);
1198: c[0].k = k-1; c[0].j = j; c[0].i = i;
1199: v[0] = - hxhydhz*(dd - gd);
1200: c[1].k = k; c[1].j = j; c[1].i = i-1;
1201: v[1] = - hyhzdhx*(dw - gw);
1202: c[2].k = k; c[2].j = j; c[2].i = i;
1203: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1204: c[3].k = k; c[3].j = j; c[3].i = i+1;
1205: v[3] = - hyhzdhx*(de + ge);
1206: c[4].k = k; c[4].j = j+1; c[4].i = i;
1207: v[4] = - hzhxdhy*(dn + gn);
1208: c[5].k = k+1; c[5].j = j; c[5].i = i;
1209: v[5] = - hxhydhz*(du + gu);
1210: ierr = MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1211: }
1213: } else if (j == my-1) {
1215: tw = x[k][j][i-1];
1216: aw = 0.5*(t0 + tw);
1217: bw = pow(aw,bm1);
1218: /* dw = bw * aw */
1219: dw = pow(aw,beta);
1220: gw = coef*bw*(t0 - tw);
1222: te = x[k][j][i+1];
1223: ae = 0.5*(t0 + te);
1224: be = pow(ae,bm1);
1225: /* de = be * ae; */
1226: de = pow(ae,beta);
1227: ge = coef*be*(te - t0);
1229: ts = x[k][j-1][i];
1230: as = 0.5*(t0 + ts);
1231: bs = pow(as,bm1);
1232: /* ds = bs * as; */
1233: ds = pow(as,beta);
1234: gs = coef*bs*(t0 - ts);
1235:
1236: /* top down interior edge */
1237: if (k == 0) {
1239: tu = x[k+1][j][i];
1240: au = 0.5*(t0 + tu);
1241: bu = pow(au,bm1);
1242: /* du = bu * au; */
1243: du = pow(au,beta);
1244: gu = coef*bu*(tu - t0);
1246: c[0].k = k; c[0].j = j-1; c[0].i = i;
1247: v[0] = - hzhxdhy*(ds - gs);
1248: c[1].k = k; c[1].j = j; c[1].i = i-1;
1249: v[1] = - hyhzdhx*(dw - gw);
1250: c[2].k = k; c[2].j = j; c[2].i = i;
1251: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1252: c[3].k = k; c[3].j = j; c[3].i = i+1;
1253: v[3] = - hyhzdhx*(de + ge);
1254: c[4].k = k+1; c[4].j = j; c[4].i = i;
1255: v[4] = - hxhydhz*(du + gu);
1256: ierr = MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1257: }
1259: /* top up interior edge */
1260: else if (k == mz-1) {
1262: td = x[k-1][j][i];
1263: ad = 0.5*(t0 + td);
1264: bd = pow(ad,bm1);
1265: /* dd = bd * ad; */
1266: dd = pow(ad,beta);
1267: gd = coef*bd*(td - t0);
1269: c[0].k = k-1; c[0].j = j; c[0].i = i;
1270: v[0] = - hxhydhz*(dd - gd);
1271: c[1].k = k; c[1].j = j-1; c[1].i = i;
1272: v[1] = - hzhxdhy*(ds - gs);
1273: c[2].k = k; c[2].j = j; c[2].i = i-1;
1274: v[2] = - hyhzdhx*(dw - gw);
1275: c[3].k = k; c[3].j = j; c[3].i = i;
1276: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1277: c[4].k = k; c[4].j = j; c[4].i = i+1;
1278: v[4] = - hyhzdhx*(de + ge);
1279: ierr = MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1280: }
1282: /* top interior plane */
1283: else {
1285: tu = x[k+1][j][i];
1286: au = 0.5*(t0 + tu);
1287: bu = pow(au,bm1);
1288: /* du = bu * au; */
1289: du = pow(au,beta);
1290: gu = coef*bu*(tu - t0);
1292: td = x[k-1][j][i];
1293: ad = 0.5*(t0 + td);
1294: bd = pow(ad,bm1);
1295: /* dd = bd * ad; */
1296: dd = pow(ad,beta);
1297: gd = coef*bd*(td - t0);
1299: c[0].k = k-1; c[0].j = j; c[0].i = i;
1300: v[0] = - hxhydhz*(dd - gd);
1301: c[1].k = k; c[1].j = j-1; c[1].i = i;
1302: v[1] = - hzhxdhy*(ds - gs);
1303: c[2].k = k; c[2].j = j; c[2].i = i-1;
1304: v[2] = - hyhzdhx*(dw - gw);
1305: c[3].k = k; c[3].j = j; c[3].i = i;
1306: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1307: c[4].k = k; c[4].j = j; c[4].i = i+1;
1308: v[4] = - hyhzdhx*(de + ge);
1309: c[5].k = k+1; c[5].j = j; c[5].i = i;
1310: v[5] = - hxhydhz*(du + gu);
1311: ierr = MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1312: }
1314: } else if (k == 0) {
1316: /* down interior plane */
1318: tw = x[k][j][i-1];
1319: aw = 0.5*(t0 + tw);
1320: bw = pow(aw,bm1);
1321: /* dw = bw * aw */
1322: dw = pow(aw,beta);
1323: gw = coef*bw*(t0 - tw);
1325: te = x[k][j][i+1];
1326: ae = 0.5*(t0 + te);
1327: be = pow(ae,bm1);
1328: /* de = be * ae; */
1329: de = pow(ae,beta);
1330: ge = coef*be*(te - t0);
1332: ts = x[k][j-1][i];
1333: as = 0.5*(t0 + ts);
1334: bs = pow(as,bm1);
1335: /* ds = bs * as; */
1336: ds = pow(as,beta);
1337: gs = coef*bs*(t0 - ts);
1338:
1339: tn = x[k][j+1][i];
1340: an = 0.5*(t0 + tn);
1341: bn = pow(an,bm1);
1342: /* dn = bn * an; */
1343: dn = pow(an,beta);
1344: gn = coef*bn*(tn - t0);
1345:
1346: tu = x[k+1][j][i];
1347: au = 0.5*(t0 + tu);
1348: bu = pow(au,bm1);
1349: /* du = bu * au; */
1350: du = pow(au,beta);
1351: gu = coef*bu*(tu - t0);
1353: c[0].k = k; c[0].j = j-1; c[0].i = i;
1354: v[0] = - hzhxdhy*(ds - gs);
1355: c[1].k = k; c[1].j = j; c[1].i = i-1;
1356: v[1] = - hyhzdhx*(dw - gw);
1357: c[2].k = k; c[2].j = j; c[2].i = i;
1358: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1359: c[3].k = k; c[3].j = j; c[3].i = i+1;
1360: v[3] = - hyhzdhx*(de + ge);
1361: c[4].k = k; c[4].j = j+1; c[4].i = i;
1362: v[4] = - hzhxdhy*(dn + gn);
1363: c[5].k = k+1; c[5].j = j; c[5].i = i;
1364: v[5] = - hxhydhz*(du + gu);
1365: ierr = MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1366: }
1367:
1368: else if (k == mz-1) {
1370: /* up interior plane */
1372: tw = x[k][j][i-1];
1373: aw = 0.5*(t0 + tw);
1374: bw = pow(aw,bm1);
1375: /* dw = bw * aw */
1376: dw = pow(aw,beta);
1377: gw = coef*bw*(t0 - tw);
1379: te = x[k][j][i+1];
1380: ae = 0.5*(t0 + te);
1381: be = pow(ae,bm1);
1382: /* de = be * ae; */
1383: de = pow(ae,beta);
1384: ge = coef*be*(te - t0);
1386: ts = x[k][j-1][i];
1387: as = 0.5*(t0 + ts);
1388: bs = pow(as,bm1);
1389: /* ds = bs * as; */
1390: ds = pow(as,beta);
1391: gs = coef*bs*(t0 - ts);
1392:
1393: tn = x[k][j+1][i];
1394: an = 0.5*(t0 + tn);
1395: bn = pow(an,bm1);
1396: /* dn = bn * an; */
1397: dn = pow(an,beta);
1398: gn = coef*bn*(tn - t0);
1399:
1400: td = x[k-1][j][i];
1401: ad = 0.5*(t0 + td);
1402: bd = pow(ad,bm1);
1403: /* dd = bd * ad; */
1404: dd = pow(ad,beta);
1405: gd = coef*bd*(t0 - td);
1406:
1407: c[0].k = k-1; c[0].j = j; c[0].i = i;
1408: v[0] = - hxhydhz*(dd - gd);
1409: c[1].k = k; c[1].j = j-1; c[1].i = i;
1410: v[1] = - hzhxdhy*(ds - gs);
1411: c[2].k = k; c[2].j = j; c[2].i = i-1;
1412: v[2] = - hyhzdhx*(dw - gw);
1413: c[3].k = k; c[3].j = j; c[3].i = i;
1414: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1415: c[4].k = k; c[4].j = j; c[4].i = i+1;
1416: v[4] = - hyhzdhx*(de + ge);
1417: c[5].k = k; c[5].j = j+1; c[5].i = i;
1418: v[5] = - hzhxdhy*(dn + gn);
1419: ierr = MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1420: }
1421: }
1422: }
1423: }
1424: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
1425: DAVecRestoreArray((DA)dmmg->dm,localX,(void**)&x);
1426: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
1427: DARestoreLocalVector((DA)dmmg->dm,&localX);
1429: PetscLogFlops((41 + 8*POWFLOP)*xm*ym);
1430: return(0);
1431: }