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: 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 = 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 5-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 David Keyes
 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,3,&user,&dmmg);

 90:   /*
 91:       Set the DA (grid structure) for the grids.
 92:   */
 93:   DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,5,5,PETSC_DECIDE,PETSC_DECIDE,1,1,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,xs,ys,xm,ym;
134:   PetscReal      tleft = user->tleft;
135:   PetscScalar    **x;


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

143:   /* Compute initial guess */
144:   for (j=ys; j<ys+ym; j++) {
145:     for (i=xs; i<xs+xm; i++) {
146:       x[j][i] = tleft;
147:     }
148:   }
149:   DAVecRestoreArray((DA)dmmg->dm,X,&x);
150:   return(0);
151: }
152: /* --------------------  Evaluate Function F(x) --------------------- */
155: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ptr)
156: {
157:   DMMG           dmmg = (DMMG)ptr;
158:   AppCtx         *user = (AppCtx*)dmmg->user;
160:   PetscInt       i,j,mx,my,xs,ys,xm,ym;
161:   PetscScalar    zero = 0.0,one = 1.0;
162:   PetscScalar    hx,hy,hxdhy,hydhx;
163:   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;
164:   PetscScalar    tleft,tright,beta;
165:   PetscScalar    **x,**f;
166:   Vec            localX;

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

183:   /* Evaluate function */
184:   for (j=ys; j<ys+ym; j++) {
185:     for (i=xs; i<xs+xm; i++) {
186:       t0 = x[j][i];

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

190:         /* general interior volume */

192:         tw = x[j][i-1];
193:         aw = 0.5*(t0 + tw);
194:         dw = PetscPowScalar(aw,beta);
195:         fw = dw*(t0 - tw);

197:         te = x[j][i+1];
198:         ae = 0.5*(t0 + te);
199:         de = PetscPowScalar(ae,beta);
200:         fe = de*(te - t0);

202:         ts = x[j-1][i];
203:         as = 0.5*(t0 + ts);
204:         ds = PetscPowScalar(as,beta);
205:         fs = ds*(t0 - ts);
206: 
207:         tn = x[j+1][i];
208:         an = 0.5*(t0 + tn);
209:         dn = PetscPowScalar(an,beta);
210:         fn = dn*(tn - t0);

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

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

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

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

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

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

245:         /* right-hand boundary */
246:         tw = x[j][i-1];
247:         aw = 0.5*(t0 + tw);
248:         dw = PetscPowScalar(aw,beta);
249:         fw = dw*(t0 - tw);
250: 
251:         te = tright;
252:         ae = 0.5*(t0 + te);
253:         de = PetscPowScalar(ae,beta);
254:         fe = de*(te - t0);
255: 
256:         if (j > 0) {
257:           ts = x[j-1][i];
258:           as = 0.5*(t0 + ts);
259:           ds = PetscPowScalar(as,beta);
260:           fs = ds*(t0 - ts);
261:         } else {
262:           fs = zero;
263:         }
264: 
265:         if (j < my-1) {
266:           tn = x[j+1][i];
267:           an = 0.5*(t0 + tn);
268:           dn = PetscPowScalar(an,beta);
269:           fn = dn*(tn - t0);
270:         } else {
271:           fn = zero;
272:         }

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

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

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

287:         fs = zero;

289:         tn = x[j+1][i];
290:         an = 0.5*(t0 + tn);
291:         dn = PetscPowScalar(an,beta);
292:         fn = dn*(tn - t0);

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

296:         /* top boundary,and i <> 0 or mx-1 */
297:         tw = x[j][i-1];
298:         aw = 0.5*(t0 + tw);
299:         dw = PetscPowScalar(aw,beta);
300:         fw = dw*(t0 - tw);

302:         te = x[j][i+1];
303:         ae = 0.5*(t0 + te);
304:         de = PetscPowScalar(ae,beta);
305:         fe = de*(te - t0);

307:         ts = x[j-1][i];
308:         as = 0.5*(t0 + ts);
309:         ds = PetscPowScalar(as,beta);
310:         fs = ds*(t0 - ts);

312:         fn = zero;

314:       }

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

318:     }
319:   }
320:   DAVecRestoreArray((DA)dmmg->dm,localX,&x);
321:   DAVecRestoreArray((DA)dmmg->dm,F,&f);
322:   DARestoreLocalVector((DA)dmmg->dm,&localX);
323:   PetscLogFlops((22 + 4*POWFLOP)*ym*xm);
324:   return(0);
325: }
326: /* --------------------  Evaluate Jacobian F(x) --------------------- */
329: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ptr)
330: {
331:   DMMG           dmmg = (DMMG)ptr;
332:   AppCtx         *user = (AppCtx*)dmmg->user;
333:   Mat            jac = *J;
335:   PetscInt       i,j,mx,my,xs,ys,xm,ym;
336:   PetscScalar    one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;
337:   PetscScalar    dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;
338:   PetscScalar    tleft,tright,beta,bm1,coef;
339:   PetscScalar    v[5],**x;
340:   Vec            localX;
341:   MatStencil     col[5],row;

344:   DAGetLocalVector((DA)dmmg->dm,&localX);
345:   *flg = SAME_NONZERO_PATTERN;
346:   DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0);
347:   hx    = one/(PetscReal)(mx-1);  hy     = one/(PetscReal)(my-1);
348:   hxdhy = hx/hy;               hydhx  = hy/hx;
349:   tleft = user->tleft;         tright = user->tright;
350:   beta  = user->beta;               bm1    = user->bm1;                coef = user->coef;

