Actual source code: ex20.c

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

 14: /*T
 15:    Concepts: SNES^solving a system of nonlinear equations
 16:    Concepts: DA^using distributed arrays
 17:    Concepts: multigrid;
 18:    Processors: n
 19: T*/

 21: /*  
 22:   
 23:     This example models the partial differential equation 
 24:    
 25:          - Div(alpha* T^beta (GRAD T)) = 0.
 26:        
 27:     where beta = 2.5 and alpha = 1.0
 28:  
 29:     BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = dT/dn_up = dT/dn_down = 0.
 30:     
 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.
 33:  
 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
 39:  
 40: */

 42:  #include petscsnes.h
 43:  #include petscda.h
 44:  #include petscmg.h
 45:  #include petscdmmg.h

 47: /* User-defined application context */

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

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


 62: int main(int argc,char **argv)
 63: {
 64:   DMMG           *dmmg;
 65:   SNES           snes;
 66:   AppCtx         user;
 68:   PetscInt       its,lits;
 69:   PetscReal      litspit;
 70:   DA             da;

 72:   PetscInitialize(&argc,&argv,PETSC_NULL,help);

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


 85:   /*
 86:       Create the multilevel DA data structure 
 87:   */
 88:   DMMGCreate(PETSC_COMM_WORLD,2,&user,&dmmg);

 90:   /*
 91:       Set the DA (grid structure) for the grids.
 92:   */
 93:   DACreate3d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,-5,-5,-5,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);
 94:   DMMGSetDM(dmmg,(DM)da);
 95:   DADestroy(da);

 97:   /*
 98:      Create the nonlinear solver, and tell the DMMG structure to use it
 99:   */
100:   DMMGSetSNES(dmmg,FormFunction,FormJacobian);
101:   DMMGSetFromOptions(dmmg);

103:   /*
104:       PreLoadBegin() means that the following section of code is run twice. The first time
105:      through the flag PreLoading is on this the nonlinear solver is only run for a single step.
106:      The second time through (the actually timed code) the maximum iterations is set to 10
107:      Preload of the executable is done to eliminate from the timing the time spent bring the 
108:      executable into memory from disk (paging in).
109:   */
110:   PreLoadBegin(PETSC_TRUE,"Solve");
111:     DMMGSetInitialGuess(dmmg,FormInitialGuess);
112:     DMMGSolve(dmmg);
113:   PreLoadEnd();
114:   snes = DMMGGetSNES(dmmg);
115:   SNESGetIterationNumber(snes,&its);
116:   SNESGetLinearSolveIterations(snes,&lits);
117:   litspit = ((PetscReal)lits)/((PetscReal)its);
118:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);
119:   PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
120:   PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / Newton = %e\n",litspit);

122:   DMMGDestroy(dmmg);
123:   PetscFinalize();

125:   return 0;
126: }
127: /* --------------------  Form initial approximation ----------------- */
130: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
131: {
132:   AppCtx         *user = (AppCtx*)dmmg->user;
133:   PetscInt       i,j,k,xs,ys,xm,ym,zs,zm;
135:   PetscReal      tleft = user->tleft;
136:   PetscScalar    ***x;


140:   /* Get ghost points */
141:   DAGetCorners((DA)dmmg->dm,&xs,&ys,&zs,&xm,&ym,&zm);
142:   DAVecGetArray((DA)dmmg->dm,X,&x);

144:   /* Compute initial guess */
145:   for (k=zs; k<zs+zm; k++) {
146:     for (j=ys; j<ys+ym; j++) {
147:       for (i=xs; i<xs+xm; i++) {
148:         x[k][j][i] = tleft;
149:       }
150:     }
151:   }
152:   DAVecRestoreArray((DA)dmmg->dm,X,&x);
153:   return(0);
154: }
155: /* --------------------  Evaluate Function F(x) --------------------- */
158: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ptr)
159: {
160:   DMMG           dmmg = (DMMG)ptr;
161:   AppCtx         *user = (AppCtx*)dmmg->user;
163:   PetscInt       i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
164:   PetscScalar    zero = 0.0,one = 1.0;
165:   PetscScalar    hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
166:   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;
167:   PetscScalar    tleft,tright,beta,td,ad,dd,fd,tu,au,du,fu;
168:   PetscScalar    ***x,***f;
169:   Vec            localX;

172:   DAGetLocalVector((DA)dmmg->dm,&localX);
173:   DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,&mz,0,0,0,0,0,0,0);
174:   hx    = one/(PetscReal)(mx-1);  hy    = one/(PetscReal)(my-1);  hz = one/(PetscReal)(mz-1);
175:   hxhydhz = hx*hy/hz;   hyhzdhx = hy*hz/hx;   hzhxdhy = hz*hx/hy;
176:   tleft = user->tleft;         tright = user->tright;
177:   beta  = user->beta;
178: 
179:   /* Get ghost points */
180:   DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
181:   DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
182:   DAGetCorners((DA)dmmg->dm,&xs,&ys,&zs,&xm,&ym,&zm);
183:   DAVecGetArray((DA)dmmg->dm,localX,&x);
184:   DAVecGetArray((DA)dmmg->dm,F,&f);

