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);

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

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

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


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

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

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

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

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

193:             /* general interior volume */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

336:           fs = zero;

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

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

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

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

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

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

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

379:           fn = zero;

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

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

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

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

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

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

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

422:            fd = zero;

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

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

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

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

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

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

457:           fu = zero;
458:         }

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

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

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

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

510:           /* general interior volume */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

746:           } else {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1018:           } else {

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1318:           /* down interior plane */

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

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

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

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

1372:           /* up interior plane */

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

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

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

1431:   PetscLogFlops((41 + 8*POWFLOP)*xm*ym);
1432:   return(0);
1433: }