352:   /* Get ghost points */
353:   DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
354:   DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
355:   DAGetCorners((DA)dmmg->dm,&xs,&ys,0,&xm,&ym,0);
356:   DAVecGetArray((DA)dmmg->dm,localX,&x);

358:   /* Evaluate Jacobian of function */
359:   for (j=ys; j<ys+ym; j++) {
360:     for (i=xs; i<xs+xm; i++) {
361:       t0 = x[j][i];

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

365:         /* general interior volume */

367:         tw = x[j][i-1];
368:         aw = 0.5*(t0 + tw);
369:         bw = PetscPowScalar(aw,bm1);
370:         /* dw = bw * aw */
371:         dw = PetscPowScalar(aw,beta);
372:         gw = coef*bw*(t0 - tw);

374:         te = x[j][i+1];
375:         ae = 0.5*(t0 + te);
376:         be = PetscPowScalar(ae,bm1);
377:         /* de = be * ae; */
378:         de = PetscPowScalar(ae,beta);
379:         ge = coef*be*(te - t0);

381:         ts = x[j-1][i];
382:         as = 0.5*(t0 + ts);
383:         bs = PetscPowScalar(as,bm1);
384:         /* ds = bs * as; */
385:         ds = PetscPowScalar(as,beta);
386:         gs = coef*bs*(t0 - ts);
387: 
388:         tn = x[j+1][i];
389:         an = 0.5*(t0 + tn);
390:         bn = PetscPowScalar(an,bm1);
391:         /* dn = bn * an; */
392:         dn = PetscPowScalar(an,beta);
393:         gn = coef*bn*(tn - t0);

395:         v[0] = - hxdhy*(ds - gs);                                      col[0].j = j-1;       col[0].i = i;
396:         v[1] = - hydhx*(dw - gw);                                      col[1].j = j;         col[1].i = i-1;
397:         v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
398:         v[3] = - hydhx*(de + ge);                                      col[3].j = j;         col[3].i = i+1;
399:         v[4] = - hxdhy*(dn + gn);                                      col[4].j = j+1;       col[4].i = i;
400:         MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);

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

404:         /* left-hand boundary */
405:         tw = tleft;
406:         aw = 0.5*(t0 + tw);
407:         bw = PetscPowScalar(aw,bm1);
408:         /* dw = bw * aw */
409:         dw = PetscPowScalar(aw,beta);
410:         gw = coef*bw*(t0 - tw);
411: 
412:         te = x[j][i + 1];
413:         ae = 0.5*(t0 + te);
414:         be = PetscPowScalar(ae,bm1);
415:         /* de = be * ae; */
416:         de = PetscPowScalar(ae,beta);
417:         ge = coef*be*(te - t0);
418: 
419:         /* left-hand bottom boundary */
420:         if (j == 0) {

422:           tn = x[j+1][i];
423:           an = 0.5*(t0 + tn);
424:           bn = PetscPowScalar(an,bm1);
425:           /* dn = bn * an; */
426:           dn = PetscPowScalar(an,beta);
427:           gn = coef*bn*(tn - t0);
428: 
429:           v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i;
430:           v[1] = - hydhx*(de + ge);                           col[1].j = j;         col[1].i = i+1;
431:           v[2] = - hxdhy*(dn + gn);                           col[2].j = j+1;       col[2].i = i;
432:           MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
433: 
434:         /* left-hand interior boundary */
435:         } else if (j < my-1) {

437:           ts = x[j-1][i];
438:           as = 0.5*(t0 + ts);
439:           bs = PetscPowScalar(as,bm1);
440:           /* ds = bs * as; */
441:           ds = PetscPowScalar(as,beta);
442:           gs = coef*bs*(t0 - ts);
443: 
444:           tn = x[j+1][i];
445:           an = 0.5*(t0 + tn);
446:           bn = PetscPowScalar(an,bm1);
447:           /* dn = bn * an; */
448:           dn = PetscPowScalar(an,beta);
449:           gn = coef*bn*(tn - t0);
450: 
451:           v[0] = - hxdhy*(ds - gs);                                      col[0].j = j-1;       col[0].i = i;
452:           v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[1].j = row.j = j; col[1].i = row.i = i;
453:           v[2] = - hydhx*(de + ge);                                      col[2].j = j;         col[2].i = i+1;
454:           v[3] = - hxdhy*(dn + gn);                                      col[3].j = j+1;       col[3].i = i;
455:           MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
456:         /* left-hand top boundary */
457:         } else {

459:           ts = x[j-1][i];
460:           as = 0.5*(t0 + ts);
461:           bs = PetscPowScalar(as,bm1);
462:           /* ds = bs * as; */
463:           ds = PetscPowScalar(as,beta);
464:           gs = coef*bs*(t0 - ts);
465: 
466:           v[0] = - hxdhy*(ds - gs);                            col[0].j = j-1;       col[0].i = i;
467:           v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[1].j = row.j = j; col[1].i = row.i = i;
468:           v[2] = - hydhx*(de + ge);                            col[2].j = j;         col[2].i = i+1;
469:           MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
470:         }

472:       } else if (i == mx-1) {
473: 
474:         /* right-hand boundary */
475:         tw = x[j][i-1];
476:         aw = 0.5*(t0 + tw);
477:         bw = PetscPowScalar(aw,bm1);
478:         /* dw = bw * aw */
479:         dw = PetscPowScalar(aw,beta);
480:         gw = coef*bw*(t0 - tw);
481: 
482:         te = tright;
483:         ae = 0.5*(t0 + te);
484:         be = PetscPowScalar(ae,bm1);
485:         /* de = be * ae; */
486:         de = PetscPowScalar(ae,beta);
487:         ge = coef*be*(te - t0);
488: 
489:         /* right-hand bottom boundary */
490:         if (j == 0) {

492:           tn = x[j+1][i];
493:           an = 0.5*(t0 + tn);
494:           bn = PetscPowScalar(an,bm1);
495:           /* dn = bn * an; */
496:           dn = PetscPowScalar(an,beta);
497:           gn = coef*bn*(tn - t0);
498: 
499:           v[0] = - hydhx*(dw - gw);                           col[0].j = j;         col[0].i = i-1;
500:           v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
501:           v[2] = - hxdhy*(dn + gn);                           col[2].j = j+1;       col[2].i = i;
502:           MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
503: 
504:         /* right-hand interior boundary */
505:         } else if (j < my-1) {

507:           ts = x[j-1][i];
508:           as = 0.5*(t0 + ts);
509:           bs = PetscPowScalar(as,bm1);
510:           /* ds = bs * as; */
511:           ds = PetscPowScalar(as,beta);
512:           gs = coef*bs*(t0 - ts);
513: 
514:           tn = x[j+1][i];
515:           an = 0.5*(t0 + tn);
516:           bn = PetscPowScalar(an,bm1);
517:           /* dn = bn * an; */
518:           dn = PetscPowScalar(an,beta);
519:           gn = coef*bn*(tn - t0);
520: 
521:           v[0] = - hxdhy*(ds - gs);                                      col[0].j = j-1;       col[0].i = i;
522:           v[1] = - hydhx*(dw - gw);                                      col[1].j = j;         col[1].i = i-1;
523:           v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
524:           v[3] = - hxdhy*(dn + gn);                                      col[3].j = j+1;       col[3].i = i;
525:           MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
526:         /* right-hand top boundary */
527:         } else {

529:           ts = x[j-1][i];
530:           as = 0.5*(t0 + ts);
531:           bs = PetscPowScalar(as,bm1);
532:           /* ds = bs * as; */
533:           ds = PetscPowScalar(as,beta);
534:           gs = coef*bs*(t0 - ts);
535: 
536:           v[0] = - hxdhy*(ds - gs);                            col[0].j = j-1;       col[0].i = i;
537:           v[1] = - hydhx*(dw - gw);                            col[1].j = j;         col[1].i = i-1;
538:           v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
539:           MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
540:         }

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

545:         tw = x[j][i-1];
546:         aw = 0.5*(t0 + tw);
547:         bw = PetscPowScalar(aw,bm1);
548:         /* dw = bw * aw */
549:         dw = PetscPowScalar(aw,beta);
550:         gw = coef*bw*(t0 - tw);

552:         te = x[j][i+1];
553:         ae = 0.5*(t0 + te);
554:         be = PetscPowScalar(ae,bm1);
555:         /* de = be * ae; */
556:         de = PetscPowScalar(ae,beta);
557:         ge = coef*be*(te - t0);

559:         tn = x[j+1][i];
560:         an = 0.5*(t0 + tn);
561:         bn = PetscPowScalar(an,bm1);
562:         /* dn = bn * an; */
563:         dn = PetscPowScalar(an,beta);
564:         gn = coef*bn*(tn - t0);
565: 
566:         v[0] = - hydhx*(dw - gw);                           col[0].j = j;         col[0].i = i-1;
567:         v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
568:         v[2] = - hydhx*(de + ge);                           col[2].j = j;         col[2].i = i+1;
569:         v[3] = - hxdhy*(dn + gn);                           col[3].j = j+1;       col[3].i = i;
570:         MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
571: 
572:       /* top boundary,and i <> 0 or mx-1 */
573:       } else if (j == my-1) {
574: 
575:         tw = x[j][i-1];
576:         aw = 0.5*(t0 + tw);
577:         bw = PetscPowScalar(aw,bm1);
578:         /* dw = bw * aw */
579:         dw = PetscPowScalar(aw,beta);
580:         gw = coef*bw*(t0 - tw);

582:         te = x[j][i+1];
583:         ae = 0.5*(t0 + te);
584:         be = PetscPowScalar(ae,bm1);
585:         /* de = be * ae; */
586:         de = PetscPowScalar(ae,beta);
587:         ge = coef*be*(te - t0);

589:         ts = x[j-1][i];
590:         as = 0.5*(t0 + ts);
591:         bs = PetscPowScalar(as,bm1);
592:          /* ds = bs * as; */
593:         ds = PetscPowScalar(as,beta);
594:         gs = coef*bs*(t0 - ts);

596:         v[0] = - hxdhy*(ds - gs);                            col[0].j = j-1;       col[0].i = i;
597:         v[1] = - hydhx*(dw - gw);                            col[1].j = j;         col[1].i = i-1;
598:         v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
599:         v[3] = - hydhx*(de + ge);                            col[3].j = j;         col[3].i = i+1;
600:         MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
601: 
602:       }
603:     }
604:   }
605:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
606:   DAVecRestoreArray((DA)dmmg->dm,localX,&x);
607:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
608:   DARestoreLocalVector((DA)dmmg->dm,&localX);

610:   PetscLogFlops((41 + 8*POWFLOP)*xm*ym);
611:   return(0);
612: }