186:   /* Evaluate function */
187:   for (k=zs; k<zs+zm; k++) {
188:     for (j=ys; j<ys+ym; j++) {
189:       for (i=xs; i<xs+xm; i++) {
190:         t0 = x[k][j][i];

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

194:             /* general interior volume */

196:           tw = x[k][j][i-1];
197:           aw = 0.5*(t0 + tw);
198:           dw = pow(aw,beta);
199:           fw = dw*(t0 - tw);

201:                  te = x[k][j][i+1];
202:           ae = 0.5*(t0 + te);
203:           de = pow(ae,beta);
204:           fe = de*(te - t0);

206:           ts = x[k][j-1][i];
207:           as = 0.5*(t0 + ts);
208:           ds = pow(as,beta);
209:           fs = ds*(t0 - ts);
210: 
211:           tn = x[k][j+1][i];
212:           an = 0.5*(t0 + tn);
213:           dn = pow(an,beta);
214:           fn = dn*(tn - t0);

216:           td = x[k-1][j][i];
217:           ad = 0.5*(t0 + td);
218:           dd = pow(ad,beta);
219:           fd = dd*(t0 - td);

221:           tu = x[k+1][j][i];
222:           au = 0.5*(t0 + tu);
223:           du = pow(au,beta);
224:           fu = du*(tu - t0);

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

228:            /* left-hand (west) boundary */
229:           tw = tleft;
230:           aw = 0.5*(t0 + tw);
231:           dw = pow(aw,beta);
232:           fw = dw*(t0 - tw);

234:           te = x[k][j][i+1];
235:           ae = 0.5*(t0 + te);
236:           de = pow(ae,beta);
237:           fe = de*(te - t0);

239:           if (j > 0) {
240:             ts = x[k][j-1][i];
241:             as = 0.5*(t0 + ts);
242:             ds = pow(as,beta);
243:             fs = ds*(t0 - ts);
244:           } else {
245:              fs = zero;
246:           }

248:           if (j < my-1) {
249:             tn = x[k][j+1][i];
250:             an = 0.5*(t0 + tn);
251:             dn = pow(an,beta);
252:             fn = dn*(tn - t0);
253:           } else {
254:             fn = zero;
255:              }

257:           if (k > 0) {
258:             td = x[k-1][j][i];
259:             ad = 0.5*(t0 + td);
260:             dd = pow(ad,beta);
261:             fd = dd*(t0 - td);
262:           } else {
263:              fd = zero;
264:           }

266:           if (k < mz-1) {
267:             tu = x[k+1][j][i];
268:             au = 0.5*(t0 + tu);
269:             du = pow(au,beta);
270:             fu = du*(tu - t0);
271:           } else {
272:              fu = zero;
273:           }

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

277:           /* right-hand (east) boundary */
278:           tw = x[k][j][i-1];
279:           aw = 0.5*(t0 + tw);
280:           dw = pow(aw,beta);
281:           fw = dw*(t0 - tw);
282: 
283:           te = tright;
284:           ae = 0.5*(t0 + te);
285:           de = pow(ae,beta);
286:           fe = de*(te - t0);
287: 
288:           if (j > 0) {
289:             ts = x[k][j-1][i];
290:             as = 0.5*(t0 + ts);
291:             ds = pow(as,beta);
292:             fs = ds*(t0 - ts);
293:           } else {
294:             fs = zero;
295:           }
296: 
297:           if (j < my-1) {
298:             tn = x[k][j+1][i];
299:             an = 0.5*(t0 + tn);
300:             dn = pow(an,beta);
301:             fn = dn*(tn - t0);
302:           } else {
303:             fn = zero;
304:           }

306:           if (k > 0) {
307:             td = x[k-1][j][i];
308:             ad = 0.5*(t0 + td);
309:             dd = pow(ad,beta);
310:             fd = dd*(t0 - td);
311:           } else {
312:              fd = zero;
313:           }

315:           if (k < mz-1) {
316:             tu = x[k+1][j][i];
317:             au = 0.5*(t0 + tu);
318:             du = pow(au,beta);
319:             fu = du*(tu - t0);
320:           } else {
321:              fu = zero;
322:           }

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

326:           /* bottom (south) boundary, and i <> 0 or mx-1 */
327:           tw = x[k][j][i-1];
328:           aw = 0.5*(t0 + tw);
329:           dw = pow(aw,beta);
330:           fw = dw*(t0 - tw);

332:           te = x[k][j][i+1];
333:           ae = 0.5*(t0 + te);
334:           de = pow(ae,beta);
335:           fe = de*(te - t0);

337:           fs = zero;

339:           tn = x[k][j+1][i];
340:           an = 0.5*(t0 + tn);
341:           dn = pow(an,beta);
342:           fn = dn*(tn - t0);

344:           if (k > 0) {
345:             td = x[k-1][j][i];
346:             ad = 0.5*(t0 + td);
347:             dd = pow(ad,beta);
348:             fd = dd*(t0 - td);
349:           } else {
350:              fd = zero;
351:           }

353:           if (k < mz-1) {
354:             tu = x[k+1][j][i];
355:             au = 0.5*(t0 + tu);
356:             du = pow(au,beta);
357:             fu = du*(tu - t0);
358:           } else {
359:              fu = zero;
360:           }

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

364:           /* top (north) boundary, and i <> 0 or mx-1 */
365:           tw = x[k][j][i-1];
366:           aw = 0.5*(t0 + tw);
367:           dw = pow(aw,beta);
368:           fw = dw*(t0 - tw);

370:           te = x[k][j][i+1];
371:           ae = 0.5*(t0 + te);
372:           de = pow(ae,beta);
373:           fe = de*(te - t0);

375:           ts = x[k][j-1][i];
376:           as = 0.5*(t0 + ts);
377:           ds = pow(as,beta);
378:           fs = ds*(t0 - ts);

380:           fn = zero;

382:           if (k > 0) {
383:             td = x[k-1][j][i];
384:             ad = 0.5*(t0 + td);
385:             dd = pow(ad,beta);
386:             fd = dd*(t0 - td);
387:           } else {
388:              fd = zero;
389:           }

391:           if (k < mz-1) {
392:             tu = x[k+1][j][i];
393:             au = 0.5*(t0 + tu);
394:             du = pow(au,beta);
395:             fu = du*(tu - t0);
396:           } else {
397:              fu = zero;
398:           }

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

402:           /* down boundary (interior only) */
403:           tw = x[k][j][i-1];
404:           aw = 0.5*(t0 + tw);
405:           dw = pow(aw,beta);
406:           fw = dw*(t0 - tw);

408:           te = x[k][j][i+1];
409:           ae = 0.5*(t0 + te);
410:           de = pow(ae,beta);
411:           fe = de*(te - t0);

413:           ts = x[k][j-1][i];
414:           as = 0.5*(t0 + ts);
415:           ds = pow(as,beta);
416:           fs = ds*(t0 - ts);

418:           tn = x[k][j+1][i];
419:           an = 0.5*(t0 + tn);
420:           dn = pow(an,beta);
421:           fn = dn*(tn - t0);

423:            fd = zero;

425:           tu = x[k+1][j][i];
426:           au = 0.5*(t0 + tu);
427:           du = pow(au,beta);
428:           fu = du*(tu - t0);
429: 
430:         } else if (k == mz-1) {

432:           /* up boundary (interior only) */
433:           tw = x[k][j][i-1];
434:           aw = 0.5*(t0 + tw);
435:           dw = pow(aw,beta);
436:           fw = dw*(t0 - tw);

438:           te = x[k][j][i+1];
439:           ae = 0.5*(t0 + te);
440:           de = pow(ae,beta);
441:           fe = de*(te - t0);

443:           ts = x[k][j-1][i];
444:           as = 0.5*(t0 + ts);
445:           ds = pow(as,beta);
446:           fs = ds*(t0 - ts);

448:           tn = x[k][j+1][i];
449:           an = 0.5*(t0 + tn);
450:           dn = pow(an,beta);
451:           fn = dn*(tn - t0);

453:           td = x[k-1][j][i];
454:           ad = 0.5*(t0 + td);
455:           dd = pow(ad,beta);
456:           fd = dd*(t0 - td);

458:           fu = zero;
459:         }

461:         f[k][j][i] = - hyhzdhx*(fe-fw) - hzhxdhy*(fn-fs) - hxhydhz*(fu-fd);
462:       }
463:     }
464:   }
465:   DAVecRestoreArray((DA)dmmg->dm,localX,&x);
466:   DAVecRestoreArray((DA)dmmg->dm,F,&f);
467:   DARestoreLocalVector((DA)dmmg->dm,&localX);
468:   PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
469:   return(0);
470: }
471: /* --------------------  Evaluate Jacobian F(x) --------------------- */
474: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ptr)
475: {
476:   DMMG           dmmg = (DMMG)ptr;
477:   AppCtx         *user = (AppCtx*)dmmg->user;
479:   PetscInt       i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
480:   PetscScalar    one = 1.0;
481:   PetscScalar    hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
482:   PetscScalar    t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw;
483:   PetscScalar    tleft,tright,beta,td,ad,dd,tu,au,du,v[7],bm1,coef;
484:   PetscScalar    ***x,bn,bs,be,bw,bu,bd,gn,gs,ge,gw,gu,gd;
485:   Vec            localX;
486:   MatStencil     c[7],row;
487:   Mat            jac = *B;

490:   DAGetLocalVector((DA)dmmg->dm,&localX);
491:   DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,&mz,0,0,0,0,0,0,0);
492:   hx    = one/(PetscReal)(mx-1);  hy    = one/(PetscReal)(my-1);  hz = one/(PetscReal)(mz-1);
493:   hxhydhz = hx*hy/hz;   hyhzdhx = hy*hz/hx;   hzhxdhy = hz*hx/hy;
494:   tleft = user->tleft;         tright = user->tright;
495:   beta  = user->beta;               bm1    = user->bm1;                coef = user->coef;

