Actual source code: ex18.c

petsc-3.13.6 2020-09-29
Report Typos and Errors

  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 Section 1.5 Writing Application Codes with PETSc 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*/