Actual source code: ex20.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  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*);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

180:           /* general interior volume */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

323:           fs = zero;

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

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

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

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

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

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

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

366:           fn = zero;

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

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

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

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

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

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

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

409:           fd = zero;

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

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

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

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

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

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

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

444:           fu = zero;
445:         }

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

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

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

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

497:           /* general interior volume */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

733:           } else {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1000:           } else {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1288:           /* down interior plane */

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

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

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

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

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

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

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

1341:           /* up interior plane */

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

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

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

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

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

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

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