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

 21: /*  
 22:   
 23:     This example models the partial differential equation 
 24:    
 25:          - Div(alpha* T^beta (GRAD T)) = 0.
 26:        
 27:     where beta = 2.5 and alpha = 1.0
 28:  
 29:     BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0.
 30:     
 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.
 33:  
 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
 39:  
 40: */

 42:  #include petscsnes.h
 43:  #include petscda.h
 44:  #include petscmg.h
 45:  #include petscdmmg.h

 47: /* User-defined application context */

 49: typedef struct {
 50:    PetscReal  tleft,tright;  /* Dirichlet boundary conditions */
 51:    PetscReal  beta,bm1,coef; /* nonlinear diffusivity parameterizations */
 52: } AppCtx;

 54: #define POWFLOP 5 /* assume a pow() takes five flops */


 62: int main(int argc,char **argv)
 63: {
 64:   DMMG           *dmmg;
 65:   SNES           snes;
 66:   AppCtx         user;
 68:   PetscInt       its,lits;
 69:   PetscReal      litspit;
 70:   DA             da;

 72:   PetscInitialize(&argc,&argv,PETSC_NULL,help);

 74:   /* set problem parameters */
 75:   user.tleft  = 1.0;
 76:   user.tright = 0.1;
 77:   user.beta   = 2.5;
 78:   PetscOptionsGetReal(PETSC_NULL,"-tleft",&user.tleft,PETSC_NULL);
 79:   PetscOptionsGetReal(PETSC_NULL,"-tright",&user.tright,PETSC_NULL);
 80:   PetscOptionsGetReal(PETSC_NULL,"-beta",&user.beta,PETSC_NULL);
 81:   user.bm1  = user.beta - 1.0;
 82:   user.coef = user.beta/2.0;


 85:   /*
 86:       Create the multilevel DA data structure 
 87:   */
 88:   DMMGCreate(PETSC_COMM_WORLD,3,&user,&dmmg);

 90:   /*
 91:       Set the DA (grid structure) for the grids.
 92:   */
 93:   DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,5,5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
 94:   DMMGSetDM(dmmg,(DM)da);
 95:   DADestroy(da);

 97:   /*
 98:      Create the nonlinear solver, and tell the DMMG structure to use it
 99:   */
100:   DMMGSetSNES(dmmg,FormFunction,FormJacobian);
101:   DMMGSetFromOptions(dmmg);

103:   /*
104:       PreLoadBegin() means that the following section of code is run twice. The first time
105:      through the flag PreLoading is on this the nonlinear solver is only run for a single step.
106:      The second time through (the actually timed code) the maximum iterations is set to 10
107:      Preload of the executable is done to eliminate from the timing the time spent bring the 
108:      executable into memory from disk (paging in).
109:   */
110:   PreLoadBegin(PETSC_TRUE,"Solve");
111:     DMMGSetInitialGuess(dmmg,FormInitialGuess);
112:     DMMGSolve(dmmg);
113:   PreLoadEnd();
114:   snes = DMMGGetSNES(dmmg);
115:   SNESGetIterationNumber(snes,&its);
116:   SNESGetLinearSolveIterations(snes,&lits);
117:   litspit = ((PetscReal)lits)/((PetscReal)its);
118:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);
119:   PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
120:   PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / Newton = %e\n",litspit);

122:   DMMGDestroy(dmmg);
123:   PetscFinalize();

