Actual source code: ex20.c


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

 14: /*

 16:     This example models the partial differential equation

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

 20:     where beta = 2.5 and alpha = 1.0

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

 24:     in the unit square, which is uniformly discretized in each of x and
 25:     y in this simple encoding.  The degrees of freedom are cell centered.

 27:     A finite volume approximation with the usual 7-point stencil
 28:     is used to discretize the boundary value problem to obtain a
 29:     nonlinear system of equations.

 31:     This code was contributed by Nickolas Jovanovic based on ex18.c

 33: */

 35: #include <petscsnes.h>
 36: #include <petscdm.h>
 37: #include <petscdmda.h>

 39: /* User-defined application context */

 41: typedef struct {
 42:   PetscReal tleft,tright;     /* Dirichlet boundary conditions */
 43:   PetscReal beta,bm1,coef;    /* nonlinear diffusivity parameterizations */
 44: } AppCtx;

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

 48: extern PetscErrorCode FormInitialGuess(SNES,Vec,void*);
 49: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 50: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);

 52: int main(int argc,char **argv)
 53: {
 54:   SNES           snes;
 55:   AppCtx         user;
 56:   PetscInt       its,lits;
 57:   PetscReal      litspit;
 58:   DM             da;

 60:   PetscInitialize(&argc,&argv,NULL,help);
 61:   /* set problem parameters */
 62:   user.tleft  = 1.0;
 63:   user.tright = 0.1;
 64:   user.beta   = 2.5;
 65:   PetscOptionsGetReal(NULL,NULL,"-tleft",&user.tleft,NULL);
 66:   PetscOptionsGetReal(NULL,NULL,"-tright",&user.tright,NULL);
 67:   PetscOptionsGetReal(NULL,NULL,"-beta",&user.beta,NULL);
 68:   user.bm1    = user.beta - 1.0;
 69:   user.coef   = user.beta/2.0;

 71:   /*
 72:       Set the DMDA (grid structure) for the grids.
 73:   */
 74:   DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,5,5,5,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);
 75:   DMSetFromOptions(da);
 76:   DMSetUp(da);
 77:   DMSetApplicationContext(da,&user);

 79:   /*
 80:      Create the nonlinear solver
 81:   */
 82:   SNESCreate(PETSC_COMM_WORLD,&snes);
 83:   SNESSetDM(snes,da);
 84:   SNESSetFunction(snes,NULL,FormFunction,&user);
 85:   SNESSetJacobian(snes,NULL,NULL,FormJacobian,&user);
 86:   SNESSetFromOptions(snes);
 87:   SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);

 89:   SNESSolve(snes,NULL,NULL);
 90:   SNESGetIterationNumber(snes,&its);
 91:   SNESGetLinearSolveIterations(snes,&lits);
 92:   litspit = ((PetscReal)lits)/((PetscReal)its);
 93:   PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
 94:   PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
 95:   PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",(double)litspit);

 97:   SNESDestroy(&snes);
 98:   DMDestroy(&da);
 99:   PetscFinalize();
