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: DMDA^using distributed arrays
 17:    Concepts: multigrid;
 18:    Processors: n
 19: T*/

 21: /*

 23:     This example models the partial differential equation

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

 27:     where beta = 2.5 and alpha = 1.0

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

 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.

 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

 40: */

 42: #include <petscsnes.h>
 43: #include <petscdm.h>
 44: #include <petscdmda.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(SNES,Vec,void*);
 56: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 57: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);

 59: int main(int argc,char **argv)
 60: {
 61:   SNES           snes;
 62:   AppCtx         user;
 63:   PetscInt       its,lits;
 64:   PetscReal      litspit;
 65:   DM             da;

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

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

 79:   /*
 80:       Create the multilevel DM data structure
 81:   */
 82:   SNESCreate(PETSC_COMM_WORLD,&snes);

 84:   /*
 85:       Set the DMDA (grid structure) for the grids.
 86:   */
 87:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,5,5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
 88:   DMSetFromOptions(da);
 89:   DMSetUp(da);
 90:   DMSetApplicationContext(da,&user);
 91:   SNESSetDM(snes,(DM)da);

 93:   /*
 94:      Create the nonlinear solver, and tell it the functions to use
 95:   */
 96:   SNESSetFunction(snes,NULL,FormFunction,&user);
 97:   SNESSetJacobian(snes,NULL,NULL,FormJacobian,&user);
 98:   SNESSetFromOptions(snes);
 99:   SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);

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

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

124:   SNESGetDM(snes,&da);
125:   DMGetApplicationContext(da,&user);
126:   tleft = user->tleft;
127:   /* Get ghost points */
128:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
129:   DMDAVecGetArray(da,X,&x);

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

154:   SNESGetDM(snes,&da);
155:   DMGetLocalVector(da,&localX);
156:   DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
157:   hx    = one/(PetscReal)(mx-1);  hy    = one/(PetscReal)(my-1);
158:   hxdhy = hx/hy;               hydhx = hy/hx;
159:   tleft = user->tleft;         tright = user->tright;
160:   beta  = user->beta;

162:   /* Get ghost points */
163:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
164:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
165:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
166:   DMDAVecGetArray(da,localX,&x);
167:   DMDAVecGetArray(da,F,&f);

