Actual source code: ex20.c

petsc-3.4.4 2014-03-13
  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: /*

 23:     This example models the partial differential equation

 25:          - Div(alpha* T^beta (GRAD T)) = 0.

 27:     where beta = 2.5 and alpha = 1.0

 29:     BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = dT/dn_up = dT/dn_down = 0.

 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.

 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

 40: */

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

 45: /* User-defined application context */

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

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

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

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

 69:   PetscInitialize(&argc,&argv,NULL,help);

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

 81:   /*
 82:       Set the DMDA (grid structure) for the grids.
 83:   */
 84:   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);
 85:   DMSetApplicationContext(da,&user);

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

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

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

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

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

128:   /* Compute initial guess */
129:   for (k=zs; k<zs+zm; k++) {
130:     for (j=ys; j<ys+ym; j++) {
131:       for (i=xs; i<xs+xm; i++) {
132:         x[k][j][i] = user->tleft;
133:       }
134:     }
135:   }
136:   DMDAVecRestoreArray(da,X,&x);
137:   return(0);
138: }
139: /* --------------------  Evaluate Function F(x) --------------------- */
142: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
143: {
144:   AppCtx         *user = (AppCtx*)ptr;
146:   PetscInt       i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
147:   PetscScalar    zero = 0.0,one = 1.0;
148:   PetscScalar    hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
149:   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;
150:   PetscScalar    tleft,tright,beta,td,ad,dd,fd=0.0,tu,au,du=0.0,fu=0.0;
151:   PetscScalar    ***x,***f;
152:   Vec            localX;
153:   DM             da;

156:   SNESGetDM(snes,&da);
157:   DMGetLocalVector(da,&localX);
158:   DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
159:   hx      = one/(PetscReal)(mx-1);  hy    = one/(PetscReal)(my-1);  hz = one/(PetscReal)(mz-1);
160:   hxhydhz = hx*hy/hz;   hyhzdhx = hy*hz/hx;   hzhxdhy = hz*hx/hy;
161:   tleft   = user->tleft;         tright = user->tright;
162:   beta    = user->beta;

164:   /* Get ghost points */
165:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
166:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
167:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
168:   DMDAVecGetArray(da,localX,&x);
169:   DMDAVecGetArray(da,F,&f);

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

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

179:           /* general interior volume */

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

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

191:           ts = x[k][j-1][i];
192:           as = 0.5*(t0 + ts);
193:           ds = PetscPowScalar(as,beta);
194:           fs = ds*(t0 - ts);

196:           tn = x[k][j+1][i];
197:           an = 0.5*(t0 + tn);
198:           dn = PetscPowScalar(an,beta);
199:           fn = dn*(tn - t0);

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

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

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

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

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

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

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

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

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

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

262:           /* right-hand (east) boundary */
263:           tw = x[k][j][i-1];
264:           aw = 0.5*(t0 + tw);
265:           dw = PetscPowScalar(aw,beta);
266:           fw = dw*(t0 - tw);

268:           te = tright;
269:           ae = 0.5*(t0 + te);
270:           de = PetscPowScalar(ae,beta);
271:           fe = de*(te - t0);

273:           if (j > 0) {
274:             ts = x[k][j-1][i];
275:             as = 0.5*(t0 + ts);
276:             ds = PetscPowScalar(as,beta);
277:             fs = ds*(t0 - ts);
278:           } else {
279:             fs = zero;
280:           }

282:           if (j < my-1) {
283:             tn = x[k][j+1][i];
284:             an = 0.5*(t0 + tn);
285:             dn = PetscPowScalar(an,beta);
286:             fn = dn*(tn - t0);
287:           } else {
288:             fn = zero;
289:           }

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

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

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

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

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

322:           fs = zero;

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

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

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

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

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

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

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

365:           fn = zero;

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

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

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

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

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

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

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

408:           fd = zero;

410:           tu = x[k+1][j][i];
411:           au = 0.5*(t0 + tu);
412:           du = PetscPowScalar(au,beta);
413:           fu = du*(tu - t0);

415:         } else if (k == mz-1) {

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

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

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

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

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

443:           fu = zero;
444:         }

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

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

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

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

497:           /* general interior volume */

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

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

513:           ts = x[k][j-1][i];
514:           as = 0.5*(t0 + ts);
515:           bs = PetscPowScalar(as,bm1);
516:           /* ds = bs * as; */
517:           ds = PetscPowScalar(as,beta);
518:           gs = coef*bs*(t0 - ts);

520:           tn = x[k][j+1][i];
521:           an = 0.5*(t0 + tn);
522:           bn = PetscPowScalar(an,bm1);
523:           /* dn = bn * an; */
524:           dn = PetscPowScalar(an,beta);
525:           gn = coef*bn*(tn - t0);

