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: /*

16:     This example models the partial differential equation

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

20:     where beta = 2.5 and alpha = 1.0

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

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

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

31:     This code was contributed by David Keyes

33: */

35: #include <petscsnes.h>
36: #include <petscdm.h>
37: #include <petscdmda.h>

39: /* User-defined application context */

41: typedef struct {
42:   PetscReal tleft, tright;   /* Dirichlet boundary conditions */
43:   PetscReal beta, bm1, coef; /* nonlinear diffusivity parameterizations */
44: } AppCtx;

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

48: extern PetscErrorCode FormInitialGuess(SNES, Vec, void *);
49: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
50: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);

52: int main(int argc, char **argv)
53: {
54:   SNES      snes;
55:   AppCtx    user;
56:   PetscInt  its, lits;
57:   PetscReal litspit;
58:   DM        da;

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

63:   /* set problem parameters */
64:   user.tleft  = 1.0;
65:   user.tright = 0.1;
66:   user.beta   = 2.5;
67:   PetscOptionsGetReal(NULL, NULL, "-tleft", &user.tleft, NULL);
68:   PetscOptionsGetReal(NULL, NULL, "-tright", &user.tright, NULL);
69:   PetscOptionsGetReal(NULL, NULL, "-beta", &user.beta, NULL);
70:   user.bm1  = user.beta - 1.0;
71:   user.coef = user.beta / 2.0;

73:   /*
74:       Create the multilevel DM data structure
75:   */
76:   SNESCreate(PETSC_COMM_WORLD, &snes);

78:   /*
79:       Set the DMDA (grid structure) for the grids.
80:   */
81:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 5, 5, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, &da);
82:   DMSetFromOptions(da);
83:   DMSetUp(da);
84:   DMSetApplicationContext(da, &user);
85:   SNESSetDM(snes, (DM)da);

87:   /*
88:      Create the nonlinear solver, and tell it the functions to use
89:   */
90:   SNESSetFunction(snes, NULL, FormFunction, &user);
91:   SNESSetJacobian(snes, NULL, NULL, FormJacobian, &user);
92:   SNESSetFromOptions(snes);
93:   SNESSetComputeInitialGuess(snes, FormInitialGuess, NULL);

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

103:   DMDestroy(&da);
104:   SNESDestroy(&snes);
105:   PetscFinalize();
106:   return 0;
107: }
108: /* --------------------  Form initial approximation ----------------- */
109: PetscErrorCode FormInitialGuess(SNES snes, Vec X, void *ctx)
110: {
111:   AppCtx       *user;
112:   PetscInt      i, j, xs, ys, xm, ym;
113:   PetscReal     tleft;
114:   PetscScalar **x;
115:   DM            da;

118:   SNESGetDM(snes, &da);
119:   DMGetApplicationContext(da, &user);
120:   tleft = user->tleft;
121:   /* Get ghost points */
122:   DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0);
123:   DMDAVecGetArray(da, X, &x);

