Actual source code: ex18.c

      2: static char help[] ="Nonlinear Radiative Transport PDE with multigrid in 2d.\n\
      3: Uses 2-dimensional distributed arrays.\n\
      4: A 2-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: DMDA^using distributed arrays
     17:    Concepts: multigrid;
     18:    Processors: n
     19: T*/

     23: /*

     25:     This example models the partial differential equation

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

     29:     where beta = 2.5 and alpha = 1.0

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

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

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

     40:     This code was contributed by David Keyes

     42: */

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

     48: /* User-defined application context */

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

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

     57: extern PetscErrorCode FormInitialGuess(SNES,Vec,void*);
     58: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
     59: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);

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

     70:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;

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

     82:   /*
     83:       Create the multilevel DM data structure
     84:   */
     85:   SNESCreate(PETSC_COMM_WORLD,&snes);

     87:   /*
     88:       Set the DMDA (grid structure) for the grids.
     89:   */
     90:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,5,5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
     91:   DMSetFromOptions(da);
     92:   DMSetUp(da);
     93:   DMSetApplicationContext(da,&user);
     94:   SNESSetDM(snes,(DM)da);

     96:   /*
     97:      Create the nonlinear solver, and tell it the functions to use
     98:   */
     99:   SNESSetFunction(snes,NULL,FormFunction,&user);
    100:   SNESSetJacobian(snes,NULL,NULL,FormJacobian,&user);
    101:   SNESSetFromOptions(snes);
    102:   SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);

    104:   SNESSolve(snes,NULL,NULL);
    105:   SNESGetIterationNumber(snes,&its);
    106:   SNESGetLinearSolveIterations(snes,&lits);
    107:   litspit = ((PetscReal)lits)/((PetscReal)its);
    108:   PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
    109:   PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
    110:   PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",(double)litspit);

    112:   DMDestroy(&da);
    113:   SNESDestroy(&snes);
    114:   PetscFinalize();
    115:   return ierr;
    116: }
    117: /* --------------------  Form initial approximation ----------------- */
    118: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
    119: {
    120:   AppCtx         *user;
    121:   PetscInt       i,j,xs,ys,xm,ym;
    123:   PetscReal      tleft;
    124:   PetscScalar    **x;
    125:   DM             da;

    128:   SNESGetDM(snes,&da);
    129:   DMGetApplicationContext(da,&user);
    130:   tleft = user->tleft;
    131:   /* Get ghost points */
    132:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
    133:   DMDAVecGetArray(da,X,&x);

    135:   /* Compute initial guess */
    136:   for (j=ys; j<ys+ym; j++) {
    137:     for (i=xs; i<xs+xm; i++) {
    138:       x[j][i] = tleft;
    139:     }
    140:   }
    141:   DMDAVecRestoreArray(da,X,&x);
    142:   return(0);
    143: }
    144: /* --------------------  Evaluate Function F(x) --------------------- */
    145: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
    146: {
    147:   AppCtx         *user = (AppCtx*)ptr;
    149:   PetscInt       i,j,mx,my,xs,ys,xm,ym;
    150:   PetscScalar    zero = 0.0,one = 1.0;
    151:   PetscScalar    hx,hy,hxdhy,hydhx;
    152:   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;
    153:   PetscScalar    tleft,tright,beta;
    154:   PetscScalar    **x,**f;
    155:   Vec            localX;
    156:   DM             da;

    159:   SNESGetDM(snes,&da);
    160:   DMGetLocalVector(da,&localX);
    161:   DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
    162:   hx    = one/(PetscReal)(mx-1);  hy    = one/(PetscReal)(my-1);
    163:   hxdhy = hx/hy;               hydhx = hy/hx;
    164:   tleft = user->tleft;         tright = user->tright;
    165:   beta  = user->beta;

    167:   /* Get ghost points */
    168:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
    169:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
    170:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
    171:   DMDAVecGetArray(da,localX,&x);
    172:   DMDAVecGetArray(da,F,&f);

    174:   /* Evaluate function */
    175:   for (j=ys; j<ys+ym; j++) {
    176:     for (i=xs; i<xs+xm; i++) {
    177:       t0 = x[j][i];

    179:       if (i > 0 && i < mx-1 && j > 0 && j < my-1) {

    181:         /* general interior volume */

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

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

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

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

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

    205:         /* left-hand boundary */
    206:         tw = tleft;
    207:         aw = 0.5*(t0 + tw);
    208:         dw = PetscPowScalar(aw,beta);
    209:         fw = dw*(t0 - tw);

    211:         te = x[j][i+1];
    212:         ae = 0.5*(t0 + te);
    213:         de = PetscPowScalar(ae,beta);
    214:         fe = de*(te - t0);

    216:         if (j > 0) {
    217:           ts = x[j-1][i];
    218:           as = 0.5*(t0 + ts);
    219:           ds = PetscPowScalar(as,beta);
    220:           fs = ds*(t0 - ts);
    221:         } else fs = zero;

    223:         if (j < my-1) {
    224:           tn = x[j+1][i];
    225:           an = 0.5*(t0 + tn);
    226:           dn = PetscPowScalar(an,beta);
    227:           fn = dn*(tn - t0);
    228:         } else fn = zero;

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

    232:         /* right-hand boundary */
    233:         tw = x[j][i-1];
    234:         aw = 0.5*(t0 + tw);
    235:         dw = PetscPowScalar(aw,beta);
    236:         fw = dw*(t0 - tw);

    238:         te = tright;
    239:         ae = 0.5*(t0 + te);
    240:         de = PetscPowScalar(ae,beta);
    241:         fe = de*(te - t0);

    243:         if (j > 0) {
    244:           ts = x[j-1][i];
    245:           as = 0.5*(t0 + ts);
    246:           ds = PetscPowScalar(as,beta);
    247:           fs = ds*(t0 - ts);
    248:         } else fs = zero;

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

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

    259:         /* bottom boundary,and i <> 0 or mx-1 */
    260:         tw = x[j][i-1];
    261:         aw = 0.5*(t0 + tw);
    262:         dw = PetscPowScalar(aw,beta);
    263:         fw = dw*(t0 - tw);

    265:         te = x[j][i+1];
    266:         ae = 0.5*(t0 + te);
    267:         de = PetscPowScalar(ae,beta);
    268:         fe = de*(te - t0);

    270:         fs = zero;

    272:         tn = x[j+1][i];
    273:         an = 0.5*(t0 + tn);
    274:         dn = PetscPowScalar(an,beta);
    275:         fn = dn*(tn - t0);

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

    279:         /* top boundary,and i <> 0 or mx-1 */
    280:         tw = x[j][i-1];
    281:         aw = 0.5*(t0 + tw);
    282:         dw = PetscPowScalar(aw,beta);
    283:         fw = dw*(t0 - tw);

    285:         te = x[j][i+1];
    286:         ae = 0.5*(t0 + te);
    287:         de = PetscPowScalar(ae,beta);
    288:         fe = de*(te - t0);

    290:         ts = x[j-1][i];
    291:         as = 0.5*(t0 + ts);
    292:         ds = PetscPowScalar(as,beta);
    293:         fs = ds*(t0 - ts);

    295:         fn = zero;

    297:       }

    299:       f[j][i] = -hydhx*(fe-fw) - hxdhy*(fn-fs);

    301:     }
    302:   }
    303:   DMDAVecRestoreArray(da,localX,&x);
    304:   DMDAVecRestoreArray(da,F,&f);
    305:   DMRestoreLocalVector(da,&localX);
    306:   PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
    307:   return(0);
    308: }
    309: /* --------------------  Evaluate Jacobian F(x) --------------------- */
    310: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat jac,Mat B,void *ptr)
    311: {
    312:   AppCtx         *user = (AppCtx*)ptr;
    314:   PetscInt       i,j,mx,my,xs,ys,xm,ym;
    315:   PetscScalar    one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;
    316:   PetscScalar    dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;
    317:   PetscScalar    tleft,tright,beta,bm1,coef;
    318:   PetscScalar    v[5],**x;
    319:   Vec            localX;
    320:   MatStencil     col[5],row;
    321:   DM             da;

    324:   SNESGetDM(snes,&da);
    325:   DMGetLocalVector(da,&localX);
    326:   DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
    327:   hx    = one/(PetscReal)(mx-1);  hy     = one/(PetscReal)(my-1);
    328:   hxdhy = hx/hy;               hydhx  = hy/hx;
    329:   tleft = user->tleft;         tright = user->tright;
    330:   beta  = user->beta;          bm1    = user->bm1;          coef = user->coef;

    332:   /* Get ghost points */
    333:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
    334:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
    335:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
    336:   DMDAVecGetArray(da,localX,&x);

    338:   /* Evaluate Jacobian of function */
    339:   for (j=ys; j<ys+ym; j++) {
    340:     for (i=xs; i<xs+xm; i++) {
    341:       t0 = x[j][i];

    343:       if (i > 0 && i < mx-1 && j > 0 && j < my-1) {

    345:         /* general interior volume */

    347:         tw = x[j][i-1];
    348:         aw = 0.5*(t0 + tw);
    349:         bw = PetscPowScalar(aw,bm1);
    350:         /* dw = bw * aw */
    351:         dw = PetscPowScalar(aw,beta);
    352:         gw = coef*bw*(t0 - tw);

    354:         te = x[j][i+1];
    355:         ae = 0.5*(t0 + te);
    356:         be = PetscPowScalar(ae,bm1);
    357:         /* de = be * ae; */
    358:         de = PetscPowScalar(ae,beta);
    359:         ge = coef*be*(te - t0);

    361:         ts = x[j-1][i];
    362:         as = 0.5*(t0 + ts);
    363:         bs = PetscPowScalar(as,bm1);
    364:         /* ds = bs * as; */
    365:         ds = PetscPowScalar(as,beta);
    366:         gs = coef*bs*(t0 - ts);

    368:         tn = x[j+1][i];
    369:         an = 0.5*(t0 + tn);
    370:         bn = PetscPowScalar(an,bm1);
    371:         /* dn = bn * an; */
    372:         dn = PetscPowScalar(an,beta);
    373:         gn = coef*bn*(tn - t0);

    375:         v[0] = -hxdhy*(ds - gs);                                       col[0].j = j-1;       col[0].i = i;
    376:         v[1] = -hydhx*(dw - gw);                                       col[1].j = j;         col[1].i = i-1;
    377:         v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
    378:         v[3] = -hydhx*(de + ge);                                       col[3].j = j;         col[3].i = i+1;
    379:         v[4] = -hxdhy*(dn + gn);                                       col[4].j = j+1;       col[4].i = i;
    380:         MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);

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

    384:         /* left-hand boundary */
    385:         tw = tleft;
    386:         aw = 0.5*(t0 + tw);
    387:         bw = PetscPowScalar(aw,bm1);
    388:         /* dw = bw * aw */
    389:         dw = PetscPowScalar(aw,beta);
    390:         gw = coef*bw*(t0 - tw);

    392:         te = x[j][i + 1];
    393:         ae = 0.5*(t0 + te);
    394:         be = PetscPowScalar(ae,bm1);
    395:         /* de = be * ae; */
    396:         de = PetscPowScalar(ae,beta);
    397:         ge = coef*be*(te - t0);

    399:         /* left-hand bottom boundary */
    400:         if (j == 0) {

    402:           tn = x[j+1][i];
    403:           an = 0.5*(t0 + tn);
    404:           bn = PetscPowScalar(an,bm1);
    405:           /* dn = bn * an; */
    406:           dn = PetscPowScalar(an,beta);
    407:           gn = coef*bn*(tn - t0);

    409:           v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i;
    410:           v[1] = -hydhx*(de + ge);                            col[1].j = j;         col[1].i = i+1;
    411:           v[2] = -hxdhy*(dn + gn);                            col[2].j = j+1;       col[2].i = i;
    412:           MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);

    414:           /* left-hand interior boundary */
    415:         } else if (j < my-1) {

    417:           ts = x[j-1][i];
    418:           as = 0.5*(t0 + ts);
    419:           bs = PetscPowScalar(as,bm1);
    420:           /* ds = bs * as; */
    421:           ds = PetscPowScalar(as,beta);
    422:           gs = coef*bs*(t0 - ts);

    424:           tn = x[j+1][i];
    425:           an = 0.5*(t0 + tn);
    426:           bn = PetscPowScalar(an,bm1);
    427:           /* dn = bn * an; */
    428:           dn = PetscPowScalar(an,beta);
    429:           gn = coef*bn*(tn - t0);

    431:           v[0] = -hxdhy*(ds - gs);                                       col[0].j = j-1;       col[0].i = i;
    432:           v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[1].j = row.j = j; col[1].i = row.i = i;
    433:           v[2] = -hydhx*(de + ge);                                       col[2].j = j;         col[2].i = i+1;
    434:           v[3] = -hxdhy*(dn + gn);                                       col[3].j = j+1;       col[3].i = i;
    435:           MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
    436:           /* left-hand top boundary */
    437:         } else {

    439:           ts = x[j-1][i];
    440:           as = 0.5*(t0 + ts);
    441:           bs = PetscPowScalar(as,bm1);
    442:           /* ds = bs * as; */
    443:           ds = PetscPowScalar(as,beta);
    444:           gs = coef*bs*(t0 - ts);

    446:           v[0] = -hxdhy*(ds - gs);                             col[0].j = j-1;       col[0].i = i;
    447:           v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[1].j = row.j = j; col[1].i = row.i = i;
    448:           v[2] = -hydhx*(de + ge);                             col[2].j = j;         col[2].i = i+1;
    449:           MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
    450:         }

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

    454:         /* right-hand boundary */
    455:         tw = x[j][i-1];
    456:         aw = 0.5*(t0 + tw);
    457:         bw = PetscPowScalar(aw,bm1);
    458:         /* dw = bw * aw */
    459:         dw = PetscPowScalar(aw,beta);
    460:         gw = coef*bw*(t0 - tw);

    462:         te = tright;
    463:         ae = 0.5*(t0 + te);
    464:         be = PetscPowScalar(ae,bm1);
    465:         /* de = be * ae; */
    466:         de = PetscPowScalar(ae,beta);
    467:         ge = coef*be*(te - t0);

    469:         /* right-hand bottom boundary */
    470:         if (j == 0) {

    472:           tn = x[j+1][i];
    473:           an = 0.5*(t0 + tn);
    474:           bn = PetscPowScalar(an,bm1);
    475:           /* dn = bn * an; */
    476:           dn = PetscPowScalar(an,beta);
    477:           gn = coef*bn*(tn - t0);

    479:           v[0] = -hydhx*(dw - gw);                            col[0].j = j;         col[0].i = i-1;
    480:           v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
    481:           v[2] = -hxdhy*(dn + gn);                            col[2].j = j+1;       col[2].i = i;
    482:           MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);

    484:           /* right-hand interior boundary */
    485:         } else if (j < my-1) {

    487:           ts = x[j-1][i];
    488:           as = 0.5*(t0 + ts);
    489:           bs = PetscPowScalar(as,bm1);
    490:           /* ds = bs * as; */
    491:           ds = PetscPowScalar(as,beta);
    492:           gs = coef*bs*(t0 - ts);

    494:           tn = x[j+1][i];
    495:           an = 0.5*(t0 + tn);
    496:           bn = PetscPowScalar(an,bm1);
    497:           /* dn = bn * an; */
    498:           dn = PetscPowScalar(an,beta);
    499:           gn = coef*bn*(tn - t0);

    501:           v[0] = -hxdhy*(ds - gs);                                       col[0].j = j-1;       col[0].i = i;
    502:           v[1] = -hydhx*(dw - gw);                                       col[1].j = j;         col[1].i = i-1;
    503:           v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
    504:           v[3] = -hxdhy*(dn + gn);                                       col[3].j = j+1;       col[3].i = i;
    505:           MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
    506:         /* right-hand top boundary */
    507:         } else {

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

    516:           v[0] = -hxdhy*(ds - gs);                             col[0].j = j-1;       col[0].i = i;
    517:           v[1] = -hydhx*(dw - gw);                             col[1].j = j;         col[1].i = i-1;
    518:           v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
    519:           MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
    520:         }

    522:         /* bottom boundary,and i <> 0 or mx-1 */
    523:       } else if (j == 0) {

    525:         tw = x[j][i-1];
    526:         aw = 0.5*(t0 + tw);
    527:         bw = PetscPowScalar(aw,bm1);
    528:         /* dw = bw * aw */
    529:         dw = PetscPowScalar(aw,beta);
    530:         gw = coef*bw*(t0 - tw);

    532:         te = x[j][i+1];
    533:         ae = 0.5*(t0 + te);
    534:         be = PetscPowScalar(ae,bm1);
    535:         /* de = be * ae; */
    536:         de = PetscPowScalar(ae,beta);
    537:         ge = coef*be*(te - t0);

    539:         tn = x[j+1][i];
    540:         an = 0.5*(t0 + tn);
    541:         bn = PetscPowScalar(an,bm1);
    542:         /* dn = bn * an; */
    543:         dn = PetscPowScalar(an,beta);
    544:         gn = coef*bn*(tn - t0);

    546:         v[0] = -hydhx*(dw - gw);                            col[0].j = j;         col[0].i = i-1;
    547:         v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
    548:         v[2] = -hydhx*(de + ge);                            col[2].j = j;         col[2].i = i+1;
    549:         v[3] = -hxdhy*(dn + gn);                            col[3].j = j+1;       col[3].i = i;
    550:         MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);

    552:         /* top boundary,and i <> 0 or mx-1 */
    553:       } else if (j == my-1) {

    555:         tw = x[j][i-1];
    556:         aw = 0.5*(t0 + tw);
    557:         bw = PetscPowScalar(aw,bm1);
    558:         /* dw = bw * aw */
    559:         dw = PetscPowScalar(aw,beta);
    560:         gw = coef*bw*(t0 - tw);

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

    569:         ts = x[j-1][i];
    570:         as = 0.5*(t0 + ts);
    571:         bs = PetscPowScalar(as,bm1);
    572:         /* ds = bs * as; */
    573:         ds = PetscPowScalar(as,beta);
    574:         gs = coef*bs*(t0 - ts);

    576:         v[0] = -hxdhy*(ds - gs);                             col[0].j = j-1;       col[0].i = i;
    577:         v[1] = -hydhx*(dw - gw);                             col[1].j = j;         col[1].i = i-1;
    578:         v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
    579:         v[3] = -hydhx*(de + ge);                             col[3].j = j;         col[3].i = i+1;
    580:         MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);

    582:       }
    583:     }
    584:   }
    585:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    586:   DMDAVecRestoreArray(da,localX,&x);
    587:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    588:   DMRestoreLocalVector(da,&localX);
    589:   if (jac != B) {
    590:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
    591:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
    592:   }

    594:   PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
    595:   return(0);
    596: }

    600: /*TEST

    602:    test:
    603:       args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view
    604:       requires: !single

    606: TEST*/