527:           td = x[k-1][j][i];
528:           ad = 0.5*(t0 + td);
529:           bd = PetscPowScalar(ad,bm1);
530:           /* dd = bd * ad; */
531:           dd = PetscPowScalar(ad,beta);
532:           gd = coef*bd*(t0 - td);

534:           tu = x[k+1][j][i];
535:           au = 0.5*(t0 + tu);
536:           bu = PetscPowScalar(au,bm1);
537:           /* du = bu * au; */
538:           du = PetscPowScalar(au,beta);
539:           gu = coef*bu*(tu - t0);

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

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

558:           /* left-hand plane boundary */
559:           tw = tleft;
560:           aw = 0.5*(t0 + tw);
561:           bw = PetscPowScalar(aw,bm1);
562:           /* dw = bw * aw */
563:           dw = PetscPowScalar(aw,beta);
564:           gw = coef*bw*(t0 - tw);

566:           te = x[k][j][i+1];
567:           ae = 0.5*(t0 + te);
568:           be = PetscPowScalar(ae,bm1);
569:           /* de = be * ae; */
570:           de = PetscPowScalar(ae,beta);
571:           ge = coef*be*(te - t0);

573:           /* left-hand bottom edge */
574:           if (j == 0) {

576:             tn = x[k][j+1][i];
577:             an = 0.5*(t0 + tn);
578:             bn = PetscPowScalar(an,bm1);
579:             /* dn = bn * an; */
580:             dn = PetscPowScalar(an,beta);
581:             gn = coef*bn*(tn - t0);

583:             /* left-hand bottom down corner */
584:             if (k == 0) {

586:               tu = x[k+1][j][i];
587:               au = 0.5*(t0 + tu);
588:               bu = PetscPowScalar(au,bm1);
589:               /* du = bu * au; */
590:               du = PetscPowScalar(au,beta);
591:               gu = coef*bu*(tu - t0);

593:               c[0].k = k; c[0].j = j; c[0].i = i;
594:               v[0]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
595:               c[1].k = k; c[1].j = j; c[1].i = i+1;
596:               v[1]   = -hyhzdhx*(de + ge);
597:               c[2].k = k; c[2].j = j+1; c[2].i = i;
598:               v[2]   = -hzhxdhy*(dn + gn);
599:               c[3].k = k+1; c[3].j = j; c[3].i = i;
600:               v[3]   = -hxhydhz*(du + gu);
601:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

606:               tu = x[k+1][j][i];
607:               au = 0.5*(t0 + tu);
608:               bu = PetscPowScalar(au,bm1);
609:               /* du = bu * au; */
610:               du = PetscPowScalar(au,beta);
611:               gu = coef*bu*(tu - t0);

613:               td = x[k-1][j][i];
614:               ad = 0.5*(t0 + td);
615:               bd = PetscPowScalar(ad,bm1);
616:               /* dd = bd * ad; */
617:               dd = PetscPowScalar(ad,beta);
618:               gd = coef*bd*(td - t0);

620:               c[0].k = k-1; c[0].j = j; c[0].i = i;
621:               v[0]   = -hxhydhz*(dd - gd);
622:               c[1].k = k; c[1].j = j; c[1].i = i;
623:               v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
624:               c[2].k = k; c[2].j = j; c[2].i = i+1;
625:               v[2]   = -hyhzdhx*(de + ge);
626:               c[3].k = k; c[3].j = j+1; c[3].i = i;
627:               v[3]   = -hzhxdhy*(dn + gn);
628:               c[4].k = k+1; c[4].j = j; c[4].i = i;
629:               v[4]   = -hxhydhz*(du + gu);
630:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

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

635:               td = x[k-1][j][i];
636:               ad = 0.5*(t0 + td);
637:               bd = PetscPowScalar(ad,bm1);
638:               /* dd = bd * ad; */
639:               dd = PetscPowScalar(ad,beta);
640:               gd = coef*bd*(td - t0);

642:               c[0].k = k-1; c[0].j = j; c[0].i = i;
643:               v[0]   = -hxhydhz*(dd - gd);
644:               c[1].k = k; c[1].j = j; c[1].i = i;
645:               v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
646:               c[2].k = k; c[2].j = j; c[2].i = i+1;
647:               v[2]   = -hyhzdhx*(de + ge);
648:               c[3].k = k; c[3].j = j+1; c[3].i = i;
649:               v[3]   = -hzhxdhy*(dn + gn);
650:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
651:             }

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

656:             ts = x[k][j-1][i];
657:             as = 0.5*(t0 + ts);
658:             bs = PetscPowScalar(as,bm1);
659:             /* ds = bs * as; */
660:             ds = PetscPowScalar(as,beta);
661:             gs = coef*bs*(ts - t0);

663:             /* left-hand top down corner */
664:             if (k == 0) {

666:               tu = x[k+1][j][i];
667:               au = 0.5*(t0 + tu);
668:               bu = PetscPowScalar(au,bm1);
669:               /* du = bu * au; */
670:               du = PetscPowScalar(au,beta);
671:               gu = coef*bu*(tu - t0);

673:               c[0].k = k; c[0].j = j-1; c[0].i = i;
674:               v[0]   = -hzhxdhy*(ds - gs);
675:               c[1].k = k; c[1].j = j; c[1].i = i;
676:               v[1]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
677:               c[2].k = k; c[2].j = j; c[2].i = i+1;
678:               v[2]   = -hyhzdhx*(de + ge);
679:               c[3].k = k+1; c[3].j = j; c[3].i = i;
680:               v[3]   = -hxhydhz*(du + gu);
681:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

686:               tu = x[k+1][j][i];
687:               au = 0.5*(t0 + tu);
688:               bu = PetscPowScalar(au,bm1);
689:               /* du = bu * au; */
690:               du = PetscPowScalar(au,beta);
691:               gu = coef*bu*(tu - t0);

693:               td = x[k-1][j][i];
694:               ad = 0.5*(t0 + td);
695:               bd = PetscPowScalar(ad,bm1);
696:               /* dd = bd * ad; */
697:               dd = PetscPowScalar(ad,beta);
698:               gd = coef*bd*(td - t0);

700:               c[0].k = k-1; c[0].j = j; c[0].i = i;
701:               v[0]   = -hxhydhz*(dd - gd);
702:               c[1].k = k; c[1].j = j-1; c[1].i = i;
703:               v[1]   = -hzhxdhy*(ds - gs);
704:               c[2].k = k; c[2].j = j; c[2].i = i;
705:               v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
706:               c[3].k = k; c[3].j = j; c[3].i = i+1;
707:               v[3]   = -hyhzdhx*(de + ge);
708:               c[4].k = k+1; c[4].j = j; c[4].i = i;
709:               v[4]   = -hxhydhz*(du + gu);
710:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

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

715:               td = x[k-1][j][i];
716:               ad = 0.5*(t0 + td);
717:               bd = PetscPowScalar(ad,bm1);
718:               /* dd = bd * ad; */
719:               dd = PetscPowScalar(ad,beta);
720:               gd = coef*bd*(td - t0);

722:               c[0].k = k-1; c[0].j = j; c[0].i = i;
723:               v[0]   = -hxhydhz*(dd - gd);
724:               c[1].k = k; c[1].j = j-1; c[1].i = i;
725:               v[1]   = -hzhxdhy*(ds - gs);
726:               c[2].k = k; c[2].j = j; c[2].i = i;
727:               v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
728:               c[3].k = k; c[3].j = j; c[3].i = i+1;
729:               v[3]   = -hyhzdhx*(de + ge);
730:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
731:             }

733:           } else {

735:             ts = x[k][j-1][i];
736:             as = 0.5*(t0 + ts);
737:             bs = PetscPowScalar(as,bm1);
738:             /* ds = bs * as; */
739:             ds = PetscPowScalar(as,beta);
740:             gs = coef*bs*(t0 - ts);

742:             tn = x[k][j+1][i];
743:             an = 0.5*(t0 + tn);
744:             bn = PetscPowScalar(an,bm1);
745:             /* dn = bn * an; */
746:             dn = PetscPowScalar(an,beta);
747:             gn = coef*bn*(tn - t0);

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

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

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

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

773:               td = x[k-1][j][i];
774:               ad = 0.5*(t0 + td);
775:               bd = PetscPowScalar(ad,bm1);
776:               /* dd = bd * ad; */
777:               dd = PetscPowScalar(ad,beta);
778:               gd = coef*bd*(t0 - td);

780:               c[0].k = k-1; c[0].j = j; c[0].i = i;
781:               v[0]   = -hxhydhz*(dd - gd);
782:               c[1].k = k; c[1].j = j-1; c[1].i = i;
783:               v[1]   = -hzhxdhy*(ds - gs);
784:               c[2].k = k; c[2].j = j; c[2].i = i;
785:               v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
786:               c[3].k = k; c[3].j = j; c[3].i = i+1;
787:               v[3]   = -hyhzdhx*(de + ge);
788:               c[4].k = k; c[4].j = j+1; c[4].i = i;
789:               v[4]   = -hzhxdhy*(dn + gn);
790:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
791:             } else { /* left-hand interior plane */

793:               td = x[k-1][j][i];
794:               ad = 0.5*(t0 + td);
795:               bd = PetscPowScalar(ad,bm1);
796:               /* dd = bd * ad; */
797:               dd = PetscPowScalar(ad,beta);
798:               gd = coef*bd*(t0 - td);

800:               tu = x[k+1][j][i];
801:               au = 0.5*(t0 + tu);
802:               bu = PetscPowScalar(au,bm1);
803:               /* du = bu * au; */
804:               du = PetscPowScalar(au,beta);
805:               gu = coef*bu*(tu - t0);

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

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

825:           /* right-hand plane boundary */
826:           tw = x[k][j][i-1];
827:           aw = 0.5*(t0 + tw);
828:           bw = PetscPowScalar(aw,bm1);
829:           /* dw = bw * aw */
830:           dw = PetscPowScalar(aw,beta);
831:           gw = coef*bw*(t0 - tw);

833:           te = tright;
834:           ae = 0.5*(t0 + te);
835:           be = PetscPowScalar(ae,bm1);
836:           /* de = be * ae; */
837:           de = PetscPowScalar(ae,beta);
838:           ge = coef*be*(te - t0);

840:           /* right-hand bottom edge */
841:           if (j == 0) {

843:             tn = x[k][j+1][i];
844:             an = 0.5*(t0 + tn);
845:             bn = PetscPowScalar(an,bm1);
846:             /* dn = bn * an; */
847:             dn = PetscPowScalar(an,beta);
848:             gn = coef*bn*(tn - t0);

850:             /* right-hand bottom down corner */
851:             if (k == 0) {

853:               tu = x[k+1][j][i];
854:               au = 0.5*(t0 + tu);
855:               bu = PetscPowScalar(au,bm1);
856:               /* du = bu * au; */
857:               du = PetscPowScalar(au,beta);
858:               gu = coef*bu*(tu - t0);

860:               c[0].k = k; c[0].j = j; c[0].i = i-1;
861:               v[0]   = -hyhzdhx*(dw - gw);
862:               c[1].k = k; c[1].j = j; c[1].i = i;
863:               v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
864:               c[2].k = k; c[2].j = j+1; c[2].i = i;
865:               v[2]   = -hzhxdhy*(dn + gn);
866:               c[3].k = k+1; c[3].j = j; c[3].i = i;
867:               v[3]   = -hxhydhz*(du + gu);
868:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

873:               tu = x[k+1][j][i];
874:               au = 0.5*(t0 + tu);
875:               bu = PetscPowScalar(au,bm1);
876:               /* du = bu * au; */
877:               du = PetscPowScalar(au,beta);
878:               gu = coef*bu*(tu - t0);

880:               td = x[k-1][j][i];
881:               ad = 0.5*(t0 + td);
882:               bd = PetscPowScalar(ad,bm1);
883:               /* dd = bd * ad; */
884:               dd = PetscPowScalar(ad,beta);
885:               gd = coef*bd*(td - t0);

887:               c[0].k = k-1; c[0].j = j; c[0].i = i;
888:               v[0]   = -hxhydhz*(dd - gd);
889:               c[1].k = k; c[1].j = j; c[1].i = i-1;
890:               v[1]   = -hyhzdhx*(dw - gw);
891:               c[2].k = k; c[2].j = j; c[2].i = i;
892:               v[2]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
893:               c[3].k = k; c[3].j = j+1; c[3].i = i;
894:               v[3]   = -hzhxdhy*(dn + gn);
895:               c[4].k = k+1; c[4].j = j; c[4].i = i;
896:               v[4]   = -hxhydhz*(du + gu);
897:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

899:               /* right-hand bottom up corner */
900:             } else {

902:               td = x[k-1][j][i];
903:               ad = 0.5*(t0 + td);
904:               bd = PetscPowScalar(ad,bm1);
905:               /* dd = bd * ad; */
906:               dd = PetscPowScalar(ad,beta);
907:               gd = coef*bd*(td - t0);

909:               c[0].k = k-1; c[0].j = j; c[0].i = i;
910:               v[0]   = -hxhydhz*(dd - gd);
911:               c[1].k = k; c[1].j = j; c[1].i = i-1;
912:               v[1]   = -hyhzdhx*(dw - gw);
913:               c[2].k = k; c[2].j = j; c[2].i = i;
914:               v[2]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
915:               c[3].k = k; c[3].j = j+1; c[3].i = i;
916:               v[3]   = -hzhxdhy*(dn + gn);
917:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
918:             }

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

923:             ts = x[k][j-1][i];
924:             as = 0.5*(t0 + ts);
925:             bs = PetscPowScalar(as,bm1);
926:             /* ds = bs * as; */
927:             ds = PetscPowScalar(as,beta);
928:             gs = coef*bs*(ts - t0);

930:             /* right-hand top down corner */
931:             if (k == 0) {

933:               tu = x[k+1][j][i];
934:               au = 0.5*(t0 + tu);
935:               bu = PetscPowScalar(au,bm1);
936:               /* du = bu * au; */
937:               du = PetscPowScalar(au,beta);
938:               gu = coef*bu*(tu - t0);

940:               c[0].k = k; c[0].j = j-1; c[0].i = i;
941:               v[0]   = -hzhxdhy*(ds - gs);
942:               c[1].k = k; c[1].j = j; c[1].i = i-1;
943:               v[1]   = -hyhzdhx*(dw - gw);
944:               c[2].k = k; c[2].j = j; c[2].i = i;
945:               v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
946:               c[3].k = k+1; c[3].j = j; c[3].i = i;
947:               v[3]   = -hxhydhz*(du + gu);
948:               MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

953:               tu = x[k+1][j][i];
954:               au = 0.5*(t0 + tu);
955:               bu = PetscPowScalar(au,bm1);
956:               /* du = bu * au; */
957:               du = PetscPowScalar(au,beta);
958:               gu = coef*bu*(tu - t0);

960:               td = x[k-1][j][i];
961:               ad = 0.5*(t0 + td);
962:               bd = PetscPowScalar(ad,bm1);
963:               /* dd = bd * ad; */
964:               dd = PetscPowScalar(ad,beta);
965:               gd = coef*bd*(td - t0);

967:               c[0].k = k-1; c[0].j = j; c[0].i = i;
968:               v[0]   = -hxhydhz*(dd - gd);
969:               c[1].k = k; c[1].j = j-1; c[1].i = i;
970:               v[1]   = -hzhxdhy*(ds - gs);
971:               c[2].k = k; c[2].j = j; c[2].i = i-1;
972:               v[2]   = -hyhzdhx*(dw - gw);
973:               c[3].k = k; c[3].j = j; c[3].i = i;
974:               v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
975:               c[4].k = k+1; c[4].j = j; c[4].i = i;
976:               v[4]   = -hxhydhz*(du + gu);
977:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

979:               /* right-hand top up corner */
980:             } else {

982:               td = x[k-1][j][i];
983:               ad = 0.5*(t0 + td);
984:               bd = PetscPowScalar(ad,bm1);
985:               /* dd = bd * ad; */
986:               dd = PetscPowScalar(ad,beta);
987:               gd = coef*bd*(td - t0);

989:               c[0].k = k-1; c[0].j = j; c[0].i = i;
990:               v[0]   = -hxhydhz*(dd - gd);
991:               c[1].k = k; c[1].j = j-1; c[1].i = i;
992:               v[1]   = -hzhxdhy*(ds - gs);
993:               c[2].k = k; c[2].j = j; c[2].i = i-1;
994:               v[2]   = -hyhzdhx*(dw - gw);
995:               c[3].k = k; c[3].j = j; c[3].i = i;
996:               v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
997:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
998:             }

1000:           } else {

1002:             ts = x[k][j-1][i];
1003:             as = 0.5*(t0 + ts);
1004:             bs = PetscPowScalar(as,bm1);
1005:             /* ds = bs * as; */
1006:             ds = PetscPowScalar(as,beta);
1007:             gs = coef*bs*(t0 - ts);

1009:             tn = x[k][j+1][i];
1010:             an = 0.5*(t0 + tn);
1011:             bn = PetscPowScalar(an,bm1);
1012:             /* dn = bn * an; */
1013:             dn = PetscPowScalar(an,beta);
1014:             gn = coef*bn*(tn - t0);

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

1019:               tu = x[k+1][j][i];
1020:               au = 0.5*(t0 + tu);
1021:               bu = PetscPowScalar(au,bm1);
1022:               /* du = bu * au; */
1023:               du = PetscPowScalar(au,beta);
1024:               gu = coef*bu*(tu - t0);

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

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

1040:               td = x[k-1][j][i];
1041:               ad = 0.5*(t0 + td);
1042:               bd = PetscPowScalar(ad,bm1);
1043:               /* dd = bd * ad; */
1044:               dd = PetscPowScalar(ad,beta);
1045:               gd = coef*bd*(t0 - td);

1047:               c[0].k = k-1; c[0].j = j; c[0].i = i;
1048:               v[0]   = -hxhydhz*(dd - gd);
1049:               c[1].k = k; c[1].j = j-1; c[1].i = i;
1050:               v[1]   = -hzhxdhy*(ds - gs);
1051:               c[2].k = k; c[2].j = j; c[2].i = i-1;
1052:               v[2]   = -hyhzdhx*(dw - gw);
1053:               c[3].k = k; c[3].j = j; c[3].i = i;
1054:               v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1055:               c[4].k = k; c[4].j = j+1; c[4].i = i;
1056:               v[4]   = -hzhxdhy*(dn + gn);
1057:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

1059:             } else { /* right-hand interior plane */

1061:               td = x[k-1][j][i];
1062:               ad = 0.5*(t0 + td);
1063:               bd = PetscPowScalar(ad,bm1);
1064:               /* dd = bd * ad; */
1065:               dd = PetscPowScalar(ad,beta);
1066:               gd = coef*bd*(t0 - td);

1068:               tu = x[k+1][j][i];
1069:               au = 0.5*(t0 + tu);
1070:               bu = PetscPowScalar(au,bm1);
1071:               /* du = bu * au; */
1072:               du = PetscPowScalar(au,beta);
1073:               gu = coef*bu*(tu - t0);

1075:               c[0].k = k-1; c[0].j = j; c[0].i = i;
1076:               v[0]   = -hxhydhz*(dd - gd);
1077:               c[1].k = k; c[1].j = j-1; c[1].i = i;
1078:               v[1]   = -hzhxdhy*(ds - gs);
1079:               c[2].k = k; c[2].j = j; c[2].i = i-1;
1080:               v[2]   = -hyhzdhx*(dw - gw);
1081:               c[3].k = k; c[3].j = j; c[3].i = i;
1082:               v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1083:               c[4].k = k; c[4].j = j+1; c[4].i = i;
1084:               v[4]   = -hzhxdhy*(dn + gn);
1085:               c[5].k = k+1; c[5].j = j; c[5].i = i;
1086:               v[5]   = -hxhydhz*(du + gu);
1087:               MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1088:             }
1089:           }

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