125:   /* Compute initial guess */
126:   for (j = ys; j < ys + ym; j++) {
127:     for (i = xs; i < xs + xm; i++) x[j][i] = tleft;
128:   }
129:   DMDAVecRestoreArray(da, X, &x);
130:   return 0;
131: }
132: /* --------------------  Evaluate Function F(x) --------------------- */
133: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr)
134: {
135:   AppCtx       *user = (AppCtx *)ptr;
136:   PetscInt      i, j, mx, my, xs, ys, xm, ym;
137:   PetscScalar   zero = 0.0, one = 1.0;
138:   PetscScalar   hx, hy, hxdhy, hydhx;
139:   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;
140:   PetscScalar   tleft, tright, beta;
141:   PetscScalar **x, **f;
142:   Vec           localX;
143:   DM            da;

146:   SNESGetDM(snes, &da);
147:   DMGetLocalVector(da, &localX);
148:   DMDAGetInfo(da, NULL, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
149:   hx     = one / (PetscReal)(mx - 1);
150:   hy     = one / (PetscReal)(my - 1);
151:   hxdhy  = hx / hy;
152:   hydhx  = hy / hx;
153:   tleft  = user->tleft;
154:   tright = user->tright;
155:   beta   = user->beta;

157:   /* Get ghost points */
158:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
159:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
160:   DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0);
161:   DMDAVecGetArray(da, localX, &x);
162:   DMDAVecGetArray(da, F, &f);

164:   /* Evaluate function */
165:   for (j = ys; j < ys + ym; j++) {
166:     for (i = xs; i < xs + xm; i++) {
167:       t0 = x[j][i];

169:       if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) {
170:         /* general interior volume */

172:         tw = x[j][i - 1];
173:         aw = 0.5 * (t0 + tw);
174:         dw = PetscPowScalar(aw, beta);
175:         fw = dw * (t0 - tw);

177:         te = x[j][i + 1];
178:         ae = 0.5 * (t0 + te);
179:         de = PetscPowScalar(ae, beta);
180:         fe = de * (te - t0);

182:         ts = x[j - 1][i];
183:         as = 0.5 * (t0 + ts);
184:         ds = PetscPowScalar(as, beta);
185:         fs = ds * (t0 - ts);

187:         tn = x[j + 1][i];
188:         an = 0.5 * (t0 + tn);
189:         dn = PetscPowScalar(an, beta);
190:         fn = dn * (tn - t0);

192:       } else if (i == 0) {
193:         /* left-hand boundary */
194:         tw = tleft;
195:         aw = 0.5 * (t0 + tw);
196:         dw = PetscPowScalar(aw, beta);
197:         fw = dw * (t0 - tw);

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

204:         if (j > 0) {
205:           ts = x[j - 1][i];
206:           as = 0.5 * (t0 + ts);
207:           ds = PetscPowScalar(as, beta);
208:           fs = ds * (t0 - ts);
209:         } else fs = zero;

211:         if (j < my - 1) {
212:           tn = x[j + 1][i];
213:           an = 0.5 * (t0 + tn);
214:           dn = PetscPowScalar(an, beta);
215:           fn = dn * (tn - t0);
216:         } else fn = zero;

218:       } else if (i == mx - 1) {
219:         /* right-hand boundary */
220:         tw = x[j][i - 1];
221:         aw = 0.5 * (t0 + tw);
222:         dw = PetscPowScalar(aw, beta);
223:         fw = dw * (t0 - tw);

225:         te = tright;
226:         ae = 0.5 * (t0 + te);
227:         de = PetscPowScalar(ae, beta);
228:         fe = de * (te - t0);

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

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

244:       } else if (j == 0) {
245:         /* bottom boundary,and i <> 0 or mx-1 */
246:         tw = x[j][i - 1];
247:         aw = 0.5 * (t0 + tw);
248:         dw = PetscPowScalar(aw, beta);
249:         fw = dw * (t0 - tw);

251:         te = x[j][i + 1];
252:         ae = 0.5 * (t0 + te);
253:         de = PetscPowScalar(ae, beta);
254:         fe = de * (te - t0);

256:         fs = zero;

258:         tn = x[j + 1][i];
259:         an = 0.5 * (t0 + tn);
260:         dn = PetscPowScalar(an, beta);
261:         fn = dn * (tn - t0);

263:       } else if (j == my - 1) {
264:         /* top boundary,and i <> 0 or mx-1 */
265:         tw = x[j][i - 1];
266:         aw = 0.5 * (t0 + tw);
267:         dw = PetscPowScalar(aw, beta);
268:         fw = dw * (t0 - tw);

270:         te = x[j][i + 1];
271:         ae = 0.5 * (t0 + te);
272:         de = PetscPowScalar(ae, beta);
273:         fe = de * (te - t0);

275:         ts = x[j - 1][i];
276:         as = 0.5 * (t0 + ts);
277:         ds = PetscPowScalar(as, beta);
278:         fs = ds * (t0 - ts);

280:         fn = zero;
281:       }

283:       f[j][i] = -hydhx * (fe - fw) - hxdhy * (fn - fs);
284:     }
285:   }
286:   DMDAVecRestoreArray(da, localX, &x);
287:   DMDAVecRestoreArray(da, F, &f);
288:   DMRestoreLocalVector(da, &localX);
289:   PetscLogFlops((22.0 + 4.0 * POWFLOP) * ym * xm);
290:   return 0;
291: }
292: /* --------------------  Evaluate Jacobian F(x) --------------------- */
293: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat jac, Mat B, void *ptr)
294: {
295:   AppCtx     *user = (AppCtx *)ptr;
296:   PetscInt    i, j, mx, my, xs, ys, xm, ym;
297:   PetscScalar one = 1.0, hx, hy, hxdhy, hydhx, t0, tn, ts, te, tw;
298:   PetscScalar dn, ds, de, dw, an, as, ae, aw, bn, bs, be, bw, gn, gs, ge, gw;
299:   PetscScalar tleft, tright, beta, bm1, coef;
300:   PetscScalar v[5], **x;
301:   Vec         localX;
302:   MatStencil  col[5], row;
303:   DM          da;

306:   SNESGetDM(snes, &da);
307:   DMGetLocalVector(da, &localX);
308:   DMDAGetInfo(da, NULL, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
309:   hx     = one / (PetscReal)(mx - 1);
310:   hy     = one / (PetscReal)(my - 1);
311:   hxdhy  = hx / hy;
312:   hydhx  = hy / hx;
313:   tleft  = user->tleft;
314:   tright = user->tright;
315:   beta   = user->beta;
316:   bm1    = user->bm1;
317:   coef   = user->coef;

319:   /* Get ghost points */
320:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
321:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
322:   DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0);
323:   DMDAVecGetArray(da, localX, &x);

325:   /* Evaluate Jacobian of function */
326:   for (j = ys; j < ys + ym; j++) {
327:     for (i = xs; i < xs + xm; i++) {
328:       t0 = x[j][i];

330:       if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) {
331:         /* general interior volume */

333:         tw = x[j][i - 1];
334:         aw = 0.5 * (t0 + tw);
335:         bw = PetscPowScalar(aw, bm1);
336:         /* dw = bw * aw */
337:         dw = PetscPowScalar(aw, beta);
338:         gw = coef * bw * (t0 - tw);

340:         te = x[j][i + 1];
341:         ae = 0.5 * (t0 + te);
342:         be = PetscPowScalar(ae, bm1);
343:         /* de = be * ae; */
344:         de = PetscPowScalar(ae, beta);
345:         ge = coef * be * (te - t0);

347:         ts = x[j - 1][i];
348:         as = 0.5 * (t0 + ts);
349:         bs = PetscPowScalar(as, bm1);
350:         /* ds = bs * as; */
351:         ds = PetscPowScalar(as, beta);
352:         gs = coef * bs * (t0 - ts);

354:         tn = x[j + 1][i];
355:         an = 0.5 * (t0 + tn);
356:         bn = PetscPowScalar(an, bm1);
357:         /* dn = bn * an; */
358:         dn = PetscPowScalar(an, beta);
359:         gn = coef * bn * (tn - t0);

361:         v[0]     = -hxdhy * (ds - gs);
362:         col[0].j = j - 1;
363:         col[0].i = i;
364:         v[1]     = -hydhx * (dw - gw);
365:         col[1].j = j;
366:         col[1].i = i - 1;
367:         v[2]     = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge);
368:         col[2].j = row.j = j;
369:         col[2].i = row.i = i;
370:         v[3]             = -hydhx * (de + ge);
371:         col[3].j         = j;
372:         col[3].i         = i + 1;
373:         v[4]             = -hxdhy * (dn + gn);
374:         col[4].j         = j + 1;
375:         col[4].i         = i;
376:         MatSetValuesStencil(B, 1, &row, 5, col, v, INSERT_VALUES);

378:       } else if (i == 0) {
379:         /* left-hand boundary */
380:         tw = tleft;
381:         aw = 0.5 * (t0 + tw);
382:         bw = PetscPowScalar(aw, bm1);
383:         /* dw = bw * aw */
384:         dw = PetscPowScalar(aw, beta);
385:         gw = coef * bw * (t0 - tw);

387:         te = x[j][i + 1];
388:         ae = 0.5 * (t0 + te);
389:         be = PetscPowScalar(ae, bm1);
390:         /* de = be * ae; */
391:         de = PetscPowScalar(ae, beta);
392:         ge = coef * be * (te - t0);

394:         /* left-hand bottom boundary */
395:         if (j == 0) {
396:           tn = x[j + 1][i];
397:           an = 0.5 * (t0 + tn);
398:           bn = PetscPowScalar(an, bm1);
399:           /* dn = bn * an; */
400:           dn = PetscPowScalar(an, beta);
401:           gn = coef * bn * (tn - t0);

403:           v[0]     = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge);
404:           col[0].j = row.j = j;
405:           col[0].i = row.i = i;
406:           v[1]             = -hydhx * (de + ge);
407:           col[1].j         = j;
408:           col[1].i         = i + 1;
409:           v[2]             = -hxdhy * (dn + gn);
410:           col[2].j         = j + 1;
411:           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) {
416:           ts = x[j - 1][i];
417:           as = 0.5 * (t0 + ts);
418:           bs = PetscPowScalar(as, bm1);
419:           /* ds = bs * as; */
420:           ds = PetscPowScalar(as, beta);
421:           gs = coef * bs * (t0 - ts);

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

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

452:           v[0]     = -hxdhy * (ds - gs);
453:           col[0].j = j - 1;
454:           col[0].i = i;
455:           v[1]     = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge);
456:           col[1].j = row.j = j;
457:           col[1].i = row.i = i;
458:           v[2]             = -hydhx * (de + ge);
459:           col[2].j         = j;
460:           col[2].i         = i + 1;
461:           MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES);
462:         }

