Actual source code: ex18.c

petsc-3.3-p7 2013-05-11
```  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*/

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 <petscdmda.h>
44: #include <petscpcmg.h>

46: /* User-defined application context */

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

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

55: extern PetscErrorCode FormInitialGuess(DM,Vec);
56: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
57: extern PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,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,PETSC_NULL,help);

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

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

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

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

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

111:   DMDestroy(&da);
112:   SNESDestroy(&snes);
113:   PetscFinalize();

115:   return 0;
116: }
117: /* --------------------  Form initial approximation ----------------- */
120: PetscErrorCode FormInitialGuess(DM da,Vec X)
121: {
122:   AppCtx         *user;
123:   PetscInt       i,j,xs,ys,xm,ym;
125:   PetscReal      tleft;
126:   PetscScalar    **x;

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) --------------------- */
147: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ptr)
148: {
149:   AppCtx         *user = (AppCtx*)ptr;
151:   PetscInt       i,j,mx,my,xs,ys,xm,ym;
152:   PetscScalar    zero = 0.0,one = 1.0;
153:   PetscScalar    hx,hy,hxdhy,hydhx;
154:   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;
155:   PetscScalar    tleft,tright,beta;
156:   PetscScalar    **x,**f;
157:   Vec            localX;
158:   DM             da;

161:   SNESGetDM(snes,&da);
162:   DMGetLocalVector(da,&localX);
163:   DMDAGetInfo(da,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
164:   hx    = one/(PetscReal)(mx-1);  hy    = one/(PetscReal)(my-1);
165:   hxdhy = hx/hy;               hydhx = hy/hx;
166:   tleft = user->tleft;         tright = user->tright;
167:   beta  = user->beta;
168:
169:   /* Get ghost points */
170:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
171:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
172:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
173:   DMDAVecGetArray(da,localX,&x);
174:   DMDAVecGetArray(da,F,&f);

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

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

183:         /* general interior volume */

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

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

195:         ts = x[j-1][i];
196:         as = 0.5*(t0 + ts);
197:         ds = PetscPowScalar(as,beta);
198:         fs = ds*(t0 - ts);
199:
200:         tn = x[j+1][i];
201:         an = 0.5*(t0 + tn);
202:         dn = PetscPowScalar(an,beta);
203:         fn = dn*(tn - t0);

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

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

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

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

227:         if (j < my-1) {
228:           tn = x[j+1][i];
229:           an = 0.5*(t0 + tn);
230:           dn = PetscPowScalar(an,beta);
231:           fn = dn*(tn - t0);
232:         } else {
233:           fn = zero;
234:         }

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

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

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

269:         /* bottom boundary,and i <> 0 or mx-1 */
270:         tw = x[j][i-1];
271:         aw = 0.5*(t0 + tw);
272:         dw = PetscPowScalar(aw,beta);
273:         fw = dw*(t0 - tw);

275:         te = x[j][i+1];
276:         ae = 0.5*(t0 + te);
277:         de = PetscPowScalar(ae,beta);
278:         fe = de*(te - t0);

280:         fs = zero;

282:         tn = x[j+1][i];
283:         an = 0.5*(t0 + tn);
284:         dn = PetscPowScalar(an,beta);
285:         fn = dn*(tn - t0);

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

289:         /* top boundary,and i <> 0 or mx-1 */
290:         tw = x[j][i-1];
291:         aw = 0.5*(t0 + tw);
292:         dw = PetscPowScalar(aw,beta);
293:         fw = dw*(t0 - tw);

295:         te = x[j][i+1];
296:         ae = 0.5*(t0 + te);
297:         de = PetscPowScalar(ae,beta);
298:         fe = de*(te - t0);

300:         ts = x[j-1][i];
301:         as = 0.5*(t0 + ts);
302:         ds = PetscPowScalar(as,beta);
303:         fs = ds*(t0 - ts);

305:         fn = zero;

307:       }

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

311:     }
312:   }
313:   DMDAVecRestoreArray(da,localX,&x);
314:   DMDAVecRestoreArray(da,F,&f);
315:   DMRestoreLocalVector(da,&localX);
316:   PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
317:   return(0);
318: }
319: /* --------------------  Evaluate Jacobian F(x) --------------------- */
322: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ptr)
323: {
324:   AppCtx         *user = (AppCtx*)ptr;
325:   Mat            jac = *J;
327:   PetscInt       i,j,mx,my,xs,ys,xm,ym;
328:   PetscScalar    one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;
329:   PetscScalar    dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;
330:   PetscScalar    tleft,tright,beta,bm1,coef;
331:   PetscScalar    v[5],**x;
332:   Vec            localX;
333:   MatStencil     col[5],row;
334:   DM             da;

337:   SNESGetDM(snes,&da);
338:   DMGetLocalVector(da,&localX);
339:   *flg = SAME_NONZERO_PATTERN;
340:   DMDAGetInfo(da,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
341:   hx    = one/(PetscReal)(mx-1);  hy     = one/(PetscReal)(my-1);
342:   hxdhy = hx/hy;               hydhx  = hy/hx;
343:   tleft = user->tleft;         tright = user->tright;
344:   beta  = user->beta;               bm1    = user->bm1;                coef = user->coef;

346:   /* Get ghost points */
347:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
348:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
349:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
350:   DMDAVecGetArray(da,localX,&x);

352:   /* Evaluate Jacobian of function */
353:   for (j=ys; j<ys+ym; j++) {
354:     for (i=xs; i<xs+xm; i++) {
355:       t0 = x[j][i];

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

359:         /* general interior volume */

361:         tw = x[j][i-1];
362:         aw = 0.5*(t0 + tw);
363:         bw = PetscPowScalar(aw,bm1);
364:         /* dw = bw * aw */
365:         dw = PetscPowScalar(aw,beta);
366:         gw = coef*bw*(t0 - tw);

368:         te = x[j][i+1];
369:         ae = 0.5*(t0 + te);
370:         be = PetscPowScalar(ae,bm1);
371:         /* de = be * ae; */
372:         de = PetscPowScalar(ae,beta);
373:         ge = coef*be*(te - t0);

375:         ts = x[j-1][i];
376:         as = 0.5*(t0 + ts);
377:         bs = PetscPowScalar(as,bm1);
378:         /* ds = bs * as; */
379:         ds = PetscPowScalar(as,beta);
380:         gs = coef*bs*(t0 - ts);
381:
382:         tn = x[j+1][i];
383:         an = 0.5*(t0 + tn);
384:         bn = PetscPowScalar(an,bm1);
385:         /* dn = bn * an; */
386:         dn = PetscPowScalar(an,beta);
387:         gn = coef*bn*(tn - t0);

389:         v[0] = - hxdhy*(ds - gs);                                      col[0].j = j-1;       col[0].i = i;
390:         v[1] = - hydhx*(dw - gw);                                      col[1].j = j;         col[1].i = i-1;
391:         v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
392:         v[3] = - hydhx*(de + ge);                                      col[3].j = j;         col[3].i = i+1;
393:         v[4] = - hxdhy*(dn + gn);                                      col[4].j = j+1;       col[4].i = i;
394:         MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);

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

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

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

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

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

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

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

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

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

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

539:         tw = x[j][i-1];
540:         aw = 0.5*(t0 + tw);
541:         bw = PetscPowScalar(aw,bm1);
542:         /* dw = bw * aw */
543:         dw = PetscPowScalar(aw,beta);
544:         gw = coef*bw*(t0 - tw);

546:         te = x[j][i+1];
547:         ae = 0.5*(t0 + te);
548:         be = PetscPowScalar(ae,bm1);
549:         /* de = be * ae; */
550:         de = PetscPowScalar(ae,beta);
551:         ge = coef*be*(te - t0);

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

576:         te = x[j][i+1];
577:         ae = 0.5*(t0 + te);
578:         be = PetscPowScalar(ae,bm1);
579:         /* de = be * ae; */
580:         de = PetscPowScalar(ae,beta);
581:         ge = coef*be*(te - t0);

583:         ts = x[j-1][i];
584:         as = 0.5*(t0 + ts);
585:         bs = PetscPowScalar(as,bm1);
586:          /* ds = bs * as; */
587:         ds = PetscPowScalar(as,beta);
588:         gs = coef*bs*(t0 - ts);

590:         v[0] = - hxdhy*(ds - gs);                            col[0].j = j-1;       col[0].i = i;
591:         v[1] = - hydhx*(dw - gw);                            col[1].j = j;         col[1].i = i-1;
592:         v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
593:         v[3] = - hydhx*(de + ge);                            col[3].j = j;         col[3].i = i+1;
594:         MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
595:
596:       }
597:     }
598:   }
599:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
600:   DMDAVecRestoreArray(da,localX,&x);
601:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
602:   DMRestoreLocalVector(da,&localX);

604:   PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
605:   return(0);
606: }

```