Actual source code: ex20.c

petsc-master 2017-03-26
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*);

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

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

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

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

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

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

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

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

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

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

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

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

174:           /* general interior volume */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

263:           te = tright;
264:           ae = 0.5*(t0 + te);
265:           de = PetscPowScalar(ae,beta);
266:           fe = de*(te - t0);

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

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

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

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

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

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

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

317:           fs = zero;

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

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

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

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

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

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

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

360:           fn = zero;

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

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

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

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

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

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

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

403:           fd = zero;

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

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

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

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

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

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

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

438:           fu = zero;
439:         }

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

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

475:   /* Get ghost points */
476:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
477:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
478:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
479:   DMDAVecGetArray(da,localX,&x);

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

489:           /* general interior volume */

491:           tw = x[k][j][i-1];
492:           aw = 0.5*(t0 + tw);
493:           bw = PetscPowScalar(aw,bm1);
494:           /* dw = bw * aw */
495:           dw = PetscPowScalar(aw,beta);
496:           gw = coef*bw*(t0 - tw);

498:           te = x[k][j][i+1];
499:           ae = 0.5*(t0 + te);
500:           be = PetscPowScalar(ae,bm1);
501:           /* de = be * ae; */
502:           de = PetscPowScalar(ae,beta);
503:           ge = coef*be*(te - t0);

505:           ts = x[k][j-1][i];
506:           as = 0.5*(t0 + ts);
507:           bs = PetscPowScalar(as,bm1);
508:           /* ds = bs * as; */
509:           ds = PetscPowScalar(as,beta);
510:           gs = coef*bs*(t0 - ts);

512:           tn = x[k][j+1][i];
513:           an = 0.5*(t0 + tn);
514:           bn = PetscPowScalar(an,bm1);
515:           /* dn = bn * an; */
516:           dn = PetscPowScalar(an,beta);
517:           gn = coef*bn*(tn - t0);

519:           td = x[k-1][j][i];
520:           ad = 0.5*(t0 + td);
521:           bd = PetscPowScalar(ad,bm1);
522:           /* dd = bd * ad; */
523:           dd = PetscPowScalar(ad,beta);
524:           gd = coef*bd*(t0 - td);

526:           tu = x[k+1][j][i];
527:           au = 0.5*(t0 + tu);
528:           bu = PetscPowScalar(au,bm1);
529:           /* du = bu * au; */
530:           du = PetscPowScalar(au,beta);
531:           gu = coef*bu*(tu - t0);

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

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

550:           /* left-hand plane boundary */
551:           tw = tleft;
552:           aw = 0.5*(t0 + tw);
553:           bw = PetscPowScalar(aw,bm1);
554:           /* dw = bw * aw */
555:           dw = PetscPowScalar(aw,beta);
556:           gw = coef*bw*(t0 - tw);

558:           te = x[k][j][i+1];
559:           ae = 0.5*(t0 + te);
560:           be = PetscPowScalar(ae,bm1);
561:           /* de = be * ae; */
562:           de = PetscPowScalar(ae,beta);
563:           ge = coef*be*(te - t0);