497:   /* Get ghost points */
498:   DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
499:   DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
500:   DAGetCorners((DA)dmmg->dm,&xs,&ys,&zs,&xm,&ym,&zm);
501:   DAVecGetArray((DA)dmmg->dm,localX,&x);

503:   /* Evaluate Jacobian of function */
504:   for (k=zs; k<zs+zm; k++) {
505:     for (j=ys; j<ys+ym; j++) {
506:       for (i=xs; i<xs+xm; i++) {
507:         t0    = x[k][j][i];
508:         row.k = k; row.j = j; row.i = i;
509:         if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {

511:           /* general interior volume */

513:           tw = x[k][j][i-1];
514:           aw = 0.5*(t0 + tw);
515:           bw = pow(aw,bm1);
516:           /* dw = bw * aw */
517:           dw = pow(aw,beta);
518:           gw = coef*bw*(t0 - tw);

520:           te = x[k][j][i+1];
521:           ae = 0.5*(t0 + te);
522:           be = pow(ae,bm1);
523:           /* de = be * ae; */
524:           de = pow(ae,beta);
525:           ge = coef*be*(te - t0);

527:           ts = x[k][j-1][i];
528:           as = 0.5*(t0 + ts);
529:           bs = pow(as,bm1);
530:           /* ds = bs * as; */
531:           ds = pow(as,beta);
532:           gs = coef*bs*(t0 - ts);
533: 
534:           tn = x[k][j+1][i];
535:           an = 0.5*(t0 + tn);
536:           bn = pow(an,bm1);
537:           /* dn = bn * an; */
538:           dn = pow(an,beta);
539:           gn = coef*bn*(tn - t0);

541:           td = x[k-1][j][i];
542:           ad = 0.5*(t0 + td);
543:           bd = pow(ad,bm1);
544:           /* dd = bd * ad; */
545:           dd = pow(ad,beta);
546:           gd = coef*bd*(t0 - td);
547: 
548:           tu = x[k+1][j][i];
549:           au = 0.5*(t0 + tu);
550:           bu = pow(au,bm1);
551:           /* du = bu * au; */
552:           du = pow(au,beta);
553:           gu = coef*bu*(tu - t0);

555:           c[0].k = k-1; c[0].j = j; c[0].i = i; v[0]   = - hxhydhz*(dd - gd);
556:           c[1].k = k; c[1].j = j-1; c[1].i = i;
557:           v[1]   = - hzhxdhy*(ds - gs);
558:           c[2].k = k; c[2].j = j; c[2].i = i-1;
559:           v[2]   = - hyhzdhx*(dw - gw);
560:           c[3].k = k; c[3].j = j; c[3].i = i;
561:           v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
562:           c[4].k = k; c[4].j = j; c[4].i = i+1;
563:           v[4]   = - hyhzdhx*(de + ge);
564:           c[5].k = k; c[5].j = j+1; c[5].i = i;
565:           v[5]   = - hzhxdhy*(dn + gn);
566:           c[6].k = k+1; c[6].j = j; c[6].i = i;
567:           v[6]   = - hxhydhz*(du + gu);
568:             MatSetValuesStencil(jac,1,&row,7,c,v,INSERT_VALUES);

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

572:           /* left-hand plane boundary */
573:           tw = tleft;
574:           aw = 0.5*(t0 + tw);
575:           bw = pow(aw,bm1);
576:           /* dw = bw * aw */
577:           dw = pow(aw,beta);
578:           gw = coef*bw*(t0 - tw);
579: 
580:           te = x[k][j][i+1];
581:           ae = 0.5*(t0 + te);
582:           be = pow(ae,bm1);
583:           /* de = be * ae; */
584:           de = pow(ae,beta);
585:           ge = coef*be*(te - t0);
586: 
587:           /* left-hand bottom edge */
588:           if (j == 0) {

590:             tn = x[k][j+1][i];
591:             an = 0.5*(t0 + tn);
592:             bn = pow(an,bm1);
593:             /* dn = bn * an; */
594:             dn = pow(an,beta);
595:             gn = coef*bn*(tn - t0);
596: 
597:             /* left-hand bottom down corner */
598:             if (k == 0) {

600:               tu = x[k+1][j][i];
601:               au = 0.5*(t0 + tu);
602:               bu = pow(au,bm1);
603:               /* du = bu * au; */
604:               du = pow(au,beta);
605:               gu = coef*bu*(tu - t0);
606: 
607:               c[0].k = k; c[0].j = j; c[0].i = i;
608:               v[0]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
609:               c[1].k = k; c[1].j = j; c[1].i = i+1;
610:               v[1]   = - hyhzdhx*(de + ge);
611:               c[2].k = k; c[2].j = j+1; c[2].i = i;
612:               v[2]   = - hzhxdhy*(dn + gn);
613:               c[3].k = k+1; c[3].j = j; c[3].i = i;
614:               v[3]   = - hxhydhz*(du + gu);
615:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

620:               tu = x[k+1][j][i];
621:               au = 0.5*(t0 + tu);
622:               bu = pow(au,bm1);
623:               /* du = bu * au; */
624:               du = pow(au,beta);
625:               gu = coef*bu*(tu - t0);
626: 
627:               td = x[k-1][j][i];
628:               ad = 0.5*(t0 + td);
629:               bd = pow(ad,bm1);
630:               /* dd = bd * ad; */
631:               dd = pow(ad,beta);
632:               gd = coef*bd*(td - t0);
633: 
634:               c[0].k = k-1; c[0].j = j; c[0].i = i;
635:               v[0]   = - hxhydhz*(dd - gd);
636:               c[1].k = k; c[1].j = j; c[1].i = i;
637:               v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
638:               c[2].k = k; c[2].j = j; c[2].i = i+1;
639:               v[2]   = - hyhzdhx*(de + ge);
640:               c[3].k = k; c[3].j = j+1; c[3].i = i;
641:               v[3]   = - hzhxdhy*(dn + gn);
642:               c[4].k = k+1; c[4].j = j; c[4].i = i;
643:               v[4]   = - hxhydhz*(du + gu);
644:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

646:             /* left-hand bottom up corner */
647:             } else {

649:               td = x[k-1][j][i];
650:               ad = 0.5*(t0 + td);
651:               bd = pow(ad,bm1);
652:               /* dd = bd * ad; */
653:               dd = pow(ad,beta);
654:               gd = coef*bd*(td - t0);
655: 
656:               c[0].k = k-1; c[0].j = j; c[0].i = i;
657:               v[0]   = - hxhydhz*(dd - gd);
658:               c[1].k = k; c[1].j = j; c[1].i = i;
659:               v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
660:               c[2].k = k; c[2].j = j; c[2].i = i+1;
661:               v[2]   = - hyhzdhx*(de + ge);
662:               c[3].k = k; c[3].j = j+1; c[3].i = i;
663:               v[3]   = - hzhxdhy*(dn + gn);
664:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
665:             }

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

670:             ts = x[k][j-1][i];
671:             as = 0.5*(t0 + ts);
672:             bs = pow(as,bm1);
673:             /* ds = bs * as; */
674:             ds = pow(as,beta);
675:             gs = coef*bs*(ts - t0);
676: 
677:             /* left-hand top down corner */
678:             if (k == 0) {

680:               tu = x[k+1][j][i];
681:               au = 0.5*(t0 + tu);
682:               bu = pow(au,bm1);
683:               /* du = bu * au; */
684:               du = pow(au,beta);
685:               gu = coef*bu*(tu - t0);
686: 
687:               c[0].k = k; c[0].j = j-1; c[0].i = i;
688:               v[0]   = - hzhxdhy*(ds - gs);
689:               c[1].k = k; c[1].j = j; c[1].i = i;
690:               v[1]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
691:               c[2].k = k; c[2].j = j; c[2].i = i+1;
692:               v[2]   = - hyhzdhx*(de + ge);
693:               c[3].k = k+1; c[3].j = j; c[3].i = i;
694:               v[3]   = - hxhydhz*(du + gu);
695:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

700:               tu = x[k+1][j][i];
701:               au = 0.5*(t0 + tu);
702:               bu = pow(au,bm1);
703:               /* du = bu * au; */
704:               du = pow(au,beta);
705:               gu = coef*bu*(tu - t0);
706: 
707:               td = x[k-1][j][i];
708:               ad = 0.5*(t0 + td);
709:               bd = pow(ad,bm1);
710:               /* dd = bd * ad; */
711:               dd = pow(ad,beta);
712:               gd = coef*bd*(td - t0);
713: 
714:               c[0].k = k-1; c[0].j = j; c[0].i = i;
715:               v[0]   = - hxhydhz*(dd - gd);
716:               c[1].k = k; c[1].j = j-1; c[1].i = i;
717:               v[1]   = - hzhxdhy*(ds - gs);
718:               c[2].k = k; c[2].j = j; c[2].i = i;
719:               v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
720:               c[3].k = k; c[3].j = j; c[3].i = i+1;
721:               v[3]   = - hyhzdhx*(de + ge);
722:               c[4].k = k+1; c[4].j = j; c[4].i = i;
723:               v[4]   = - hxhydhz*(du + gu);
724:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

726:             /* left-hand top up corner */
727:             } else {

729:               td = x[k-1][j][i];
730:               ad = 0.5*(t0 + td);
731:               bd = pow(ad,bm1);
732:               /* dd = bd * ad; */
733:               dd = pow(ad,beta);
734:               gd = coef*bd*(td - t0);
735: 
736:               c[0].k = k-1; c[0].j = j; c[0].i = i;
737:               v[0]   = - hxhydhz*(dd - gd);
738:               c[1].k = k; c[1].j = j-1; c[1].i = i;
739:               v[1]   = - hzhxdhy*(ds - gs);
740:               c[2].k = k; c[2].j = j; c[2].i = i;
741:               v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
742:               c[3].k = k; c[3].j = j; c[3].i = i+1;
743:               v[3]   = - hyhzdhx*(de + ge);
744:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
745:             }

747:           } else {

749:             ts = x[k][j-1][i];
750:             as = 0.5*(t0 + ts);
751:             bs = pow(as,bm1);
752:             /* ds = bs * as; */
753:             ds = pow(as,beta);
754:             gs = coef*bs*(t0 - ts);
755: 
756:             tn = x[k][j+1][i];
757:             an = 0.5*(t0 + tn);
758:             bn = pow(an,bm1);
759:             /* dn = bn * an; */
760:             dn = pow(an,beta);
761:             gn = coef*bn*(tn - t0);

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

766:               tu = x[k+1][j][i];
767:               au = 0.5*(t0 + tu);
768:               bu = pow(au,bm1);
769:               /* du = bu * au; */
770:               du = pow(au,beta);
771:               gu = coef*bu*(tu - t0);

773:               c[0].k = k; c[0].j = j-1; c[0].i = i;
774:               v[0]   = - hzhxdhy*(ds - gs);
775:               c[1].k = k; c[1].j = j; c[1].i = i;
776:               v[1]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
777:               c[2].k = k; c[2].j = j; c[2].i = i+1;
778:               v[2]   = - hyhzdhx*(de + ge);
779:               c[3].k = k; c[3].j = j+1; c[3].i = i;
780:               v[3]   = - hzhxdhy*(dn + gn);
781:               c[4].k = k+1; c[4].j = j; c[4].i = i;
782:               v[4]   = - hxhydhz*(du + gu);
783:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
784:             }

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

789:               td = x[k-1][j][i];
790:               ad = 0.5*(t0 + td);
791:               bd = pow(ad,bm1);
792:               /* dd = bd * ad; */
793:               dd = pow(ad,beta);
794:               gd = coef*bd*(t0 - td);
795: 
796:               c[0].k = k-1; c[0].j = j; c[0].i = i;
797:               v[0]   = - hxhydhz*(dd - gd);
798:               c[1].k = k; c[1].j = j-1; c[1].i = i;
799:               v[1]   = - hzhxdhy*(ds - gs);
800:               c[2].k = k; c[2].j = j; c[2].i = i;
801:               v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
802:               c[3].k = k; c[3].j = j; c[3].i = i+1;
803:               v[3]   = - hyhzdhx*(de + ge);
804:               c[4].k = k; c[4].j = j+1; c[4].i = i;
805:               v[4]   = - hzhxdhy*(dn + gn);
806:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
807:             }

809:             /* left-hand interior plane */
810:             else {

812:               td = x[k-1][j][i];
813:               ad = 0.5*(t0 + td);
814:               bd = pow(ad,bm1);
815:               /* dd = bd * ad; */
816:               dd = pow(ad,beta);
817:               gd = coef*bd*(t0 - td);
818: 
819:               tu = x[k+1][j][i];
820:               au = 0.5*(t0 + tu);
821:               bu = pow(au,bm1);
822:               /* du = bu * au; */
823:               du = pow(au,beta);
824:               gu = coef*bu*(tu - t0);

826:               c[0].k = k-1; c[0].j = j; c[0].i = i;
827:               v[0]   = - hxhydhz*(dd - gd);
828:               c[1].k = k; c[1].j = j-1; c[1].i = i;
829:               v[1]   = - hzhxdhy*(ds - gs);
830:               c[2].k = k; c[2].j = j; c[2].i = i;
831:               v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
832:               c[3].k = k; c[3].j = j; c[3].i = i+1;
833:               v[3]   = - hyhzdhx*(de + ge);
834:               c[4].k = k; c[4].j = j+1; c[4].i = i;
835:               v[4]   = - hzhxdhy*(dn + gn);
836:               c[5].k = k+1; c[5].j = j; c[5].i = i;
837:               v[5]   = - hxhydhz*(du + gu);
838:                 MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
839:             }
840:           }

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

844:           /* right-hand plane boundary */
845:           tw = x[k][j][i-1];
846:           aw = 0.5*(t0 + tw);
847:           bw = pow(aw,bm1);
848:           /* dw = bw * aw */
849:           dw = pow(aw,beta);
850:           gw = coef*bw*(t0 - tw);
851: 
852:           te = tright;
853:           ae = 0.5*(t0 + te);
854:           be = pow(ae,bm1);
855:           /* de = be * ae; */
856:           de = pow(ae,beta);
857:           ge = coef*be*(te - t0);
858: 
859:           /* right-hand bottom edge */
860:           if (j == 0) {

862:             tn = x[k][j+1][i];
863:             an = 0.5*(t0 + tn);
864:             bn = pow(an,bm1);
865:             /* dn = bn * an; */
866:             dn = pow(an,beta);
867:             gn = coef*bn*(tn - t0);
868: 
869:             /* right-hand bottom down corner */
870:             if (k == 0) {

872:               tu = x[k+1][j][i];
873:               au = 0.5*(t0 + tu);
874:               bu = pow(au,bm1);
875:               /* du = bu * au; */
876:               du = pow(au,beta);
877:               gu = coef*bu*(tu - t0);
878: 
879:               c[0].k = k; c[0].j = j; c[0].i = i-1;
880:               v[0]   = - hyhzdhx*(dw - gw);
881:               c[1].k = k; c[1].j = j; c[1].i = i;
882:               v[1]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
883:               c[2].k = k; c[2].j = j+1; c[2].i = i;
884:               v[2]   = - hzhxdhy*(dn + gn);
885:               c[3].k = k+1; c[3].j = j; c[3].i = i;
886:               v[3]   = - hxhydhz*(du + gu);
887:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

892:               tu = x[k+1][j][i];
893:               au = 0.5*(t0 + tu);
894:               bu = pow(au,bm1);
895:               /* du = bu * au; */
896:               du = pow(au,beta);
897:               gu = coef*bu*(tu - t0);
898: 
899:               td = x[k-1][j][i];
900:               ad = 0.5*(t0 + td);
901:               bd = pow(ad,bm1);
902:               /* dd = bd * ad; */
903:               dd = pow(ad,beta);
904:               gd = coef*bd*(td - t0);
905: 
906:               c[0].k = k-1; c[0].j = j; c[0].i = i;
907:               v[0]   = - hxhydhz*(dd - gd);
908:               c[1].k = k; c[1].j = j; c[1].i = i-1;
909:               v[1]   = - hyhzdhx*(dw - gw);
910:               c[2].k = k; c[2].j = j; c[2].i = i;
911:               v[2]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
912:               c[3].k = k; c[3].j = j+1; c[3].i = i;
913:               v[3]   = - hzhxdhy*(dn + gn);
914:               c[4].k = k+1; c[4].j = j; c[4].i = i;
915:               v[4]   = - hxhydhz*(du + gu);
916:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

918:             /* right-hand bottom up corner */
919:             } else {

921:               td = x[k-1][j][i];
922:               ad = 0.5*(t0 + td);
923:               bd = pow(ad,bm1);
924:               /* dd = bd * ad; */
925:               dd = pow(ad,beta);
926:               gd = coef*bd*(td - t0);
927: 
928:               c[0].k = k-1; c[0].j = j; c[0].i = i;
929:               v[0]   = - hxhydhz*(dd - gd);
930:               c[1].k = k; c[1].j = j; c[1].i = i-1;
931:               v[1]   = - hyhzdhx*(dw - gw);
932:               c[2].k = k; c[2].j = j; c[2].i = i;
933:               v[2]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
934:               c[3].k = k; c[3].j = j+1; c[3].i = i;
935:               v[3]   = - hzhxdhy*(dn + gn);
936:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
937:             }

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

942:             ts = x[k][j-1][i];
943:             as = 0.5*(t0 + ts);
944:             bs = pow(as,bm1);
945:             /* ds = bs * as; */
946:             ds = pow(as,beta);
947:             gs = coef*bs*(ts - t0);
948: 
949:             /* right-hand top down corner */
950:             if (k == 0) {

952:               tu = x[k+1][j][i];
953:               au = 0.5*(t0 + tu);
954:               bu = pow(au,bm1);
955:               /* du = bu * au; */
956:               du = pow(au,beta);
957:               gu = coef*bu*(tu - t0);
958: 
959:               c[0].k = k; c[0].j = j-1; c[0].i = i;
960:               v[0]   = - hzhxdhy*(ds - gs);
961:               c[1].k = k; c[1].j = j; c[1].i = i-1;
962:               v[1]   = - hyhzdhx*(dw - gw);
963:               c[2].k = k; c[2].j = j; c[2].i = i;
964:               v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
965:               c[3].k = k+1; c[3].j = j; c[3].i = i;
966:               v[3]   = - hxhydhz*(du + gu);
967:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);

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

972:               tu = x[k+1][j][i];
973:               au = 0.5*(t0 + tu);
974:               bu = pow(au,bm1);
975:               /* du = bu * au; */
976:               du = pow(au,beta);
977:               gu = coef*bu*(tu - t0);
978: 
979:               td = x[k-1][j][i];
980:               ad = 0.5*(t0 + td);
981:               bd = pow(ad,bm1);
982:               /* dd = bd * ad; */
983:               dd = pow(ad,beta);
984:               gd = coef*bd*(td - t0);
985: 
986:               c[0].k = k-1; c[0].j = j; c[0].i = i;
987:               v[0]   = - hxhydhz*(dd - gd);
988:               c[1].k = k; c[1].j = j-1; c[1].i = i;
989:               v[1]   = - hzhxdhy*(ds - gs);
990:               c[2].k = k; c[2].j = j; c[2].i = i-1;
991:               v[2]   = - hyhzdhx*(dw - gw);
992:               c[3].k = k; c[3].j = j; c[3].i = i;
993:               v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
994:               c[4].k = k+1; c[4].j = j; c[4].i = i;
995:               v[4]   = - hxhydhz*(du + gu);
996:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);

998:             /* right-hand top up corner */
999:             } else {

1001:               td = x[k-1][j][i];
1002:               ad = 0.5*(t0 + td);
1003:               bd = pow(ad,bm1);
1004:               /* dd = bd * ad; */
1005:               dd = pow(ad,beta);
1006:               gd = coef*bd*(td - t0);
1007: 
1008:               c[0].k = k-1; c[0].j = j; c[0].i = i;
1009:               v[0]   = - hxhydhz*(dd - gd);
1010:               c[1].k = k; c[1].j = j-1; c[1].i = i;
1011:               v[1]   = - hzhxdhy*(ds - gs);
1012:               c[2].k = k; c[2].j = j; c[2].i = i-1;
1013:               v[2]   = - hyhzdhx*(dw - gw);
1014:               c[3].k = k; c[3].j = j; c[3].i = i;
1015:               v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1016:                 MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
1017:             }

1019:           } else {

1021:             ts = x[k][j-1][i];
1022:             as = 0.5*(t0 + ts);
1023:             bs = pow(as,bm1);
1024:             /* ds = bs * as; */
1025:             ds = pow(as,beta);
1026:             gs = coef*bs*(t0 - ts);
1027: 
1028:             tn = x[k][j+1][i];
1029:             an = 0.5*(t0 + tn);
1030:             bn = pow(an,bm1);
1031:             /* dn = bn * an; */
1032:             dn = pow(an,beta);
1033:             gn = coef*bn*(tn - t0);

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

1038:               tu = x[k+1][j][i];
1039:               au = 0.5*(t0 + tu);
1040:               bu = pow(au,bm1);
1041:               /* du = bu * au; */
1042:               du = pow(au,beta);
1043:               gu = coef*bu*(tu - t0);

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

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

1061:               td = x[k-1][j][i];
1062:               ad = 0.5*(t0 + td);
1063:               bd = pow(ad,bm1);
1064:               /* dd = bd * ad; */
1065:               dd = pow(ad,beta);
1066:               gd = coef*bd*(t0 - td);
1067: 
1068:               c[0].k = k-1; c[0].j = j; c[0].i = i;
1069:               v[0]   = - hxhydhz*(dd - gd);
1070:               c[1].k = k; c[1].j = j-1; c[1].i = i;
1071:               v[1]   = - hzhxdhy*(ds - gs);
1072:               c[2].k = k; c[2].j = j; c[2].i = i-1;
1073:               v[2]   = - hyhzdhx*(dw - gw);
1074:               c[3].k = k; c[3].j = j; c[3].i = i;
1075:               v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1076:               c[4].k = k; c[4].j = j+1; c[4].i = i;
1077:               v[4]   = - hzhxdhy*(dn + gn);
1078:                 MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1079:             }