464:       } else if (i == mx - 1) {
465:         /* right-hand boundary */
466:         tw = x[j][i - 1];
467:         aw = 0.5 * (t0 + tw);
468:         bw = PetscPowScalar(aw, bm1);
469:         /* dw = bw * aw */
470:         dw = PetscPowScalar(aw, beta);
471:         gw = coef * bw * (t0 - tw);

473:         te = tright;
474:         ae = 0.5 * (t0 + te);
475:         be = PetscPowScalar(ae, bm1);
476:         /* de = be * ae; */
477:         de = PetscPowScalar(ae, beta);
478:         ge = coef * be * (te - t0);

480:         /* right-hand bottom boundary */
481:         if (j == 0) {
482:           tn = x[j + 1][i];
483:           an = 0.5 * (t0 + tn);
484:           bn = PetscPowScalar(an, bm1);
485:           /* dn = bn * an; */
486:           dn = PetscPowScalar(an, beta);
487:           gn = coef * bn * (tn - t0);

489:           v[0]     = -hydhx * (dw - gw);
490:           col[0].j = j;
491:           col[0].i = i - 1;
492:           v[1]     = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge);
493:           col[1].j = row.j = j;
494:           col[1].i = row.i = i;
495:           v[2]             = -hxdhy * (dn + gn);
496:           col[2].j         = j + 1;
497:           col[2].i         = i;
498:           MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES);

