Actual source code: ex18.c

petsc-master 2016-07-24
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*);

 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:   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 ----------------- */
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) --------------------- */
147: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
148: {
149:   AppCtx         *user = (AppCtx*)ptr;
151:   PetscInt       i,j,mx,my,xs,ys,xm,ym;
152:   PetscScalar    zero = 0.0,one = 1.0;
153:   PetscScalar    hx,hy,hxdhy,hydhx;
154:   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;
155:   PetscScalar    tleft,tright,beta;
156:   PetscScalar    **x,**f;
157:   Vec            localX;
158:   DM             da;

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

169:   /* Get ghost points */
170:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
171:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
172:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
173:   DMDAVecGetArray(da,localX,&x);
174:   DMDAVecGetArray(da,F,&f);

176:   /* Evaluate function */
177:   for (j=ys; j<ys+ym; j++) {
178:     for (i=xs; i<xs+xm; i++) {
179:       t0 = x[j][i];

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

183:         /* general interior volume */

185:         tw = x[j][i-1];
186:         aw = 0.5*(t0 + tw);
187:         dw = PetscPowScalar(aw,beta);
188:         fw = dw*(t0 - tw);

190:         te = x[j][i+1];
191:         ae = 0.5*(t0 + te);
192:         de = PetscPowScalar(ae,beta);
193:         fe = de*(te - t0);

195:         ts = x[j-1][i];
196:         as = 0.5*(t0 + ts);
197:         ds = PetscPowScalar(as,beta);
198:         fs = ds*(t0 - ts);

200:         tn = x[j+1][i];
201:         an = 0.5*(t0 + tn);
202:         dn = PetscPowScalar(an,beta);
203:         fn = dn*(tn - t0);

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

207:         /* left-hand boundary */
208:         tw = tleft;
209:         aw = 0.5*(t0 + tw);
210:         dw = PetscPowScalar(aw,beta);
211:         fw = dw*(t0 - tw);

213:         te = x[j][i+1];
214:         ae = 0.5*(t0 + te);
215:         de = PetscPowScalar(ae,beta);
216:         fe = de*(te - t0);

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

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

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

234:         /* right-hand boundary */
235:         tw = x[j][i-1];
236:         aw = 0.5*(t0 + tw);
237:         dw = PetscPowScalar(aw,beta);
238:         fw = dw*(t0 - tw);

240:         te = tright;
241:         ae = 0.5*(t0 + te);
242:         de = PetscPowScalar(ae,beta);
243:         fe = de*(te - t0);

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

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

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

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

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

272:         fs = zero;

274:         tn = x[j+1][i];
275:         an = 0.5*(t0 + tn);
276:         dn = PetscPowScalar(an,beta);
277:         fn = dn*(tn - t0);

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

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

287:         te = x[j][i+1];
288:         ae = 0.5*(t0 + te);
289:         de = PetscPowScalar(ae,beta);
290:         fe = de*(te - t0);

292:         ts = x[j-1][i];
293:         as = 0.5*(t0 + ts);
294:         ds = PetscPowScalar(as,beta);
295:         fs = ds*(t0 - ts);

297:         fn = zero;

299:       }

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

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

328:   SNESGetDM(snes,&da);
329:   DMGetLocalVector(da,&localX);
330:   DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
331:   hx    = one/(PetscReal)(mx-1);  hy     = one/(PetscReal)(my-1);
332:   hxdhy = hx/hy;               hydhx  = hy/hx;
333:   tleft = user->tleft;         tright = user->tright;
334:   beta  = user->beta;          bm1    = user->bm1;          coef = user->coef;

336:   /* Get ghost points */
337:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
338:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
339:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
340:   DMDAVecGetArray(da,localX,&x);

342:   /* Evaluate Jacobian of function */
343:   for (j=ys; j<ys+ym; j++) {
344:     for (i=xs; i<xs+xm; i++) {
345:       t0 = x[j][i];

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

349:         /* general interior volume */

351:         tw = x[j][i-1];
352:         aw = 0.5*(t0 + tw);
353:         bw = PetscPowScalar(aw,bm1);
354:         /* dw = bw * aw */
355:         dw = PetscPowScalar(aw,beta);
356:         gw = coef*bw*(t0 - tw);

358:         te = x[j][i+1];
359:         ae = 0.5*(t0 + te);
360:         be = PetscPowScalar(ae,bm1);
361:         /* de = be * ae; */
362:         de = PetscPowScalar(ae,beta);
363:         ge = coef*be*(te - t0);

365:         ts = x[j-1][i];
366:         as = 0.5*(t0 + ts);
367:         bs = PetscPowScalar(as,bm1);
368:         /* ds = bs * as; */
369:         ds = PetscPowScalar(as,beta);
370:         gs = coef*bs*(t0 - ts);

372:         tn = x[j+1][i];
373:         an = 0.5*(t0 + tn);
374:         bn = PetscPowScalar(an,bm1);
375:         /* dn = bn * an; */
376:         dn = PetscPowScalar(an,beta);
377:         gn = coef*bn*(tn - t0);

379:         v[0] = -hxdhy*(ds - gs);                                       col[0].j = j-1;       col[0].i = i;
380:         v[1] = -hydhx*(dw - gw);                                       col[1].j = j;         col[1].i = i-1;
381:         v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
382:         v[3] = -hydhx*(de + ge);                                       col[3].j = j;         col[3].i = i+1;
383:         v[4] = -hxdhy*(dn + gn);                                       col[4].j = j+1;       col[4].i = i;
384:         MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);

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

388:         /* left-hand boundary */
389:         tw = tleft;
390:         aw = 0.5*(t0 + tw);
391:         bw = PetscPowScalar(aw,bm1);
392:         /* dw = bw * aw */
393:         dw = PetscPowScalar(aw,beta);
394:         gw = coef*bw*(t0 - tw);

396:         te = x[j][i + 1];
397:         ae = 0.5*(t0 + te);
398:         be = PetscPowScalar(ae,bm1);
399:         /* de = be * ae; */
400:         de = PetscPowScalar(ae,beta);
401:         ge = coef*be*(te - t0);

403:         /* left-hand bottom boundary */
404:         if (j == 0) {

406:           tn = x[j+1][i];
407:           an = 0.5*(t0 + tn);
408:           bn = PetscPowScalar(an,bm1);
409:           /* dn = bn * an; */
410:           dn = PetscPowScalar(an,beta);
411:           gn = coef*bn*(tn - t0);

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

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

421:           ts = x[j-1][i];
422:           as = 0.5*(t0 + ts);
423:           bs = PetscPowScalar(as,bm1);
424:           /* ds = bs * as; */
425:           ds = PetscPowScalar(as,beta);
426:           gs = coef*bs*(t0 - ts);

428:           tn = x[j+1][i];
429:           an = 0.5*(t0 + tn);
430:           bn = PetscPowScalar(an,bm1);
431:           /* dn = bn * an; */
432:           dn = PetscPowScalar(an,beta);
433:           gn = coef*bn*(tn - t0);

435:           v[0] = -hxdhy*(ds - gs);                                       col[0].j = j-1;       col[0].i = i;
436:           v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[1].j = row.j = j; col[1].i = row.i = i;
437:           v[2] = -hydhx*(de + ge);                                       col[2].j = j;         col[2].i = i+1;
438:           v[3] = -hxdhy*(dn + gn);                                       col[3].j = j+1;       col[3].i = i;
439:           MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
440:           /* left-hand top boundary */
441:         } else {

443:           ts = x[j-1][i];
444:           as = 0.5*(t0 + ts);
445:           bs = PetscPowScalar(as,bm1);
446:           /* ds = bs * as; */
447:           ds = PetscPowScalar(as,beta);
448:           gs = coef*bs*(t0 - ts);

450:           v[0] = -hxdhy*(ds - gs);                             col[0].j = j-1;       col[0].i = i;
451:           v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[1].j = row.j = j; col[1].i = row.i = i;
452:           v[2] = -hydhx*(de + ge);                             col[2].j = j;         col[2].i = i+1;
453:           MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
454:         }

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

458:         /* right-hand boundary */
459:         tw = x[j][i-1];
460:         aw = 0.5*(t0 + tw);
461:         bw = PetscPowScalar(aw,bm1);
462:         /* dw = bw * aw */
463:         dw = PetscPowScalar(aw,beta);
464:         gw = coef*bw*(t0 - tw);

466:         te = tright;
467:         ae = 0.5*(t0 + te);
468:         be = PetscPowScalar(ae,bm1);
469:         /* de = be * ae; */
470:         de = PetscPowScalar(ae,beta);
471:         ge = coef*be*(te - t0);

473:         /* right-hand bottom boundary */
474:         if (j == 0) {

476:           tn = x[j+1][i];
477:           an = 0.5*(t0 + tn);
478:           bn = PetscPowScalar(an,bm1);
479:           /* dn = bn * an; */
480:           dn = PetscPowScalar(an,beta);
481:           gn = coef*bn*(tn - t0);

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

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

491:           ts = x[j-1][i];
492:           as = 0.5*(t0 + ts);
493:           bs = PetscPowScalar(as,bm1);
494:           /* ds = bs * as; */
495:           ds = PetscPowScalar(as,beta);
496:           gs = coef*bs*(t0 - ts);

498:           tn = x[j+1][i];
499:           an = 0.5*(t0 + tn);
500:           bn = PetscPowScalar(an,bm1);
501:           /* dn = bn * an; */
502:           dn = PetscPowScalar(an,beta);
503:           gn = coef*bn*(tn - t0);

505:           v[0] = -hxdhy*(ds - gs);                                       col[0].j = j-1;       col[0].i = i;
506:           v[1] = -hydhx*(dw - gw);                                       col[1].j = j;         col[1].i = i-1;
507:           v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
508:           v[3] = -hxdhy*(dn + gn);                                       col[3].j = j+1;       col[3].i = i;
509:           MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
510:         /* right-hand top boundary */
511:         } else {

513:           ts = x[j-1][i];
514:           as = 0.5*(t0 + ts);
515:           bs = PetscPowScalar(as,bm1);
516:           /* ds = bs * as; */
517:           ds = PetscPowScalar(as,beta);
518:           gs = coef*bs*(t0 - ts);

520:           v[0] = -hxdhy*(ds - gs);                             col[0].j = j-1;       col[0].i = i;
521:           v[1] = -hydhx*(dw - gw);                             col[1].j = j;         col[1].i = i-1;
522:           v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
523:           MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
524:         }

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

529:         tw = x[j][i-1];
530:         aw = 0.5*(t0 + tw);
531:         bw = PetscPowScalar(aw,bm1);
532:         /* dw = bw * aw */
533:         dw = PetscPowScalar(aw,beta);
534:         gw = coef*bw*(t0 - tw);

536:         te = x[j][i+1];
537:         ae = 0.5*(t0 + te);
538:         be = PetscPowScalar(ae,bm1);
539:         /* de = be * ae; */
540:         de = PetscPowScalar(ae,beta);
541:         ge = coef*be*(te - t0);

543:         tn = x[j+1][i];
544:         an = 0.5*(t0 + tn);
545:         bn = PetscPowScalar(an,bm1);
546:         /* dn = bn * an; */
547:         dn = PetscPowScalar(an,beta);
548:         gn = coef*bn*(tn - t0);

550:         v[0] = -hydhx*(dw - gw);                            col[0].j = j;         col[0].i = i-1;
551:         v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
552:         v[2] = -hydhx*(de + ge);                            col[2].j = j;         col[2].i = i+1;
553:         v[3] = -hxdhy*(dn + gn);                            col[3].j = j+1;       col[3].i = i;
554:         MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);

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

559:         tw = x[j][i-1];
560:         aw = 0.5*(t0 + tw);
561:         bw = PetscPowScalar(aw,bm1);
562:         /* dw = bw * aw */
563:         dw = PetscPowScalar(aw,beta);
564:         gw = coef*bw*(t0 - tw);

566:         te = x[j][i+1];
567:         ae = 0.5*(t0 + te);
568:         be = PetscPowScalar(ae,bm1);
569:         /* de = be * ae; */
570:         de = PetscPowScalar(ae,beta);
571:         ge = coef*be*(te - t0);

573:         ts = x[j-1][i];
574:         as = 0.5*(t0 + ts);
575:         bs = PetscPowScalar(as,bm1);
576:         /* ds = bs * as; */
577:         ds = PetscPowScalar(as,beta);
578:         gs = coef*bs*(t0 - ts);

580:         v[0] = -hxdhy*(ds - gs);                             col[0].j = j-1;       col[0].i = i;
581:         v[1] = -hydhx*(dw - gw);                             col[1].j = j;         col[1].i = i-1;
582:         v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
583:         v[3] = -hydhx*(de + ge);                             col[3].j = j;         col[3].i = i+1;
584:         MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);

586:       }
587:     }
588:   }
589:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
590:   DMDAVecRestoreArray(da,localX,&x);
591:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
592:   DMRestoreLocalVector(da,&localX);
593:   if (jac != B) {
594:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
595:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
596:   }

598:   PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
599:   return(0);
600: }