100:   return 0;
101: }
102: /* --------------------  Form initial approximation ----------------- */
103: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
104: {
105:   AppCtx         *user;
106:   PetscInt       i,j,k,xs,ys,xm,ym,zs,zm;
107:   PetscScalar    ***x;
108:   DM             da;

111:   SNESGetDM(snes,&da);
112:   DMGetApplicationContext(da,&user);
113:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
114:   DMDAVecGetArray(da,X,&x);

116:   /* Compute initial guess */
117:   for (k=zs; k<zs+zm; k++) {
118:     for (j=ys; j<ys+ym; j++) {
119:       for (i=xs; i<xs+xm; i++) {
120:         x[k][j][i] = user->tleft;
121:       }
122:     }
123:   }
124:   DMDAVecRestoreArray(da,X,&x);
125:   return 0;
126: }
127: /* --------------------  Evaluate Function F(x) --------------------- */
128: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
129: {
130:   AppCtx         *user = (AppCtx*)ptr;
131:   PetscInt       i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
132:   PetscScalar    zero = 0.0,one = 1.0;
133:   PetscScalar    hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
134:   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;
135:   PetscScalar    tleft,tright,beta,td,ad,dd,fd=0.0,tu,au,du=0.0,fu=0.0;
136:   PetscScalar    ***x,***f;
137:   Vec            localX;
138:   DM             da;

141:   SNESGetDM(snes,&da);
142:   DMGetLocalVector(da,&localX);
143:   DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
144:   hx      = one/(PetscReal)(mx-1);  hy    = one/(PetscReal)(my-1);  hz = one/(PetscReal)(mz-1);
145:   hxhydhz = hx*hy/hz;   hyhzdhx = hy*hz/hx;   hzhxdhy = hz*hx/hy;
146:   tleft   = user->tleft;         tright = user->tright;
147:   beta    = user->beta;

149:   /* Get ghost points */
150:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
151:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
152:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
153:   DMDAVecGetArray(da,localX,&x);
154:   DMDAVecGetArray(da,F,&f);

156:   /* Evaluate function */
157:   for (k=zs; k<zs+zm; k++) {
158:     for (j=ys; j<ys+ym; j++) {
159:       for (i=xs; i<xs+xm; i++) {
160:         t0 = x[k][j][i];

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

164:           /* general interior volume */

166:           tw = x[k][j][i-1];
167:           aw = 0.5*(t0 + tw);
168:           dw = PetscPowScalar(aw,beta);
169:           fw = dw*(t0 - tw);

171:           te = x[k][j][i+1];
172:           ae = 0.5*(t0 + te);
173:           de = PetscPowScalar(ae,beta);
174:           fe = de*(te - t0);

176:           ts = x[k][j-1][i];
177:           as = 0.5*(t0 + ts);
178:           ds = PetscPowScalar(as,beta);
179:           fs = ds*(t0 - ts);

181:           tn = x[k][j+1][i];
182:           an = 0.5*(t0 + tn);
183:           dn = PetscPowScalar(an,beta);
184:           fn = dn*(tn - t0);

186:           td = x[k-1][j][i];
187:           ad = 0.5*(t0 + td);
188:           dd = PetscPowScalar(ad,beta);
189:           fd = dd*(t0 - td);

191:           tu = x[k+1][j][i];
192:           au = 0.5*(t0 + tu);
193:           du = PetscPowScalar(au,beta);
194:           fu = du*(tu - t0);

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

198:           /* left-hand (west) boundary */
199:           tw = tleft;
200:           aw = 0.5*(t0 + tw);
201:           dw = PetscPowScalar(aw,beta);
202:           fw = dw*(t0 - tw);

204:           te = x[k][j][i+1];
205:           ae = 0.5*(t0 + te);
206:           de = PetscPowScalar(ae,beta);
207:           fe = de*(te - t0);

209:           if (j > 0) {
210:             ts = x[k][j-1][i];
211:             as = 0.5*(t0 + ts);
212:             ds = PetscPowScalar(as,beta);
213:             fs = ds*(t0 - ts);
214:           } else {
215:             fs = zero;
216:           }

218:           if (j < my-1) {
219:             tn = x[k][j+1][i];
220:             an = 0.5*(t0 + tn);
221:             dn = PetscPowScalar(an,beta);
222:             fn = dn*(tn - t0);
223:           } else {
224:             fn = zero;
225:           }

227:           if (k > 0) {
228:             td = x[k-1][j][i];
229:             ad = 0.5*(t0 + td);
230:             dd = PetscPowScalar(ad,beta);
231:             fd = dd*(t0 - td);
232:           } else {
233:             fd = zero;
234:           }

236:           if (k < mz-1) {
237:             tu = x[k+1][j][i];
238:             au = 0.5*(t0 + tu);
239:             du = PetscPowScalar(au,beta);
240:             fu = du*(tu - t0);
241:           } else {
242:             fu = zero;
243:           }

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

247:           /* right-hand (east) boundary */
248:           tw = x[k][j][i-1];
249:           aw = 0.5*(t0 + tw);
250:           dw = PetscPowScalar(aw,beta);
251:           fw = dw*(t0 - tw);

253:           te = tright;
254:           ae = 0.5*(t0 + te);
255:           de = PetscPowScalar(ae,beta);
256:           fe = de*(te - t0);

258:           if (j > 0) {
259:             ts = x[k][j-1][i];
260:             as = 0.5*(t0 + ts);
261:             ds = PetscPowScalar(as,beta);
262:             fs = ds*(t0 - ts);
263:           } else {
264:             fs = zero;
265:           }

267:           if (j < my-1) {
268:             tn = x[k][j+1][i];
269:             an = 0.5*(t0 + tn);
270:             dn = PetscPowScalar(an,beta);
271:             fn = dn*(tn - t0);
272:           } else {
273:             fn = zero;
274:           }

276:           if (k > 0) {
277:             td = x[k-1][j][i];
278:             ad = 0.5*(t0 + td);
279:             dd = PetscPowScalar(ad,beta);
280:             fd = dd*(t0 - td);
281:           } else {
282:             fd = zero;
283:           }

285:           if (k < mz-1) {
286:             tu = x[k+1][j][i];
287:             au = 0.5*(t0 + tu);
288:             du = PetscPowScalar(au,beta);
289:             fu = du*(tu - t0);
290:           } else {
291:             fu = zero;
292:           }

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

296:           /* bottom (south) boundary, and i <> 0 or mx-1 */
297:           tw = x[k][j][i-1];
298:           aw = 0.5*(t0 + tw);
299:           dw = PetscPowScalar(aw,beta);
300:           fw = dw*(t0 - tw);

302:           te = x[k][j][i+1];
303:           ae = 0.5*(t0 + te);
304:           de = PetscPowScalar(ae,beta);
305:           fe = de*(te - t0);

307:           fs = zero;

309:           tn = x[k][j+1][i];
310:           an = 0.5*(t0 + tn);
311:           dn = PetscPowScalar(an,beta);
312:           fn = dn*(tn - t0);

314:           if (k > 0) {
315:             td = x[k-1][j][i];
316:             ad = 0.5*(t0 + td);
317:             dd = PetscPowScalar(ad,beta);
318:             fd = dd*(t0 - td);
319:           } else {
320:             fd = zero;
321:           }

323:           if (k < mz-1) {
324:             tu = x[k+1][j][i];
325:             au = 0.5*(t0 + tu);
326:             du = PetscPowScalar(au,beta);
327:             fu = du*(tu - t0);
328:           } else {
329:             fu = zero;
330:           }

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

334:           /* top (north) boundary, and i <> 0 or mx-1 */
335:           tw = x[k][j][i-1];
336:           aw = 0.5*(t0 + tw);
337:           dw = PetscPowScalar(aw,beta);
338:           fw = dw*(t0 - tw);

340:           te = x[k][j][i+1];
341:           ae = 0.5*(t0 + te);
342:           de = PetscPowScalar(ae,beta);
343:           fe = de*(te - t0);

345:           ts = x[k][j-1][i];
346:           as = 0.5*(t0 + ts);
347:           ds = PetscPowScalar(as,beta);
348:           fs = ds*(t0 - ts);

350:           fn = zero;

352:           if (k > 0) {
353:             td = x[k-1][j][i];
354:             ad = 0.5*(t0 + td);
355:             dd = PetscPowScalar(ad,beta);
356:             fd = dd*(t0 - td);
357:           } else {
358:             fd = zero;
359:           }

361:           if (k < mz-1) {
362:             tu = x[k+1][j][i];
363:             au = 0.5*(t0 + tu);
364:             du = PetscPowScalar(au,beta);
365:             fu = du*(tu - t0);
366:           } else {
367:             fu = zero;
368:           }

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

372:           /* down boundary (interior only) */
373:           tw = x[k][j][i-1];
374:           aw = 0.5*(t0 + tw);
375:           dw = PetscPowScalar(aw,beta);
376:           fw = dw*(t0 - tw);

378:           te = x[k][j][i+1];
379:           ae = 0.5*(t0 + te);
380:           de = PetscPowScalar(ae,beta);
381:           fe = de*(te - t0);

383:           ts = x[k][j-1][i];
384:           as = 0.5*(t0 + ts);
385:           ds = PetscPowScalar(as,beta);
386:           fs = ds*(t0 - ts);

388:           tn = x[k][j+1][i];
389:           an = 0.5*(t0 + tn);
390:           dn = PetscPowScalar(an,beta);
391:           fn = dn*(tn - t0);

393:           fd = zero;

395:           tu = x[k+1][j][i];
396:           au = 0.5*(t0 + tu);
397:           du = PetscPowScalar(au,beta);
398:           fu = du*(tu - t0);

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

402:           /* up boundary (interior only) */
403:           tw = x[k][j][i-1];
404:           aw = 0.5*(t0 + tw);
405:           dw = PetscPowScalar(aw,beta);
406:           fw = dw*(t0 - tw);

408:           te = x[k][j][i+1];
409:           ae = 0.5*(t0 + te);
410:           de = PetscPowScalar(ae,beta);
411:           fe = de*(te - t0);

413:           ts = x[k][j-1][i];
414:           as = 0.5*(t0 + ts);
415:           ds = PetscPowScalar(as,beta);
416:           fs = ds*(t0 - ts);

418:           tn = x[k][j+1][i];
419:           an = 0.5*(t0 + tn);
420:           dn = PetscPowScalar(an,beta);
421:           fn = dn*(tn - t0);

423:           td = x[k-1][j][i];
424:           ad = 0.5*(t0 + td);
425:           dd = PetscPowScalar(ad,beta);
426:           fd = dd*(t0 - td);

428:           fu = zero;
429:         }

431:         f[k][j][i] = -hyhzdhx*(fe-fw) - hzhxdhy*(fn-fs) - hxhydhz*(fu-fd);
432:       }
433:     }
434:   }
435:   DMDAVecRestoreArray(da,localX,&x);
436:   DMDAVecRestoreArray(da,F,&f);
437:   DMRestoreLocalVector(da,&localX);
438:   PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
439:   return 0;
440: }
441: /* --------------------  Evaluate Jacobian F(x) --------------------- */
442: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
443: {
444:   AppCtx         *user = (AppCtx*)ptr;
445:   PetscInt       i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
446:   PetscScalar    one = 1.0;
447:   PetscScalar    hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
448:   PetscScalar    t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw;
449:   PetscScalar    tleft,tright,beta,td,ad,dd,tu,au,du,v[7],bm1,coef;
450:   PetscScalar    ***x,bn,bs,be,bw,bu,bd,gn,gs,ge,gw,gu,gd;
451:   Vec            localX;
452:   MatStencil     c[7],row;
453:   DM             da;

456:   SNESGetDM(snes,&da);
457:   DMGetLocalVector(da,&localX);
458:   DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
459:   hx      = one/(PetscReal)(mx-1);  hy    = one/(PetscReal)(my-1);  hz = one/(PetscReal)(mz-1);
460:   hxhydhz = hx*hy/hz;   hyhzdhx = hy*hz/hx;   hzhxdhy = hz*hx/hy;
461:   tleft   = user->tleft;         tright = user->tright;
462:   beta    = user->beta;          bm1    = user->bm1;              coef = user->coef;

464:   /* Get ghost points */
465:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
466:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
467:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
468:   DMDAVecGetArray(da,localX,&x);

470:   /* Evaluate Jacobian of function */
471:   for (k=zs; k<zs+zm; k++) {
472:     for (j=ys; j<ys+ym; j++) {
473:       for (i=xs; i<xs+xm; i++) {
474:         t0    = x[k][j][i];
475:         row.k = k; row.j = j; row.i = i;
476:         if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {

478:           /* general interior volume */

480:           tw = x[k][j][i-1];
481:           aw = 0.5*(t0 + tw);
482:           bw = PetscPowScalar(aw,bm1);
483:           /* dw = bw * aw */
484:           dw = PetscPowScalar(aw,beta);
485:           gw = coef*bw*(t0 - tw);

487:           te = x[k][j][i+1];
488:           ae = 0.5*(t0 + te);
489:           be = PetscPowScalar(ae,bm1);
490:           /* de = be * ae; */
491:           de = PetscPowScalar(ae,beta);
492:           ge = coef*be*(te - t0);

494:           ts = x[k][j-1][i];
495:           as = 0.5*(t0 + ts);
496:           bs = PetscPowScalar(as,bm1);
497:           /* ds = bs * as; */
498:           ds = PetscPowScalar(as,beta);
499:           gs = coef*bs*(t0 - ts);

501:           tn = x[k][j+1][i];
502:           an = 0.5*(t0 + tn);
503:           bn = PetscPowScalar(an,bm1);
504:           /* dn = bn * an; */
505:           dn = PetscPowScalar(an,beta);
506:           gn = coef*bn*(tn - t0);

508:           td = x[k-1][j][i];
509:           ad = 0.5*(t0 + td);
510:           bd = PetscPowScalar(ad,bm1);
511:           /* dd = bd * ad; */
512:           dd = PetscPowScalar(ad,beta);
513:           gd = coef*bd*(t0 - td);

515:           tu = x[k+1][j][i];
516:           au = 0.5*(t0 + tu);
517:           bu = PetscPowScalar(au,bm1);
518:           /* du = bu * au; */
519:           du = PetscPowScalar(au,beta);
520:           gu = coef*bu*(tu - t0);

522:           c[0].k = k-1; c[0].j = j; c[0].i = i; v[0]   = -hxhydhz*(dd - gd);
523:           c[1].k = k; c[1].j = j-1; c[1].i = i;
524:           v[1]   = -hzhxdhy*(ds - gs);
525:           c[2].k = k; c[2].j = j; c[2].i = i-1;
526:           v[2]   = -hyhzdhx*(dw - gw);
527:           c[3].k = k; c[3].j = j; c[3].i = i;
528:           v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
529:           c[4].k = k; c[4].j = j; c[4].i = i+1;
530:           v[4]   = -hyhzdhx*(de + ge);
531:           c[5].k = k; c[5].j = j+1; c[5].i = i;
532:           v[5]   = -hzhxdhy*(dn + gn);
533:           c[6].k = k+1; c[6].j = j; c[6].i = i;
534:           v[6]   = -hxhydhz*(du + gu);
535:           MatSetValuesStencil(jac,1,&row,7,c,v,INSERT_VALUES);

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

539:           /* left-hand plane boundary */
540:           tw = tleft;
541:           aw = 0.5*(t0 + tw);
542:           bw = PetscPowScalar(aw,bm1);
543:           /* dw = bw * aw */
544:           dw = PetscPowScalar(aw,beta);
545:           gw = coef*bw*(t0 - tw);

547:           te = x[k][j][i+1];
548:           ae = 0.5*(t0 + te);
549:           be = PetscPowScalar(ae,bm1);
550:           /* de = be * ae; */
551:           de = PetscPowScalar(ae,beta);
552:           ge = coef*be*(te - t0);

554:           /* left-hand bottom edge */
555:           if (j == 0) {

557:             tn = x[k][j+1][i];
558:             an = 0.5*(t0 + tn);
559:             bn = PetscPowScalar(an,bm1);
560:             /* dn = bn * an; */
561:             dn = PetscPowScalar(an,beta);
562:             gn = coef*bn*(tn - t0);

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

567:               tu = x[k+1][j][i];
568:               au = 0.5*(t0 + tu);
569:               bu = PetscPowScalar(au,bm1);
570:               /* du = bu * au; */
571:               du = PetscPowScalar(au,beta);
572:               gu = coef*bu*(tu - t0);

574:               c[0].k = k; c[0].j = j; c[0].i = i;
575:               v[0]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
576:               c[1].k = k; c[1].j = j; c[1].i = i+1;
577:               v[1]   = -hyhzdhx*(de + ge);
578:               c[2].k = k; c[2].j = j+1; c[2].i = i;
579:               v[2]   = -hzhxdhy*(dn + gn);
580:               c[3].k = k+1; c[3].j = j; c[3].i = i;
581:               v[3]   = -hxhydhz*(du + gu);
582:               MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

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

594:               td = x[k-1][j][i];
595:               ad = 0.5*(t0 + td);
596:               bd = PetscPowScalar(ad,bm1);
597:               /* dd = bd * ad; */
598:               dd = PetscPowScalar(ad,beta);
599:               gd = coef*bd*(td - t0);

601:               c[0].k = k-1; c[0].j = j; c[0].i = i;
602:               v[0]   = -hxhydhz*(dd - gd);
603:               c[1].k = k; c[1].j = j; c[1].i = i;
604:               v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
605:               c[2].k = k; c[2].j = j; c[2].i = i+1;
606:               v[2]   = -hyhzdhx*(de + ge);
607:               c[3].k = k; c[3].j = j+1; c[3].i = i;
608:               v[3]   = -hzhxdhy*(dn + gn);
609:               c[4].k = k+1; c[4].j = j; c[4].i = i;
610:               v[4]   = -hxhydhz*(du + gu);
611:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

613:               /* left-hand bottom up corner */
614:             } else {

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

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

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

637:             ts = x[k][j-1][i];
638:             as = 0.5*(t0 + ts);
639:             bs = PetscPowScalar(as,bm1);
640:             /* ds = bs * as; */
641:             ds = PetscPowScalar(as,beta);
642:             gs = coef*bs*(ts - t0);

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

647:               tu = x[k+1][j][i];
648:               au = 0.5*(t0 + tu);
649:               bu = PetscPowScalar(au,bm1);
650:               /* du = bu * au; */
651:               du = PetscPowScalar(au,beta);
652:               gu = coef*bu*(tu - t0);

654:               c[0].k = k; c[0].j = j-1; c[0].i = i;
655:               v[0]   = -hzhxdhy*(ds - gs);
656:               c[1].k = k; c[1].j = j; c[1].i = i;
657:               v[1]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
658:               c[2].k = k; c[2].j = j; c[2].i = i+1;
659:               v[2]   = -hyhzdhx*(de + ge);
660:               c[3].k = k+1; c[3].j = j; c[3].i = i;
661:               v[3]   = -hxhydhz*(du + gu);
662:               MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

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

674:               td = x[k-1][j][i];
675:               ad = 0.5*(t0 + td);
676:               bd = PetscPowScalar(ad,bm1);
677:               /* dd = bd * ad; */
678:               dd = PetscPowScalar(ad,beta);
679:               gd = coef*bd*(td - t0);

681:               c[0].k = k-1; c[0].j = j; c[0].i = i;
682:               v[0]   = -hxhydhz*(dd - gd);
683:               c[1].k = k; c[1].j = j-1; c[1].i = i;
684:               v[1]   = -hzhxdhy*(ds - gs);
685:               c[2].k = k; c[2].j = j; c[2].i = i;
686:               v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
687:               c[3].k = k; c[3].j = j; c[3].i = i+1;
688:               v[3]   = -hyhzdhx*(de + ge);
689:               c[4].k = k+1; c[4].j = j; c[4].i = i;
690:               v[4]   = -hxhydhz*(du + gu);
691:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

693:               /* left-hand top up corner */
694:             } else {

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

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

714:           } else {

716:             ts = x[k][j-1][i];
717:             as = 0.5*(t0 + ts);
718:             bs = PetscPowScalar(as,bm1);
719:             /* ds = bs * as; */
720:             ds = PetscPowScalar(as,beta);
721:             gs = coef*bs*(t0 - ts);

723:             tn = x[k][j+1][i];
724:             an = 0.5*(t0 + tn);
725:             bn = PetscPowScalar(an,bm1);
726:             /* dn = bn * an; */
727:             dn = PetscPowScalar(an,beta);
728:             gn = coef*bn*(tn - t0);

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

733:               tu = x[k+1][j][i];
734:               au = 0.5*(t0 + tu);
735:               bu = PetscPowScalar(au,bm1);
736:               /* du = bu * au; */
737:               du = PetscPowScalar(au,beta);
738:               gu = coef*bu*(tu - t0);

740:               c[0].k = k; c[0].j = j-1; c[0].i = i;
741:               v[0]   = -hzhxdhy*(ds - gs);
742:               c[1].k = k; c[1].j = j; c[1].i = i;
743:               v[1]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
744:               c[2].k = k; c[2].j = j; c[2].i = i+1;
745:               v[2]   = -hyhzdhx*(de + ge);
746:               c[3].k = k; c[3].j = j+1; c[3].i = i;
747:               v[3]   = -hzhxdhy*(dn + gn);
748:               c[4].k = k+1; c[4].j = j; c[4].i = i;
749:               v[4]   = -hxhydhz*(du + gu);
750:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

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

754:               td = x[k-1][j][i];
755:               ad = 0.5*(t0 + td);
756:               bd = PetscPowScalar(ad,bm1);
757:               /* dd = bd * ad; */
758:               dd = PetscPowScalar(ad,beta);
759:               gd = coef*bd*(t0 - td);

761:               c[0].k = k-1; c[0].j = j; c[0].i = i;
762:               v[0]   = -hxhydhz*(dd - gd);
763:               c[1].k = k; c[1].j = j-1; c[1].i = i;
764:               v[1]   = -hzhxdhy*(ds - gs);
765:               c[2].k = k; c[2].j = j; c[2].i = i;
766:               v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
767:               c[3].k = k; c[3].j = j; c[3].i = i+1;
768:               v[3]   = -hyhzdhx*(de + ge);
769:               c[4].k = k; c[4].j = j+1; c[4].i = i;
770:               v[4]   = -hzhxdhy*(dn + gn);
771:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
772:             } else { /* left-hand interior plane */

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

781:               tu = x[k+1][j][i];
782:               au = 0.5*(t0 + tu);
783:               bu = PetscPowScalar(au,bm1);
784:               /* du = bu * au; */
785:               du = PetscPowScalar(au,beta);
786:               gu = coef*bu*(tu - t0);

788:               c[0].k = k-1; c[0].j = j; c[0].i = i;
789:               v[0]   = -hxhydhz*(dd - gd);
790:               c[1].k = k; c[1].j = j-1; c[1].i = i;
791:               v[1]   = -hzhxdhy*(ds - gs);
792:               c[2].k = k; c[2].j = j; c[2].i = i;
793:               v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
794:               c[3].k = k; c[3].j = j; c[3].i = i+1;
795:               v[3]   = -hyhzdhx*(de + ge);
796:               c[4].k = k; c[4].j = j+1; c[4].i = i;
797:               v[4]   = -hzhxdhy*(dn + gn);
798:               c[5].k = k+1; c[5].j = j; c[5].i = i;
799:               v[5]   = -hxhydhz*(du + gu);
800:               MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
801:             }
802:           }

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

806:           /* right-hand plane boundary */
807:           tw = x[k][j][i-1];
808:           aw = 0.5*(t0 + tw);
809:           bw = PetscPowScalar(aw,bm1);
810:           /* dw = bw * aw */
811:           dw = PetscPowScalar(aw,beta);
812:           gw = coef*bw*(t0 - tw);

814:           te = tright;
815:           ae = 0.5*(t0 + te);
816:           be = PetscPowScalar(ae,bm1);
817:           /* de = be * ae; */
818:           de = PetscPowScalar(ae,beta);
819:           ge = coef*be*(te - t0);

821:           /* right-hand bottom edge */
822:           if (j == 0) {

824:             tn = x[k][j+1][i];
825:             an = 0.5*(t0 + tn);
826:             bn = PetscPowScalar(an,bm1);
827:             /* dn = bn * an; */
828:             dn = PetscPowScalar(an,beta);
829:             gn = coef*bn*(tn - t0);

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

834:               tu = x[k+1][j][i];
835:               au = 0.5*(t0 + tu);
836:               bu = PetscPowScalar(au,bm1);
837:               /* du = bu * au; */
838:               du = PetscPowScalar(au,beta);
839:               gu = coef*bu*(tu - t0);

841:               c[0].k = k; c[0].j = j; c[0].i = i-1;
842:               v[0]   = -hyhzdhx*(dw - gw);
843:               c[1].k = k; c[1].j = j; c[1].i = i;
844:               v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
845:               c[2].k = k; c[2].j = j+1; c[2].i = i;
846:               v[2]   = -hzhxdhy*(dn + gn);
847:               c[3].k = k+1; c[3].j = j; c[3].i = i;
848:               v[3]   = -hxhydhz*(du + gu);
849:               MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

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

861:               td = x[k-1][j][i];
862:               ad = 0.5*(t0 + td);
863:               bd = PetscPowScalar(ad,bm1);
864:               /* dd = bd * ad; */
865:               dd = PetscPowScalar(ad,beta);
866:               gd = coef*bd*(td - t0);

868:               c[0].k = k-1; c[0].j = j; c[0].i = i;
869:               v[0]   = -hxhydhz*(dd - gd);
870:               c[1].k = k; c[1].j = j; c[1].i = i-1;
871:               v[1]   = -hyhzdhx*(dw - gw);
872:               c[2].k = k; c[2].j = j; c[2].i = i;
873:               v[2]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
874:               c[3].k = k; c[3].j = j+1; c[3].i = i;
875:               v[3]   = -hzhxdhy*(dn + gn);
876:               c[4].k = k+1; c[4].j = j; c[4].i = i;
877:               v[4]   = -hxhydhz*(du + gu);
878:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

880:               /* right-hand bottom up corner */
881:             } else {

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

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

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

904:             ts = x[k][j-1][i];
905:             as = 0.5*(t0 + ts);
906:             bs = PetscPowScalar(as,bm1);
907:             /* ds = bs * as; */
908:             ds = PetscPowScalar(as,beta);
909:             gs = coef*bs*(ts - t0);

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

914:               tu = x[k+1][j][i];
915:               au = 0.5*(t0 + tu);
916:               bu = PetscPowScalar(au,bm1);
917:               /* du = bu * au; */
918:               du = PetscPowScalar(au,beta);
919:               gu = coef*bu*(tu - t0);

921:               c[0].k = k; c[0].j = j-1; c[0].i = i;
922:               v[0]   = -hzhxdhy*(ds - gs);
923:               c[1].k = k; c[1].j = j; c[1].i = i-1;
924:               v[1]   = -hyhzdhx*(dw - gw);
925:               c[2].k = k; c[2].j = j; c[2].i = i;
926:               v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
927:               c[3].k = k+1; c[3].j = j; c[3].i = i;
928:               v[3]   = -hxhydhz*(du + gu);
929:               MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

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

941:               td = x[k-1][j][i];
942:               ad = 0.5*(t0 + td);
943:               bd = PetscPowScalar(ad,bm1);
944:               /* dd = bd * ad; */
945:               dd = PetscPowScalar(ad,beta);
946:               gd = coef*bd*(td - t0);

948:               c[0].k = k-1; c[0].j = j; c[0].i = i;
949:               v[0]   = -hxhydhz*(dd - gd);
950:               c[1].k = k; c[1].j = j-1; c[1].i = i;
951:               v[1]   = -hzhxdhy*(ds - gs);
952:               c[2].k = k; c[2].j = j; c[2].i = i-1;
953:               v[2]   = -hyhzdhx*(dw - gw);
954:               c[3].k = k; c[3].j = j; c[3].i = i;
955:               v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
956:               c[4].k = k+1; c[4].j = j; c[4].i = i;
957:               v[4]   = -hxhydhz*(du + gu);
958:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

960:               /* right-hand top up corner */
961:             } else {

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

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

981:           } else {

983:             ts = x[k][j-1][i];
984:             as = 0.5*(t0 + ts);
985:             bs = PetscPowScalar(as,bm1);
986:             /* ds = bs * as; */
987:             ds = PetscPowScalar(as,beta);
988:             gs = coef*bs*(t0 - ts);

990:             tn = x[k][j+1][i];
991:             an = 0.5*(t0 + tn);
992:             bn = PetscPowScalar(an,bm1);
993:             /* dn = bn * an; */
994:             dn = PetscPowScalar(an,beta);
995:             gn = coef*bn*(tn - t0);

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

1000:               tu = x[k+1][j][i];
1001:               au = 0.5*(t0 + tu);
1002:               bu = PetscPowScalar(au,bm1);
1003:               /* du = bu * au; */
1004:               du = PetscPowScalar(au,beta);
1005:               gu = coef*bu*(tu - t0);

1007:               c[0].k = k; c[0].j = j-1; c[0].i = i;
1008:               v[0]   = -hzhxdhy*(ds - gs);
1009:               c[1].k = k; c[1].j = j; c[1].i = i-1;
1010:               v[1]   = -hyhzdhx*(dw - gw);
1011:               c[2].k = k; c[2].j = j; c[2].i = i;
1012:               v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1013:               c[3].k = k; c[3].j = j+1; c[3].i = i;
1014:               v[3]   = -hzhxdhy*(dn + gn);
1015:               c[4].k = k+1; c[4].j = j; c[4].i = i;
1016:               v[4]   = -hxhydhz*(du + gu);
1017:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

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

1021:               td = x[k-1][j][i];
1022:               ad = 0.5*(t0 + td);
1023:               bd = PetscPowScalar(ad,bm1);
1024:               /* dd = bd * ad; */
1025:               dd = PetscPowScalar(ad,beta);
1026:               gd = coef*bd*(t0 - td);

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

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

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

1049:               tu = x[k+1][j][i];
1050:               au = 0.5*(t0 + tu);
1051:               bu = PetscPowScalar(au,bm1);
1052:               /* du = bu * au; */
1053:               du = PetscPowScalar(au,beta);
1054:               gu = coef*bu*(tu - t0);

1056:               c[0].k = k-1; c[0].j = j; c[0].i = i;
1057:               v[0]   = -hxhydhz*(dd - gd);
1058:               c[1].k = k; c[1].j = j-1; c[1].i = i;
1059:               v[1]   = -hzhxdhy*(ds - gs);
1060:               c[2].k = k; c[2].j = j; c[2].i = i-1;
1061:               v[2]   = -hyhzdhx*(dw - gw);
1062:               c[3].k = k; c[3].j = j; c[3].i = i;
1063:               v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1064:               c[4].k = k; c[4].j = j+1; c[4].i = i;
1065:               v[4]   = -hzhxdhy*(dn + gn);
1066:               c[5].k = k+1; c[5].j = j; c[5].i = i;
1067:               v[5]   = -hxhydhz*(du + gu);
1068:               MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1069:             }
1070:           }

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

1074:           tw = x[k][j][i-1];
1075:           aw = 0.5*(t0 + tw);
1076:           bw = PetscPowScalar(aw,bm1);
1077:           /* dw = bw * aw */
1078:           dw = PetscPowScalar(aw,beta);
1079:           gw = coef*bw*(t0 - tw);

1081:           te = x[k][j][i+1];
1082:           ae = 0.5*(t0 + te);
1083:           be = PetscPowScalar(ae,bm1);
1084:           /* de = be * ae; */
1085:           de = PetscPowScalar(ae,beta);
1086:           ge = coef*be*(te - t0);

1088:           tn = x[k][j+1][i];
1089:           an = 0.5*(t0 + tn);
1090:           bn = PetscPowScalar(an,bm1);
1091:           /* dn = bn * an; */
1092:           dn = PetscPowScalar(an,beta);
1093:           gn = coef*bn*(tn - t0);

1095:           /* bottom down interior edge */
1096:           if (k == 0) {

1098:             tu = x[k+1][j][i];
1099:             au = 0.5*(t0 + tu);
1100:             bu = PetscPowScalar(au,bm1);
1101:             /* du = bu * au; */
1102:             du = PetscPowScalar(au,beta);
1103:             gu = coef*bu*(tu - t0);

1105:             c[0].k = k; c[0].j = j; c[0].i = i-1;
1106:             v[0]   = -hyhzdhx*(dw - gw);
1107:             c[1].k = k; c[1].j = j; c[1].i = i;
1108:             v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1109:             c[2].k = k; c[2].j = j; c[2].i = i+1;
1110:             v[2]   = -hyhzdhx*(de + ge);
1111:             c[3].k = k; c[3].j = j+1; c[3].i = i;
1112:             v[3]   = -hzhxdhy*(dn + gn);
1113:             c[4].k = k+1; c[4].j = j; c[4].i = i;
1114:             v[4]   = -hxhydhz*(du + gu);
1115:             MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

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

1119:             td = x[k-1][j][i];
1120:             ad = 0.5*(t0 + td);
1121:             bd = PetscPowScalar(ad,bm1);
1122:             /* dd = bd * ad; */
1123:             dd = PetscPowScalar(ad,beta);
1124:             gd = coef*bd*(td - t0);

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

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

1140:             tu = x[k+1][j][i];
1141:             au = 0.5*(t0 + tu);
1142:             bu = PetscPowScalar(au,bm1);
1143:             /* du = bu * au; */
1144:             du = PetscPowScalar(au,beta);
1145:             gu = coef*bu*(tu - t0);

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

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

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

1171:           tw = x[k][j][i-1];
1172:           aw = 0.5*(t0 + tw);
1173:           bw = PetscPowScalar(aw,bm1);
1174:           /* dw = bw * aw */
1175:           dw = PetscPowScalar(aw,beta);
1176:           gw = coef*bw*(t0 - tw);

1178:           te = x[k][j][i+1];
1179:           ae = 0.5*(t0 + te);
1180:           be = PetscPowScalar(ae,bm1);
1181:           /* de = be * ae; */
1182:           de = PetscPowScalar(ae,beta);
1183:           ge = coef*be*(te - t0);

1185:           ts = x[k][j-1][i];
1186:           as = 0.5*(t0 + ts);
1187:           bs = PetscPowScalar(as,bm1);
1188:           /* ds = bs * as; */
1189:           ds = PetscPowScalar(as,beta);
1190:           gs = coef*bs*(t0 - ts);

1192:           /* top down interior edge */
1193:           if (k == 0) {

1195:             tu = x[k+1][j][i];
1196:             au = 0.5*(t0 + tu);
1197:             bu = PetscPowScalar(au,bm1);
1198:             /* du = bu * au; */
1199:             du = PetscPowScalar(au,beta);
1200:             gu = coef*bu*(tu - t0);

1202:             c[0].k = k; c[0].j = j-1; c[0].i = i;
1203:             v[0]   = -hzhxdhy*(ds - gs);
1204:             c[1].k = k; c[1].j = j; c[1].i = i-1;
1205:             v[1]   = -hyhzdhx*(dw - gw);
1206:             c[2].k = k; c[2].j = j; c[2].i = i;
1207:             v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1208:             c[3].k = k; c[3].j = j; c[3].i = i+1;
1209:             v[3]   = -hyhzdhx*(de + ge);
1210:             c[4].k = k+1; c[4].j = j; c[4].i = i;
1211:             v[4]   = -hxhydhz*(du + gu);
1212:             MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

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

1216:             td = x[k-1][j][i];
1217:             ad = 0.5*(t0 + td);
1218:             bd = PetscPowScalar(ad,bm1);
1219:             /* dd = bd * ad; */
1220:             dd = PetscPowScalar(ad,beta);
1221:             gd = coef*bd*(td - t0);

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

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

1237:             tu = x[k+1][j][i];
1238:             au = 0.5*(t0 + tu);
1239:             bu = PetscPowScalar(au,bm1);
1240:             /* du = bu * au; */
1241:             du = PetscPowScalar(au,beta);
1242:             gu = coef*bu*(tu - t0);

1244:             td = x[k-1][j][i];
1245:             ad = 0.5*(t0 + td);
1246:             bd = PetscPowScalar(ad,bm1);
1247:             /* dd = bd * ad; */
1248:             dd = PetscPowScalar(ad,beta);
1249:             gd = coef*bd*(td - t0);

1251:             c[0].k = k-1; c[0].j = j; c[0].i = i;
1252:             v[0]   = -hxhydhz*(dd - gd);
1253:             c[1].k = k; c[1].j = j-1; c[1].i = i;
1254:             v[1]   = -hzhxdhy*(ds - gs);
1255:             c[2].k = k; c[2].j = j; c[2].i = i-1;
1256:             v[2]   = -hyhzdhx*(dw - gw);
1257:             c[3].k = k; c[3].j = j; c[3].i = i;
1258:             v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1259:             c[4].k = k; c[4].j = j; c[4].i = i+1;
1260:             v[4]   = -hyhzdhx*(de + ge);
1261:             c[5].k = k+1; c[5].j = j; c[5].i = i;
1262:             v[5]   = -hxhydhz*(du + gu);
1263:             MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1264:           }

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

1268:           /* down interior plane */

1270:           tw = x[k][j][i-1];
1271:           aw = 0.5*(t0 + tw);
1272:           bw = PetscPowScalar(aw,bm1);
1273:           /* dw = bw * aw */
1274:           dw = PetscPowScalar(aw,beta);
1275:           gw = coef*bw*(t0 - tw);

1277:           te = x[k][j][i+1];
1278:           ae = 0.5*(t0 + te);
1279:           be = PetscPowScalar(ae,bm1);
1280:           /* de = be * ae; */
1281:           de = PetscPowScalar(ae,beta);
1282:           ge = coef*be*(te - t0);

1284:           ts = x[k][j-1][i];
1285:           as = 0.5*(t0 + ts);
1286:           bs = PetscPowScalar(as,bm1);
1287:           /* ds = bs * as; */
1288:           ds = PetscPowScalar(as,beta);
1289:           gs = coef*bs*(t0 - ts);

1291:           tn = x[k][j+1][i];
1292:           an = 0.5*(t0 + tn);
1293:           bn = PetscPowScalar(an,bm1);
1294:           /* dn = bn * an; */
1295:           dn = PetscPowScalar(an,beta);
1296:           gn = coef*bn*(tn - t0);

1298:           tu = x[k+1][j][i];
1299:           au = 0.5*(t0 + tu);
1300:           bu = PetscPowScalar(au,bm1);
1301:           /* du = bu * au; */
1302:           du = PetscPowScalar(au,beta);
1303:           gu = coef*bu*(tu - t0);

1305:           c[0].k = k; c[0].j = j-1; c[0].i = i;
1306:           v[0]   = -hzhxdhy*(ds - gs);
1307:           c[1].k = k; c[1].j = j; c[1].i = i-1;
1308:           v[1]   = -hyhzdhx*(dw - gw);
1309:           c[2].k = k; c[2].j = j; c[2].i = i;
1310:           v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1311:           c[3].k = k; c[3].j = j; c[3].i = i+1;
1312:           v[3]   = -hyhzdhx*(de + ge);
1313:           c[4].k = k; c[4].j = j+1; c[4].i = i;
1314:           v[4]   = -hzhxdhy*(dn + gn);
1315:           c[5].k = k+1; c[5].j = j; c[5].i = i;
1316:           v[5]   = -hxhydhz*(du + gu);
1317:           MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);

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

1321:           /* up interior plane */

1323:           tw = x[k][j][i-1];
1324:           aw = 0.5*(t0 + tw);
1325:           bw = PetscPowScalar(aw,bm1);
1326:           /* dw = bw * aw */
1327:           dw = PetscPowScalar(aw,beta);
1328:           gw = coef*bw*(t0 - tw);

1330:           te = x[k][j][i+1];
1331:           ae = 0.5*(t0 + te);
1332:           be = PetscPowScalar(ae,bm1);
1333:           /* de = be * ae; */
1334:           de = PetscPowScalar(ae,beta);
1335:           ge = coef*be*(te - t0);

1337:           ts = x[k][j-1][i];
1338:           as = 0.5*(t0 + ts);
1339:           bs = PetscPowScalar(as,bm1);
1340:           /* ds = bs * as; */
1341:           ds = PetscPowScalar(as,beta);
1342:           gs = coef*bs*(t0 - ts);

1344:           tn = x[k][j+1][i];
1345:           an = 0.5*(t0 + tn);
1346:           bn = PetscPowScalar(an,bm1);
1347:           /* dn = bn * an; */
1348:           dn = PetscPowScalar(an,beta);
1349:           gn = coef*bn*(tn - t0);

1351:           td = x[k-1][j][i];
1352:           ad = 0.5*(t0 + td);
1353:           bd = PetscPowScalar(ad,bm1);
1354:           /* dd = bd * ad; */
1355:           dd = PetscPowScalar(ad,beta);
1356:           gd = coef*bd*(t0 - td);

1358:           c[0].k = k-1; c[0].j = j; c[0].i = i;
1359:           v[0]   = -hxhydhz*(dd - gd);
1360:           c[1].k = k; c[1].j = j-1; c[1].i = i;
1361:           v[1]   = -hzhxdhy*(ds - gs);
1362:           c[2].k = k; c[2].j = j; c[2].i = i-1;
1363:           v[2]   = -hyhzdhx*(dw - gw);
1364:           c[3].k = k; c[3].j = j; c[3].i = i;
1365:           v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1366:           c[4].k = k; c[4].j = j; c[4].i = i+1;
1367:           v[4]   = -hyhzdhx*(de + ge);
1368:           c[5].k = k; c[5].j = j+1; c[5].i = i;
1369:           v[5]   = -hzhxdhy*(dn + gn);
1370:           MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1371:         }
1372:       }
1373:     }
1374:   }
1375:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
1376:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
1377:   if (jac != J) {
1378:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1379:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1380:   }
1381:   DMDAVecRestoreArray(da,localX,&x);
1382:   DMRestoreLocalVector(da,&localX);

1384:   PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
1385:   return 0;
1386: }

1388: /*TEST

1390:    test:
1391:       nsize: 4
1392:       args: -snes_monitor_short -pc_mg_type full -ksp_type fgmres -pc_type mg -snes_view -pc_mg_levels 2 -pc_mg_galerkin pmat
1393:       requires: !single

1395: TEST*/