Actual source code: ex20.c
petsc-3.4.2 2013-07-02

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

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

     21: /*

     23:     This example models the partial differential equation

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

     27:     where beta = 2.5 and alpha = 1.0

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

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

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

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

     40: */

     42: #include <petscsnes.h>
     43: #include <petscdmda.h>

     45: /* User-defined application context */

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

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

     54: extern PetscErrorCode FormInitialGuess(SNES,Vec,void*);
     55: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
     56: extern PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);

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

     69:   PetscInitialize(&argc,&argv,NULL,help);

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

     81:   /*
     82:       Set the DMDA (grid structure) for the grids.
     83:   */
     84:   DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-5,-5,-5,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);
     85:   DMSetApplicationContext(da,&user);

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

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

    105:   SNESDestroy(&snes);
    106:   DMDestroy(&da);
    107:   PetscFinalize();

    109:   return 0;
    110: }
    111: /* --------------------  Form initial approximation ----------------- */
    114: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
    115: {
    116:   AppCtx         *user;
    117:   PetscInt       i,j,k,xs,ys,xm,ym,zs,zm;
    119:   PetscScalar    ***x;
    120:   DM             da;

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

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

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

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

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

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

    179:           /* general interior volume */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    322:           fs = zero;

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

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

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

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

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

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

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

    365:           fn = zero;

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

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

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

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

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

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

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

    408:           fd = zero;

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

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

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

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

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

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

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

    443:           fu = zero;
    444:         }

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

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

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

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

    497:           /* general interior volume */

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

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

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

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

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

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

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

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

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

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

    573:           /* left-hand bottom edge */
    574:           if (j == 0) {

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

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

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

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

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

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

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

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

    632:               /* left-hand bottom up corner */
    633:             } else {

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

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

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

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

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

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

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

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

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

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

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

    712:               /* left-hand top up corner */
    713:             } else {

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

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

    733:           } else {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    840:           /* right-hand bottom edge */
    841:           if (j == 0) {

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

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

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

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

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

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

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

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

    899:               /* right-hand bottom up corner */
    900:             } else {

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

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

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

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

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

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

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

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

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

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

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

    979:               /* right-hand top up corner */
    980:             } else {

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

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

    1000:           } else {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    1115:           /* bottom down interior edge */
    1116:           if (k == 0) {

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

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

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

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

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

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

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

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

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

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

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

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

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

    1212:           /* top down interior edge */
    1213:           if (k == 0) {

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

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

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

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

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

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

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

    1264:             td = x[k-1][j][i];
    1265:             ad = 0.5*(t0 + td);
    1266:             bd = PetscPowScalar(ad,bm1);
    1267:             /* dd = bd * ad; */
    1268:             dd = PetscPowScalar(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 + du + gd - gu);
    1279:             c[4].k = k; c[4].j = j; c[4].i = i+1;
    1280:             v[4]   = -hyhzdhx*(de + ge);
    1281:             c[5].k = k+1; c[5].j = j; c[5].i = i;
    1282:             v[5]   = -hxhydhz*(du + gu);
    1283:               MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
    1284:           }

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

    1288:           /* down interior plane */

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

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

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

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

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

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

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

    1341:           /* up interior plane */

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

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

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

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

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

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

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