1081:             /* right-hand interior plane */
1082:             else {

1084:               td = x[k-1][j][i];
1085:               ad = 0.5*(t0 + td);
1086:               bd = pow(ad,bm1);
1087:               /* dd = bd * ad; */
1088:               dd = pow(ad,beta);
1089:               gd = coef*bd*(t0 - td);
1090: 
1091:               tu = x[k+1][j][i];
1092:               au = 0.5*(t0 + tu);
1093:               bu = pow(au,bm1);
1094:               /* du = bu * au; */
1095:               du = pow(au,beta);
1096:               gu = coef*bu*(tu - t0);

1098:               c[0].k = k-1; c[0].j = j; c[0].i = i;
1099:               v[0]   = - hxhydhz*(dd - gd);
1100:               c[1].k = k; c[1].j = j-1; c[1].i = i;
1101:               v[1]   = - hzhxdhy*(ds - gs);
1102:               c[2].k = k; c[2].j = j; c[2].i = i-1;
1103:               v[2]   = - hyhzdhx*(dw - gw);
1104:               c[3].k = k; c[3].j = j; c[3].i = i;
1105:               v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1106:               c[4].k = k; c[4].j = j+1; c[4].i = i;
1107:               v[4]   = - hzhxdhy*(dn + gn);
1108:               c[5].k = k+1; c[5].j = j; c[5].i = i;
1109:               v[5]   = - hxhydhz*(du + gu);
1110:                 MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1111:             }
1112:           }

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