169:   /* Evaluate function */
170:   for (j=ys; j<ys+ym; j++) {
171:     for (i=xs; i<xs+xm; i++) {
172:       t0 = x[j][i];

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

176:         /* general interior volume */

178:         tw = x[j][i-1];
179:         aw = 0.5*(t0 + tw);
180:         dw = PetscPowScalar(aw,beta);
181:         fw = dw*(t0 - tw);

183:         te = x[j][i+1];
184:         ae = 0.5*(t0 + te);
185:         de = PetscPowScalar(ae,beta);
186:         fe = de*(te - t0);

188:         ts = x[j-1][i];
189:         as = 0.5*(t0 + ts);
190:         ds = PetscPowScalar(as,beta);
191:         fs = ds*(t0 - ts);

193:         tn = x[j+1][i];
194:         an = 0.5*(t0 + tn);
195:         dn = PetscPowScalar(an,beta);
196:         fn = dn*(tn - t0);

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

200:         /* left-hand boundary */
201:         tw = tleft;
202:         aw = 0.5*(t0 + tw);
203:         dw = PetscPowScalar(aw,beta);
204:         fw = dw*(t0 - tw);

206:         te = x[j][i+1];
207:         ae = 0.5*(t0 + te);
208:         de = PetscPowScalar(ae,beta);
209:         fe = de*(te - t0);

211:         if (j > 0) {
212:           ts = x[j-1][i];
213:           as = 0.5*(t0 + ts);
214:           ds = PetscPowScalar(as,beta);
215:           fs = ds*(t0 - ts);
216:         } else fs = zero;

218:         if (j < my-1) {
219:           tn = x[j+1][i];
220:           an = 0.5*(t0 + tn);
221:           dn = PetscPowScalar(an,beta);
222:           fn = dn*(tn - t0);
223:         } else fn = zero;

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

227:         /* right-hand boundary */
228:         tw = x[j][i-1];
229:         aw = 0.5*(t0 + tw);
230:         dw = PetscPowScalar(aw,beta);
231:         fw = dw*(t0 - tw);

233:         te = tright;
234:         ae = 0.5*(t0 + te);
235:         de = PetscPowScalar(ae,beta);
236:         fe = de*(te - t0);

238:         if (j > 0) {
239:           ts = x[j-1][i];
240:           as = 0.5*(t0 + ts);
241:           ds = PetscPowScalar(as,beta);
242:           fs = ds*(t0 - ts);
243:         } else fs = zero;

245:         if (j < my-1) {
246:           tn = x[j+1][i];
247:           an = 0.5*(t0 + tn);
248:           dn = PetscPowScalar(an,beta);
249:           fn = dn*(tn - t0);
250:         } else fn = zero;

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

254:         /* bottom boundary,and i <> 0 or mx-1 */
255:         tw = x[j][i-1];
256:         aw = 0.5*(t0 + tw);
257:         dw = PetscPowScalar(aw,beta);
258:         fw = dw*(t0 - tw);

260:         te = x[j][i+1];
261:         ae = 0.5*(t0 + te);
262:         de = PetscPowScalar(ae,beta);
263:         fe = de*(te - t0);

265:         fs = zero;

267:         tn = x[j+1][i];
268:         an = 0.5*(t0 + tn);
269:         dn = PetscPowScalar(an,beta);
270:         fn = dn*(tn - t0);

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

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

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

285:         ts = x[j-1][i];
286:         as = 0.5*(t0 + ts);
287:         ds = PetscPowScalar(as,beta);
288:         fs = ds*(t0 - ts);

290:         fn = zero;

292:       }

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

296:     }
297:   }
298:   DMDAVecRestoreArray(da,localX,&x);
299:   DMDAVecRestoreArray(da,F,&f);
300:   DMRestoreLocalVector(da,&localX);
301:   PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
302:   return 0;
303: }
304: /* --------------------  Evaluate Jacobian F(x) --------------------- */
305: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat jac,Mat B,void *ptr)
306: {
307:   AppCtx         *user = (AppCtx*)ptr;
308:   PetscInt       i,j,mx,my,xs,ys,xm,ym;
309:   PetscScalar    one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;
310:   PetscScalar    dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;
311:   PetscScalar    tleft,tright,beta,bm1,coef;
312:   PetscScalar    v[5],**x;
313:   Vec            localX;
314:   MatStencil     col[5],row;
315:   DM             da;

318:   SNESGetDM(snes,&da);
319:   DMGetLocalVector(da,&localX);
320:   DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
321:   hx    = one/(PetscReal)(mx-1);  hy     = one/(PetscReal)(my-1);
322:   hxdhy = hx/hy;               hydhx  = hy/hx;
323:   tleft = user->tleft;         tright = user->tright;
324:   beta  = user->beta;          bm1    = user->bm1;          coef = user->coef;

326:   /* Get ghost points */
327:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
328:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
329:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
330:   DMDAVecGetArray(da,localX,&x);

332:   /* Evaluate Jacobian of function */
333:   for (j=ys; j<ys+ym; j++) {
334:     for (i=xs; i<xs+xm; i++) {
335:       t0 = x[j][i];

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

339:         /* general interior volume */

341:         tw = x[j][i-1];
342:         aw = 0.5*(t0 + tw);
343:         bw = PetscPowScalar(aw,bm1);
344:         /* dw = bw * aw */
345:         dw = PetscPowScalar(aw,beta);
346:         gw = coef*bw*(t0 - tw);

348:         te = x[j][i+1];
349:         ae = 0.5*(t0 + te);
350:         be = PetscPowScalar(ae,bm1);
351:         /* de = be * ae; */
352:         de = PetscPowScalar(ae,beta);
353:         ge = coef*be*(te - t0);

355:         ts = x[j-1][i];
356:         as = 0.5*(t0 + ts);
357:         bs = PetscPowScalar(as,bm1);
358:         /* ds = bs * as; */
359:         ds = PetscPowScalar(as,beta);
360:         gs = coef*bs*(t0 - ts);

362:         tn = x[j+1][i];
363:         an = 0.5*(t0 + tn);
364:         bn = PetscPowScalar(an,bm1);
365:         /* dn = bn * an; */
366:         dn = PetscPowScalar(an,beta);
367:         gn = coef*bn*(tn - t0);

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

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

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

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

393:         /* left-hand bottom boundary */
394:         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); col[0].j = row.j = j; col[0].i = row.i = i;
404:           v[1] = -hydhx*(de + ge);                            col[1].j = j;         col[1].i = i+1;
405:           v[2] = -hxdhy*(dn + gn);                            col[2].j = j+1;       col[2].i = i;
406:           MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);

408:           /* left-hand interior boundary */
409:         } else if (j < my-1) {

411:           ts = x[j-1][i];
412:           as = 0.5*(t0 + ts);
413:           bs = PetscPowScalar(as,bm1);
414:           /* ds = bs * as; */
415:           ds = PetscPowScalar(as,beta);
416:           gs = coef*bs*(t0 - ts);

418:           tn = x[j+1][i];
419:           an = 0.5*(t0 + tn);
420:           bn = PetscPowScalar(an,bm1);
421:           /* dn = bn * an; */
422:           dn = PetscPowScalar(an,beta);
423:           gn = coef*bn*(tn - t0);

425:           v[0] = -hxdhy*(ds - gs);                                       col[0].j = j-1;       col[0].i = i;
426:           v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[1].j = row.j = j; col[1].i = row.i = i;
427:           v[2] = -hydhx*(de + ge);                                       col[2].j = j;         col[2].i = i+1;
428:           v[3] = -hxdhy*(dn + gn);                                       col[3].j = j+1;       col[3].i = i;
429:           MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
430:           /* left-hand top boundary */
431:         } else {

433:           ts = x[j-1][i];
434:           as = 0.5*(t0 + ts);
435:           bs = PetscPowScalar(as,bm1);
436:           /* ds = bs * as; */
437:           ds = PetscPowScalar(as,beta);
438:           gs = coef*bs*(t0 - ts);

440:           v[0] = -hxdhy*(ds - gs);                             col[0].j = j-1;       col[0].i = i;
441:           v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[1].j = row.j = j; col[1].i = row.i = i;
442:           v[2] = -hydhx*(de + ge);                             col[2].j = j;         col[2].i = i+1;
443:           MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
444:         }

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