1093:           tw = x[k][j][i-1];
1094:           aw = 0.5*(t0 + tw);
1095:           bw = PetscPowScalar(aw,bm1);
1096:           /* dw = bw * aw */
1097:           dw = PetscPowScalar(aw,beta);
1098:           gw = coef*bw*(t0 - tw);

1100:           te = x[k][j][i+1];
1101:           ae = 0.5*(t0 + te);
1102:           be = PetscPowScalar(ae,bm1);
1103:           /* de = be * ae; */
1104:           de = PetscPowScalar(ae,beta);
1105:           ge = coef*be*(te - t0);

1107:           tn = x[k][j+1][i];
1108:           an = 0.5*(t0 + tn);
1109:           bn = PetscPowScalar(an,bm1);
1110:           /* dn = bn * an; */
1111:           dn = PetscPowScalar(an,beta);
1112:           gn = coef*bn*(tn - t0);


1115:           /* bottom down interior edge */
1116:           if (k == 0) {

1118:             tu = x[k+1][j][i];
1119:             au = 0.5*(t0 + tu);
1120:             bu = PetscPowScalar(au,bm1);
1121:             /* du = bu * au; */
1122:             du = PetscPowScalar(au,beta);
1123:             gu = coef*bu*(tu - t0);

1125:             c[0].k = k; c[0].j = j; c[0].i = i-1;
1126:             v[0]   = -hyhzdhx*(dw - gw);
1127:             c[1].k = k; c[1].j = j; c[1].i = i;
1128:             v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1129:             c[2].k = k; c[2].j = j; c[2].i = i+1;
1130:             v[2]   = -hyhzdhx*(de + ge);
1131:             c[3].k = k; c[3].j = j+1; c[3].i = i;
1132:             v[3]   = -hzhxdhy*(dn + gn);
1133:             c[4].k = k+1; c[4].j = j; c[4].i = i;
1134:             v[4]   = -hxhydhz*(du + gu);
1135:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

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

1139:             td = x[k-1][j][i];
1140:             ad = 0.5*(t0 + td);
1141:             bd = PetscPowScalar(ad,bm1);
1142:             /* dd = bd * ad; */
1143:             dd = PetscPowScalar(ad,beta);
1144:             gd = coef*bd*(td - t0);

1146:             c[0].k = k-1; c[0].j = j; c[0].i = i;
1147:             v[0]   = -hxhydhz*(dd - gd);
1148:             c[1].k = k; c[1].j = j; c[1].i = i-1;
1149:             v[1]   = -hyhzdhx*(dw - gw);
1150:             c[2].k = k; c[2].j = j; c[2].i = i;
1151:             v[2]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1152:             c[3].k = k; c[3].j = j; c[3].i = i+1;
1153:             v[3]   = -hyhzdhx*(de + ge);
1154:             c[4].k = k; c[4].j = j+1; c[4].i = i;
1155:             v[4]   = -hzhxdhy*(dn + gn);
1156:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

1158:           } else { /* bottom interior plane */

1160:             tu = x[k+1][j][i];
1161:             au = 0.5*(t0 + tu);
1162:             bu = PetscPowScalar(au,bm1);
1163:             /* du = bu * au; */
1164:             du = PetscPowScalar(au,beta);
1165:             gu = coef*bu*(tu - t0);

1167:             td = x[k-1][j][i];
1168:             ad = 0.5*(t0 + td);
1169:             bd = PetscPowScalar(ad,bm1);
1170:             /* dd = bd * ad; */
1171:             dd = PetscPowScalar(ad,beta);
1172:             gd = coef*bd*(td - t0);

1174:             c[0].k = k-1; c[0].j = j; c[0].i = i;
1175:             v[0]   = -hxhydhz*(dd - gd);
1176:             c[1].k = k; c[1].j = j; c[1].i = i-1;
1177:             v[1]   = -hyhzdhx*(dw - gw);
1178:             c[2].k = k; c[2].j = j; c[2].i = i;
1179:             v[2]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1180:             c[3].k = k; c[3].j = j; c[3].i = i+1;
1181:             v[3]   = -hyhzdhx*(de + ge);
1182:             c[4].k = k; c[4].j = j+1; c[4].i = i;
1183:             v[4]   = -hzhxdhy*(dn + gn);
1184:             c[5].k = k+1; c[5].j = j; c[5].i = i;
1185:             v[5]   = -hxhydhz*(du + gu);
1186:             MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1187:           }

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

