Actual source code: ex18.c

petsc-3.7.3 2016-08-01
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);

 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;


 83:   /*
 84:       Create the multilevel DM data structure
 85:   */
 86:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

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

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

111:   DMDestroy(&da);
112:   SNESDestroy(&snes);
113:   PetscFinalize();

115:   return 0;
116: }
117: /* --------------------  Form initial approximation ----------------- */
120: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
121: {
122:   AppCtx         *user;
123:   PetscInt       i,j,xs,ys,xm,ym;
125:   PetscReal      tleft;
126:   PetscScalar    **x;
127:   DM             da;

130:   SNESGetDM(snes,&da);
131:   DMGetApplicationContext(da,&user);
132:   tleft = user->tleft;
133:   /* Get ghost points */
134:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
135:   DMDAVecGetArray(da,X,&x);

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

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

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

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

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

185:         /* general interior volume */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

274:         fs = zero;

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

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

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

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

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

299:         fn = zero;

301:       }

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

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

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

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

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

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

351:         /* general interior volume */

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

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

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

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

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

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

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

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

405:         /* left-hand bottom boundary */
406:         if (j == 0) {

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

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

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

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

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

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

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

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

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

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

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

475:         /* right-hand bottom boundary */
476:         if (j == 0) {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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