Actual source code: ex18.c

petsc-master 2017-02-22
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*/

 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;
 64:   PetscInt       its,lits;
 65:   PetscReal      litspit;
 66:   DM             da;

 68:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;

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

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

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

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

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

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

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

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

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

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

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

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

179:         /* general interior volume */

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

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

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

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

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

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

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

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

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

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

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

236:         te = tright;
237:         ae = 0.5*(t0 + te);
238:         de = PetscPowScalar(ae,beta);
239:         fe = de*(te - t0);

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

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

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

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

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

268:         fs = zero;

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

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

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

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

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

293:         fn = zero;

295:       }

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

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

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

330:   /* Get ghost points */
331:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
332:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
333:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
334:   DMDAVecGetArray(da,localX,&x);

336:   /* Evaluate Jacobian of function */
337:   for (j=ys; j<ys+ym; j++) {
338:     for (i=xs; i<xs+xm; i++) {
339:       t0 = x[j][i];

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

343:         /* general interior volume */

345:         tw = x[j][i-1];
346:         aw = 0.5*(t0 + tw);
347:         bw = PetscPowScalar(aw,bm1);
348:         /* dw = bw * aw */
349:         dw = PetscPowScalar(aw,beta);
350:         gw = coef*bw*(t0 - tw);

352:         te = x[j][i+1];
353:         ae = 0.5*(t0 + te);
354:         be = PetscPowScalar(ae,bm1);
355:         /* de = be * ae; */
356:         de = PetscPowScalar(ae,beta);
357:         ge = coef*be*(te - t0);

359:         ts = x[j-1][i];
360:         as = 0.5*(t0 + ts);
361:         bs = PetscPowScalar(as,bm1);
362:         /* ds = bs * as; */
363:         ds = PetscPowScalar(as,beta);
364:         gs = coef*bs*(t0 - ts);

366:         tn = x[j+1][i];
367:         an = 0.5*(t0 + tn);
368:         bn = PetscPowScalar(an,bm1);
369:         /* dn = bn * an; */
370:         dn = PetscPowScalar(an,beta);
371:         gn = coef*bn*(tn - t0);

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

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

382:         /* left-hand boundary */
383:         tw = tleft;
384:         aw = 0.5*(t0 + tw);
385:         bw = PetscPowScalar(aw,bm1);
386:         /* dw = bw * aw */
387:         dw = PetscPowScalar(aw,beta);
388:         gw = coef*bw*(t0 - tw);

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

397:         /* left-hand bottom boundary */
398:         if (j == 0) {

400:           tn = x[j+1][i];
401:           an = 0.5*(t0 + tn);
402:           bn = PetscPowScalar(an,bm1);
403:           /* dn = bn * an; */
404:           dn = PetscPowScalar(an,beta);
405:           gn = coef*bn*(tn - t0);

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

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

415:           ts = x[j-1][i];
416:           as = 0.5*(t0 + ts);
417:           bs = PetscPowScalar(as,bm1);
418:           /* ds = bs * as; */
419:           ds = PetscPowScalar(as,beta);
420:           gs = coef*bs*(t0 - ts);

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

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

437:           ts = x[j-1][i];
438:           as = 0.5*(t0 + ts);
439:           bs = PetscPowScalar(as,bm1);
440:           /* ds = bs * as; */
441:           ds = PetscPowScalar(as,beta);
442:           gs = coef*bs*(t0 - ts);

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

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

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

460:         te = tright;
461:         ae = 0.5*(t0 + te);
462:         be = PetscPowScalar(ae,bm1);
463:         /* de = be * ae; */
464:         de = PetscPowScalar(ae,beta);
465:         ge = coef*be*(te - t0);

467:         /* right-hand bottom boundary */
468:         if (j == 0) {

470:           tn = x[j+1][i];
471:           an = 0.5*(t0 + tn);
472:           bn = PetscPowScalar(an,bm1);
473:           /* dn = bn * an; */
474:           dn = PetscPowScalar(an,beta);
475:           gn = coef*bn*(tn - t0);

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

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

485:           ts = x[j-1][i];
486:           as = 0.5*(t0 + ts);
487:           bs = PetscPowScalar(as,bm1);
488:           /* ds = bs * as; */
489:           ds = PetscPowScalar(as,beta);
490:           gs = coef*bs*(t0 - ts);

492:           tn = x[j+1][i];
493:           an = 0.5*(t0 + tn);
494:           bn = PetscPowScalar(an,bm1);
495:           /* dn = bn * an; */
496:           dn = PetscPowScalar(an,beta);
497:           gn = coef*bn*(tn - t0);

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

507:           ts = x[j-1][i];
508:           as = 0.5*(t0 + ts);
509:           bs = PetscPowScalar(as,bm1);
510:           /* ds = bs * as; */
511:           ds = PetscPowScalar(as,beta);
512:           gs = coef*bs*(t0 - ts);

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

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

523:         tw = x[j][i-1];
524:         aw = 0.5*(t0 + tw);
525:         bw = PetscPowScalar(aw,bm1);
526:         /* dw = bw * aw */
527:         dw = PetscPowScalar(aw,beta);
528:         gw = coef*bw*(t0 - tw);

530:         te = x[j][i+1];
531:         ae = 0.5*(t0 + te);
532:         be = PetscPowScalar(ae,bm1);
533:         /* de = be * ae; */
534:         de = PetscPowScalar(ae,beta);
535:         ge = coef*be*(te - t0);

537:         tn = x[j+1][i];
538:         an = 0.5*(t0 + tn);
539:         bn = PetscPowScalar(an,bm1);
540:         /* dn = bn * an; */
541:         dn = PetscPowScalar(an,beta);
542:         gn = coef*bn*(tn - t0);

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

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

553:         tw = x[j][i-1];
554:         aw = 0.5*(t0 + tw);
555:         bw = PetscPowScalar(aw,bm1);
556:         /* dw = bw * aw */
557:         dw = PetscPowScalar(aw,beta);
558:         gw = coef*bw*(t0 - tw);

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

567:         ts = x[j-1][i];
568:         as = 0.5*(t0 + ts);
569:         bs = PetscPowScalar(as,bm1);
570:         /* ds = bs * as; */
571:         ds = PetscPowScalar(as,beta);
572:         gs = coef*bs*(t0 - ts);

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

580:       }
581:     }
582:   }
583:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
584:   DMDAVecRestoreArray(da,localX,&x);
585:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
586:   DMRestoreLocalVector(da,&localX);
587:   if (jac != B) {
588:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
589:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
590:   }

592:   PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
593:   return(0);
594: }