500:           /* right-hand interior boundary */
501:         } else if (j < my - 1) {
502:           ts = x[j - 1][i];
503:           as = 0.5 * (t0 + ts);
504:           bs = PetscPowScalar(as, bm1);
505:           /* ds = bs * as; */
506:           ds = PetscPowScalar(as, beta);
507:           gs = coef * bs * (t0 - ts);

509:           tn = x[j + 1][i];
510:           an = 0.5 * (t0 + tn);
511:           bn = PetscPowScalar(an, bm1);
512:           /* dn = bn * an; */
513:           dn = PetscPowScalar(an, beta);
514:           gn = coef * bn * (tn - t0);

516:           v[0]     = -hxdhy * (ds - gs);
517:           col[0].j = j - 1;
518:           col[0].i = i;
519:           v[1]     = -hydhx * (dw - gw);
520:           col[1].j = j;
521:           col[1].i = i - 1;
522:           v[2]     = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge);
523:           col[2].j = row.j = j;
524:           col[2].i = row.i = i;
525:           v[3]             = -hxdhy * (dn + gn);
526:           col[3].j         = j + 1;
527:           col[3].i         = i;
528:           MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES);
529:           /* right-hand top boundary */
530:         } else {
531:           ts = x[j - 1][i];
532:           as = 0.5 * (t0 + ts);
533:           bs = PetscPowScalar(as, bm1);
534:           /* ds = bs * as; */
535:           ds = PetscPowScalar(as, beta);
536:           gs = coef * bs * (t0 - ts);

538:           v[0]     = -hxdhy * (ds - gs);
539:           col[0].j = j - 1;
540:           col[0].i = i;
541:           v[1]     = -hydhx * (dw - gw);
542:           col[1].j = j;
543:           col[1].i = i - 1;
544:           v[2]     = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge);
545:           col[2].j = row.j = j;
546:           col[2].i = row.i = i;
547:           MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES);
548:         }