125:   return 0;
126: }
127: /* --------------------  Form initial approximation ----------------- */
130: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
131: {
132:   AppCtx         *user = (AppCtx*)dmmg->user;
133:   PetscInt       i,j,xs,ys,xm,ym;
135:   PetscReal      tleft = user->tleft;
136:   PetscScalar    **x;


140:   /* Get ghost points */
141:   DAGetCorners((DA)dmmg->dm,&xs,&ys,0,&xm,&ym,0);
142:   DAVecGetArray((DA)dmmg->dm,X,&x);

144:   /* Compute initial guess */
145:   for (j=ys; j<ys+ym; j++) {
146:     for (i=xs; i<xs+xm; i++) {
147:       x[j][i] = tleft;
148:     }
149:   }
150:   DAVecRestoreArray((DA)dmmg->dm,X,&x);
151:   return(0);
152: }
153: /* --------------------  Evaluate Function F(x) --------------------- */
156: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ptr)
157: {
158:   DMMG           dmmg = (DMMG)ptr;
159:   AppCtx         *user = (AppCtx*)dmmg->user;
161:   PetscInt       i,j,mx,my,xs,ys,xm,ym;
162:   PetscScalar    zero = 0.0,one = 1.0;
163:   PetscScalar    hx,hy,hxdhy,hydhx;
164:   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;
165:   PetscScalar    tleft,tright,beta;
166:   PetscScalar    **x,**f;
167:   Vec            localX;

170:   DAGetLocalVector((DA)dmmg->dm,&localX);
171:   DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0);
172:   hx    = one/(PetscReal)(mx-1);  hy    = one/(PetscReal)(my-1);
173:   hxdhy = hx/hy;               hydhx = hy/hx;
174:   tleft = user->tleft;         tright = user->tright;
175:   beta  = user->beta;
176: 
177:   /* Get ghost points */
178:   DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
179:   DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
180:   DAGetCorners((DA)dmmg->dm,&xs,&ys,0,&xm,&ym,0);
181:   DAVecGetArray((DA)dmmg->dm,localX,&x);
182:   DAVecGetArray((DA)dmmg->dm,F,&f);

184:   /* Evaluate function */
185:   for (j=ys; j<ys+ym; j++) {
186:     for (i=xs; i<xs+xm; i++) {
187:       t0 = x[j][i];

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

191:         /* general interior volume */

193:         tw = x[j][i-1];
194:         aw = 0.5*(t0 + tw);
195:         dw = PetscPowScalar(aw,beta);
196:         fw = dw*(t0 - tw);

198:         te = x[j][i+1];
199:         ae = 0.5*(t0 + te);
200:         de = PetscPowScalar(ae,beta);
201:         fe = de*(te - t0);

203:         ts = x[j-1][i];
204:         as = 0.5*(t0 + ts);
205:         ds = PetscPowScalar(as,beta);
206:         fs = ds*(t0 - ts);
207: 
208:         tn = x[j+1][i];
209:         an = 0.5*(t0 + tn);
210:         dn = PetscPowScalar(an,beta);
211:         fn = dn*(tn - t0);

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

215:         /* left-hand boundary */
216:         tw = tleft;
217:         aw = 0.5*(t0 + tw);
218:         dw = PetscPowScalar(aw,beta);
219:         fw = dw*(t0 - tw);

221:         te = x[j][i+1];
222:         ae = 0.5*(t0 + te);
223:         de = PetscPowScalar(ae,beta);
224:         fe = de*(te - t0);

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

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

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

246:         /* right-hand boundary */
247:         tw = x[j][i-1];
248:         aw = 0.5*(t0 + tw);
249:         dw = PetscPowScalar(aw,beta);
250:         fw = dw*(t0 - tw);
251: 
252:         te = tright;
253:         ae = 0.5*(t0 + te);
254:         de = PetscPowScalar(ae,beta);
255:         fe = de*(te - t0);
256: 
257:         if (j > 0) {
258:           ts = x[j-1][i];
259:           as = 0.5*(t0 + ts);
260:           ds = PetscPowScalar(as,beta);
261:           fs = ds*(t0 - ts);
262:         } else {
263:           fs = zero;
264:         }
265: 
266:         if (j < my-1) {
267:           tn = x[j+1][i];
268:           an = 0.5*(t0 + tn);
269:           dn = PetscPowScalar(an,beta);
270:           fn = dn*(tn - t0);
271:         } else {
272:           fn = zero;
273:         }

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

277:         /* bottom 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:         fs = zero;

290:         tn = x[j+1][i];
291:         an = 0.5*(t0 + tn);
292:         dn = PetscPowScalar(an,beta);
293:         fn = dn*(tn - t0);

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

297:         /* top boundary,and i <> 0 or mx-1 */
298:         tw = x[j][i-1];
299:         aw = 0.5*(t0 + tw);
300:         dw = PetscPowScalar(aw,beta);
301:         fw = dw*(t0 - tw);

303:         te = x[j][i+1];
304:         ae = 0.5*(t0 + te);
305:         de = PetscPowScalar(ae,beta);
306:         fe = de*(te - t0);

308:         ts = x[j-1][i];
309:         as = 0.5*(t0 + ts);
310:         ds = PetscPowScalar(as,beta);
311:         fs = ds*(t0 - ts);

313:         fn = zero;

315:       }

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

319:     }
320:   }
321:   DAVecRestoreArray((DA)dmmg->dm,localX,&x);
322:   DAVecRestoreArray((DA)dmmg->dm,F,&f);
323:   DARestoreLocalVector((DA)dmmg->dm,&localX);
324:   PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
325:   return(0);
326: }
327: /* --------------------  Evaluate Jacobian F(x) --------------------- */
330: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ptr)
331: {
332:   DMMG           dmmg = (DMMG)ptr;
333:   AppCtx         *user = (AppCtx*)dmmg->user;
334:   Mat            jac = *J;
336:   PetscInt       i,j,mx,my,xs,ys,xm,ym;
337:   PetscScalar    one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;
338:   PetscScalar    dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;
339:   PetscScalar    tleft,tright,beta,bm1,coef;
340:   PetscScalar    v[5],**x;
341:   Vec            localX;
342:   MatStencil     col[5],row;

345:   DAGetLocalVector((DA)dmmg->dm,&localX);
346:   *flg = SAME_NONZERO_PATTERN;
347:   DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0);
348:   hx    = one/(PetscReal)(mx-1);  hy     = one/(PetscReal)(my-1);
349:   hxdhy = hx/hy;               hydhx  = hy/hx;
350:   tleft = user->tleft;         tright = user->tright;
351:   beta  = user->beta;               bm1    = user->bm1;                coef = user->coef;