1191:           tw = x[k][j][i-1];
1192:           aw = 0.5*(t0 + tw);
1193:           bw = PetscPowScalar(aw,bm1);
1194:           /* dw = bw * aw */
1195:           dw = PetscPowScalar(aw,beta);
1196:           gw = coef*bw*(t0 - tw);

1198:           te = x[k][j][i+1];
1199:           ae = 0.5*(t0 + te);
1200:           be = PetscPowScalar(ae,bm1);
1201:           /* de = be * ae; */
1202:           de = PetscPowScalar(ae,beta);
1203:           ge = coef*be*(te - t0);

1205:           ts = x[k][j-1][i];
1206:           as = 0.5*(t0 + ts);
1207:           bs = PetscPowScalar(as,bm1);
1208:           /* ds = bs * as; */
1209:           ds = PetscPowScalar(as,beta);
1210:           gs = coef*bs*(t0 - ts);

1212:           /* top down interior edge */
1213:           if (k == 0) {

1215:             tu = x[k+1][j][i];
1216:             au = 0.5*(t0 + tu);
1217:             bu = PetscPowScalar(au,bm1);
1218:             /* du = bu * au; */
1219:             du = PetscPowScalar(au,beta);
1220:             gu = coef*bu*(tu - t0);

1222:             c[0].k = k; c[0].j = j-1; c[0].i = i;
1223:             v[0]   = -hzhxdhy*(ds - gs);
1224:             c[1].k = k; c[1].j = j; c[1].i = i-1;
1225:             v[1]   = -hyhzdhx*(dw - gw);
1226:             c[2].k = k; c[2].j = j; c[2].i = i;
1227:             v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1228:             c[3].k = k; c[3].j = j; c[3].i = i+1;
1229:             v[3]   = -hyhzdhx*(de + ge);
1230:             c[4].k = k+1; c[4].j = j; c[4].i = i;
1231:             v[4]   = -hxhydhz*(du + gu);
1232:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

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

1236:             td = x[k-1][j][i];
1237:             ad = 0.5*(t0 + td);
1238:             bd = PetscPowScalar(ad,bm1);
1239:             /* dd = bd * ad; */
1240:             dd = PetscPowScalar(ad,beta);
1241:             gd = coef*bd*(td - t0);

1243:             c[0].k = k-1; c[0].j = j; c[0].i = i;
1244:             v[0]   = -hxhydhz*(dd - gd);
1245:             c[1].k = k; c[1].j = j-1; c[1].i = i;
1246:             v[1]   = -hzhxdhy*(ds - gs);
1247:             c[2].k = k; c[2].j = j; c[2].i = i-1;
1248:             v[2]   = -hyhzdhx*(dw - gw);
1249:             c[3].k = k; c[3].j = j; c[3].i = i;
1250:             v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1251:             c[4].k = k; c[4].j = j; c[4].i = i+1;
1252:             v[4]   = -hyhzdhx*(de + ge);
1253:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

1255:           } else { /* top interior plane */

1257:             tu = x[k+1][j][i];
1258:             au = 0.5*(t0 + tu);
1259:             bu = PetscPowScalar(au,bm1);
1260:             /* du = bu * au; */
1261:             du = PetscPowScalar(au,beta);
1262:             gu = coef*bu*(tu - t0);

1264:             td = x[k-1][j][i];
1265:             ad = 0.5*(t0 + td);
1266:             bd = PetscPowScalar(ad,bm1);
1267:             /* dd = bd * ad; */
1268:             dd = PetscPowScalar(ad,beta);
1269:             gd = coef*bd*(td - t0);

1271:             c[0].k = k-1; c[0].j = j; c[0].i = i;
1272:             v[0]   = -hxhydhz*(dd - gd);
1273:             c[1].k = k; c[1].j = j-1; c[1].i = i;
1274:             v[1]   = -hzhxdhy*(ds - gs);
1275:             c[2].k = k; c[2].j = j; c[2].i = i-1;
1276:             v[2]   = -hyhzdhx*(dw - gw);
1277:             c[3].k = k; c[3].j = j; c[3].i = i;
1278:             v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1279:             c[4].k = k; c[4].j = j; c[4].i = i+1;
1280:             v[4]   = -hyhzdhx*(de + ge);
1281:             c[5].k = k+1; c[5].j = j; c[5].i = i;
1282:             v[5]   = -hxhydhz*(du + gu);
1283:               MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1284:           }

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