550:         /* bottom boundary,and i <> 0 or mx-1 */
551:       } else if (j == 0) {
552:         tw = x[j][i - 1];
553:         aw = 0.5 * (t0 + tw);
554:         bw = PetscPowScalar(aw, bm1);
555:         /* dw = bw * aw */
556:         dw = PetscPowScalar(aw, beta);
557:         gw = coef * bw * (t0 - tw);

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

566:         tn = x[j + 1][i];
567:         an = 0.5 * (t0 + tn);
568:         bn = PetscPowScalar(an, bm1);
569:         /* dn = bn * an; */
570:         dn = PetscPowScalar(an, beta);
571:         gn = coef * bn * (tn - t0);

573:         v[0]     = -hydhx * (dw - gw);
574:         col[0].j = j;
575:         col[0].i = i - 1;
576:         v[1]     = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge);
577:         col[1].j = row.j = j;
578:         col[1].i = row.i = i;
579:         v[2]             = -hydhx * (de + ge);
580:         col[2].j         = j;
581:         col[2].i         = i + 1;
582:         v[3]             = -hxdhy * (dn + gn);
583:         col[3].j         = j + 1;
584:         col[3].i         = i;
585:         MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES);

587:         /* top boundary,and i <> 0 or mx-1 */
588:       } else if (j == my - 1) {
589:         tw = x[j][i - 1];
590:         aw = 0.5 * (t0 + tw);
591:         bw = PetscPowScalar(aw, bm1);
592:         /* dw = bw * aw */
593:         dw = PetscPowScalar(aw, beta);
594:         gw = coef * bw * (t0 - tw);

596:         te = x[j][i + 1];
597:         ae = 0.5 * (t0 + te);
598:         be = PetscPowScalar(ae, bm1);
599:         /* de = be * ae; */
600:         de = PetscPowScalar(ae, beta);
601:         ge = coef * be * (te - t0);

603:         ts = x[j - 1][i];
604:         as = 0.5 * (t0 + ts);
605:         bs = PetscPowScalar(as, bm1);
606:         /* ds = bs * as; */
607:         ds = PetscPowScalar(as, beta);
608:         gs = coef * bs * (t0 - ts);

610:         v[0]     = -hxdhy * (ds - gs);
611:         col[0].j = j - 1;
612:         col[0].i = i;
613:         v[1]     = -hydhx * (dw - gw);
614:         col[1].j = j;
615:         col[1].i = i - 1;
616:         v[2]     = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge);
617:         col[2].j = row.j = j;
618:         col[2].i = row.i = i;
619:         v[3]             = -hydhx * (de + ge);
620:         col[3].j         = j;
621:         col[3].i         = i + 1;
622:         MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES);
623:       }
624:     }
625:   }
626:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
627:   DMDAVecRestoreArray(da, localX, &x);
628:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
629:   DMRestoreLocalVector(da, &localX);
630:   if (jac != B) {
631:     MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
632:     MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
633:   }

635:   PetscLogFlops((41.0 + 8.0 * POWFLOP) * xm * ym);
636:   return 0;
637: }

639: /*TEST

641:    test:
642:       suffix: 1
643:       args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view
644:       requires: !single

646:    test:
647:       suffix: 2
648:       args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view -snes_type newtontrdc
649:       requires: !single

651: TEST*/
```