Actual source code: ex20.c


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

 14: /*T
 15:    Concepts: SNES^solving a system of nonlinear equations
 16:    Concepts: 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 <petscdm.h>
 44: #include <petscdmda.h>

 46: /* User-defined application context */

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

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

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

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

 67:   PetscInitialize(&argc,&argv,NULL,help);
 68:   /* set problem parameters */
 69:   user.tleft  = 1.0;
 70:   user.tright = 0.1;
 71:   user.beta   = 2.5;
 72:   PetscOptionsGetReal(NULL,NULL,"-tleft",&user.tleft,NULL);
 73:   PetscOptionsGetReal(NULL,NULL,"-tright",&user.tright,NULL);
 74:   PetscOptionsGetReal(NULL,NULL,"-beta",&user.beta,NULL);
 75:   user.bm1    = user.beta - 1.0;
 76:   user.coef   = user.beta/2.0;

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

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

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

104:   SNESDestroy(&snes);
105:   DMDestroy(&da);
106:   PetscFinalize();
107:   return 0;
108: }
109: /* --------------------  Form initial approximation ----------------- */
110: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
111: {
112:   AppCtx         *user;
113:   PetscInt       i,j,k,xs,ys,xm,ym,zs,zm;
114:   PetscScalar    ***x;
115:   DM             da;

118:   SNESGetDM(snes,&da);
119:   DMGetApplicationContext(da,&user);
120:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
121:   DMDAVecGetArray(da,X,&x);

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

148:   SNESGetDM(snes,&da);
149:   DMGetLocalVector(da,&localX);
150:   DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
151:   hx      = one/(PetscReal)(mx-1);  hy    = one/(PetscReal)(my-1);  hz = one/(PetscReal)(mz-1);
152:   hxhydhz = hx*hy/hz;   hyhzdhx = hy*hz/hx;   hzhxdhy = hz*hx/hy;
153:   tleft   = user->tleft;         tright = user->tright;
154:   beta    = user->beta;

156:   /* Get ghost points */
157:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
158:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
159:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
160:   DMDAVecGetArray(da,localX,&x);
161:   DMDAVecGetArray(da,F,&f);

163:   /* Evaluate function */
164:   for (k=zs; k<zs+zm; k++) {
165:     for (j=ys; j<ys+ym; j++) {
166:       for (i=xs; i<xs+xm; i++) {
167:         t0 = x[k][j][i];

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

171:           /* general interior volume */

173:           tw = x[k][j][i-1];
174:           aw = 0.5*(t0 + tw);
175:           dw = PetscPowScalar(aw,beta);
176:           fw = dw*(t0 - tw);

178:           te = x[k][j][i+1];
179:           ae = 0.5*(t0 + te);
180:           de = PetscPowScalar(ae,beta);
181:           fe = de*(te - t0);

183:           ts = x[k][j-1][i];
184:           as = 0.5*(t0 + ts);
185:           ds = PetscPowScalar(as,beta);
186:           fs = ds*(t0 - ts);

188:           tn = x[k][j+1][i];
189:           an = 0.5*(t0 + tn);
190:           dn = PetscPowScalar(an,beta);
191:           fn = dn*(tn - t0);

193:           td = x[k-1][j][i];
194:           ad = 0.5*(t0 + td);
195:           dd = PetscPowScalar(ad,beta);
196:           fd = dd*(t0 - td);

198:           tu = x[k+1][j][i];
199:           au = 0.5*(t0 + tu);
200:           du = PetscPowScalar(au,beta);
201:           fu = du*(tu - t0);

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

205:           /* left-hand (west) boundary */
206:           tw = tleft;
207:           aw = 0.5*(t0 + tw);
208:           dw = PetscPowScalar(aw,beta);
209:           fw = dw*(t0 - tw);

211:           te = x[k][j][i+1];
212:           ae = 0.5*(t0 + te);
213:           de = PetscPowScalar(ae,beta);
214:           fe = de*(te - t0);

216:           if (j > 0) {
217:             ts = x[k][j-1][i];
218:             as = 0.5*(t0 + ts);
219:             ds = PetscPowScalar(as,beta);
220:             fs = ds*(t0 - ts);
221:           } else {
222:             fs = zero;
223:           }

225:           if (j < my-1) {
226:             tn = x[k][j+1][i];
227:             an = 0.5*(t0 + tn);
228:             dn = PetscPowScalar(an,beta);
229:             fn = dn*(tn - t0);
230:           } else {
231:             fn = zero;
232:           }

234:           if (k > 0) {
235:             td = x[k-1][j][i];
236:             ad = 0.5*(t0 + td);
237:             dd = PetscPowScalar(ad,beta);
238:             fd = dd*(t0 - td);
239:           } else {
240:             fd = zero;
241:           }

243:           if (k < mz-1) {
244:             tu = x[k+1][j][i];
245:             au = 0.5*(t0 + tu);
246:             du = PetscPowScalar(au,beta);
247:             fu = du*(tu - t0);
248:           } else {
249:             fu = zero;
250:           }

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

254:           /* right-hand (east) boundary */
255:           tw = x[k][j][i-1];
256:           aw = 0.5*(t0 + tw);
257:           dw = PetscPowScalar(aw,beta);
258:           fw = dw*(t0 - tw);

260:           te = tright;
261:           ae = 0.5*(t0 + te);
262:           de = PetscPowScalar(ae,beta);
263:           fe = de*(te - t0);

265:           if (j > 0) {
266:             ts = x[k][j-1][i];
267:             as = 0.5*(t0 + ts);
268:             ds = PetscPowScalar(as,beta);
269:             fs = ds*(t0 - ts);
270:           } else {
271:             fs = zero;
272:           }

274:           if (j < my-1) {
275:             tn = x[k][j+1][i];
276:             an = 0.5*(t0 + tn);
277:             dn = PetscPowScalar(an,beta);
278:             fn = dn*(tn - t0);
279:           } else {
280:             fn = zero;
281:           }

283:           if (k > 0) {
284:             td = x[k-1][j][i];
285:             ad = 0.5*(t0 + td);
286:             dd = PetscPowScalar(ad,beta);
287:             fd = dd*(t0 - td);
288:           } else {
289:             fd = zero;
290:           }

292:           if (k < mz-1) {
293:             tu = x[k+1][j][i];
294:             au = 0.5*(t0 + tu);
295:             du = PetscPowScalar(au,beta);
296:             fu = du*(tu - t0);
297:           } else {
298:             fu = zero;
299:           }

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

303:           /* bottom (south) boundary, and i <> 0 or mx-1 */
304:           tw = x[k][j][i-1];
305:           aw = 0.5*(t0 + tw);
306:           dw = PetscPowScalar(aw,beta);
307:           fw = dw*(t0 - tw);

309:           te = x[k][j][i+1];
310:           ae = 0.5*(t0 + te);
311:           de = PetscPowScalar(ae,beta);
312:           fe = de*(te - t0);

314:           fs = zero;

316:           tn = x[k][j+1][i];
317:           an = 0.5*(t0 + tn);
318:           dn = PetscPowScalar(an,beta);
319:           fn = dn*(tn - t0);

321:           if (k > 0) {
322:             td = x[k-1][j][i];
323:             ad = 0.5*(t0 + td);
324:             dd = PetscPowScalar(ad,beta);
325:             fd = dd*(t0 - td);
326:           } else {
327:             fd = zero;
328:           }

330:           if (k < mz-1) {
331:             tu = x[k+1][j][i];
332:             au = 0.5*(t0 + tu);
333:             du = PetscPowScalar(au,beta);
334:             fu = du*(tu - t0);
335:           } else {
336:             fu = zero;
337:           }

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

341:           /* top (north) boundary, and i <> 0 or mx-1 */
342:           tw = x[k][j][i-1];
343:           aw = 0.5*(t0 + tw);
344:           dw = PetscPowScalar(aw,beta);
345:           fw = dw*(t0 - tw);

347:           te = x[k][j][i+1];
348:           ae = 0.5*(t0 + te);
349:           de = PetscPowScalar(ae,beta);
350:           fe = de*(te - t0);

352:           ts = x[k][j-1][i];
353:           as = 0.5*(t0 + ts);
354:           ds = PetscPowScalar(as,beta);
355:           fs = ds*(t0 - ts);

357:           fn = zero;

359:           if (k > 0) {
360:             td = x[k-1][j][i];
361:             ad = 0.5*(t0 + td);
362:             dd = PetscPowScalar(ad,beta);
363:             fd = dd*(t0 - td);
364:           } else {
365:             fd = zero;
366:           }

368:           if (k < mz-1) {
369:             tu = x[k+1][j][i];
370:             au = 0.5*(t0 + tu);
371:             du = PetscPowScalar(au,beta);
372:             fu = du*(tu - t0);
373:           } else {
374:             fu = zero;
375:           }

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

379:           /* down boundary (interior only) */
380:           tw = x[k][j][i-1];
381:           aw = 0.5*(t0 + tw);
382:           dw = PetscPowScalar(aw,beta);
383:           fw = dw*(t0 - tw);

385:           te = x[k][j][i+1];
386:           ae = 0.5*(t0 + te);
387:           de = PetscPowScalar(ae,beta);
388:           fe = de*(te - t0);

390:           ts = x[k][j-1][i];
391:           as = 0.5*(t0 + ts);
392:           ds = PetscPowScalar(as,beta);
393:           fs = ds*(t0 - ts);

395:           tn = x[k][j+1][i];
396:           an = 0.5*(t0 + tn);
397:           dn = PetscPowScalar(an,beta);
398:           fn = dn*(tn - t0);

400:           fd = zero;

402:           tu = x[k+1][j][i];
403:           au = 0.5*(t0 + tu);
404:           du = PetscPowScalar(au,beta);
405:           fu = du*(tu - t0);

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

409:           /* up boundary (interior only) */
410:           tw = x[k][j][i-1];
411:           aw = 0.5*(t0 + tw);
412:           dw = PetscPowScalar(aw,beta);
413:           fw = dw*(t0 - tw);

415:           te = x[k][j][i+1];
416:           ae = 0.5*(t0 + te);
417:           de = PetscPowScalar(ae,beta);
418:           fe = de*(te - t0);

420:           ts = x[k][j-1][i];
421:           as = 0.5*(t0 + ts);
422:           ds = PetscPowScalar(as,beta);
423:           fs = ds*(t0 - ts);

425:           tn = x[k][j+1][i];
426:           an = 0.5*(t0 + tn);
427:           dn = PetscPowScalar(an,beta);
428:           fn = dn*(tn - t0);

430:           td = x[k-1][j][i];
431:           ad = 0.5*(t0 + td);
432:           dd = PetscPowScalar(ad,beta);
433:           fd = dd*(t0 - td);

435:           fu = zero;
436:         }

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

463:   SNESGetDM(snes,&da);
464:   DMGetLocalVector(da,&localX);
465:   DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
466:   hx      = one/(PetscReal)(mx-1);  hy    = one/(PetscReal)(my-1);  hz = one/(PetscReal)(mz-1);
467:   hxhydhz = hx*hy/hz;   hyhzdhx = hy*hz/hx;   hzhxdhy = hz*hx/hy;
468:   tleft   = user->tleft;         tright = user->tright;
469:   beta    = user->beta;          bm1    = user->bm1;              coef = user->coef;

471:   /* Get ghost points */
472:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
473:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
474:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
475:   DMDAVecGetArray(da,localX,&x);

477:   /* Evaluate Jacobian of function */
478:   for (k=zs; k<zs+zm; k++) {
479:     for (j=ys; j<ys+ym; j++) {
480:       for (i=xs; i<xs+xm; i++) {
481:         t0    = x[k][j][i];
482:         row.k = k; row.j = j; row.i = i;
483:         if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {

485:           /* general interior volume */

487:           tw = x[k][j][i-1];
488:           aw = 0.5*(t0 + tw);
489:           bw = PetscPowScalar(aw,bm1);
490:           /* dw = bw * aw */
491:           dw = PetscPowScalar(aw,beta);
492:           gw = coef*bw*(t0 - tw);

494:           te = x[k][j][i+1];
495:           ae = 0.5*(t0 + te);
496:           be = PetscPowScalar(ae,bm1);
497:           /* de = be * ae; */
498:           de = PetscPowScalar(ae,beta);
499:           ge = coef*be*(te - t0);

501:           ts = x[k][j-1][i];
502:           as = 0.5*(t0 + ts);
503:           bs = PetscPowScalar(as,bm1);
504:           /* ds = bs * as; */
505:           ds = PetscPowScalar(as,beta);
506:           gs = coef*bs*(t0 - ts);

508:           tn = x[k][j+1][i];
509:           an = 0.5*(t0 + tn);
510:           bn = PetscPowScalar(an,bm1);
511:           /* dn = bn * an; */
512:           dn = PetscPowScalar(an,beta);
513:           gn = coef*bn*(tn - t0);

515:           td = x[k-1][j][i];
516:           ad = 0.5*(t0 + td);
517:           bd = PetscPowScalar(ad,bm1);
518:           /* dd = bd * ad; */
519:           dd = PetscPowScalar(ad,beta);
520:           gd = coef*bd*(t0 - td);

522:           tu = x[k+1][j][i];
523:           au = 0.5*(t0 + tu);
524:           bu = PetscPowScalar(au,bm1);
525:           /* du = bu * au; */
526:           du = PetscPowScalar(au,beta);
527:           gu = coef*bu*(tu - t0);

529:           c[0].k = k-1; c[0].j = j; c[0].i = i; v[0]   = -hxhydhz*(dd - gd);
530:           c[1].k = k; c[1].j = j-1; c[1].i = i;
531:           v[1]   = -hzhxdhy*(ds - gs);
532:           c[2].k = k; c[2].j = j; c[2].i = i-1;
533:           v[2]   = -hyhzdhx*(dw - gw);
534:           c[3].k = k; c[3].j = j; c[3].i = i;
535:           v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
536:           c[4].k = k; c[4].j = j; c[4].i = i+1;
537:           v[4]   = -hyhzdhx*(de + ge);
538:           c[5].k = k; c[5].j = j+1; c[5].i = i;
539:           v[5]   = -hzhxdhy*(dn + gn);
540:           c[6].k = k+1; c[6].j = j; c[6].i = i;
541:           v[6]   = -hxhydhz*(du + gu);
542:           MatSetValuesStencil(jac,1,&row,7,c,v,INSERT_VALUES);

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

546:           /* left-hand plane boundary */
547:           tw = tleft;
548:           aw = 0.5*(t0 + tw);
549:           bw = PetscPowScalar(aw,bm1);
550:           /* dw = bw * aw */
551:           dw = PetscPowScalar(aw,beta);
552:           gw = coef*bw*(t0 - tw);

554:           te = x[k][j][i+1];
555:           ae = 0.5*(t0 + te);
556:           be = PetscPowScalar(ae,bm1);
557:           /* de = be * ae; */
558:           de = PetscPowScalar(ae,beta);
559:           ge = coef*be*(te - t0);

561:           /* left-hand bottom edge */
562:           if (j == 0) {

564:             tn = x[k][j+1][i];
565:             an = 0.5*(t0 + tn);
566:             bn = PetscPowScalar(an,bm1);
567:             /* dn = bn * an; */
568:             dn = PetscPowScalar(an,beta);
569:             gn = coef*bn*(tn - t0);

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

574:               tu = x[k+1][j][i];
575:               au = 0.5*(t0 + tu);
576:               bu = PetscPowScalar(au,bm1);
577:               /* du = bu * au; */
578:               du = PetscPowScalar(au,beta);
579:               gu = coef*bu*(tu - t0);

581:               c[0].k = k; c[0].j = j; c[0].i = i;
582:               v[0]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
583:               c[1].k = k; c[1].j = j; c[1].i = i+1;
584:               v[1]   = -hyhzdhx*(de + ge);
585:               c[2].k = k; c[2].j = j+1; c[2].i = i;
586:               v[2]   = -hzhxdhy*(dn + gn);
587:               c[3].k = k+1; c[3].j = j; c[3].i = i;
588:               v[3]   = -hxhydhz*(du + gu);
589:               MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

594:               tu = x[k+1][j][i];
595:               au = 0.5*(t0 + tu);
596:               bu = PetscPowScalar(au,bm1);
597:               /* du = bu * au; */
598:               du = PetscPowScalar(au,beta);
599:               gu = coef*bu*(tu - t0);

601:               td = x[k-1][j][i];
602:               ad = 0.5*(t0 + td);
603:               bd = PetscPowScalar(ad,bm1);
604:               /* dd = bd * ad; */
605:               dd = PetscPowScalar(ad,beta);
606:               gd = coef*bd*(td - t0);

608:               c[0].k = k-1; c[0].j = j; c[0].i = i;
609:               v[0]   = -hxhydhz*(dd - gd);
610:               c[1].k = k; c[1].j = j; c[1].i = i;
611:               v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
612:               c[2].k = k; c[2].j = j; c[2].i = i+1;
613:               v[2]   = -hyhzdhx*(de + ge);
614:               c[3].k = k; c[3].j = j+1; c[3].i = i;
615:               v[3]   = -hzhxdhy*(dn + gn);
616:               c[4].k = k+1; c[4].j = j; c[4].i = i;
617:               v[4]   = -hxhydhz*(du + gu);
618:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

620:               /* left-hand bottom up corner */
621:             } else {

623:               td = x[k-1][j][i];
624:               ad = 0.5*(t0 + td);
625:               bd = PetscPowScalar(ad,bm1);
626:               /* dd = bd * ad; */
627:               dd = PetscPowScalar(ad,beta);
628:               gd = coef*bd*(td - t0);

630:               c[0].k = k-1; c[0].j = j; c[0].i = i;
631:               v[0]   = -hxhydhz*(dd - gd);
632:               c[1].k = k; c[1].j = j; c[1].i = i;
633:               v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
634:               c[2].k = k; c[2].j = j; c[2].i = i+1;
635:               v[2]   = -hyhzdhx*(de + ge);
636:               c[3].k = k; c[3].j = j+1; c[3].i = i;
637:               v[3]   = -hzhxdhy*(dn + gn);
638:               MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
639:             }

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

644:             ts = x[k][j-1][i];
645:             as = 0.5*(t0 + ts);
646:             bs = PetscPowScalar(as,bm1);
647:             /* ds = bs * as; */
648:             ds = PetscPowScalar(as,beta);
649:             gs = coef*bs*(ts - t0);

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

654:               tu = x[k+1][j][i];
655:               au = 0.5*(t0 + tu);
656:               bu = PetscPowScalar(au,bm1);
657:               /* du = bu * au; */
658:               du = PetscPowScalar(au,beta);
659:               gu = coef*bu*(tu - t0);

661:               c[0].k = k; c[0].j = j-1; c[0].i = i;
662:               v[0]   = -hzhxdhy*(ds - gs);
663:               c[1].k = k; c[1].j = j; c[1].i = i;
664:               v[1]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
665:               c[2].k = k; c[2].j = j; c[2].i = i+1;
666:               v[2]   = -hyhzdhx*(de + ge);
667:               c[3].k = k+1; c[3].j = j; c[3].i = i;
668:               v[3]   = -hxhydhz*(du + gu);
669:               MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

674:               tu = x[k+1][j][i];
675:               au = 0.5*(t0 + tu);
676:               bu = PetscPowScalar(au,bm1);
677:               /* du = bu * au; */
678:               du = PetscPowScalar(au,beta);
679:               gu = coef*bu*(tu - t0);

681:               td = x[k-1][j][i];
682:               ad = 0.5*(t0 + td);
683:               bd = PetscPowScalar(ad,bm1);
684:               /* dd = bd * ad; */
685:               dd = PetscPowScalar(ad,beta);
686:               gd = coef*bd*(td - t0);

688:               c[0].k = k-1; c[0].j = j; c[0].i = i;
689:               v[0]   = -hxhydhz*(dd - gd);
690:               c[1].k = k; c[1].j = j-1; c[1].i = i;
691:               v[1]   = -hzhxdhy*(ds - gs);
692:               c[2].k = k; c[2].j = j; c[2].i = i;
693:               v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
694:               c[3].k = k; c[3].j = j; c[3].i = i+1;
695:               v[3]   = -hyhzdhx*(de + ge);
696:               c[4].k = k+1; c[4].j = j; c[4].i = i;
697:               v[4]   = -hxhydhz*(du + gu);
698:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

700:               /* left-hand top up corner */
701:             } else {

703:               td = x[k-1][j][i];
704:               ad = 0.5*(t0 + td);
705:               bd = PetscPowScalar(ad,bm1);
706:               /* dd = bd * ad; */
707:               dd = PetscPowScalar(ad,beta);
708:               gd = coef*bd*(td - t0);

710:               c[0].k = k-1; c[0].j = j; c[0].i = i;
711:               v[0]   = -hxhydhz*(dd - gd);
712:               c[1].k = k; c[1].j = j-1; c[1].i = i;
713:               v[1]   = -hzhxdhy*(ds - gs);
714:               c[2].k = k; c[2].j = j; c[2].i = i;
715:               v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
716:               c[3].k = k; c[3].j = j; c[3].i = i+1;
717:               v[3]   = -hyhzdhx*(de + ge);
718:               MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
719:             }

721:           } else {

723:             ts = x[k][j-1][i];
724:             as = 0.5*(t0 + ts);
725:             bs = PetscPowScalar(as,bm1);
726:             /* ds = bs * as; */
727:             ds = PetscPowScalar(as,beta);
728:             gs = coef*bs*(t0 - ts);

730:             tn = x[k][j+1][i];
731:             an = 0.5*(t0 + tn);
732:             bn = PetscPowScalar(an,bm1);
733:             /* dn = bn * an; */
734:             dn = PetscPowScalar(an,beta);
735:             gn = coef*bn*(tn - t0);

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

740:               tu = x[k+1][j][i];
741:               au = 0.5*(t0 + tu);
742:               bu = PetscPowScalar(au,bm1);
743:               /* du = bu * au; */
744:               du = PetscPowScalar(au,beta);
745:               gu = coef*bu*(tu - t0);

747:               c[0].k = k; c[0].j = j-1; c[0].i = i;
748:               v[0]   = -hzhxdhy*(ds - gs);
749:               c[1].k = k; c[1].j = j; c[1].i = i;
750:               v[1]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
751:               c[2].k = k; c[2].j = j; c[2].i = i+1;
752:               v[2]   = -hyhzdhx*(de + ge);
753:               c[3].k = k; c[3].j = j+1; c[3].i = i;
754:               v[3]   = -hzhxdhy*(dn + gn);
755:               c[4].k = k+1; c[4].j = j; c[4].i = i;
756:               v[4]   = -hxhydhz*(du + gu);
757:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

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

761:               td = x[k-1][j][i];
762:               ad = 0.5*(t0 + td);
763:               bd = PetscPowScalar(ad,bm1);
764:               /* dd = bd * ad; */
765:               dd = PetscPowScalar(ad,beta);
766:               gd = coef*bd*(t0 - td);

768:               c[0].k = k-1; c[0].j = j; c[0].i = i;
769:               v[0]   = -hxhydhz*(dd - gd);
770:               c[1].k = k; c[1].j = j-1; c[1].i = i;
771:               v[1]   = -hzhxdhy*(ds - gs);
772:               c[2].k = k; c[2].j = j; c[2].i = i;
773:               v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
774:               c[3].k = k; c[3].j = j; c[3].i = i+1;
775:               v[3]   = -hyhzdhx*(de + ge);
776:               c[4].k = k; c[4].j = j+1; c[4].i = i;
777:               v[4]   = -hzhxdhy*(dn + gn);
778:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
779:             } else { /* left-hand interior plane */

781:               td = x[k-1][j][i];
782:               ad = 0.5*(t0 + td);
783:               bd = PetscPowScalar(ad,bm1);
784:               /* dd = bd * ad; */
785:               dd = PetscPowScalar(ad,beta);
786:               gd = coef*bd*(t0 - td);

788:               tu = x[k+1][j][i];
789:               au = 0.5*(t0 + tu);
790:               bu = PetscPowScalar(au,bm1);
791:               /* du = bu * au; */
792:               du = PetscPowScalar(au,beta);
793:               gu = coef*bu*(tu - t0);

795:               c[0].k = k-1; c[0].j = j; c[0].i = i;
796:               v[0]   = -hxhydhz*(dd - gd);
797:               c[1].k = k; c[1].j = j-1; c[1].i = i;
798:               v[1]   = -hzhxdhy*(ds - gs);
799:               c[2].k = k; c[2].j = j; c[2].i = i;
800:               v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
801:               c[3].k = k; c[3].j = j; c[3].i = i+1;
802:               v[3]   = -hyhzdhx*(de + ge);
803:               c[4].k = k; c[4].j = j+1; c[4].i = i;
804:               v[4]   = -hzhxdhy*(dn + gn);
805:               c[5].k = k+1; c[5].j = j; c[5].i = i;
806:               v[5]   = -hxhydhz*(du + gu);
807:               MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
808:             }
809:           }

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

813:           /* right-hand plane boundary */
814:           tw = x[k][j][i-1];
815:           aw = 0.5*(t0 + tw);
816:           bw = PetscPowScalar(aw,bm1);
817:           /* dw = bw * aw */
818:           dw = PetscPowScalar(aw,beta);
819:           gw = coef*bw*(t0 - tw);

821:           te = tright;
822:           ae = 0.5*(t0 + te);
823:           be = PetscPowScalar(ae,bm1);
824:           /* de = be * ae; */
825:           de = PetscPowScalar(ae,beta);
826:           ge = coef*be*(te - t0);

828:           /* right-hand bottom edge */
829:           if (j == 0) {

831:             tn = x[k][j+1][i];
832:             an = 0.5*(t0 + tn);
833:             bn = PetscPowScalar(an,bm1);
834:             /* dn = bn * an; */
835:             dn = PetscPowScalar(an,beta);
836:             gn = coef*bn*(tn - t0);

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

841:               tu = x[k+1][j][i];
842:               au = 0.5*(t0 + tu);
843:               bu = PetscPowScalar(au,bm1);
844:               /* du = bu * au; */
845:               du = PetscPowScalar(au,beta);
846:               gu = coef*bu*(tu - t0);

848:               c[0].k = k; c[0].j = j; c[0].i = i-1;
849:               v[0]   = -hyhzdhx*(dw - gw);
850:               c[1].k = k; c[1].j = j; c[1].i = i;
851:               v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
852:               c[2].k = k; c[2].j = j+1; c[2].i = i;
853:               v[2]   = -hzhxdhy*(dn + gn);
854:               c[3].k = k+1; c[3].j = j; c[3].i = i;
855:               v[3]   = -hxhydhz*(du + gu);
856:               MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

861:               tu = x[k+1][j][i];
862:               au = 0.5*(t0 + tu);
863:               bu = PetscPowScalar(au,bm1);
864:               /* du = bu * au; */
865:               du = PetscPowScalar(au,beta);
866:               gu = coef*bu*(tu - t0);

868:               td = x[k-1][j][i];
869:               ad = 0.5*(t0 + td);
870:               bd = PetscPowScalar(ad,bm1);
871:               /* dd = bd * ad; */
872:               dd = PetscPowScalar(ad,beta);
873:               gd = coef*bd*(td - t0);

875:               c[0].k = k-1; c[0].j = j; c[0].i = i;
876:               v[0]   = -hxhydhz*(dd - gd);
877:               c[1].k = k; c[1].j = j; c[1].i = i-1;
878:               v[1]   = -hyhzdhx*(dw - gw);
879:               c[2].k = k; c[2].j = j; c[2].i = i;
880:               v[2]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
881:               c[3].k = k; c[3].j = j+1; c[3].i = i;
882:               v[3]   = -hzhxdhy*(dn + gn);
883:               c[4].k = k+1; c[4].j = j; c[4].i = i;
884:               v[4]   = -hxhydhz*(du + gu);
885:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

887:               /* right-hand bottom up corner */
888:             } else {

890:               td = x[k-1][j][i];
891:               ad = 0.5*(t0 + td);
892:               bd = PetscPowScalar(ad,bm1);
893:               /* dd = bd * ad; */
894:               dd = PetscPowScalar(ad,beta);
895:               gd = coef*bd*(td - t0);

897:               c[0].k = k-1; c[0].j = j; c[0].i = i;
898:               v[0]   = -hxhydhz*(dd - gd);
899:               c[1].k = k; c[1].j = j; c[1].i = i-1;
900:               v[1]   = -hyhzdhx*(dw - gw);
901:               c[2].k = k; c[2].j = j; c[2].i = i;
902:               v[2]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
903:               c[3].k = k; c[3].j = j+1; c[3].i = i;
904:               v[3]   = -hzhxdhy*(dn + gn);
905:               MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
906:             }

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

911:             ts = x[k][j-1][i];
912:             as = 0.5*(t0 + ts);
913:             bs = PetscPowScalar(as,bm1);
914:             /* ds = bs * as; */
915:             ds = PetscPowScalar(as,beta);
916:             gs = coef*bs*(ts - t0);

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

921:               tu = x[k+1][j][i];
922:               au = 0.5*(t0 + tu);
923:               bu = PetscPowScalar(au,bm1);
924:               /* du = bu * au; */
925:               du = PetscPowScalar(au,beta);
926:               gu = coef*bu*(tu - t0);

928:               c[0].k = k; c[0].j = j-1; c[0].i = i;
929:               v[0]   = -hzhxdhy*(ds - gs);
930:               c[1].k = k; c[1].j = j; c[1].i = i-1;
931:               v[1]   = -hyhzdhx*(dw - gw);
932:               c[2].k = k; c[2].j = j; c[2].i = i;
933:               v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
934:               c[3].k = k+1; c[3].j = j; c[3].i = i;
935:               v[3]   = -hxhydhz*(du + gu);
936:               MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

941:               tu = x[k+1][j][i];
942:               au = 0.5*(t0 + tu);
943:               bu = PetscPowScalar(au,bm1);
944:               /* du = bu * au; */
945:               du = PetscPowScalar(au,beta);
946:               gu = coef*bu*(tu - t0);

948:               td = x[k-1][j][i];
949:               ad = 0.5*(t0 + td);
950:               bd = PetscPowScalar(ad,bm1);
951:               /* dd = bd * ad; */
952:               dd = PetscPowScalar(ad,beta);
953:               gd = coef*bd*(td - t0);

955:               c[0].k = k-1; c[0].j = j; c[0].i = i;
956:               v[0]   = -hxhydhz*(dd - gd);
957:               c[1].k = k; c[1].j = j-1; c[1].i = i;
958:               v[1]   = -hzhxdhy*(ds - gs);
959:               c[2].k = k; c[2].j = j; c[2].i = i-1;
960:               v[2]   = -hyhzdhx*(dw - gw);
961:               c[3].k = k; c[3].j = j; c[3].i = i;
962:               v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
963:               c[4].k = k+1; c[4].j = j; c[4].i = i;
964:               v[4]   = -hxhydhz*(du + gu);
965:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

967:               /* right-hand top up corner */
968:             } else {

970:               td = x[k-1][j][i];
971:               ad = 0.5*(t0 + td);
972:               bd = PetscPowScalar(ad,bm1);
973:               /* dd = bd * ad; */
974:               dd = PetscPowScalar(ad,beta);
975:               gd = coef*bd*(td - t0);

977:               c[0].k = k-1; c[0].j = j; c[0].i = i;
978:               v[0]   = -hxhydhz*(dd - gd);
979:               c[1].k = k; c[1].j = j-1; c[1].i = i;
980:               v[1]   = -hzhxdhy*(ds - gs);
981:               c[2].k = k; c[2].j = j; c[2].i = i-1;
982:               v[2]   = -hyhzdhx*(dw - gw);
983:               c[3].k = k; c[3].j = j; c[3].i = i;
984:               v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
985:               MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
986:             }

988:           } else {

990:             ts = x[k][j-1][i];
991:             as = 0.5*(t0 + ts);
992:             bs = PetscPowScalar(as,bm1);
993:             /* ds = bs * as; */
994:             ds = PetscPowScalar(as,beta);
995:             gs = coef*bs*(t0 - ts);

997:             tn = x[k][j+1][i];
998:             an = 0.5*(t0 + tn);
999:             bn = PetscPowScalar(an,bm1);
1000:             /* dn = bn * an; */
1001:             dn = PetscPowScalar(an,beta);
1002:             gn = coef*bn*(tn - t0);

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

1007:               tu = x[k+1][j][i];
1008:               au = 0.5*(t0 + tu);
1009:               bu = PetscPowScalar(au,bm1);
1010:               /* du = bu * au; */
1011:               du = PetscPowScalar(au,beta);
1012:               gu = coef*bu*(tu - t0);

1014:               c[0].k = k; c[0].j = j-1; c[0].i = i;
1015:               v[0]   = -hzhxdhy*(ds - gs);
1016:               c[1].k = k; c[1].j = j; c[1].i = i-1;
1017:               v[1]   = -hyhzdhx*(dw - gw);
1018:               c[2].k = k; c[2].j = j; c[2].i = i;
1019:               v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1020:               c[3].k = k; c[3].j = j+1; c[3].i = i;
1021:               v[3]   = -hzhxdhy*(dn + gn);
1022:               c[4].k = k+1; c[4].j = j; c[4].i = i;
1023:               v[4]   = -hxhydhz*(du + gu);
1024:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

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

1028:               td = x[k-1][j][i];
1029:               ad = 0.5*(t0 + td);
1030:               bd = PetscPowScalar(ad,bm1);
1031:               /* dd = bd * ad; */
1032:               dd = PetscPowScalar(ad,beta);
1033:               gd = coef*bd*(t0 - td);

1035:               c[0].k = k-1; c[0].j = j; c[0].i = i;
1036:               v[0]   = -hxhydhz*(dd - gd);
1037:               c[1].k = k; c[1].j = j-1; c[1].i = i;
1038:               v[1]   = -hzhxdhy*(ds - gs);
1039:               c[2].k = k; c[2].j = j; c[2].i = i-1;
1040:               v[2]   = -hyhzdhx*(dw - gw);
1041:               c[3].k = k; c[3].j = j; c[3].i = i;
1042:               v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1043:               c[4].k = k; c[4].j = j+1; c[4].i = i;
1044:               v[4]   = -hzhxdhy*(dn + gn);
1045:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

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

1049:               td = x[k-1][j][i];
1050:               ad = 0.5*(t0 + td);
1051:               bd = PetscPowScalar(ad,bm1);
1052:               /* dd = bd * ad; */
1053:               dd = PetscPowScalar(ad,beta);
1054:               gd = coef*bd*(t0 - td);

1056:               tu = x[k+1][j][i];
1057:               au = 0.5*(t0 + tu);
1058:               bu = PetscPowScalar(au,bm1);
1059:               /* du = bu * au; */
1060:               du = PetscPowScalar(au,beta);
1061:               gu = coef*bu*(tu - t0);

1063:               c[0].k = k-1; c[0].j = j; c[0].i = i;
1064:               v[0]   = -hxhydhz*(dd - gd);
1065:               c[1].k = k; c[1].j = j-1; c[1].i = i;
1066:               v[1]   = -hzhxdhy*(ds - gs);
1067:               c[2].k = k; c[2].j = j; c[2].i = i-1;
1068:               v[2]   = -hyhzdhx*(dw - gw);
1069:               c[3].k = k; c[3].j = j; c[3].i = i;
1070:               v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1071:               c[4].k = k; c[4].j = j+1; c[4].i = i;
1072:               v[4]   = -hzhxdhy*(dn + gn);
1073:               c[5].k = k+1; c[5].j = j; c[5].i = i;
1074:               v[5]   = -hxhydhz*(du + gu);
1075:               MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1076:             }
1077:           }

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

1081:           tw = x[k][j][i-1];
1082:           aw = 0.5*(t0 + tw);
1083:           bw = PetscPowScalar(aw,bm1);
1084:           /* dw = bw * aw */
1085:           dw = PetscPowScalar(aw,beta);
1086:           gw = coef*bw*(t0 - tw);

1088:           te = x[k][j][i+1];
1089:           ae = 0.5*(t0 + te);
1090:           be = PetscPowScalar(ae,bm1);
1091:           /* de = be * ae; */
1092:           de = PetscPowScalar(ae,beta);
1093:           ge = coef*be*(te - t0);

1095:           tn = x[k][j+1][i];
1096:           an = 0.5*(t0 + tn);
1097:           bn = PetscPowScalar(an,bm1);
1098:           /* dn = bn * an; */
1099:           dn = PetscPowScalar(an,beta);
1100:           gn = coef*bn*(tn - t0);

1102:           /* bottom down interior edge */
1103:           if (k == 0) {

1105:             tu = x[k+1][j][i];
1106:             au = 0.5*(t0 + tu);
1107:             bu = PetscPowScalar(au,bm1);
1108:             /* du = bu * au; */
1109:             du = PetscPowScalar(au,beta);
1110:             gu = coef*bu*(tu - t0);

1112:             c[0].k = k; c[0].j = j; c[0].i = i-1;
1113:             v[0]   = -hyhzdhx*(dw - gw);
1114:             c[1].k = k; c[1].j = j; c[1].i = i;
1115:             v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1116:             c[2].k = k; c[2].j = j; c[2].i = i+1;
1117:             v[2]   = -hyhzdhx*(de + ge);
1118:             c[3].k = k; c[3].j = j+1; c[3].i = i;
1119:             v[3]   = -hzhxdhy*(dn + gn);
1120:             c[4].k = k+1; c[4].j = j; c[4].i = i;
1121:             v[4]   = -hxhydhz*(du + gu);
1122:             MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

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

1126:             td = x[k-1][j][i];
1127:             ad = 0.5*(t0 + td);
1128:             bd = PetscPowScalar(ad,bm1);
1129:             /* dd = bd * ad; */
1130:             dd = PetscPowScalar(ad,beta);
1131:             gd = coef*bd*(td - t0);

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

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

1147:             tu = x[k+1][j][i];
1148:             au = 0.5*(t0 + tu);
1149:             bu = PetscPowScalar(au,bm1);
1150:             /* du = bu * au; */
1151:             du = PetscPowScalar(au,beta);
1152:             gu = coef*bu*(tu - t0);

1154:             td = x[k-1][j][i];
1155:             ad = 0.5*(t0 + td);
1156:             bd = PetscPowScalar(ad,bm1);
1157:             /* dd = bd * ad; */
1158:             dd = PetscPowScalar(ad,beta);
1159:             gd = coef*bd*(td - t0);

1161:             c[0].k = k-1; c[0].j = j; c[0].i = i;
1162:             v[0]   = -hxhydhz*(dd - gd);
1163:             c[1].k = k; c[1].j = j; c[1].i = i-1;
1164:             v[1]   = -hyhzdhx*(dw - gw);
1165:             c[2].k = k; c[2].j = j; c[2].i = i;
1166:             v[2]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1167:             c[3].k = k; c[3].j = j; c[3].i = i+1;
1168:             v[3]   = -hyhzdhx*(de + ge);
1169:             c[4].k = k; c[4].j = j+1; c[4].i = i;
1170:             v[4]   = -hzhxdhy*(dn + gn);
1171:             c[5].k = k+1; c[5].j = j; c[5].i = i;
1172:             v[5]   = -hxhydhz*(du + gu);
1173:             MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1174:           }

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

1178:           tw = x[k][j][i-1];
1179:           aw = 0.5*(t0 + tw);
1180:           bw = PetscPowScalar(aw,bm1);
1181:           /* dw = bw * aw */
1182:           dw = PetscPowScalar(aw,beta);
1183:           gw = coef*bw*(t0 - tw);

1185:           te = x[k][j][i+1];
1186:           ae = 0.5*(t0 + te);
1187:           be = PetscPowScalar(ae,bm1);
1188:           /* de = be * ae; */
1189:           de = PetscPowScalar(ae,beta);
1190:           ge = coef*be*(te - t0);

1192:           ts = x[k][j-1][i];
1193:           as = 0.5*(t0 + ts);
1194:           bs = PetscPowScalar(as,bm1);
1195:           /* ds = bs * as; */
1196:           ds = PetscPowScalar(as,beta);
1197:           gs = coef*bs*(t0 - ts);

1199:           /* top down interior edge */
1200:           if (k == 0) {

1202:             tu = x[k+1][j][i];
1203:             au = 0.5*(t0 + tu);
1204:             bu = PetscPowScalar(au,bm1);
1205:             /* du = bu * au; */
1206:             du = PetscPowScalar(au,beta);
1207:             gu = coef*bu*(tu - t0);

1209:             c[0].k = k; c[0].j = j-1; c[0].i = i;
1210:             v[0]   = -hzhxdhy*(ds - gs);
1211:             c[1].k = k; c[1].j = j; c[1].i = i-1;
1212:             v[1]   = -hyhzdhx*(dw - gw);
1213:             c[2].k = k; c[2].j = j; c[2].i = i;
1214:             v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1215:             c[3].k = k; c[3].j = j; c[3].i = i+1;
1216:             v[3]   = -hyhzdhx*(de + ge);
1217:             c[4].k = k+1; c[4].j = j; c[4].i = i;
1218:             v[4]   = -hxhydhz*(du + gu);
1219:             MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

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

1223:             td = x[k-1][j][i];
1224:             ad = 0.5*(t0 + td);
1225:             bd = PetscPowScalar(ad,bm1);
1226:             /* dd = bd * ad; */
1227:             dd = PetscPowScalar(ad,beta);
1228:             gd = coef*bd*(td - t0);

1230:             c[0].k = k-1; c[0].j = j; c[0].i = i;
1231:             v[0]   = -hxhydhz*(dd - gd);
1232:             c[1].k = k; c[1].j = j-1; c[1].i = i;
1233:             v[1]   = -hzhxdhy*(ds - gs);
1234:             c[2].k = k; c[2].j = j; c[2].i = i-1;
1235:             v[2]   = -hyhzdhx*(dw - gw);
1236:             c[3].k = k; c[3].j = j; c[3].i = i;
1237:             v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1238:             c[4].k = k; c[4].j = j; c[4].i = i+1;
1239:             v[4]   = -hyhzdhx*(de + ge);
1240:             MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

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

1244:             tu = x[k+1][j][i];
1245:             au = 0.5*(t0 + tu);
1246:             bu = PetscPowScalar(au,bm1);
1247:             /* du = bu * au; */
1248:             du = PetscPowScalar(au,beta);
1249:             gu = coef*bu*(tu - t0);

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

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

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

1275:           /* down interior plane */

1277:           tw = x[k][j][i-1];
1278:           aw = 0.5*(t0 + tw);
1279:           bw = PetscPowScalar(aw,bm1);
1280:           /* dw = bw * aw */
1281:           dw = PetscPowScalar(aw,beta);
1282:           gw = coef*bw*(t0 - tw);

1284:           te = x[k][j][i+1];
1285:           ae = 0.5*(t0 + te);
1286:           be = PetscPowScalar(ae,bm1);
1287:           /* de = be * ae; */
1288:           de = PetscPowScalar(ae,beta);
1289:           ge = coef*be*(te - t0);

1291:           ts = x[k][j-1][i];
1292:           as = 0.5*(t0 + ts);
1293:           bs = PetscPowScalar(as,bm1);
1294:           /* ds = bs * as; */
1295:           ds = PetscPowScalar(as,beta);
1296:           gs = coef*bs*(t0 - ts);

1298:           tn = x[k][j+1][i];
1299:           an = 0.5*(t0 + tn);
1300:           bn = PetscPowScalar(an,bm1);
1301:           /* dn = bn * an; */
1302:           dn = PetscPowScalar(an,beta);
1303:           gn = coef*bn*(tn - t0);

1305:           tu = x[k+1][j][i];
1306:           au = 0.5*(t0 + tu);
1307:           bu = PetscPowScalar(au,bm1);
1308:           /* du = bu * au; */
1309:           du = PetscPowScalar(au,beta);
1310:           gu = coef*bu*(tu - t0);

1312:           c[0].k = k; c[0].j = j-1; c[0].i = i;
1313:           v[0]   = -hzhxdhy*(ds - gs);
1314:           c[1].k = k; c[1].j = j; c[1].i = i-1;
1315:           v[1]   = -hyhzdhx*(dw - gw);
1316:           c[2].k = k; c[2].j = j; c[2].i = i;
1317:           v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1318:           c[3].k = k; c[3].j = j; c[3].i = i+1;
1319:           v[3]   = -hyhzdhx*(de + ge);
1320:           c[4].k = k; c[4].j = j+1; c[4].i = i;
1321:           v[4]   = -hzhxdhy*(dn + gn);
1322:           c[5].k = k+1; c[5].j = j; c[5].i = i;
1323:           v[5]   = -hxhydhz*(du + gu);
1324:           MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);

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

1328:           /* up interior plane */

1330:           tw = x[k][j][i-1];
1331:           aw = 0.5*(t0 + tw);
1332:           bw = PetscPowScalar(aw,bm1);
1333:           /* dw = bw * aw */
1334:           dw = PetscPowScalar(aw,beta);
1335:           gw = coef*bw*(t0 - tw);

1337:           te = x[k][j][i+1];
1338:           ae = 0.5*(t0 + te);
1339:           be = PetscPowScalar(ae,bm1);
1340:           /* de = be * ae; */
1341:           de = PetscPowScalar(ae,beta);
1342:           ge = coef*be*(te - t0);

1344:           ts = x[k][j-1][i];
1345:           as = 0.5*(t0 + ts);
1346:           bs = PetscPowScalar(as,bm1);
1347:           /* ds = bs * as; */
1348:           ds = PetscPowScalar(as,beta);
1349:           gs = coef*bs*(t0 - ts);

1351:           tn = x[k][j+1][i];
1352:           an = 0.5*(t0 + tn);
1353:           bn = PetscPowScalar(an,bm1);
1354:           /* dn = bn * an; */
1355:           dn = PetscPowScalar(an,beta);
1356:           gn = coef*bn*(tn - t0);

1358:           td = x[k-1][j][i];
1359:           ad = 0.5*(t0 + td);
1360:           bd = PetscPowScalar(ad,bm1);
1361:           /* dd = bd * ad; */
1362:           dd = PetscPowScalar(ad,beta);
1363:           gd = coef*bd*(t0 - td);

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

1391:   PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
1392:   return 0;
1393: }

1395: /*TEST

1397:    test:
1398:       nsize: 4
1399:       args: -snes_monitor_short -pc_mg_type full -ksp_type fgmres -pc_type mg -snes_view -pc_mg_levels 2 -pc_mg_galerkin pmat
1400:       requires: !single

1402: TEST*/