448:         /* right-hand boundary */
449:         tw = x[j][i-1];
450:         aw = 0.5*(t0 + tw);
451:         bw = PetscPowScalar(aw,bm1);
452:         /* dw = bw * aw */
453:         dw = PetscPowScalar(aw,beta);
454:         gw = coef*bw*(t0 - tw);

456:         te = tright;
457:         ae = 0.5*(t0 + te);
458:         be = PetscPowScalar(ae,bm1);
459:         /* de = be * ae; */
460:         de = PetscPowScalar(ae,beta);
461:         ge = coef*be*(te - t0);

463:         /* right-hand bottom boundary */
464:         if (j == 0) {

466:           tn = x[j+1][i];
467:           an = 0.5*(t0 + tn);
468:           bn = PetscPowScalar(an,bm1);
469:           /* dn = bn * an; */
470:           dn = PetscPowScalar(an,beta);
471:           gn = coef*bn*(tn - t0);

473:           v[0] = -hydhx*(dw - gw);                            col[0].j = j;         col[0].i = i-1;
474:           v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
475:           v[2] = -hxdhy*(dn + gn);                            col[2].j = j+1;       col[2].i = i;
476:           MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);

478:           /* right-hand interior boundary */
479:         } else if (j < my-1) {

481:           ts = x[j-1][i];
482:           as = 0.5*(t0 + ts);
483:           bs = PetscPowScalar(as,bm1);
484:           /* ds = bs * as; */
485:           ds = PetscPowScalar(as,beta);
486:           gs = coef*bs*(t0 - ts);

488:           tn = x[j+1][i];
489:           an = 0.5*(t0 + tn);
490:           bn = PetscPowScalar(an,bm1);
491:           /* dn = bn * an; */
492:           dn = PetscPowScalar(an,beta);
493:           gn = coef*bn*(tn - t0);

495:           v[0] = -hxdhy*(ds - gs);                                       col[0].j = j-1;       col[0].i = i;
496:           v[1] = -hydhx*(dw - gw);                                       col[1].j = j;         col[1].i = i-1;
497:           v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
498:           v[3] = -hxdhy*(dn + gn);                                       col[3].j = j+1;       col[3].i = i;
499:           MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
500:         /* right-hand top boundary */
501:         } else {

503:           ts = x[j-1][i];
504:           as = 0.5*(t0 + ts);
505:           bs = PetscPowScalar(as,bm1);
506:           /* ds = bs * as; */
507:           ds = PetscPowScalar(as,beta);
508:           gs = coef*bs*(t0 - ts);

510:           v[0] = -hxdhy*(ds - gs);                             col[0].j = j-1;       col[0].i = i;
511:           v[1] = -hydhx*(dw - gw);                             col[1].j = j;         col[1].i = i-1;
512:           v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
513:           MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
514:         }

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