1288:           /* down interior plane */

1290:           tw = x[k][j][i-1];
1291:           aw = 0.5*(t0 + tw);
1292:           bw = PetscPowScalar(aw,bm1);
1293:           /* dw = bw * aw */
1294:           dw = PetscPowScalar(aw,beta);
1295:           gw = coef*bw*(t0 - tw);

1297:           te = x[k][j][i+1];
1298:           ae = 0.5*(t0 + te);
1299:           be = PetscPowScalar(ae,bm1);
1300:           /* de = be * ae; */
1301:           de = PetscPowScalar(ae,beta);
1302:           ge = coef*be*(te - t0);

1304:           ts = x[k][j-1][i];
1305:           as = 0.5*(t0 + ts);
1306:           bs = PetscPowScalar(as,bm1);
1307:           /* ds = bs * as; */
1308:           ds = PetscPowScalar(as,beta);
1309:           gs = coef*bs*(t0 - ts);

1311:           tn = x[k][j+1][i];
1312:           an = 0.5*(t0 + tn);
1313:           bn = PetscPowScalar(an,bm1);
1314:           /* dn = bn * an; */
1315:           dn = PetscPowScalar(an,beta);
1316:           gn = coef*bn*(tn - t0);

1318:           tu = x[k+1][j][i];
1319:           au = 0.5*(t0 + tu);
1320:           bu = PetscPowScalar(au,bm1);
1321:           /* du = bu * au; */
1322:           du = PetscPowScalar(au,beta);
1323:           gu = coef*bu*(tu - t0);

