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

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:   */
110:     DMMGSetInitialGuess(dmmg,FormInitialGuess);
111:     DMMGSolve(dmmg);
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: }

```