Actual source code: ex20.c

petsc-3.3-p7 2013-05-11
  2: static char help[] ="Nonlinear Radiative Transport PDE with multigrid in 3d.\n\
  3: Uses 3-dimensional distributed arrays.\n\
  4: A 3-dim simplified Radiative Transport test problem is used, with analytic Jacobian. \n\
  5: \n\
  6:   Solves the linear systems via multilevel methods \n\
  7: \n\
  8: The command line\n\
  9: options are:\n\
 10:   -tleft <tl>, where <tl> indicates the left Diriclet BC \n\
 11:   -tright <tr>, where <tr> indicates the right Diriclet BC \n\
 12:   -beta <beta>, where <beta> indicates the exponent in T \n\n";

 14: /*T
 15:    Concepts: SNES^solving a system of nonlinear equations
 16:    Concepts: DM^using distributed arrays
 17:    Concepts: multigrid;
 18:    Processors: n
 19: T*/

 21: /*  
 22:   
 23:     This example models the partial differential equation 
 24:    
 25:          - Div(alpha* T^beta (GRAD T)) = 0.
 26:        
 27:     where beta = 2.5 and alpha = 1.0
 28:  
 29:     BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = dT/dn_up = dT/dn_down = 0.
 30:     
 31:     in the unit square, which is uniformly discretized in each of x and 
 32:     y in this simple encoding.  The degrees of freedom are cell centered.
 33:  
 34:     A finite volume approximation with the usual 7-point stencil 
 35:     is used to discretize the boundary value problem to obtain a 
 36:     nonlinear system of equations. 

 38:     This code was contributed by Nickolas Jovanovic based on ex18.c
 39:  
 40: */

 42: #include <petscsnes.h>
 43: #include <petscdmda.h>
 44: #include <petscpcmg.h>

 46: /* User-defined application context */

 48: typedef struct {
 49:    PetscReal   tleft,tright;  /* Dirichlet boundary conditions */
 50:    PetscReal   beta,bm1,coef; /* nonlinear diffusivity parameterizations */
 51: } AppCtx;

 53: #define POWFLOP 5 /* assume a pow() takes five flops */

 55: extern PetscErrorCode FormInitialGuess(DM,Vec);
 56: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 57: extern PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);

 61: int main(int argc,char **argv)
 62: {
 63:   SNES           snes;
 64:   AppCtx         user;
 66:   PetscInt       its,lits;
 67:   PetscReal      litspit;
 68:   DM             da;

 70:   PetscInitialize(&argc,&argv,PETSC_NULL,help);

 72:   /* set problem parameters */
 73:   user.tleft  = 1.0;
 74:   user.tright = 0.1;
 75:   user.beta   = 2.5;
 76:   PetscOptionsGetReal(PETSC_NULL,"-tleft",&user.tleft,PETSC_NULL);
 77:   PetscOptionsGetReal(PETSC_NULL,"-tright",&user.tright,PETSC_NULL);
 78:   PetscOptionsGetReal(PETSC_NULL,"-beta",&user.beta,PETSC_NULL);
 79:   user.bm1  = user.beta - 1.0;
 80:   user.coef = user.beta/2.0;

 82:   /*
 83:       Set the DMDA (grid structure) for the grids.
 84:   */
 85:   DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-5,-5,-5,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);
 86:   DMSetApplicationContext(da,&user);

 88:   /*
 89:      Create the nonlinear solver
 90:   */
 91:   SNESCreate(PETSC_COMM_WORLD,&snes);
 92:   SNESSetDM(snes,da);
 93:   SNESSetFunction(snes,PETSC_NULL,FormFunction,&user);
 94:   SNESSetJacobian(snes,PETSC_NULL,PETSC_NULL,FormJacobian,&user);
 95:   SNESSetFromOptions(snes);
 96:   DMSetInitialGuess(da,FormInitialGuess);

 98:   SNESSolve(snes,PETSC_NULL,PETSC_NULL);
 99:   SNESGetIterationNumber(snes,&its);
100:   SNESGetLinearSolveIterations(snes,&lits);
101:   litspit = ((PetscReal)lits)/((PetscReal)its);
102:   PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
103:   PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
104:   PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",litspit);

106:   SNESDestroy(&snes);
107:   DMDestroy(&da);
108:   PetscFinalize();

