Actual source code: ex20.c

petsc-3.12.5 2020-03-29
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*/



 23: /*

 25:     This example models the partial differential equation

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

 29:     where beta = 2.5 and alpha = 1.0

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

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

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

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

 42: */

 44:  #include <petscsnes.h>
 45:  #include <petscdm.h>
 46:  #include <petscdmda.h>

 48: /* User-defined application context */

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

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

 57: extern PetscErrorCode FormInitialGuess(SNES,Vec,void*);
 58: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 59: 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:   DMSetFromOptions(da);
 86:   DMSetUp(da);
 87:   DMSetApplicationContext(da,&user);

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

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

107:   SNESDestroy(&snes);
108:   DMDestroy(&da);
109:   PetscFinalize();
110:   return ierr;
111: }
112: /* --------------------  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) --------------------- */
139: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
140: {
141:   AppCtx         *user = (AppCtx*)ptr;
143:   PetscInt       i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
144:   PetscScalar    zero = 0.0,one = 1.0;
145:   PetscScalar    hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
146:   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;
147:   PetscScalar    tleft,tright,beta,td,ad,dd,fd=0.0,tu,au,du=0.0,fu=0.0;
148:   PetscScalar    ***x,***f;
149:   Vec            localX;
150:   DM             da;

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

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

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

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

176:           /* general interior volume */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

319:           fs = zero;

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

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

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

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

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

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

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

362:           fn = zero;

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

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

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

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

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

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

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

405:           fd = zero;

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

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

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

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

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

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

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

440:           fu = zero;
441:         }

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

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

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

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

491:           /* general interior volume */

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

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

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

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

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

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

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

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

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

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

567:           /* left-hand bottom edge */
568:           if (j == 0) {

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

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

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

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

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

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

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

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

626:               /* left-hand bottom up corner */
627:             } else {

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

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

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

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

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

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

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

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

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

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

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

706:               /* left-hand top up corner */
707:             } else {

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

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

727:           } else {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

834:           /* right-hand bottom edge */
835:           if (j == 0) {

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

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

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

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

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

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

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

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

893:               /* right-hand bottom up corner */
894:             } else {

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

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

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

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

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

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

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

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

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

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

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

973:               /* right-hand top up corner */
974:             } else {

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

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

994:           } else {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


1109:           /* bottom down interior edge */
1110:           if (k == 0) {

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

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

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

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

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

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

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

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

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

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

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

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

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

1206:           /* top down interior edge */
1207:           if (k == 0) {

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

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

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

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

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

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

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

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

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

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

1282:           /* down interior plane */

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

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

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

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

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

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

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

1335:           /* up interior plane */

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

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

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

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

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

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

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


1403: /*TEST

1405:    test:
1406:       nsize: 4
1407:       args: -snes_monitor_short -pc_mg_type full -ksp_type fgmres -pc_type mg -snes_view -pc_mg_levels 2 -pc_mg_galerkin pmat
1408:       requires: !single

1410: TEST*/