1325:           c[0].k = k; c[0].j = j-1; c[0].i = i;
1326:           v[0]   = -hzhxdhy*(ds - gs);
1327:           c[1].k = k; c[1].j = j; c[1].i = i-1;
1328:           v[1]   = -hyhzdhx*(dw - gw);
1329:           c[2].k = k; c[2].j = j; c[2].i = i;
1330:           v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1331:           c[3].k = k; c[3].j = j; c[3].i = i+1;
1332:           v[3]   = -hyhzdhx*(de + ge);
1333:           c[4].k = k; c[4].j = j+1; c[4].i = i;
1334:           v[4]   = -hzhxdhy*(dn + gn);
1335:           c[5].k = k+1; c[5].j = j; c[5].i = i;
1336:           v[5]   = -hxhydhz*(du + gu);
1337:             MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);

1339:         } else if (k == mz-1) {

1341:           /* up interior plane */

1343:           tw = x[k][j][i-1];
1344:           aw = 0.5*(t0 + tw);
1345:           bw = PetscPowScalar(aw,bm1);
1346:           /* dw = bw * aw */
1347:           dw = PetscPowScalar(aw,beta);
1348:           gw = coef*bw*(t0 - tw);

1350:           te = x[k][j][i+1];
1351:           ae = 0.5*(t0 + te);
1352:           be = PetscPowScalar(ae,bm1);
1353:           /* de = be * ae; */
1354:           de = PetscPowScalar(ae,beta);
1355:           ge = coef*be*(te - t0);

1357:           ts = x[k][j-1][i];
1358:           as = 0.5*(t0 + ts);
1359:           bs = PetscPowScalar(as,bm1);
1360:           /* ds = bs * as; */
1361:           ds = PetscPowScalar(as,beta);
1362:           gs = coef*bs*(t0 - ts);

1364:           tn = x[k][j+1][i];
1365:           an = 0.5*(t0 + tn);
1366:           bn = PetscPowScalar(an,bm1);
1367:           /* dn = bn * an; */
1368:           dn = PetscPowScalar(an,beta);
1369:           gn = coef*bn*(tn - t0);

1371:           td = x[k-1][j][i];
1372:           ad = 0.5*(t0 + td);
1373:           bd = PetscPowScalar(ad,bm1);
1374:           /* dd = bd * ad; */
1375:           dd = PetscPowScalar(ad,beta);
1376:           gd = coef*bd*(t0 - td);

1378:           c[0].k = k-1; c[0].j = j; c[0].i = i;
1379:           v[0]   = -hxhydhz*(dd - gd);
1380:           c[1].k = k; c[1].j = j-1; c[1].i = i;
1381:           v[1]   = -hzhxdhy*(ds - gs);
1382:           c[2].k = k; c[2].j = j; c[2].i = i-1;
1383:           v[2]   = -hyhzdhx*(dw - gw);
1384:           c[3].k = k; c[3].j = j; c[3].i = i;
1385:           v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1386:           c[4].k = k; c[4].j = j; c[4].i = i+1;
1387:           v[4]   = -hyhzdhx*(de + ge);
1388:           c[5].k = k; c[5].j = j+1; c[5].i = i;
1389:           v[5]   = -hzhxdhy*(dn + gn);
1390:             MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1391:         }
1392:       }
1393:     }
1394:   }
1395:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
1396:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
1397:   if (jac != *J) {
1398:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1399:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1400:   }
1401:   DMDAVecRestoreArray(da,localX,&x);
1402:   DMRestoreLocalVector(da,&localX);

1404:   PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
1405:   return(0);
1406: }