110:   return 0;
111: }
112: /* --------------------  Form initial approximation ----------------- */
115: PetscErrorCode FormInitialGuess(DM da,Vec X)
116: {
117:   AppCtx         *user;
118:   PetscInt       i,j,k,xs,ys,xm,ym,zs,zm;
120:   PetscScalar    ***x;

123:   DMGetApplicationContext(da,&user);
124:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
125:   DMDAVecGetArray(da,X,&x);

127:   /* Compute initial guess */
128:   for (k=zs; k<zs+zm; k++) {
129:     for (j=ys; j<ys+ym; j++) {
130:       for (i=xs; i<xs+xm; i++) {
131:         x[k][j][i] = user->tleft;
132:       }
133:     }
134:   }
135:   DMDAVecRestoreArray(da,X,&x);
136:   return(0);
137: }
138: /* --------------------  Evaluate Function F(x) --------------------- */
141: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ptr)
142: {
143:   AppCtx         *user = (AppCtx*)ptr;
145:   PetscInt       i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
146:   PetscScalar    zero = 0.0,one = 1.0;
147:   PetscScalar    hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
148:   PetscScalar    t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw,fn = 0.0,fs = 0.0,fe =0.0,fw = 0.0;
149:   PetscScalar    tleft,tright,beta,td,ad,dd,fd=0.0,tu,au,du=0.0,fu=0.0;
150:   PetscScalar    ***x,***f;
151:   Vec            localX;
152:   DM             da;

155:   SNESGetDM(snes,&da);
156:   DMGetLocalVector(da,&localX);
157:   DMDAGetInfo(da,PETSC_NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
158:   hx    = one/(PetscReal)(mx-1);  hy    = one/(PetscReal)(my-1);  hz = one/(PetscReal)(mz-1);
159:   hxhydhz = hx*hy/hz;   hyhzdhx = hy*hz/hx;   hzhxdhy = hz*hx/hy;
160:   tleft = user->tleft;         tright = user->tright;
161:   beta  = user->beta;
162: 
163:   /* Get ghost points */
164:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
165:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
166:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
167:   DMDAVecGetArray(da,localX,&x);
168:   DMDAVecGetArray(da,F,&f);

170:   /* Evaluate function */
171:   for (k=zs; k<zs+zm; k++) {
172:     for (j=ys; j<ys+ym; j++) {
173:       for (i=xs; i<xs+xm; i++) {
174:         t0 = x[k][j][i];

176:         if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {

178:             /* general interior volume */

180:           tw = x[k][j][i-1];
181:           aw = 0.5*(t0 + tw);
182:           dw = pow(aw,beta);
183:           fw = dw*(t0 - tw);

185:                  te = x[k][j][i+1];
186:           ae = 0.5*(t0 + te);
187:           de = pow(ae,beta);
188:           fe = de*(te - t0);

190:           ts = x[k][j-1][i];
191:           as = 0.5*(t0 + ts);
192:           ds = pow(as,beta);
193:           fs = ds*(t0 - ts);
194: 
195:           tn = x[k][j+1][i];
196:           an = 0.5*(t0 + tn);
197:           dn = pow(an,beta);
198:           fn = dn*(tn - t0);

200:           td = x[k-1][j][i];
201:           ad = 0.5*(t0 + td);
202:           dd = pow(ad,beta);
203:           fd = dd*(t0 - td);

205:           tu = x[k+1][j][i];
206:           au = 0.5*(t0 + tu);
207:           du = pow(au,beta);
208:           fu = du*(tu - t0);

210:         } else if (i == 0) {

212:            /* left-hand (west) boundary */
213:           tw = tleft;
214:           aw = 0.5*(t0 + tw);
215:           dw = pow(aw,beta);
216:           fw = dw*(t0 - tw);

218:           te = x[k][j][i+1];
219:           ae = 0.5*(t0 + te);
220:           de = pow(ae,beta);
221:           fe = de*(te - t0);

223:           if (j > 0) {
224:             ts = x[k][j-1][i];
225:             as = 0.5*(t0 + ts);
226:             ds = pow(as,beta);
227:             fs = ds*(t0 - ts);
228:           } else {
229:              fs = zero;
230:           }

232:           if (j < my-1) {
233:             tn = x[k][j+1][i];
234:             an = 0.5*(t0 + tn);
235:             dn = pow(an,beta);
236:             fn = dn*(tn - t0);
237:           } else {
238:             fn = zero;
239:              }

241:           if (k > 0) {
242:             td = x[k-1][j][i];
243:             ad = 0.5*(t0 + td);
244:             dd = pow(ad,beta);
245:             fd = dd*(t0 - td);
246:           } else {
247:              fd = zero;
248:           }

250:           if (k < mz-1) {
251:             tu = x[k+1][j][i];
252:             au = 0.5*(t0 + tu);
253:             du = pow(au,beta);
254:             fu = du*(tu - t0);
255:           } else {
256:              fu = zero;
257:           }

259:         } else if (i == mx-1) {

261:           /* right-hand (east) boundary */
262:           tw = x[k][j][i-1];
263:           aw = 0.5*(t0 + tw);
264:           dw = pow(aw,beta);
265:           fw = dw*(t0 - tw);
266: 
267:           te = tright;
268:           ae = 0.5*(t0 + te);
269:           de = pow(ae,beta);
270:           fe = de*(te - t0);
271: 
272:           if (j > 0) {
273:             ts = x[k][j-1][i];
274:             as = 0.5*(t0 + ts);
275:             ds = pow(as,beta);
276:             fs = ds*(t0 - ts);
277:           } else {
278:             fs = zero;
279:           }
280: 
281:           if (j < my-1) {
282:             tn = x[k][j+1][i];
283:             an = 0.5*(t0 + tn);
284:             dn = pow(an,beta);
285:             fn = dn*(tn - t0);
286:           } else {
287:             fn = zero;
288:           }

290:           if (k > 0) {
291:             td = x[k-1][j][i];
292:             ad = 0.5*(t0 + td);
293:             dd = pow(ad,beta);
294:             fd = dd*(t0 - td);
295:           } else {
296:              fd = zero;
297:           }

299:           if (k < mz-1) {
300:             tu = x[k+1][j][i];
301:             au = 0.5*(t0 + tu);
302:             du = pow(au,beta);
303:             fu = du*(tu - t0);
304:           } else {
305:              fu = zero;
306:           }

308:         } else if (j == 0) {

310:           /* bottom (south) boundary, and i <> 0 or mx-1 */
311:           tw = x[k][j][i-1];
312:           aw = 0.5*(t0 + tw);
313:           dw = pow(aw,beta);
314:           fw = dw*(t0 - tw);

316:           te = x[k][j][i+1];
317:           ae = 0.5*(t0 + te);
318:           de = pow(ae,beta);
319:           fe = de*(te - t0);

321:           fs = zero;

323:           tn = x[k][j+1][i];
324:           an = 0.5*(t0 + tn);
325:           dn = pow(an,beta);
326:           fn = dn*(tn - t0);

328:           if (k > 0) {
329:             td = x[k-1][j][i];
330:             ad = 0.5*(t0 + td);
331:             dd = pow(ad,beta);
332:             fd = dd*(t0 - td);
333:           } else {
334:              fd = zero;
335:           }

337:           if (k < mz-1) {
338:             tu = x[k+1][j][i];
339:             au = 0.5*(t0 + tu);
340:             du = pow(au,beta);
341:             fu = du*(tu - t0);
342:           } else {
343:              fu = zero;
344:           }

346:         } else if (j == my-1) {

348:           /* top (north) boundary, and i <> 0 or mx-1 */
349:           tw = x[k][j][i-1];
350:           aw = 0.5*(t0 + tw);
351:           dw = pow(aw,beta);
352:           fw = dw*(t0 - tw);

354:           te = x[k][j][i+1];
355:           ae = 0.5*(t0 + te);
356:           de = pow(ae,beta);
357:           fe = de*(te - t0);

359:           ts = x[k][j-1][i];
360:           as = 0.5*(t0 + ts);
361:           ds = pow(as,beta);
362:           fs = ds*(t0 - ts);

364:           fn = zero;

366:           if (k > 0) {
367:             td = x[k-1][j][i];
368:             ad = 0.5*(t0 + td);
369:             dd = pow(ad,beta);
370:             fd = dd*(t0 - td);
371:           } else {
372:              fd = zero;
373:           }

375:           if (k < mz-1) {
376:             tu = x[k+1][j][i];
377:             au = 0.5*(t0 + tu);
378:             du = pow(au,beta);
379:             fu = du*(tu - t0);
380:           } else {
381:              fu = zero;
382:           }

384:         } else if (k == 0) {

386:           /* down boundary (interior only) */
387:           tw = x[k][j][i-1];
388:           aw = 0.5*(t0 + tw);
389:           dw = pow(aw,beta);
390:           fw = dw*(t0 - tw);

392:           te = x[k][j][i+1];
393:           ae = 0.5*(t0 + te);
394:           de = pow(ae,beta);
395:           fe = de*(te - t0);

397:           ts = x[k][j-1][i];
398:           as = 0.5*(t0 + ts);
399:           ds = pow(as,beta);
400:           fs = ds*(t0 - ts);

402:           tn = x[k][j+1][i];
403:           an = 0.5*(t0 + tn);
404:           dn = pow(an,beta);
405:           fn = dn*(tn - t0);

407:            fd = zero;

409:           tu = x[k+1][j][i];
410:           au = 0.5*(t0 + tu);
411:           du = pow(au,beta);
412:           fu = du*(tu - t0);
413: 
414:         } else if (k == mz-1) {

416:           /* up boundary (interior only) */
417:           tw = x[k][j][i-1];
418:           aw = 0.5*(t0 + tw);
419:           dw = pow(aw,beta);
420:           fw = dw*(t0 - tw);

422:           te = x[k][j][i+1];
423:           ae = 0.5*(t0 + te);
424:           de = pow(ae,beta);
425:           fe = de*(te - t0);

427:           ts = x[k][j-1][i];
428:           as = 0.5*(t0 + ts);
429:           ds = pow(as,beta);
430:           fs = ds*(t0 - ts);

432:           tn = x[k][j+1][i];
433:           an = 0.5*(t0 + tn);
434:           dn = pow(an,beta);
435:           fn = dn*(tn - t0);

437:           td = x[k-1][j][i];
438:           ad = 0.5*(t0 + td);
439:           dd = pow(ad,beta);
440:           fd = dd*(t0 - td);

442:           fu = zero;
443:         }

445:         f[k][j][i] = - hyhzdhx*(fe-fw) - hzhxdhy*(fn-fs) - hxhydhz*(fu-fd);
446:       }
447:     }
448:   }
449:   DMDAVecRestoreArray(da,localX,&x);
450:   DMDAVecRestoreArray(da,F,&f);
451:   DMRestoreLocalVector(da,&localX);
452:   PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
453:   return(0);
454: }
455: /* --------------------  Evaluate Jacobian F(x) --------------------- */
458: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ptr)
459: {
460:   AppCtx         *user = (AppCtx*)ptr;
462:   PetscInt       i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
463:   PetscScalar    one = 1.0;
464:   PetscScalar    hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
465:   PetscScalar    t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw;
466:   PetscScalar    tleft,tright,beta,td,ad,dd,tu,au,du,v[7],bm1,coef;
467:   PetscScalar    ***x,bn,bs,be,bw,bu,bd,gn,gs,ge,gw,gu,gd;
468:   Vec            localX;
469:   MatStencil     c[7],row;
470:   Mat            jac = *B;
471:   DM             da;

474:   SNESGetDM(snes,&da);
475:   DMGetLocalVector(da,&localX);
476:   DMDAGetInfo(da,PETSC_NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
477:   hx    = one/(PetscReal)(mx-1);  hy    = one/(PetscReal)(my-1);  hz = one/(PetscReal)(mz-1);
478:   hxhydhz = hx*hy/hz;   hyhzdhx = hy*hz/hx;   hzhxdhy = hz*hx/hy;
479:   tleft = user->tleft;         tright = user->tright;
480:   beta  = user->beta;               bm1    = user->bm1;                coef = user->coef;

482:   /* Get ghost points */
483:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
484:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
485:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
486:   DMDAVecGetArray(da,localX,&x);

488:   /* Evaluate Jacobian of function */
489:   for (k=zs; k<zs+zm; k++) {
490:     for (j=ys; j<ys+ym; j++) {
491:       for (i=xs; i<xs+xm; i++) {
492:         t0    = x[k][j][i];
493:         row.k = k; row.j = j; row.i = i;
494:         if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {

496:           /* general interior volume */

498:           tw = x[k][j][i-1];
499:           aw = 0.5*(t0 + tw);
500:           bw = pow(aw,bm1);
501:           /* dw = bw * aw */
502:           dw = pow(aw,beta);
503:           gw = coef*bw*(t0 - tw);

505:           te = x[k][j][i+1];
506:           ae = 0.5*(t0 + te);
507:           be = pow(ae,bm1);
508:           /* de = be * ae; */
509:           de = pow(ae,beta);
510:           ge = coef*be*(te - t0);

512:           ts = x[k][j-1][i];
513:           as = 0.5*(t0 + ts);
514:           bs = pow(as,bm1);
515:           /* ds = bs * as; */
516:           ds = pow(as,beta);
517:           gs = coef*bs*(t0 - ts);
518: 
519:           tn = x[k][j+1][i];
520:           an = 0.5*(t0 + tn);
521:           bn = pow(an,bm1);
522:           /* dn = bn * an; */
523:           dn = pow(an,beta);
524:           gn = coef*bn*(tn - t0);

526:           td = x[k-1][j][i];
527:           ad = 0.5*(t0 + td);
528:           bd = pow(ad,bm1);
529:           /* dd = bd * ad; */
530:           dd = pow(ad,beta);
531:           gd = coef*bd*(t0 - td);
532: 
533:           tu = x[k+1][j][i];
534:           au = 0.5*(t0 + tu);
535:           bu = pow(au,bm1);
536:           /* du = bu * au; */
537:           du = pow(au,beta);
538:           gu = coef*bu*(tu - t0);

540:           c[0].k = k-1; c[0].j = j; c[0].i = i; v[0]   = - hxhydhz*(dd - gd);
541:           c[1].k = k; c[1].j = j-1; c[1].i = i;
542:           v[1]   = - hzhxdhy*(ds - gs);
543:           c[2].k = k; c[2].j = j; c[2].i = i-1;
544:           v[2]   = - hyhzdhx*(dw - gw);
545:           c[3].k = k; c[3].j = j; c[3].i = i;
546:           v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
547:           c[4].k = k; c[4].j = j; c[4].i = i+1;
548:           v[4]   = - hyhzdhx*(de + ge);
549:           c[5].k = k; c[5].j = j+1; c[5].i = i;
550:           v[5]   = - hzhxdhy*(dn + gn);
551:           c[6].k = k+1; c[6].j = j; c[6].i = i;
552:           v[6]   = - hxhydhz*(du + gu);
553:             MatSetValuesStencil(jac,1,&row,7,c,v,INSERT_VALUES);

555:         } else if (i == 0) {

557:           /* left-hand plane boundary */
558:           tw = tleft;
559:           aw = 0.5*(t0 + tw);
560:           bw = pow(aw,bm1);
561:           /* dw = bw * aw */
562:           dw = pow(aw,beta);
563:           gw = coef*bw*(t0 - tw);
564: 
565:           te = x[k][j][i+1];
566:           ae = 0.5*(t0 + te);
567:           be = pow(ae,bm1);
568:           /* de = be * ae; */
569:           de = pow(ae,beta);
570:           ge = coef*be*(te - t0);
571: 
572:           /* left-hand bottom edge */
573:           if (j == 0) {

575:             tn = x[k][j+1][i];
576:             an = 0.5*(t0 + tn);
577:             bn = pow(an,bm1);
578:             /* dn = bn * an; */
579:             dn = pow(an,beta);
580:             gn = coef*bn*(tn - t0);
581: 
582:             /* left-hand bottom down corner */
583:             if (k == 0) {

585:               tu = x[k+1][j][i];
586:               au = 0.5*(t0 + tu);
587:               bu = pow(au,bm1);
588:               /* du = bu * au; */
589:               du = pow(au,beta);
590:               gu = coef*bu*(tu - t0);
591: 
592:               c[0].k = k; c[0].j = j; c[0].i = i;
593:               v[0]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
594:               c[1].k = k; c[1].j = j; c[1].i = i+1;
595:               v[1]   = - hyhzdhx*(de + ge);
596:               c[2].k = k; c[2].j = j+1; c[2].i = i;
597:               v[2]   = - hzhxdhy*(dn + gn);
598:               c[3].k = k+1; c[3].j = j; c[3].i = i;
599:               v[3]   = - hxhydhz*(du + gu);
600:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

602:             /* left-hand bottom interior edge */
603:             } else if (k < mz-1) {

605:               tu = x[k+1][j][i];
606:               au = 0.5*(t0 + tu);
607:               bu = pow(au,bm1);
608:               /* du = bu * au; */
609:               du = pow(au,beta);
610:               gu = coef*bu*(tu - t0);
611: 
612:               td = x[k-1][j][i];
613:               ad = 0.5*(t0 + td);
614:               bd = pow(ad,bm1);
615:               /* dd = bd * ad; */
616:               dd = pow(ad,beta);
617:               gd = coef*bd*(td - t0);
618: 
619:               c[0].k = k-1; c[0].j = j; c[0].i = i;
620:               v[0]   = - hxhydhz*(dd - gd);
621:               c[1].k = k; c[1].j = j; c[1].i = i;
622:               v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
623:               c[2].k = k; c[2].j = j; c[2].i = i+1;
624:               v[2]   = - hyhzdhx*(de + ge);
625:               c[3].k = k; c[3].j = j+1; c[3].i = i;
626:               v[3]   = - hzhxdhy*(dn + gn);
627:               c[4].k = k+1; c[4].j = j; c[4].i = i;
628:               v[4]   = - hxhydhz*(du + gu);
629:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

631:             /* left-hand bottom up corner */
632:             } else {

634:               td = x[k-1][j][i];
635:               ad = 0.5*(t0 + td);
636:               bd = pow(ad,bm1);
637:               /* dd = bd * ad; */
638:               dd = pow(ad,beta);
639:               gd = coef*bd*(td - t0);
640: 
641:               c[0].k = k-1; c[0].j = j; c[0].i = i;
642:               v[0]   = - hxhydhz*(dd - gd);
643:               c[1].k = k; c[1].j = j; c[1].i = i;
644:               v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
645:               c[2].k = k; c[2].j = j; c[2].i = i+1;
646:               v[2]   = - hyhzdhx*(de + ge);
647:               c[3].k = k; c[3].j = j+1; c[3].i = i;
648:               v[3]   = - hzhxdhy*(dn + gn);
649:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
650:             }

652:           /* left-hand top edge */
653:           } else if (j == my-1) {

655:             ts = x[k][j-1][i];
656:             as = 0.5*(t0 + ts);
657:             bs = pow(as,bm1);
658:             /* ds = bs * as; */
659:             ds = pow(as,beta);
660:             gs = coef*bs*(ts - t0);
661: 
662:             /* left-hand top down corner */
663:             if (k == 0) {

665:               tu = x[k+1][j][i];
666:               au = 0.5*(t0 + tu);
667:               bu = pow(au,bm1);
668:               /* du = bu * au; */
669:               du = pow(au,beta);
670:               gu = coef*bu*(tu - t0);
671: 
672:               c[0].k = k; c[0].j = j-1; c[0].i = i;
673:               v[0]   = - hzhxdhy*(ds - gs);
674:               c[1].k = k; c[1].j = j; c[1].i = i;
675:               v[1]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
676:               c[2].k = k; c[2].j = j; c[2].i = i+1;
677:               v[2]   = - hyhzdhx*(de + ge);
678:               c[3].k = k+1; c[3].j = j; c[3].i = i;
679:               v[3]   = - hxhydhz*(du + gu);
680:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

682:             /* left-hand top interior edge */
683:             } else if (k < mz-1) {

685:               tu = x[k+1][j][i];
686:               au = 0.5*(t0 + tu);
687:               bu = pow(au,bm1);
688:               /* du = bu * au; */
689:               du = pow(au,beta);
690:               gu = coef*bu*(tu - t0);
691: 
692:               td = x[k-1][j][i];
693:               ad = 0.5*(t0 + td);
694:               bd = pow(ad,bm1);
695:               /* dd = bd * ad; */
696:               dd = pow(ad,beta);
697:               gd = coef*bd*(td - t0);
698: 
699:               c[0].k = k-1; c[0].j = j; c[0].i = i;
700:               v[0]   = - hxhydhz*(dd - gd);
701:               c[1].k = k; c[1].j = j-1; c[1].i = i;
702:               v[1]   = - hzhxdhy*(ds - gs);
703:               c[2].k = k; c[2].j = j; c[2].i = i;
704:               v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
705:               c[3].k = k; c[3].j = j; c[3].i = i+1;
706:               v[3]   = - hyhzdhx*(de + ge);
707:               c[4].k = k+1; c[4].j = j; c[4].i = i;
708:               v[4]   = - hxhydhz*(du + gu);
709:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

711:             /* left-hand top up corner */
712:             } else {

714:               td = x[k-1][j][i];
715:               ad = 0.5*(t0 + td);
716:               bd = pow(ad,bm1);
717:               /* dd = bd * ad; */
718:               dd = pow(ad,beta);
719:               gd = coef*bd*(td - t0);
720: 
721:               c[0].k = k-1; c[0].j = j; c[0].i = i;
722:               v[0]   = - hxhydhz*(dd - gd);
723:               c[1].k = k; c[1].j = j-1; c[1].i = i;
724:               v[1]   = - hzhxdhy*(ds - gs);
725:               c[2].k = k; c[2].j = j; c[2].i = i;
726:               v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
727:               c[3].k = k; c[3].j = j; c[3].i = i+1;
728:               v[3]   = - hyhzdhx*(de + ge);
729:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
730:             }

732:           } else {

734:             ts = x[k][j-1][i];
735:             as = 0.5*(t0 + ts);
736:             bs = pow(as,bm1);
737:             /* ds = bs * as; */
738:             ds = pow(as,beta);
739:             gs = coef*bs*(t0 - ts);
740: 
741:             tn = x[k][j+1][i];
742:             an = 0.5*(t0 + tn);
743:             bn = pow(an,bm1);
744:             /* dn = bn * an; */
745:             dn = pow(an,beta);
746:             gn = coef*bn*(tn - t0);

748:             /* left-hand down interior edge */
749:             if (k == 0) {

751:               tu = x[k+1][j][i];
752:               au = 0.5*(t0 + tu);
753:               bu = pow(au,bm1);
754:               /* du = bu * au; */
755:               du = pow(au,beta);
756:               gu = coef*bu*(tu - t0);

758:               c[0].k = k; c[0].j = j-1; c[0].i = i;
759:               v[0]   = - hzhxdhy*(ds - gs);
760:               c[1].k = k; c[1].j = j; c[1].i = i;
761:               v[1]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
762:               c[2].k = k; c[2].j = j; c[2].i = i+1;
763:               v[2]   = - hyhzdhx*(de + ge);
764:               c[3].k = k; c[3].j = j+1; c[3].i = i;
765:               v[3]   = - hzhxdhy*(dn + gn);
766:               c[4].k = k+1; c[4].j = j; c[4].i = i;
767:               v[4]   = - hxhydhz*(du + gu);
768:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
769:             }

771:             /* left-hand up interior edge */
772:             else if (k == mz-1) {

774:               td = x[k-1][j][i];
775:               ad = 0.5*(t0 + td);
776:               bd = pow(ad,bm1);
777:               /* dd = bd * ad; */
778:               dd = pow(ad,beta);
779:               gd = coef*bd*(t0 - td);
780: 
781:               c[0].k = k-1; c[0].j = j; c[0].i = i;
782:               v[0]   = - hxhydhz*(dd - gd);
783:               c[1].k = k; c[1].j = j-1; c[1].i = i;
784:               v[1]   = - hzhxdhy*(ds - gs);
785:               c[2].k = k; c[2].j = j; c[2].i = i;
786:               v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
787:               c[3].k = k; c[3].j = j; c[3].i = i+1;
788:               v[3]   = - hyhzdhx*(de + ge);
789:               c[4].k = k; c[4].j = j+1; c[4].i = i;
790:               v[4]   = - hzhxdhy*(dn + gn);
791:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
792:             }

794:             /* left-hand interior plane */
795:             else {

797:               td = x[k-1][j][i];
798:               ad = 0.5*(t0 + td);
799:               bd = pow(ad,bm1);
800:               /* dd = bd * ad; */
801:               dd = pow(ad,beta);
802:               gd = coef*bd*(t0 - td);
803: 
804:               tu = x[k+1][j][i];
805:               au = 0.5*(t0 + tu);
806:               bu = pow(au,bm1);
807:               /* du = bu * au; */
808:               du = pow(au,beta);
809:               gu = coef*bu*(tu - t0);

811:               c[0].k = k-1; c[0].j = j; c[0].i = i;
812:               v[0]   = - hxhydhz*(dd - gd);
813:               c[1].k = k; c[1].j = j-1; c[1].i = i;
814:               v[1]   = - hzhxdhy*(ds - gs);
815:               c[2].k = k; c[2].j = j; c[2].i = i;
816:               v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
817:               c[3].k = k; c[3].j = j; c[3].i = i+1;
818:               v[3]   = - hyhzdhx*(de + ge);
819:               c[4].k = k; c[4].j = j+1; c[4].i = i;
820:               v[4]   = - hzhxdhy*(dn + gn);
821:               c[5].k = k+1; c[5].j = j; c[5].i = i;
822:               v[5]   = - hxhydhz*(du + gu);
823:                 MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
824:             }
825:           }

827:         } else if (i == mx-1) {

829:           /* right-hand plane boundary */
830:           tw = x[k][j][i-1];
831:           aw = 0.5*(t0 + tw);
832:           bw = pow(aw,bm1);
833:           /* dw = bw * aw */
834:           dw = pow(aw,beta);
835:           gw = coef*bw*(t0 - tw);
836: 
837:           te = tright;
838:           ae = 0.5*(t0 + te);
839:           be = pow(ae,bm1);
840:           /* de = be * ae; */
841:           de = pow(ae,beta);
842:           ge = coef*be*(te - t0);
843: 
844:           /* right-hand bottom edge */
845:           if (j == 0) {

847:             tn = x[k][j+1][i];
848:             an = 0.5*(t0 + tn);
849:             bn = pow(an,bm1);
850:             /* dn = bn * an; */
851:             dn = pow(an,beta);
852:             gn = coef*bn*(tn - t0);
853: 
854:             /* right-hand bottom down corner */
855:             if (k == 0) {

857:               tu = x[k+1][j][i];
858:               au = 0.5*(t0 + tu);
859:               bu = pow(au,bm1);
860:               /* du = bu * au; */
861:               du = pow(au,beta);
862:               gu = coef*bu*(tu - t0);
863: 
864:               c[0].k = k; c[0].j = j; c[0].i = i-1;
865:               v[0]   = - hyhzdhx*(dw - gw);
866:               c[1].k = k; c[1].j = j; c[1].i = i;
867:               v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
868:               c[2].k = k; c[2].j = j+1; c[2].i = i;
869:               v[2]   = - hzhxdhy*(dn + gn);
870:               c[3].k = k+1; c[3].j = j; c[3].i = i;
871:               v[3]   = - hxhydhz*(du + gu);
872:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

874:             /* right-hand bottom interior edge */
875:             } else if (k < mz-1) {

877:               tu = x[k+1][j][i];
878:               au = 0.5*(t0 + tu);
879:               bu = pow(au,bm1);
880:               /* du = bu * au; */
881:               du = pow(au,beta);
882:               gu = coef*bu*(tu - t0);
883: 
884:               td = x[k-1][j][i];
885:               ad = 0.5*(t0 + td);
886:               bd = pow(ad,bm1);
887:               /* dd = bd * ad; */
888:               dd = pow(ad,beta);
889:               gd = coef*bd*(td - t0);
890: 
891:               c[0].k = k-1; c[0].j = j; c[0].i = i;
892:               v[0]   = - hxhydhz*(dd - gd);
893:               c[1].k = k; c[1].j = j; c[1].i = i-1;
894:               v[1]   = - hyhzdhx*(dw - gw);
895:               c[2].k = k; c[2].j = j; c[2].i = i;
896:               v[2]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
897:               c[3].k = k; c[3].j = j+1; c[3].i = i;
898:               v[3]   = - hzhxdhy*(dn + gn);
899:               c[4].k = k+1; c[4].j = j; c[4].i = i;
900:               v[4]   = - hxhydhz*(du + gu);
901:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

903:             /* right-hand bottom up corner */
904:             } else {

906:               td = x[k-1][j][i];
907:               ad = 0.5*(t0 + td);
908:               bd = pow(ad,bm1);
909:               /* dd = bd * ad; */
910:               dd = pow(ad,beta);
911:               gd = coef*bd*(td - t0);
912: 
913:               c[0].k = k-1; c[0].j = j; c[0].i = i;
914:               v[0]   = - hxhydhz*(dd - gd);
915:               c[1].k = k; c[1].j = j; c[1].i = i-1;
916:               v[1]   = - hyhzdhx*(dw - gw);
917:               c[2].k = k; c[2].j = j; c[2].i = i;
918:               v[2]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
919:               c[3].k = k; c[3].j = j+1; c[3].i = i;
920:               v[3]   = - hzhxdhy*(dn + gn);
921:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
922:             }

924:           /* right-hand top edge */
925:           } else if (j == my-1) {

927:             ts = x[k][j-1][i];
928:             as = 0.5*(t0 + ts);
929:             bs = pow(as,bm1);
930:             /* ds = bs * as; */
931:             ds = pow(as,beta);
932:             gs = coef*bs*(ts - t0);
933: 
934:             /* right-hand top down corner */
935:             if (k == 0) {

937:               tu = x[k+1][j][i];
938:               au = 0.5*(t0 + tu);
939:               bu = pow(au,bm1);
940:               /* du = bu * au; */
941:               du = pow(au,beta);
942:               gu = coef*bu*(tu - t0);
943: 
944:               c[0].k = k; c[0].j = j-1; c[0].i = i;
945:               v[0]   = - hzhxdhy*(ds - gs);
946:               c[1].k = k; c[1].j = j; c[1].i = i-1;
947:               v[1]   = - hyhzdhx*(dw - gw);
948:               c[2].k = k; c[2].j = j; c[2].i = i;
949:               v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
950:               c[3].k = k+1; c[3].j = j; c[3].i = i;
951:               v[3]   = - hxhydhz*(du + gu);
952:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

954:             /* right-hand top interior edge */
955:             } else if (k < mz-1) {

957:               tu = x[k+1][j][i];
958:               au = 0.5*(t0 + tu);
959:               bu = pow(au,bm1);
960:               /* du = bu * au; */
961:               du = pow(au,beta);
962:               gu = coef*bu*(tu - t0);
963: 
964:               td = x[k-1][j][i];
965:               ad = 0.5*(t0 + td);
966:               bd = pow(ad,bm1);
967:               /* dd = bd * ad; */
968:               dd = pow(ad,beta);
969:               gd = coef*bd*(td - t0);
970: 
971:               c[0].k = k-1; c[0].j = j; c[0].i = i;
972:               v[0]   = - hxhydhz*(dd - gd);
973:               c[1].k = k; c[1].j = j-1; c[1].i = i;
974:               v[1]   = - hzhxdhy*(ds - gs);
975:               c[2].k = k; c[2].j = j; c[2].i = i-1;
976:               v[2]   = - hyhzdhx*(dw - gw);
977:               c[3].k = k; c[3].j = j; c[3].i = i;
978:               v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
979:               c[4].k = k+1; c[4].j = j; c[4].i = i;
980:               v[4]   = - hxhydhz*(du + gu);
981:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

983:             /* right-hand top up corner */
984:             } else {

986:               td = x[k-1][j][i];
987:               ad = 0.5*(t0 + td);
988:               bd = pow(ad,bm1);
989:               /* dd = bd * ad; */
990:               dd = pow(ad,beta);
991:               gd = coef*bd*(td - t0);
992: 
993:               c[0].k = k-1; c[0].j = j; c[0].i = i;
994:               v[0]   = - hxhydhz*(dd - gd);
995:               c[1].k = k; c[1].j = j-1; c[1].i = i;
996:               v[1]   = - hzhxdhy*(ds - gs);
997:               c[2].k = k; c[2].j = j; c[2].i = i-1;
998:               v[2]   = - hyhzdhx*(dw - gw);
999:               c[3].k = k; c[3].j = j; c[3].i = i;
1000:               v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1001:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
1002:             }

1004:           } else {

1006:             ts = x[k][j-1][i];
1007:             as = 0.5*(t0 + ts);
1008:             bs = pow(as,bm1);
1009:             /* ds = bs * as; */
1010:             ds = pow(as,beta);
1011:             gs = coef*bs*(t0 - ts);
1012: 
1013:             tn = x[k][j+1][i];
1014:             an = 0.5*(t0 + tn);
1015:             bn = pow(an,bm1);
1016:             /* dn = bn * an; */
1017:             dn = pow(an,beta);
1018:             gn = coef*bn*(tn - t0);

1020:             /* right-hand down interior edge */
1021:             if (k == 0) {

1023:               tu = x[k+1][j][i];
1024:               au = 0.5*(t0 + tu);
1025:               bu = pow(au,bm1);
1026:               /* du = bu * au; */
1027:               du = pow(au,beta);
1028:               gu = coef*bu*(tu - t0);

1030:               c[0].k = k; c[0].j = j-1; c[0].i = i;
1031:               v[0]   = - hzhxdhy*(ds - gs);
1032:               c[1].k = k; c[1].j = j; c[1].i = i-1;
1033:               v[1]   = - hyhzdhx*(dw - gw);
1034:               c[2].k = k; c[2].j = j; c[2].i = i;
1035:               v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1036:               c[3].k = k; c[3].j = j+1; c[3].i = i;
1037:               v[3]   = - hzhxdhy*(dn + gn);
1038:               c[4].k = k+1; c[4].j = j; c[4].i = i;
1039:               v[4]   = - hxhydhz*(du + gu);
1040:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1041:             }

1043:             /* right-hand up interior edge */
1044:             else if (k == mz-1) {

1046:               td = x[k-1][j][i];
1047:               ad = 0.5*(t0 + td);
1048:               bd = pow(ad,bm1);
1049:               /* dd = bd * ad; */
1050:               dd = pow(ad,beta);
1051:               gd = coef*bd*(t0 - td);
1052: 
1053:               c[0].k = k-1; c[0].j = j; c[0].i = i;
1054:               v[0]   = - hxhydhz*(dd - gd);
1055:               c[1].k = k; c[1].j = j-1; c[1].i = i;
1056:               v[1]   = - hzhxdhy*(ds - gs);
1057:               c[2].k = k; c[2].j = j; c[2].i = i-1;
1058:               v[2]   = - hyhzdhx*(dw - gw);
1059:               c[3].k = k; c[3].j = j; c[3].i = i;
1060:               v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1061:               c[4].k = k; c[4].j = j+1; c[4].i = i;
1062:               v[4]   = - hzhxdhy*(dn + gn);
1063:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1064:             }

1066:             /* right-hand interior plane */
1067:             else {

1069:               td = x[k-1][j][i];
1070:               ad = 0.5*(t0 + td);
1071:               bd = pow(ad,bm1);
1072:               /* dd = bd * ad; */
1073:               dd = pow(ad,beta);
1074:               gd = coef*bd*(t0 - td);
1075: 
1076:               tu = x[k+1][j][i];
1077:               au = 0.5*(t0 + tu);
1078:               bu = pow(au,bm1);
1079:               /* du = bu * au; */
1080:               du = pow(au,beta);
1081:               gu = coef*bu*(tu - t0);

1083:               c[0].k = k-1; c[0].j = j; c[0].i = i;
1084:               v[0]   = - hxhydhz*(dd - gd);
1085:               c[1].k = k; c[1].j = j-1; c[1].i = i;
1086:               v[1]   = - hzhxdhy*(ds - gs);
1087:               c[2].k = k; c[2].j = j; c[2].i = i-1;
1088:               v[2]   = - hyhzdhx*(dw - gw);
1089:               c[3].k = k; c[3].j = j; c[3].i = i;
1090:               v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1091:               c[4].k = k; c[4].j = j+1; c[4].i = i;
1092:               v[4]   = - hzhxdhy*(dn + gn);
1093:               c[5].k = k+1; c[5].j = j; c[5].i = i;
1094:               v[5]   = - hxhydhz*(du + gu);
1095:                 MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1096:             }
1097:           }

1099:         } else if (j == 0) {

1101:           tw = x[k][j][i-1];
1102:           aw = 0.5*(t0 + tw);
1103:           bw = pow(aw,bm1);
1104:           /* dw = bw * aw */
1105:           dw = pow(aw,beta);
1106:           gw = coef*bw*(t0 - tw);

1108:           te = x[k][j][i+1];
1109:           ae = 0.5*(t0 + te);
1110:           be = pow(ae,bm1);
1111:           /* de = be * ae; */
1112:           de = pow(ae,beta);
1113:           ge = coef*be*(te - t0);

1115:           tn = x[k][j+1][i];
1116:           an = 0.5*(t0 + tn);
1117:           bn = pow(an,bm1);
1118:           /* dn = bn * an; */
1119:           dn = pow(an,beta);
1120:           gn = coef*bn*(tn - t0);


1123:           /* bottom down interior edge */
1124:           if (k == 0) {

1126:             tu = x[k+1][j][i];
1127:             au = 0.5*(t0 + tu);
1128:             bu = pow(au,bm1);
1129:             /* du = bu * au; */
1130:             du = pow(au,beta);
1131:             gu = coef*bu*(tu - t0);

1133:             c[0].k = k; c[0].j = j; c[0].i = i-1;
1134:             v[0]   = - hyhzdhx*(dw - gw);
1135:             c[1].k = k; c[1].j = j; c[1].i = i;
1136:             v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1137:             c[2].k = k; c[2].j = j; c[2].i = i+1;
1138:             v[2]   = - hyhzdhx*(de + ge);
1139:             c[3].k = k; c[3].j = j+1; c[3].i = i;
1140:             v[3]   = - hzhxdhy*(dn + gn);
1141:             c[4].k = k+1; c[4].j = j; c[4].i = i;
1142:             v[4]   = - hxhydhz*(du + gu);
1143:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1144:           }

1146:           /* bottom up interior edge */
1147:           else if (k == mz-1) {

1149:             td = x[k-1][j][i];
1150:             ad = 0.5*(t0 + td);
1151:             bd = pow(ad,bm1);
1152:             /* dd = bd * ad; */
1153:             dd = pow(ad,beta);
1154:             gd = coef*bd*(td - t0);

1156:             c[0].k = k-1; c[0].j = j; c[0].i = i;
1157:             v[0]   = - hxhydhz*(dd - gd);
1158:             c[1].k = k; c[1].j = j; c[1].i = i-1;
1159:             v[1]   = - hyhzdhx*(dw - gw);
1160:             c[2].k = k; c[2].j = j; c[2].i = i;
1161:             v[2]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1162:             c[3].k = k; c[3].j = j; c[3].i = i+1;
1163:             v[3]   = - hyhzdhx*(de + ge);
1164:             c[4].k = k; c[4].j = j+1; c[4].i = i;
1165:             v[4]   = - hzhxdhy*(dn + gn);
1166:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1167:           }

1169:           /* bottom interior plane */
1170:           else {

1172:             tu = x[k+1][j][i];
1173:             au = 0.5*(t0 + tu);
1174:             bu = pow(au,bm1);
1175:             /* du = bu * au; */
1176:             du = pow(au,beta);
1177:             gu = coef*bu*(tu - t0);

1179:             td = x[k-1][j][i];
1180:             ad = 0.5*(t0 + td);
1181:             bd = pow(ad,bm1);
1182:             /* dd = bd * ad; */
1183:             dd = pow(ad,beta);
1184:             gd = coef*bd*(td - t0);

1186:             c[0].k = k-1; c[0].j = j; c[0].i = i;
1187:             v[0]   = - hxhydhz*(dd - gd);
1188:             c[1].k = k; c[1].j = j; c[1].i = i-1;
1189:             v[1]   = - hyhzdhx*(dw - gw);
1190:             c[2].k = k; c[2].j = j; c[2].i = i;
1191:             v[2]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1192:             c[3].k = k; c[3].j = j; c[3].i = i+1;
1193:             v[3]   = - hyhzdhx*(de + ge);
1194:             c[4].k = k; c[4].j = j+1; c[4].i = i;
1195:             v[4]   = - hzhxdhy*(dn + gn);
1196:             c[5].k = k+1; c[5].j = j; c[5].i = i;
1197:             v[5]   = - hxhydhz*(du + gu);
1198:               MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1199:           }

1201:         } else if (j == my-1) {

1203:           tw = x[k][j][i-1];
1204:           aw = 0.5*(t0 + tw);
1205:           bw = pow(aw,bm1);
1206:           /* dw = bw * aw */
1207:           dw = pow(aw,beta);
1208:           gw = coef*bw*(t0 - tw);

1210:           te = x[k][j][i+1];
1211:           ae = 0.5*(t0 + te);
1212:           be = pow(ae,bm1);
1213:           /* de = be * ae; */
1214:           de = pow(ae,beta);
1215:           ge = coef*be*(te - t0);

1217:           ts = x[k][j-1][i];
1218:           as = 0.5*(t0 + ts);
1219:           bs = pow(as,bm1);
1220:           /* ds = bs * as; */
1221:           ds = pow(as,beta);
1222:           gs = coef*bs*(t0 - ts);
1223: 
1224:           /* top down interior edge */
1225:           if (k == 0) {

1227:             tu = x[k+1][j][i];
1228:             au = 0.5*(t0 + tu);
1229:             bu = pow(au,bm1);
1230:             /* du = bu * au; */
1231:             du = pow(au,beta);
1232:             gu = coef*bu*(tu - t0);

1234:             c[0].k = k; c[0].j = j-1; c[0].i = i;
1235:             v[0]   = - hzhxdhy*(ds - gs);
1236:             c[1].k = k; c[1].j = j; c[1].i = i-1;
1237:             v[1]   = - hyhzdhx*(dw - gw);
1238:             c[2].k = k; c[2].j = j; c[2].i = i;
1239:             v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1240:             c[3].k = k; c[3].j = j; c[3].i = i+1;
1241:             v[3]   = - hyhzdhx*(de + ge);
1242:             c[4].k = k+1; c[4].j = j; c[4].i = i;
1243:             v[4]   = - hxhydhz*(du + gu);
1244:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1245:           }

1247:           /* top up interior edge */
1248:           else if (k == mz-1) {

1250:             td = x[k-1][j][i];
1251:             ad = 0.5*(t0 + td);
1252:             bd = pow(ad,bm1);
1253:             /* dd = bd * ad; */
1254:             dd = pow(ad,beta);
1255:             gd = coef*bd*(td - t0);

1257:             c[0].k = k-1; c[0].j = j; c[0].i = i;
1258:             v[0]   = - hxhydhz*(dd - gd);
1259:             c[1].k = k; c[1].j = j-1; c[1].i = i;
1260:             v[1]   = - hzhxdhy*(ds - gs);
1261:             c[2].k = k; c[2].j = j; c[2].i = i-1;
1262:             v[2]   = - hyhzdhx*(dw - gw);
1263:             c[3].k = k; c[3].j = j; c[3].i = i;
1264:             v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1265:             c[4].k = k; c[4].j = j; c[4].i = i+1;
1266:             v[4]   = - hyhzdhx*(de + ge);
1267:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1268:           }

1270:           /* top interior plane */
1271:           else {

1273:             tu = x[k+1][j][i];
1274:             au = 0.5*(t0 + tu);
1275:             bu = pow(au,bm1);
1276:             /* du = bu * au; */
1277:             du = pow(au,beta);
1278:             gu = coef*bu*(tu - t0);

1280:             td = x[k-1][j][i];
1281:             ad = 0.5*(t0 + td);
1282:             bd = pow(ad,bm1);
1283:             /* dd = bd * ad; */
1284:             dd = pow(ad,beta);
1285:             gd = coef*bd*(td - t0);

1287:             c[0].k = k-1; c[0].j = j; c[0].i = i;
1288:             v[0]   = - hxhydhz*(dd - gd);
1289:             c[1].k = k; c[1].j = j-1; c[1].i = i;
1290:             v[1]   = - hzhxdhy*(ds - gs);
1291:             c[2].k = k; c[2].j = j; c[2].i = i-1;
1292:             v[2]   = - hyhzdhx*(dw - gw);
1293:             c[3].k = k; c[3].j = j; c[3].i = i;
1294:             v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1295:             c[4].k = k; c[4].j = j; c[4].i = i+1;
1296:             v[4]   = - hyhzdhx*(de + ge);
1297:             c[5].k = k+1; c[5].j = j; c[5].i = i;
1298:             v[5]   = - hxhydhz*(du + gu);
1299:               MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1300:           }

1302:         } else if (k == 0) {

1304:           /* down interior plane */

1306:           tw = x[k][j][i-1];
1307:           aw = 0.5*(t0 + tw);
1308:           bw = pow(aw,bm1);
1309:           /* dw = bw * aw */
1310:           dw = pow(aw,beta);
1311:           gw = coef*bw*(t0 - tw);

1313:           te = x[k][j][i+1];
1314:           ae = 0.5*(t0 + te);
1315:           be = pow(ae,bm1);
1316:           /* de = be * ae; */
1317:           de = pow(ae,beta);
1318:           ge = coef*be*(te - t0);

1320:           ts = x[k][j-1][i];
1321:           as = 0.5*(t0 + ts);
1322:           bs = pow(as,bm1);
1323:           /* ds = bs * as; */
1324:           ds = pow(as,beta);
1325:           gs = coef*bs*(t0 - ts);
1326: 
1327:           tn = x[k][j+1][i];
1328:           an = 0.5*(t0 + tn);
1329:           bn = pow(an,bm1);
1330:           /* dn = bn * an; */
1331:           dn = pow(an,beta);
1332:           gn = coef*bn*(tn - t0);
1333: 
1334:           tu = x[k+1][j][i];
1335:           au = 0.5*(t0 + tu);
1336:           bu = pow(au,bm1);
1337:           /* du = bu * au; */
1338:           du = pow(au,beta);
1339:           gu = coef*bu*(tu - t0);

1341:           c[0].k = k; c[0].j = j-1; c[0].i = i;
1342:           v[0]   = - hzhxdhy*(ds - gs);
1343:           c[1].k = k; c[1].j = j; c[1].i = i-1;
1344:           v[1]   = - hyhzdhx*(dw - gw);
1345:           c[2].k = k; c[2].j = j; c[2].i = i;
1346:           v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1347:           c[3].k = k; c[3].j = j; c[3].i = i+1;
1348:           v[3]   = - hyhzdhx*(de + ge);
1349:           c[4].k = k; c[4].j = j+1; c[4].i = i;
1350:           v[4]   = - hzhxdhy*(dn + gn);
1351:           c[5].k = k+1; c[5].j = j; c[5].i = i;
1352:           v[5]   = - hxhydhz*(du + gu);
1353:             MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1354:         }
1355: 
1356:         else if (k == mz-1) {

1358:           /* up interior plane */

1360:           tw = x[k][j][i-1];
1361:           aw = 0.5*(t0 + tw);
1362:           bw = pow(aw,bm1);
1363:           /* dw = bw * aw */
1364:           dw = pow(aw,beta);
1365:           gw = coef*bw*(t0 - tw);

1367:           te = x[k][j][i+1];
1368:           ae = 0.5*(t0 + te);
1369:           be = pow(ae,bm1);
1370:           /* de = be * ae; */
1371:           de = pow(ae,beta);
1372:           ge = coef*be*(te - t0);

1374:           ts = x[k][j-1][i];
1375:           as = 0.5*(t0 + ts);
1376:           bs = pow(as,bm1);
1377:           /* ds = bs * as; */
1378:           ds = pow(as,beta);
1379:           gs = coef*bs*(t0 - ts);
1380: 
1381:           tn = x[k][j+1][i];
1382:           an = 0.5*(t0 + tn);
1383:           bn = pow(an,bm1);
1384:           /* dn = bn * an; */
1385:           dn = pow(an,beta);
1386:           gn = coef*bn*(tn - t0);
1387: 
1388:           td = x[k-1][j][i];
1389:           ad = 0.5*(t0 + td);
1390:           bd = pow(ad,bm1);
1391:           /* dd = bd * ad; */
1392:           dd = pow(ad,beta);
1393:           gd = coef*bd*(t0 - td);
1394: 
1395:           c[0].k = k-1; c[0].j = j; c[0].i = i;
1396:           v[0]   = - hxhydhz*(dd - gd);
1397:           c[1].k = k; c[1].j = j-1; c[1].i = i;
1398:           v[1]   = - hzhxdhy*(ds - gs);
1399:           c[2].k = k; c[2].j = j; c[2].i = i-1;
1400:           v[2]   = - hyhzdhx*(dw - gw);
1401:           c[3].k = k; c[3].j = j; c[3].i = i;
1402:           v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1403:           c[4].k = k; c[4].j = j; c[4].i = i+1;
1404:           v[4]   = - hyhzdhx*(de + ge);
1405:           c[5].k = k; c[5].j = j+1; c[5].i = i;
1406:           v[5]   = - hzhxdhy*(dn + gn);
1407:             MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1408:         }
1409:       }
1410:     }
1411:   }
1412:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
1413:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
1414:   if (jac != *J) {
1415:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1416:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1417:   }
1418:   DMDAVecRestoreArray(da,localX,&x);
1419:   DMRestoreLocalVector(da,&localX);

1421:   PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
1422:   return(0);
1423: }