519:         tw = x[j][i-1];
520:         aw = 0.5*(t0 + tw);
521:         bw = PetscPowScalar(aw,bm1);
522:         /* dw = bw * aw */
523:         dw = PetscPowScalar(aw,beta);
524:         gw = coef*bw*(t0 - tw);

526:         te = x[j][i+1];
527:         ae = 0.5*(t0 + te);
528:         be = PetscPowScalar(ae,bm1);
529:         /* de = be * ae; */
530:         de = PetscPowScalar(ae,beta);
531:         ge = coef*be*(te - t0);

533:         tn = x[j+1][i];
534:         an = 0.5*(t0 + tn);
535:         bn = PetscPowScalar(an,bm1);
536:         /* dn = bn * an; */
537:         dn = PetscPowScalar(an,beta);
538:         gn = coef*bn*(tn - t0);

540:         v[0] = -hydhx*(dw - gw);                            col[0].j = j;         col[0].i = i-1;
541:         v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
542:         v[2] = -hydhx*(de + ge);                            col[2].j = j;         col[2].i = i+1;
543:         v[3] = -hxdhy*(dn + gn);                            col[3].j = j+1;       col[3].i = i;
544:         MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);

546:         /* top boundary,and i <> 0 or mx-1 */
547:       } else if (j == my-1) {

549:         tw = x[j][i-1];
550:         aw = 0.5*(t0 + tw);
551:         bw = PetscPowScalar(aw,bm1);
552:         /* dw = bw * aw */
553:         dw = PetscPowScalar(aw,beta);
554:         gw = coef*bw*(t0 - tw);

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

563:         ts = x[j-1][i];
564:         as = 0.5*(t0 + ts);
565:         bs = PetscPowScalar(as,bm1);
566:         /* ds = bs * as; */
567:         ds = PetscPowScalar(as,beta);
568:         gs = coef*bs*(t0 - ts);

570:         v[0] = -hxdhy*(ds - gs);                             col[0].j = j-1;       col[0].i = i;
571:         v[1] = -hydhx*(dw - gw);                             col[1].j = j;         col[1].i = i-1;
572:         v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
573:         v[3] = -hydhx*(de + ge);                             col[3].j = j;         col[3].i = i+1;
574:         MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);

576:       }
577:     }
578:   }
579:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
580:   DMDAVecRestoreArray(da,localX,&x);
581:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
582:   DMRestoreLocalVector(da,&localX);
583:   if (jac != B) {
584:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
585:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
586:   }

588:   PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
589:   return 0;
590: }

592: /*TEST

594:    test:
595:       suffix: 1
596:       args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view
597:       requires: !single

599:    test:
600:       suffix: 2
601:       args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view -snes_type newtontrdc
602:       requires: !single

604: TEST*/