1116:           tw = x[k][j][i-1];
1117:           aw = 0.5*(t0 + tw);
1118:           bw = pow(aw,bm1);
1119:           /* dw = bw * aw */
1120:           dw = pow(aw,beta);
1121:           gw = coef*bw*(t0 - tw);

1123:           te = x[k][j][i+1];
1124:           ae = 0.5*(t0 + te);
1125:           be = pow(ae,bm1);
1126:           /* de = be * ae; */
1127:           de = pow(ae,beta);
1128:           ge = coef*be*(te - t0);

1130:           tn = x[k][j+1][i];
1131:           an = 0.5*(t0 + tn);
1132:           bn = pow(an,bm1);
1133:           /* dn = bn * an; */
1134:           dn = pow(an,beta);
1135:           gn = coef*bn*(tn - t0);


1138:           /* bottom down interior edge */
1139:           if (k == 0) {

1141:             tu = x[k+1][j][i];
1142:             au = 0.5*(t0 + tu);
1143:             bu = pow(au,bm1);
1144:             /* du = bu * au; */
1145:             du = pow(au,beta);
1146:             gu = coef*bu*(tu - t0);

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

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

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

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

1184:           /* bottom interior plane */
1185:           else {

1187:             tu = x[k+1][j][i];
1188:             au = 0.5*(t0 + tu);
1189:             bu = pow(au,bm1);
1190:             /* du = bu * au; */
1191:             du = pow(au,beta);
1192:             gu = coef*bu*(tu - t0);

1194:             td = x[k-1][j][i];
1195:             ad = 0.5*(t0 + td);
1196:             bd = pow(ad,bm1);
1197:             /* dd = bd * ad; */
1198:             dd = pow(ad,beta);
1199:             gd = coef*bd*(td - t0);

1201:             c[0].k = k-1; c[0].j = j; c[0].i = i;
1202:             v[0]   = - hxhydhz*(dd - gd);
1203:             c[1].k = k; c[1].j = j; c[1].i = i-1;
1204:             v[1]   = - hyhzdhx*(dw - gw);
1205:             c[2].k = k; c[2].j = j; c[2].i = i;
1206:             v[2]   =   hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1207:             c[3].k = k; c[3].j = j; c[3].i = i+1;
1208:             v[3]   = - hyhzdhx*(de + ge);
1209:             c[4].k = k; c[4].j = j+1; c[4].i = i;
1210:             v[4]   = - hzhxdhy*(dn + gn);
1211:             c[5].k = k+1; c[5].j = j; c[5].i = i;
1212:             v[5]   = - hxhydhz*(du + gu);
1213:               MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1214:           }

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

1218:           tw = x[k][j][i-1];
1219:           aw = 0.5*(t0 + tw);
1220:           bw = pow(aw,bm1);
1221:           /* dw = bw * aw */
1222:           dw = pow(aw,beta);
1223:           gw = coef*bw*(t0 - tw);

1225:           te = x[k][j][i+1];
1226:           ae = 0.5*(t0 + te);
1227:           be = pow(ae,bm1);
1228:           /* de = be * ae; */
1229:           de = pow(ae,beta);
1230:           ge = coef*be*(te - t0);

1232:           ts = x[k][j-1][i];
1233:           as = 0.5*(t0 + ts);
1234:           bs = pow(as,bm1);
1235:           /* ds = bs * as; */
1236:           ds = pow(as,beta);
1237:           gs = coef*bs*(t0 - ts);
1238: 
1239:           /* top down interior edge */
1240:           if (k == 0) {

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

1249:             c[0].k = k; c[0].j = j-1; c[0].i = i;
1250:             v[0]   = - hzhxdhy*(ds - gs);
1251:             c[1].k = k; c[1].j = j; c[1].i = i-1;
1252:             v[1]   = - hyhzdhx*(dw - gw);
1253:             c[2].k = k; c[2].j = j; c[2].i = i;
1254:             v[2]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1255:             c[3].k = k; c[3].j = j; c[3].i = i+1;
1256:             v[3]   = - hyhzdhx*(de + ge);
1257:             c[4].k = k+1; c[4].j = j; c[4].i = i;
1258:             v[4]   = - hxhydhz*(du + gu);
1259:               MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1260:           }

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

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

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

1285:           /* top interior plane */
1286:           else {

1288:             tu = x[k+1][j][i];
1289:             au = 0.5*(t0 + tu);
1290:             bu = pow(au,bm1);
1291:             /* du = bu * au; */
1292:             du = pow(au,beta);
1293:             gu = coef*bu*(tu - t0);

1295:             td = x[k-1][j][i];
1296:             ad = 0.5*(t0 + td);
1297:             bd = pow(ad,bm1);
1298:             /* dd = bd * ad; */
1299:             dd = pow(ad,beta);
1300:             gd = coef*bd*(td - t0);

1302:             c[0].k = k-1; c[0].j = j; c[0].i = i;
1303:             v[0]   = - hxhydhz*(dd - gd);
1304:             c[1].k = k; c[1].j = j-1; c[1].i = i;
1305:             v[1]   = - hzhxdhy*(ds - gs);
1306:             c[2].k = k; c[2].j = j; c[2].i = i-1;
1307:             v[2]   = - hyhzdhx*(dw - gw);
1308:             c[3].k = k; c[3].j = j; c[3].i = i;
1309:             v[3]   =   hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1310:             c[4].k = k; c[4].j = j; c[4].i = i+1;
1311:             v[4]   = - hyhzdhx*(de + ge);
1312:             c[5].k = k+1; c[5].j = j; c[5].i = i;
1313:             v[5]   = - hxhydhz*(du + gu);
1314:               MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1315:           }

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

1319:           /* down interior plane */

1321:           tw = x[k][j][i-1];
1322:           aw = 0.5*(t0 + tw);
1323:           bw = pow(aw,bm1);
1324:           /* dw = bw * aw */
1325:           dw = pow(aw,beta);
1326:           gw = coef*bw*(t0 - tw);

1328:           te = x[k][j][i+1];
1329:           ae = 0.5*(t0 + te);
1330:           be = pow(ae,bm1);
1331:           /* de = be * ae; */
1332:           de = pow(ae,beta);
1333:           ge = coef*be*(te - t0);

1335:           ts = x[k][j-1][i];
1336:           as = 0.5*(t0 + ts);
1337:           bs = pow(as,bm1);
1338:           /* ds = bs * as; */
1339:           ds = pow(as,beta);
1340:           gs = coef*bs*(t0 - ts);
1341: 
1342:           tn = x[k][j+1][i];
1343:           an = 0.5*(t0 + tn);
1344:           bn = pow(an,bm1);
1345:           /* dn = bn * an; */
1346:           dn = pow(an,beta);
1347:           gn = coef*bn*(tn - t0);
1348: 
1349:           tu = x[k+1][j][i];
1350:           au = 0.5*(t0 + tu);
1351:           bu = pow(au,bm1);
1352:           /* du = bu * au; */
1353:           du = pow(au,beta);
1354:           gu = coef*bu*(tu - t0);

1356:           c[0].k = k; c[0].j = j-1; c[0].i = i;
1357:           v[0]   = - hzhxdhy*(ds - gs);
1358:           c[1].k = k; c[1].j = j; c[1].i = i-1;
1359:           v[1]   = - hyhzdhx*(dw - gw);
1360:           c[2].k = k; c[2].j = j; c[2].i = i;
1361:           v[2]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1362:           c[3].k = k; c[3].j = j; c[3].i = i+1;
1363:           v[3]   = - hyhzdhx*(de + ge);
1364:           c[4].k = k; c[4].j = j+1; c[4].i = i;
1365:           v[4]   = - hzhxdhy*(dn + gn);
1366:           c[5].k = k+1; c[5].j = j; c[5].i = i;
1367:           v[5]   = - hxhydhz*(du + gu);
1368:             MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1369:         }
1370: 
1371:         else if (k == mz-1) {

1373:           /* up interior plane */

1375:           tw = x[k][j][i-1];
1376:           aw = 0.5*(t0 + tw);
1377:           bw = pow(aw,bm1);
1378:           /* dw = bw * aw */
1379:           dw = pow(aw,beta);
1380:           gw = coef*bw*(t0 - tw);

1382:           te = x[k][j][i+1];
1383:           ae = 0.5*(t0 + te);
1384:           be = pow(ae,bm1);
1385:           /* de = be * ae; */
1386:           de = pow(ae,beta);
1387:           ge = coef*be*(te - t0);

1389:           ts = x[k][j-1][i];
1390:           as = 0.5*(t0 + ts);
1391:           bs = pow(as,bm1);
1392:           /* ds = bs * as; */
1393:           ds = pow(as,beta);
1394:           gs = coef*bs*(t0 - ts);
1395: 
1396:           tn = x[k][j+1][i];
1397:           an = 0.5*(t0 + tn);
1398:           bn = pow(an,bm1);
1399:           /* dn = bn * an; */
1400:           dn = pow(an,beta);
1401:           gn = coef*bn*(tn - t0);
1402: 
1403:           td = x[k-1][j][i];
1404:           ad = 0.5*(t0 + td);
1405:           bd = pow(ad,bm1);
1406:           /* dd = bd * ad; */
1407:           dd = pow(ad,beta);
1408:           gd = coef*bd*(t0 - td);
1409: 
1410:           c[0].k = k-1; c[0].j = j; c[0].i = i;
1411:           v[0]   = - hxhydhz*(dd - gd);
1412:           c[1].k = k; c[1].j = j-1; c[1].i = i;
1413:           v[1]   = - hzhxdhy*(ds - gs);
1414:           c[2].k = k; c[2].j = j; c[2].i = i-1;
1415:           v[2]   = - hyhzdhx*(dw - gw);
1416:           c[3].k = k; c[3].j = j; c[3].i = i;
1417:           v[3]   =   hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1418:           c[4].k = k; c[4].j = j; c[4].i = i+1;
1419:           v[4]   = - hyhzdhx*(de + ge);
1420:           c[5].k = k; c[5].j = j+1; c[5].i = i;
1421:           v[5]   = - hzhxdhy*(dn + gn);
1422:             MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1423:         }
1424:       }
1425:     }
1426:   }
1427:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
1428:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
1429:   if (jac != *J) {
1430:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1431:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1432:   }
1433:   DAVecRestoreArray((DA)dmmg->dm,localX,&x);
1434:   DARestoreLocalVector((DA)dmmg->dm,&localX);

1436:   PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
1437:   return(0);
1438: }