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;

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

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

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

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

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

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

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

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

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

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

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

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

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

169:         /* general interior volume */

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

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

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

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

191:       } 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) {

220:         /* right-hand boundary */
221:         tw = x[j][i-1];
222:         aw = 0.5*(t0 + tw);
223:         dw = PetscPowScalar(aw,beta);
224:         fw = dw*(t0 - tw);

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

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

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

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

247:         /* bottom boundary,and i <> 0 or mx-1 */
248:         tw = x[j][i-1];
249:         aw = 0.5*(t0 + tw);
250:         dw = PetscPowScalar(aw,beta);
251:         fw = dw*(t0 - tw);

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

258:         fs = zero;

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

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

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

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

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

283:         fn = zero;

285:       }

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

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

311:   SNESGetDM(snes,&da);
312:   DMGetLocalVector(da,&localX);
313:   DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
314:   hx    = one/(PetscReal)(mx-1);  hy     = one/(PetscReal)(my-1);
315:   hxdhy = hx/hy;               hydhx  = hy/hx;
316:   tleft = user->tleft;         tright = user->tright;
317:   beta  = user->beta;          bm1    = user->bm1;          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) {

332:         /* general interior volume */

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

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

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

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

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

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

371:         /* left-hand boundary */
372:         tw = tleft;
373:         aw = 0.5*(t0 + tw);
374:         bw = PetscPowScalar(aw,bm1);
375:         /* dw = bw * aw */
376:         dw = PetscPowScalar(aw,beta);
377:         gw = coef*bw*(t0 - tw);

379:         te = x[j][i + 1];
380:         ae = 0.5*(t0 + te);
381:         be = PetscPowScalar(ae,bm1);
382:         /* de = be * ae; */
383:         de = PetscPowScalar(ae,beta);
384:         ge = coef*be*(te - t0);

386:         /* left-hand bottom boundary */
387:         if (j == 0) {

389:           tn = x[j+1][i];
390:           an = 0.5*(t0 + tn);
391:           bn = PetscPowScalar(an,bm1);
392:           /* dn = bn * an; */
393:           dn = PetscPowScalar(an,beta);
394:           gn = coef*bn*(tn - t0);

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

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

404:           ts = x[j-1][i];
405:           as = 0.5*(t0 + ts);
406:           bs = PetscPowScalar(as,bm1);
407:           /* ds = bs * as; */
408:           ds = PetscPowScalar(as,beta);
409:           gs = coef*bs*(t0 - ts);

411:           tn = x[j+1][i];
412:           an = 0.5*(t0 + tn);
413:           bn = PetscPowScalar(an,bm1);
414:           /* dn = bn * an; */
415:           dn = PetscPowScalar(an,beta);
416:           gn = coef*bn*(tn - t0);

418:           v[0] = -hxdhy*(ds - gs);                                       col[0].j = j-1;       col[0].i = i;
419:           v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[1].j = row.j = j; col[1].i = row.i = i;
420:           v[2] = -hydhx*(de + ge);                                       col[2].j = j;         col[2].i = i+1;
421:           v[3] = -hxdhy*(dn + gn);                                       col[3].j = j+1;       col[3].i = i;
422:           MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
423:           /* left-hand top boundary */
424:         } else {

426:           ts = x[j-1][i];
427:           as = 0.5*(t0 + ts);
428:           bs = PetscPowScalar(as,bm1);
429:           /* ds = bs * as; */
430:           ds = PetscPowScalar(as,beta);
431:           gs = coef*bs*(t0 - ts);

433:           v[0] = -hxdhy*(ds - gs);                             col[0].j = j-1;       col[0].i = i;
434:           v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[1].j = row.j = j; col[1].i = row.i = i;
435:           v[2] = -hydhx*(de + ge);                             col[2].j = j;         col[2].i = i+1;
436:           MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
437:         }

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

441:         /* right-hand boundary */
442:         tw = x[j][i-1];
443:         aw = 0.5*(t0 + tw);
444:         bw = PetscPowScalar(aw,bm1);
445:         /* dw = bw * aw */
446:         dw = PetscPowScalar(aw,beta);
447:         gw = coef*bw*(t0 - tw);

449:         te = tright;
450:         ae = 0.5*(t0 + te);
451:         be = PetscPowScalar(ae,bm1);
452:         /* de = be * ae; */
453:         de = PetscPowScalar(ae,beta);
454:         ge = coef*be*(te - t0);

456:         /* right-hand bottom boundary */
457:         if (j == 0) {

459:           tn = x[j+1][i];
460:           an = 0.5*(t0 + tn);
461:           bn = PetscPowScalar(an,bm1);
462:           /* dn = bn * an; */
463:           dn = PetscPowScalar(an,beta);
464:           gn = coef*bn*(tn - t0);

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

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

474:           ts = x[j-1][i];
475:           as = 0.5*(t0 + ts);
476:           bs = PetscPowScalar(as,bm1);
477:           /* ds = bs * as; */
478:           ds = PetscPowScalar(as,beta);
479:           gs = coef*bs*(t0 - ts);

481:           tn = x[j+1][i];
482:           an = 0.5*(t0 + tn);
483:           bn = PetscPowScalar(an,bm1);
484:           /* dn = bn * an; */
485:           dn = PetscPowScalar(an,beta);
486:           gn = coef*bn*(tn - t0);

488:           v[0] = -hxdhy*(ds - gs);                                       col[0].j = j-1;       col[0].i = i;
489:           v[1] = -hydhx*(dw - gw);                                       col[1].j = j;         col[1].i = i-1;
490:           v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
491:           v[3] = -hxdhy*(dn + gn);                                       col[3].j = j+1;       col[3].i = i;
492:           MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
493:         /* right-hand top boundary */
494:         } else {

496:           ts = x[j-1][i];
497:           as = 0.5*(t0 + ts);
498:           bs = PetscPowScalar(as,bm1);
499:           /* ds = bs * as; */
500:           ds = PetscPowScalar(as,beta);
501:           gs = coef*bs*(t0 - ts);

503:           v[0] = -hxdhy*(ds - gs);                             col[0].j = j-1;       col[0].i = i;
504:           v[1] = -hydhx*(dw - gw);                             col[1].j = j;         col[1].i = i-1;
505:           v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
506:           MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
507:         }

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

512:         tw = x[j][i-1];
513:         aw = 0.5*(t0 + tw);
514:         bw = PetscPowScalar(aw,bm1);
515:         /* dw = bw * aw */
516:         dw = PetscPowScalar(aw,beta);
517:         gw = coef*bw*(t0 - tw);

519:         te = x[j][i+1];
520:         ae = 0.5*(t0 + te);
521:         be = PetscPowScalar(ae,bm1);
522:         /* de = be * ae; */
523:         de = PetscPowScalar(ae,beta);
524:         ge = coef*be*(te - t0);

526:         tn = x[j+1][i];
527:         an = 0.5*(t0 + tn);
528:         bn = PetscPowScalar(an,bm1);
529:         /* dn = bn * an; */
530:         dn = PetscPowScalar(an,beta);
531:         gn = coef*bn*(tn - t0);

533:         v[0] = -hydhx*(dw - gw);                            col[0].j = j;         col[0].i = i-1;
534:         v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
535:         v[2] = -hydhx*(de + ge);                            col[2].j = j;         col[2].i = i+1;
536:         v[3] = -hxdhy*(dn + gn);                            col[3].j = j+1;       col[3].i = i;
537:         MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);

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

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

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

556:         ts = x[j-1][i];
557:         as = 0.5*(t0 + ts);
558:         bs = PetscPowScalar(as,bm1);
559:         /* ds = bs * as; */
560:         ds = PetscPowScalar(as,beta);
561:         gs = coef*bs*(t0 - ts);

563:         v[0] = -hxdhy*(ds - gs);                             col[0].j = j-1;       col[0].i = i;
564:         v[1] = -hydhx*(dw - gw);                             col[1].j = j;         col[1].i = i-1;
565:         v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
566:         v[3] = -hydhx*(de + ge);                             col[3].j = j;         col[3].i = i+1;
567:         MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);

569:       }
570:     }
571:   }
572:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
573:   DMDAVecRestoreArray(da,localX,&x);
574:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
575:   DMRestoreLocalVector(da,&localX);
576:   if (jac != B) {
577:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
578:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
579:   }

581:   PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
582:   return 0;
583: }

585: /*TEST

587:    test:
588:       suffix: 1
589:       args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view
590:       requires: !single

592:    test:
593:       suffix: 2
594:       args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view -snes_type newtontrdc
595:       requires: !single

597: TEST*/