Actual source code: ex20.c

petsc-master 2016-08-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*);

 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);if (ierr) return ierr;
 71:   /* set problem parameters */
 72:   user.tleft  = 1.0;
 73:   user.tright = 0.1;
 74:   user.beta   = 2.5;
 75:   PetscOptionsGetReal(NULL,NULL,"-tleft",&user.tleft,NULL);
 76:   PetscOptionsGetReal(NULL,NULL,"-tright",&user.tright,NULL);
 77:   PetscOptionsGetReal(NULL,NULL,"-beta",&user.beta,NULL);
 78:   user.bm1    = user.beta - 1.0;
 79:   user.coef   = user.beta/2.0;

 81:   /*
 82:       Set the DMDA (grid structure) for the grids.
 83:   */
 84:   DMDACreate3d(PETSC_COMM_WORLD,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);
 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 ----------------- */
113: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
114: {
115:   AppCtx         *user;
116:   PetscInt       i,j,k,xs,ys,xm,ym,zs,zm;
118:   PetscScalar    ***x;
119:   DM             da;

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

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

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

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

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

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

178:           /* general interior volume */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

321:           fs = zero;

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

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

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

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

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

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

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

364:           fn = zero;

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

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

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

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

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

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

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

407:           fd = zero;

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

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

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

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

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

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

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

442:           fu = zero;
443:         }

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

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

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

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

495:           /* general interior volume */

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

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

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

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

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

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

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

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

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

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

571:           /* left-hand bottom edge */
572:           if (j == 0) {

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

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

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

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

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

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

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

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

630:               /* left-hand bottom up corner */
631:             } else {

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

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

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

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

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

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

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

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

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

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

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

710:               /* left-hand top up corner */
711:             } else {

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

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

731:           } else {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

838:           /* right-hand bottom edge */
839:           if (j == 0) {

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

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

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

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

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

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

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

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

897:               /* right-hand bottom up corner */
898:             } else {

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

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

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

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

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

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

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

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

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

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

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

977:               /* right-hand top up corner */
978:             } else {

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

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

998:           } else {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


1113:           /* bottom down interior edge */
1114:           if (k == 0) {

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

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

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

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

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

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

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

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

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

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

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

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

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

1210:           /* top down interior edge */
1211:           if (k == 0) {

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

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

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

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

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

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

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

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

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

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

1286:           /* down interior plane */

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

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

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

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

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

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

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

1339:           /* up interior plane */

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

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

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

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

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

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

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