353:   /* Get ghost points */
354:   DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
355:   DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
356:   DAGetCorners((DA)dmmg->dm,&xs,&ys,0,&xm,&ym,0);
357:   DAVecGetArray((DA)dmmg->dm,localX,&x);

359:   /* Evaluate Jacobian of function */
360:   for (j=ys; j<ys+ym; j++) {
361:     for (i=xs; i<xs+xm; i++) {
362:       t0 = x[j][i];

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

366:         /* general interior volume */

368:         tw = x[j][i-1];
369:         aw = 0.5*(t0 + tw);
370:         bw = PetscPowScalar(aw,bm1);
371:         /* dw = bw * aw */
372:         dw = PetscPowScalar(aw,beta);
373:         gw = coef*bw*(t0 - tw);

375:         te = x[j][i+1];
376:         ae = 0.5*(t0 + te);
377:         be = PetscPowScalar(ae,bm1);
378:         /* de = be * ae; */
379:         de = PetscPowScalar(ae,beta);
380:         ge = coef*be*(te - t0);

382:         ts = x[j-1][i];
383:         as = 0.5*(t0 + ts);
384:         bs = PetscPowScalar(as,bm1);
385:         /* ds = bs * as; */
386:         ds = PetscPowScalar(as,beta);
387:         gs = coef*bs*(t0 - ts);
388: 
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*(ds - gs);                                      col[0].j = j-1;       col[0].i = i;
397:         v[1] = - hydhx*(dw - gw);                                      col[1].j = j;         col[1].i = i-1;
398:         v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
399:         v[3] = - hydhx*(de + ge);                                      col[3].j = j;         col[3].i = i+1;
400:         v[4] = - hxdhy*(dn + gn);                                      col[4].j = j+1;       col[4].i = i;
401:         MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);

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

405:         /* left-hand boundary */
406:         tw = tleft;
407:         aw = 0.5*(t0 + tw);
408:         bw = PetscPowScalar(aw,bm1);
409:         /* dw = bw * aw */
410:         dw = PetscPowScalar(aw,beta);
411:         gw = coef*bw*(t0 - tw);
412: 
413:         te = x[j][i + 1];
414:         ae = 0.5*(t0 + te);
415:         be = PetscPowScalar(ae,bm1);
416:         /* de = be * ae; */
417:         de = PetscPowScalar(ae,beta);
418:         ge = coef*be*(te - t0);
419: 
420:         /* left-hand bottom boundary */
421:         if (j == 0) {

423:           tn = x[j+1][i];
424:           an = 0.5*(t0 + tn);
425:           bn = PetscPowScalar(an,bm1);
426:           /* dn = bn * an; */
427:           dn = PetscPowScalar(an,beta);
428:           gn = coef*bn*(tn - t0);
429: 
430:           v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i;
431:           v[1] = - hydhx*(de + ge);                           col[1].j = j;         col[1].i = i+1;
432:           v[2] = - hxdhy*(dn + gn);                           col[2].j = j+1;       col[2].i = i;
433:           MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
434: 
435:         /* left-hand interior boundary */
436:         } else if (j < my-1) {

438:           ts = x[j-1][i];
439:           as = 0.5*(t0 + ts);
440:           bs = PetscPowScalar(as,bm1);
441:           /* ds = bs * as; */
442:           ds = PetscPowScalar(as,beta);
443:           gs = coef*bs*(t0 - ts);
444: 
445:           tn = x[j+1][i];
446:           an = 0.5*(t0 + tn);
447:           bn = PetscPowScalar(an,bm1);
448:           /* dn = bn * an; */
449:           dn = PetscPowScalar(an,beta);
450:           gn = coef*bn*(tn - t0);
451: 
452:           v[0] = - hxdhy*(ds - gs);                                      col[0].j = j-1;       col[0].i = i;
453:           v[1] = hxdhy*(ds + dn + gs - gn) + 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:           v[3] = - hxdhy*(dn + gn);                                      col[3].j = j+1;       col[3].i = i;
456:           MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
457:         /* left-hand top boundary */
458:         } else {

460:           ts = x[j-1][i];
461:           as = 0.5*(t0 + ts);
462:           bs = PetscPowScalar(as,bm1);
463:           /* ds = bs * as; */
464:           ds = PetscPowScalar(as,beta);
465:           gs = coef*bs*(t0 - ts);
466: 
467:           v[0] = - hxdhy*(ds - gs);                            col[0].j = j-1;       col[0].i = i;
468:           v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[1].j = row.j = j; col[1].i = row.i = i;
469:           v[2] = - hydhx*(de + ge);                            col[2].j = j;         col[2].i = i+1;
470:           MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
471:         }

473:       } else if (i == mx-1) {
474: 
475:         /* right-hand boundary */
476:         tw = x[j][i-1];
477:         aw = 0.5*(t0 + tw);
478:         bw = PetscPowScalar(aw,bm1);
479:         /* dw = bw * aw */
480:         dw = PetscPowScalar(aw,beta);
481:         gw = coef*bw*(t0 - tw);
482: 
483:         te = tright;
484:         ae = 0.5*(t0 + te);
485:         be = PetscPowScalar(ae,bm1);
486:         /* de = be * ae; */
487:         de = PetscPowScalar(ae,beta);
488:         ge = coef*be*(te - t0);
489: 
490:         /* right-hand bottom boundary */
491:         if (j == 0) {

493:           tn = x[j+1][i];
494:           an = 0.5*(t0 + tn);
495:           bn = PetscPowScalar(an,bm1);
496:           /* dn = bn * an; */
497:           dn = PetscPowScalar(an,beta);
498:           gn = coef*bn*(tn - t0);
499: 
500:           v[0] = - hydhx*(dw - gw);                           col[0].j = j;         col[0].i = i-1;
501:           v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
502:           v[2] = - hxdhy*(dn + gn);                           col[2].j = j+1;       col[2].i = i;
503:           MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
504: 
505:         /* right-hand interior boundary */
506:         } else if (j < my-1) {

508:           ts = x[j-1][i];
509:           as = 0.5*(t0 + ts);
510:           bs = PetscPowScalar(as,bm1);
511:           /* ds = bs * as; */
512:           ds = PetscPowScalar(as,beta);
513:           gs = coef*bs*(t0 - ts);
514: 
515:           tn = x[j+1][i];
516:           an = 0.5*(t0 + tn);
517:           bn = PetscPowScalar(an,bm1);
518:           /* dn = bn * an; */
519:           dn = PetscPowScalar(an,beta);
520:           gn = coef*bn*(tn - t0);
521: 
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 + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
525:           v[3] = - hxdhy*(dn + gn);                                      col[3].j = j+1;       col[3].i = i;
526:           MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
527:         /* right-hand top boundary */
528:         } else {

530:           ts = x[j-1][i];
531:           as = 0.5*(t0 + ts);
532:           bs = PetscPowScalar(as,bm1);
533:           /* ds = bs * as; */
534:           ds = PetscPowScalar(as,beta);
535:           gs = coef*bs*(t0 - ts);
536: 
537:           v[0] = - hxdhy*(ds - gs);                            col[0].j = j-1;       col[0].i = i;
538:           v[1] = - hydhx*(dw - gw);                            col[1].j = j;         col[1].i = i-1;
539:           v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
540:           MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
541:         }

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

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

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

560:         tn = x[j+1][i];
561:         an = 0.5*(t0 + tn);
562:         bn = PetscPowScalar(an,bm1);
563:         /* dn = bn * an; */
564:         dn = PetscPowScalar(an,beta);
565:         gn = coef*bn*(tn - t0);
566: 
567:         v[0] = - hydhx*(dw - gw);                           col[0].j = j;         col[0].i = i-1;
568:         v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
569:         v[2] = - hydhx*(de + ge);                           col[2].j = j;         col[2].i = i+1;
570:         v[3] = - hxdhy*(dn + gn);                           col[3].j = j+1;       col[3].i = i;
571:         MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
572: 
573:       /* top boundary,and i <> 0 or mx-1 */
574:       } else if (j == my-1) {
575: 
576:         tw = x[j][i-1];
577:         aw = 0.5*(t0 + tw);
578:         bw = PetscPowScalar(aw,bm1);
579:         /* dw = bw * aw */
580:         dw = PetscPowScalar(aw,beta);
581:         gw = coef*bw*(t0 - tw);

583:         te = x[j][i+1];
584:         ae = 0.5*(t0 + te);
585:         be = PetscPowScalar(ae,bm1);
586:         /* de = be * ae; */
587:         de = PetscPowScalar(ae,beta);
588:         ge = coef*be*(te - t0);

590:         ts = x[j-1][i];
591:         as = 0.5*(t0 + ts);
592:         bs = PetscPowScalar(as,bm1);
593:          /* ds = bs * as; */
594:         ds = PetscPowScalar(as,beta);
595:         gs = coef*bs*(t0 - ts);

597:         v[0] = - hxdhy*(ds - gs);                            col[0].j = j-1;       col[0].i = i;
598:         v[1] = - hydhx*(dw - gw);                            col[1].j = j;         col[1].i = i-1;
599:         v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
600:         v[3] = - hydhx*(de + ge);                            col[3].j = j;         col[3].i = i+1;
601:         MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
602: 
603:       }
604:     }
605:   }
606:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
607:   DAVecRestoreArray((DA)dmmg->dm,localX,&x);
608:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
609:   DARestoreLocalVector((DA)dmmg->dm,&localX);

611:   PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
612:   return(0);
613: }