565:           /* left-hand bottom edge */
566:           if (j == 0) {

568:             tn = x[k][j+1][i];
569:             an = 0.5*(t0 + tn);
570:             bn = PetscPowScalar(an,bm1);
571:             /* dn = bn * an; */
572:             dn = PetscPowScalar(an,beta);
573:             gn = coef*bn*(tn - t0);

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

578:               tu = x[k+1][j][i];
579:               au = 0.5*(t0 + tu);
580:               bu = PetscPowScalar(au,bm1);
581:               /* du = bu * au; */
582:               du = PetscPowScalar(au,beta);
583:               gu = coef*bu*(tu - t0);

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

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

598:               tu = x[k+1][j][i];
599:               au = 0.5*(t0 + tu);
600:               bu = PetscPowScalar(au,bm1);
601:               /* du = bu * au; */
602:               du = PetscPowScalar(au,beta);
603:               gu = coef*bu*(tu - t0);

605:               td = x[k-1][j][i];
606:               ad = 0.5*(t0 + td);
607:               bd = PetscPowScalar(ad,bm1);
608:               /* dd = bd * ad; */
609:               dd = PetscPowScalar(ad,beta);
610:               gd = coef*bd*(td - t0);

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

624:               /* left-hand bottom up corner */
625:             } else {

627:               td = x[k-1][j][i];
628:               ad = 0.5*(t0 + td);
629:               bd = PetscPowScalar(ad,bm1);
630:               /* dd = bd * ad; */
631:               dd = PetscPowScalar(ad,beta);
632:               gd = coef*bd*(td - t0);

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

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

648:             ts = x[k][j-1][i];
649:             as = 0.5*(t0 + ts);
650:             bs = PetscPowScalar(as,bm1);
651:             /* ds = bs * as; */
652:             ds = PetscPowScalar(as,beta);
653:             gs = coef*bs*(ts - t0);

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

658:               tu = x[k+1][j][i];
659:               au = 0.5*(t0 + tu);
660:               bu = PetscPowScalar(au,bm1);
661:               /* du = bu * au; */
662:               du = PetscPowScalar(au,beta);
663:               gu = coef*bu*(tu - t0);

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

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

678:               tu = x[k+1][j][i];
679:               au = 0.5*(t0 + tu);
680:               bu = PetscPowScalar(au,bm1);
681:               /* du = bu * au; */
682:               du = PetscPowScalar(au,beta);
683:               gu = coef*bu*(tu - t0);

685:               td = x[k-1][j][i];
686:               ad = 0.5*(t0 + td);
687:               bd = PetscPowScalar(ad,bm1);
688:               /* dd = bd * ad; */
689:               dd = PetscPowScalar(ad,beta);
690:               gd = coef*bd*(td - t0);

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

704:               /* left-hand top up corner */
705:             } else {

707:               td = x[k-1][j][i];
708:               ad = 0.5*(t0 + td);
709:               bd = PetscPowScalar(ad,bm1);
710:               /* dd = bd * ad; */
711:               dd = PetscPowScalar(ad,beta);
712:               gd = coef*bd*(td - t0);

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

725:           } else {

727:             ts = x[k][j-1][i];
728:             as = 0.5*(t0 + ts);
729:             bs = PetscPowScalar(as,bm1);
730:             /* ds = bs * as; */
731:             ds = PetscPowScalar(as,beta);
732:             gs = coef*bs*(t0 - ts);

734:             tn = x[k][j+1][i];
735:             an = 0.5*(t0 + tn);
736:             bn = PetscPowScalar(an,bm1);
737:             /* dn = bn * an; */
738:             dn = PetscPowScalar(an,beta);
739:             gn = coef*bn*(tn - t0);

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

744:               tu = x[k+1][j][i];
745:               au = 0.5*(t0 + tu);
746:               bu = PetscPowScalar(au,bm1);
747:               /* du = bu * au; */
748:               du = PetscPowScalar(au,beta);
749:               gu = coef*bu*(tu - t0);

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

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

765:               td = x[k-1][j][i];
766:               ad = 0.5*(t0 + td);
767:               bd = PetscPowScalar(ad,bm1);
768:               /* dd = bd * ad; */
769:               dd = PetscPowScalar(ad,beta);
770:               gd = coef*bd*(t0 - td);

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

785:               td = x[k-1][j][i];
786:               ad = 0.5*(t0 + td);
787:               bd = PetscPowScalar(ad,bm1);
788:               /* dd = bd * ad; */
789:               dd = PetscPowScalar(ad,beta);
790:               gd = coef*bd*(t0 - td);

792:               tu = x[k+1][j][i];
793:               au = 0.5*(t0 + tu);
794:               bu = PetscPowScalar(au,bm1);
795:               /* du = bu * au; */
796:               du = PetscPowScalar(au,beta);
797:               gu = coef*bu*(tu - t0);

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

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

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

825:           te = tright;
826:           ae = 0.5*(t0 + te);
827:           be = PetscPowScalar(ae,bm1);
828:           /* de = be * ae; */
829:           de = PetscPowScalar(ae,beta);
830:           ge = coef*be*(te - t0);

832:           /* right-hand bottom edge */
833:           if (j == 0) {

835:             tn = x[k][j+1][i];
836:             an = 0.5*(t0 + tn);
837:             bn = PetscPowScalar(an,bm1);
838:             /* dn = bn * an; */
839:             dn = PetscPowScalar(an,beta);
840:             gn = coef*bn*(tn - t0);

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

845:               tu = x[k+1][j][i];
846:               au = 0.5*(t0 + tu);
847:               bu = PetscPowScalar(au,bm1);
848:               /* du = bu * au; */
849:               du = PetscPowScalar(au,beta);
850:               gu = coef*bu*(tu - t0);

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

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

865:               tu = x[k+1][j][i];
866:               au = 0.5*(t0 + tu);
867:               bu = PetscPowScalar(au,bm1);
868:               /* du = bu * au; */
869:               du = PetscPowScalar(au,beta);
870:               gu = coef*bu*(tu - t0);

872:               td = x[k-1][j][i];
873:               ad = 0.5*(t0 + td);
874:               bd = PetscPowScalar(ad,bm1);
875:               /* dd = bd * ad; */
876:               dd = PetscPowScalar(ad,beta);
877:               gd = coef*bd*(td - t0);

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

891:               /* right-hand bottom up corner */
892:             } else {

894:               td = x[k-1][j][i];
895:               ad = 0.5*(t0 + td);
896:               bd = PetscPowScalar(ad,bm1);
897:               /* dd = bd * ad; */
898:               dd = PetscPowScalar(ad,beta);
899:               gd = coef*bd*(td - t0);

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

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

915:             ts = x[k][j-1][i];
916:             as = 0.5*(t0 + ts);
917:             bs = PetscPowScalar(as,bm1);
918:             /* ds = bs * as; */
919:             ds = PetscPowScalar(as,beta);
920:             gs = coef*bs*(ts - t0);

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

925:               tu = x[k+1][j][i];
926:               au = 0.5*(t0 + tu);
927:               bu = PetscPowScalar(au,bm1);
928:               /* du = bu * au; */
929:               du = PetscPowScalar(au,beta);
930:               gu = coef*bu*(tu - t0);

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

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

945:               tu = x[k+1][j][i];
946:               au = 0.5*(t0 + tu);
947:               bu = PetscPowScalar(au,bm1);
948:               /* du = bu * au; */
949:               du = PetscPowScalar(au,beta);
950:               gu = coef*bu*(tu - t0);

952:               td = x[k-1][j][i];
953:               ad = 0.5*(t0 + td);
954:               bd = PetscPowScalar(ad,bm1);
955:               /* dd = bd * ad; */
956:               dd = PetscPowScalar(ad,beta);
957:               gd = coef*bd*(td - t0);

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

971:               /* right-hand top up corner */
972:             } else {

974:               td = x[k-1][j][i];
975:               ad = 0.5*(t0 + td);
976:               bd = PetscPowScalar(ad,bm1);
977:               /* dd = bd * ad; */
978:               dd = PetscPowScalar(ad,beta);
979:               gd = coef*bd*(td - t0);

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

992:           } else {

994:             ts = x[k][j-1][i];
995:             as = 0.5*(t0 + ts);
996:             bs = PetscPowScalar(as,bm1);
997:             /* ds = bs * as; */
998:             ds = PetscPowScalar(as,beta);
999:             gs = coef*bs*(t0 - ts);

1001:             tn = x[k][j+1][i];
1002:             an = 0.5*(t0 + tn);
1003:             bn = PetscPowScalar(an,bm1);
1004:             /* dn = bn * an; */
1005:             dn = PetscPowScalar(an,beta);
1006:             gn = coef*bn*(tn - t0);

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

1011:               tu = x[k+1][j][i];
1012:               au = 0.5*(t0 + tu);
1013:               bu = PetscPowScalar(au,bm1);
1014:               /* du = bu * au; */
1015:               du = PetscPowScalar(au,beta);
1016:               gu = coef*bu*(tu - t0);

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

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

1032:               td = x[k-1][j][i];
1033:               ad = 0.5*(t0 + td);
1034:               bd = PetscPowScalar(ad,bm1);
1035:               /* dd = bd * ad; */
1036:               dd = PetscPowScalar(ad,beta);
1037:               gd = coef*bd*(t0 - td);

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

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

1053:               td = x[k-1][j][i];
1054:               ad = 0.5*(t0 + td);
1055:               bd = PetscPowScalar(ad,bm1);
1056:               /* dd = bd * ad; */
1057:               dd = PetscPowScalar(ad,beta);
1058:               gd = coef*bd*(t0 - td);

1060:               tu = x[k+1][j][i];
1061:               au = 0.5*(t0 + tu);
1062:               bu = PetscPowScalar(au,bm1);
1063:               /* du = bu * au; */
1064:               du = PetscPowScalar(au,beta);
1065:               gu = coef*bu*(tu - t0);

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

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

1085:           tw = x[k][j][i-1];
1086:           aw = 0.5*(t0 + tw);
1087:           bw = PetscPowScalar(aw,bm1);
1088:           /* dw = bw * aw */
1089:           dw = PetscPowScalar(aw,beta);
1090:           gw = coef*bw*(t0 - tw);

1092:           te = x[k][j][i+1];
1093:           ae = 0.5*(t0 + te);
1094:           be = PetscPowScalar(ae,bm1);
1095:           /* de = be * ae; */
1096:           de = PetscPowScalar(ae,beta);
1097:           ge = coef*be*(te - t0);

1099:           tn = x[k][j+1][i];
1100:           an = 0.5*(t0 + tn);
1101:           bn = PetscPowScalar(an,bm1);
1102:           /* dn = bn * an; */
1103:           dn = PetscPowScalar(an,beta);
1104:           gn = coef*bn*(tn - t0);


1107:           /* bottom down interior edge */
1108:           if (k == 0) {

1110:             tu = x[k+1][j][i];
1111:             au = 0.5*(t0 + tu);
1112:             bu = PetscPowScalar(au,bm1);
1113:             /* du = bu * au; */
1114:             du = PetscPowScalar(au,beta);
1115:             gu = coef*bu*(tu - t0);

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

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

1131:             td = x[k-1][j][i];
1132:             ad = 0.5*(t0 + td);
1133:             bd = PetscPowScalar(ad,bm1);
1134:             /* dd = bd * ad; */
1135:             dd = PetscPowScalar(ad,beta);
1136:             gd = coef*bd*(td - t0);

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

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

1152:             tu = x[k+1][j][i];
1153:             au = 0.5*(t0 + tu);
1154:             bu = PetscPowScalar(au,bm1);
1155:             /* du = bu * au; */
1156:             du = PetscPowScalar(au,beta);
1157:             gu = coef*bu*(tu - t0);

1159:             td = x[k-1][j][i];
1160:             ad = 0.5*(t0 + td);
1161:             bd = PetscPowScalar(ad,bm1);
1162:             /* dd = bd * ad; */
1163:             dd = PetscPowScalar(ad,beta);
1164:             gd = coef*bd*(td - t0);

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

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

1183:           tw = x[k][j][i-1];
1184:           aw = 0.5*(t0 + tw);
1185:           bw = PetscPowScalar(aw,bm1);
1186:           /* dw = bw * aw */
1187:           dw = PetscPowScalar(aw,beta);
1188:           gw = coef*bw*(t0 - tw);

1190:           te = x[k][j][i+1];
1191:           ae = 0.5*(t0 + te);
1192:           be = PetscPowScalar(ae,bm1);
1193:           /* de = be * ae; */
1194:           de = PetscPowScalar(ae,beta);
1195:           ge = coef*be*(te - t0);

1197:           ts = x[k][j-1][i];
1198:           as = 0.5*(t0 + ts);
1199:           bs = PetscPowScalar(as,bm1);
1200:           /* ds = bs * as; */
1201:           ds = PetscPowScalar(as,beta);
1202:           gs = coef*bs*(t0 - ts);

1204:           /* top down interior edge */
1205:           if (k == 0) {

1207:             tu = x[k+1][j][i];
1208:             au = 0.5*(t0 + tu);
1209:             bu = PetscPowScalar(au,bm1);
1210:             /* du = bu * au; */
1211:             du = PetscPowScalar(au,beta);
1212:             gu = coef*bu*(tu - t0);

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

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

1228:             td = x[k-1][j][i];
1229:             ad = 0.5*(t0 + td);
1230:             bd = PetscPowScalar(ad,bm1);
1231:             /* dd = bd * ad; */
1232:             dd = PetscPowScalar(ad,beta);
1233:             gd = coef*bd*(td - t0);

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

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

1249:             tu = x[k+1][j][i];
1250:             au = 0.5*(t0 + tu);
1251:             bu = PetscPowScalar(au,bm1);
1252:             /* du = bu * au; */
1253:             du = PetscPowScalar(au,beta);
1254:             gu = coef*bu*(tu - t0);

1256:             td = x[k-1][j][i];
1257:             ad = 0.5*(t0 + td);
1258:             bd = PetscPowScalar(ad,bm1);
1259:             /* dd = bd * ad; */
1260:             dd = PetscPowScalar(ad,beta);
1261:             gd = coef*bd*(td - t0);

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

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

1280:           /* down interior plane */

1282:           tw = x[k][j][i-1];
1283:           aw = 0.5*(t0 + tw);
1284:           bw = PetscPowScalar(aw,bm1);
1285:           /* dw = bw * aw */
1286:           dw = PetscPowScalar(aw,beta);
1287:           gw = coef*bw*(t0 - tw);

1289:           te = x[k][j][i+1];
1290:           ae = 0.5*(t0 + te);
1291:           be = PetscPowScalar(ae,bm1);
1292:           /* de = be * ae; */
1293:           de = PetscPowScalar(ae,beta);
1294:           ge = coef*be*(te - t0);

1296:           ts = x[k][j-1][i];
1297:           as = 0.5*(t0 + ts);
1298:           bs = PetscPowScalar(as,bm1);
1299:           /* ds = bs * as; */
1300:           ds = PetscPowScalar(as,beta);
1301:           gs = coef*bs*(t0 - ts);

1303:           tn = x[k][j+1][i];
1304:           an = 0.5*(t0 + tn);
1305:           bn = PetscPowScalar(an,bm1);
1306:           /* dn = bn * an; */
1307:           dn = PetscPowScalar(an,beta);
1308:           gn = coef*bn*(tn - t0);

1310:           tu = x[k+1][j][i];
1311:           au = 0.5*(t0 + tu);
1312:           bu = PetscPowScalar(au,bm1);
1313:           /* du = bu * au; */
1314:           du = PetscPowScalar(au,beta);
1315:           gu = coef*bu*(tu - t0);

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

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

1333:           /* up interior plane */

1335:           tw = x[k][j][i-1];
1336:           aw = 0.5*(t0 + tw);
1337:           bw = PetscPowScalar(aw,bm1);
1338:           /* dw = bw * aw */
1339:           dw = PetscPowScalar(aw,beta);
1340:           gw = coef*bw*(t0 - tw);

1342:           te = x[k][j][i+1];
1343:           ae = 0.5*(t0 + te);
1344:           be = PetscPowScalar(ae,bm1);
1345:           /* de = be * ae; */
1346:           de = PetscPowScalar(ae,beta);
1347:           ge = coef*be*(te - t0);

1349:           ts = x[k][j-1][i];
1350:           as = 0.5*(t0 + ts);
1351:           bs = PetscPowScalar(as,bm1);
1352:           /* ds = bs * as; */
1353:           ds = PetscPowScalar(as,beta);
1354:           gs = coef*bs*(t0 - ts);

1356:           tn = x[k][j+1][i];
1357:           an = 0.5*(t0 + tn);
1358:           bn = PetscPowScalar(an,bm1);
1359:           /* dn = bn * an; */
1360:           dn = PetscPowScalar(an,beta);
1361:           gn = coef*bn*(tn - t0);

1363:           td = x[k-1][j][i];
1364:           ad = 0.5*(t0 + td);
1365:           bd = PetscPowScalar(ad,bm1);
1366:           /* dd = bd * ad; */
1367:           dd = PetscPowScalar(ad,beta);
1368:           gd = coef*bd*(t0 - td);

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

1396:   PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
1397:   return(0);
1398: }