Actual source code: ex4.c

  1: /*$Id: ex4.c,v 1.64 2001/08/07 21:31:12 bsmith Exp $*/

  3: /* NOTE:  THIS PROGRAM HAS NOT YET BEEN SET UP IN TUTORIAL STYLE. */

  5: static char help[] ="This program demonstrates use of the SNES package to solve systems of nonlinear equations on a single processor.\n\
  6:   Both of these examples employ\n\
  7: sparse storage of the Jacobian matrices and are taken from the MINPACK-2\n\
  8: test suite.  By default the Bratu (SFI - solid fuel ignition) test problem\n\
  9: is solved.  The command line options are:\n\
 10:    -cavity : Solve FDC (flow in a driven cavity) problem\n\
 11:    -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
 12:       problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
 13:       problem FDC:  <parameter> = Reynolds number (par > 0)\n\
 14:    -mx <xg>, where <xg> = number of grid points in the x-direction\n\
 15:    -my <yg>, where <yg> = number of grid points in the y-direction\n\n";

 17: /*  
 18:     1) Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 19:     the partial differential equation
 20:   
 21:             -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
 22:   
 23:     with boundary conditions
 24:    
 25:              u = 0  for  x = 0, x = 1, y = 0, y = 1.
 26:   
 27:     A finite difference approximation with the usual 5-point stencil
 28:     is used to discretize the boundary value problem to obtain a nonlinear 
 29:     system of equations.
 30:   
 31:     2) Flow in a Driven Cavity (FDC) problem. The problem is
 32:     formulated as a boundary value problem, which is discretized by
 33:     standard finite difference approximations to obtain a system of
 34:     nonlinear equations. 
 35: */

 37:  #include petscsnes.h

 39: typedef struct {
 40:       PetscReal   param;        /* test problem parameter */
 41:       int         mx;           /* Discretization in x-direction */
 42:       int         my;           /* Discretization in y-direction */
 43: } AppCtx;

 45: extern int  FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*),
 46:                           FormFunction1(SNES,Vec,Vec,void*),
 47:                           FormInitialGuess1(AppCtx*,Vec);
 48: extern int  FormJacobian2(SNES,Vec,Mat*,Mat*,MatStructure*,void*),
 49:                    FormFunction2(SNES,Vec,Vec,void*),
 50:                    FormInitialGuess2(AppCtx*,Vec);
 51: 
 54: int main(int argc, char **argv)
 55: {
 56:   SNES         snes;                 /* SNES context */
 57:   SNESType     method = SNESLS;      /* default nonlinear solution method */
 58:   Vec          x, r;                 /* solution, residual vectors */
 59:   Mat          J;                    /* Jacobian matrix */
 60:   AppCtx       user;                 /* user-defined application context */
 61:   PetscDraw    draw;                 /* drawing context */
 62:   int          ierr, its, N, nfails;
 63:   PetscTruth   flg,cavity;
 64:   PetscReal    bratu_lambda_max = 6.81, bratu_lambda_min = 0.;
 65:   PetscScalar  *xvalues;

 67:   PetscInitialize(&argc, &argv,(char *)0,help);
 68:   /* PetscDrawOpenX(PETSC_COMM_WORLD,0,"Solution",300,0,300,300,&draw); */
 69:   PetscDrawCreate(PETSC_COMM_WORLD,0,"Solution",300,0,300,300,&draw);
 70:   PetscDrawSetType(draw,PETSC_DRAW_X);

 72:   user.mx    = 4;
 73:   user.my    = 4;
 74:   PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
 75:   PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
 76:   PetscOptionsHasName(PETSC_NULL,"-cavity",&cavity);
 77:   if (cavity) user.param = 100.0;
 78:   else        user.param = 6.0;
 79:   PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
 80:   if (!cavity && (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min)) {
 81:     SETERRQ(1,"Lambda is out of range");
 82:   }
 83:   N = user.mx*user.my;
 84: 
 85:   /* Set up data structures */
 86:   VecCreateSeq(PETSC_COMM_SELF,N,&x);
 87:   VecDuplicate(x,&r);
 88:   MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,5,PETSC_NULL,&J);

 90:   /* Create nonlinear solver */
 91:   SNESCreate(PETSC_COMM_WORLD,&snes);
 92:   SNESSetType(snes,method);

 94:   /* Set various routines and compute initial guess for nonlinear solver */
 95:   if (cavity){
 96:     FormInitialGuess2(&user,x);
 97:     SNESSetFunction(snes,r,FormFunction2,(void *)&user);
 98:     SNESSetJacobian(snes,J,J,FormJacobian2,(void *)&user);
 99:   } else {
100:     FormInitialGuess1(&user,x);
101:     SNESSetFunction(snes,r,FormFunction1,(void *)&user);
102:     SNESSetJacobian(snes,J,J,FormJacobian1,(void *)&user);
103:   }

105:   /* Set options and solve nonlinear system */
106:   SNESSetFromOptions(snes);
107:   SNESSolve(snes,x);
108:   SNESGetIterationNumber(snes,&its);
109:   SNESGetNumberUnsuccessfulSteps(snes,&nfails);

111:   PetscPrintf(PETSC_COMM_SELF,"number of Newton iterations = %d, ",its);
112:   PetscPrintf(PETSC_COMM_SELF,"number of unsuccessful steps = %d\n\n",nfails);
113:   VecGetArray(x,&xvalues);
114:   PetscDrawTensorContour(draw,user.mx,user.my,0,0,xvalues);
115:   VecRestoreArray(x,&xvalues);

117:   /* Free data structures */
118:   VecDestroy(x);  VecDestroy(r);
119:   MatDestroy(J);  SNESDestroy(snes);
120:   PetscDrawDestroy(draw);
121:   PetscFinalize();

123:   return 0;
124: }
125: /* ------------------------------------------------------------------ */
126: /*           Bratu (Solid Fuel Ignition) Test Problem                 */
127: /* ------------------------------------------------------------------ */

129: /* --------------------  Form initial approximation ----------------- */

133: int FormInitialGuess1(AppCtx *user,Vec X)
134: {
135:   int         i, j, row, mx, my, ierr;
136:   PetscReal   lambda, temp1, temp, hx, hy, hxdhy, hydhx,sc;
137:   PetscScalar *x;

139:   mx         = user->mx;
140:   my         = user->my;
141:   lambda = user->param;

143:   hx    = 1.0 / (PetscReal)(mx-1);
144:   hy    = 1.0 / (PetscReal)(my-1);
145:   sc    = hx*hy;
146:   hxdhy = hx/hy;
147:   hydhx = hy/hx;

149:   VecGetArray(X,&x);
150:   temp1 = lambda/(lambda + 1.0);
151:   for (j=0; j<my; j++) {
152:     temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
153:     for (i=0; i<mx; i++) {
154:       row = i + j*mx;
155:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
156:         x[row] = 0.0;
157:         continue;
158:       }
159:       x[row] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
160:     }
161:   }
162:   VecRestoreArray(X,&x);
163:   return 0;
164: }
165: /* --------------------  Evaluate Function F(x) --------------------- */
166: 
169: int FormFunction1(SNES snes,Vec X,Vec F,void *ptr)
170: {
171:   AppCtx       *user = (AppCtx*)ptr;
172:   int          ierr, i, j, row, mx, my;
173:   PetscReal    two = 2.0, one = 1.0, lambda,hx, hy, hxdhy, hydhx;
174:   PetscScalar  ut, ub, ul, ur, u, uxx, uyy, sc,*x,*f;

176:   mx         = user->mx;
177:   my         = user->my;
178:   lambda = user->param;

180:   hx    = one / (PetscReal)(mx-1);
181:   hy    = one / (PetscReal)(my-1);
182:   sc    = hx*hy;
183:   hxdhy = hx/hy;
184:   hydhx = hy/hx;

186:   VecGetArray(X,&x);
187:   VecGetArray(F,&f);
188:   for (j=0; j<my; j++) {
189:     for (i=0; i<mx; i++) {
190:       row = i + j*mx;
191:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
192:         f[row] = x[row];
193:         continue;
194:       }
195:       u = x[row];
196:       ub = x[row - mx];
197:       ul = x[row - 1];
198:       ut = x[row + mx];
199:       ur = x[row + 1];
200:       uxx = (-ur + two*u - ul)*hydhx;
201:       uyy = (-ut + two*u - ub)*hxdhy;
202:       f[row] = uxx + uyy - sc*lambda*exp(u);
203:     }
204:   }
205:   VecRestoreArray(X,&x);
206:   VecRestoreArray(F,&f);
207:   return 0;
208: }
209: /* --------------------  Evaluate Jacobian F'(x) -------------------- */

213: int FormJacobian1(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
214: {
215:   AppCtx       *user = (AppCtx*)ptr;
216:   Mat          jac = *J;
217:   int          i, j, row, mx, my, col[5], ierr;
218:   PetscScalar  two = 2.0, one = 1.0, lambda, v[5],sc, *x;
219:   PetscReal    hx, hy, hxdhy, hydhx;


222:   mx         = user->mx;
223:   my         = user->my;
224:   lambda = user->param;

226:   hx    = 1.0 / (PetscReal)(mx-1);
227:   hy    = 1.0 / (PetscReal)(my-1);
228:   sc    = hx*hy;
229:   hxdhy = hx/hy;
230:   hydhx = hy/hx;

232:   VecGetArray(X,&x);
233:   for (j=0; j<my; j++) {
234:     for (i=0; i<mx; i++) {
235:       row = i + j*mx;
236:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
237:         MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
238:         continue;
239:       }
240:       v[0] = -hxdhy; col[0] = row - mx;
241:       v[1] = -hydhx; col[1] = row - 1;
242:       v[2] = two*(hydhx + hxdhy) - sc*lambda*exp(x[row]); col[2] = row;
243:       v[3] = -hydhx; col[3] = row + 1;
244:       v[4] = -hxdhy; col[4] = row + mx;
245:       MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
246:     }
247:   }
248:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
249:   VecRestoreArray(X,&x);
250:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
251:   *flag = SAME_NONZERO_PATTERN;
252:   return 0;
253: }
254: /* ------------------------------------------------------------------ */
255: /*                       Driven Cavity Test Problem                   */
256: /* ------------------------------------------------------------------ */

258: /* --------------------  Form initial approximation ----------------- */

262: int FormInitialGuess2(AppCtx *user,Vec X)
263: {
264:   int          ierr, i, j, row, mx, my;
265:   PetscScalar  xx,yy,*x;
266:   PetscReal    hx, hy;

268:   mx         = user->mx;
269:   my         = user->my;

271:   /* Test for invalid input parameters */
272:   if ((mx <= 0) || (my <= 0)) SETERRQ(1,0);

274:   hx    = 1.0 / (PetscReal)(mx-1);
275:   hy    = 1.0 / (PetscReal)(my-1);

277:   VecGetArray(X,&x);
278:   yy = 0.0;
279:   for (j=0; j<my; j++) {
280:     xx = 0.0;
281:     for (i=0; i<mx; i++) {
282:       row = i + j*mx;
283:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
284:         x[row] = 0.0;
285:       }
286:       else {
287:         x[row] = - xx*(1.0 - xx)*yy*(1.0 - yy);
288:       }
289:       xx = xx + hx;
290:     }
291:     yy = yy + hy;
292:   }
293:   VecRestoreArray(X,&x);
294:   return 0;
295: }
296: /* --------------------  Evaluate Function F(x) --------------------- */

300: int FormFunction2(SNES snes,Vec X,Vec F,void *pptr)
301: {
302:   AppCtx       *user = (AppCtx*)pptr;
303:   int          i, j, row, mx, my, ierr;
304:   PetscScalar  two = 2.0, zero = 0.0, pb, pbb,pbr, pl,pll,p,pr,prr;
305:   PetscScalar  ptl,pt,ptt,dpdy,dpdx,pblap,ptlap,rey,pbl,ptr,pllap,plap,prlap;
306:   PetscScalar  *x,*f, hx2, hy2, hxhy2;
307:   PetscReal    hx, hy;

309:   mx         = user->mx;
310:   my         = user->my;
311:   hx     = 1.0 / (PetscReal)(mx-1);
312:   hy     = 1.0 / (PetscReal)(my-1);
313:   hx2    = hx*hx;
314:   hy2    = hy*hy;
315:   hxhy2  = hx2*hy2;
316:   rey    = user->param;

318:   VecGetArray(X,&x);
319:   VecGetArray(F,&f);
320:   for (j=0; j<my; j++) {
321:     for (i=0; i<mx; i++) {
322:       row = i + j*mx;
323:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
324:         f[row] = x[row];
325:         continue;
326:       }
327:       if (i == 1 || j == 1) {
328:            pbl = zero;
329:       }
330:       else {
331:            pbl = x[row-mx-1];
332:       }
333:       if (j == 1) {
334:            pb = zero;
335:            pbb = x[row];
336:       }
337:       else if (j == 2) {
338:            pb = x[row-mx];
339:            pbb = zero;
340:       }
341:       else {
342:            pb = x[row-mx];
343:            pbb = x[row-2*mx];
344:       }
345:       if (j == 1 || i == mx-2) {
346:            pbr = zero;
347:       }
348:       else {
349:            pbr = x[row-mx+1];
350:       }
351:       if (i == 1) {
352:            pl = zero;
353:            pll = x[row];
354:       }
355:       else if (i == 2) {
356:            pl = x[row-1];
357:            pll = zero;
358:       }
359:       else {
360:            pl = x[row-1];
361:            pll = x[row-2];
362:       }
363:       p = x[row];
364:       if (i == mx-3) {
365:            pr = x[row+1];
366:            prr = zero;
367:       }
368:       else if (i == mx-2) {
369:            pr = zero;
370:            prr = x[row];
371:       }
372:       else {
373:            pr = x[row+1];
374:            prr = x[row+2];
375:       }
376:       if (j == my-2 || i == 1) {
377:            ptl = zero;
378:       }
379:       else {
380:            ptl = x[row+mx-1];
381:       }
382:       if (j == my-3) {
383:            pt = x[row+mx];
384:            ptt = zero;
385:       }
386:       else if (j == my-2) {
387:            pt = zero;
388:            ptt = x[row] + two*hy;
389:       }
390:       else {
391:            pt = x[row+mx];
392:            ptt = x[row+2*mx];
393:       }
394:       if (j == my-2 || i == mx-2) {
395:            ptr = zero;
396:       }
397:       else {
398:            ptr = x[row+mx+1];
399:       }

401:       dpdy = (pt - pb)/(two*hy);
402:       dpdx = (pr - pl)/(two*hx);

404:       /*  Laplacians at each point in the 5 point stencil */
405:       pblap = (pbr - two*pb + pbl)/hx2 + (p   - two*pb + pbb)/hy2;
406:       pllap = (p   - two*pl + pll)/hx2 + (ptl - two*pl + pbl)/hy2;
407:       plap =  (pr  - two*p  + pl)/hx2 + (pt  - two*p  + pb)/hy2;
408:       prlap = (prr - two*pr + p)/hx2 + (ptr - two*pr + pbr)/hy2;
409:       ptlap = (ptr - two*pt + ptl)/hx2 + (ptt - two*pt + p)/hy2;

411:       f[row] = hxhy2*((prlap - two*plap + pllap)/hx2
412:                         + (ptlap - two*plap + pblap)/hy2
413:                         - rey*(dpdy*(prlap - pllap)/(two*hx)
414:                         - dpdx*(ptlap - pblap)/(two*hy)));
415:     }
416:   }
417:   VecRestoreArray(X,&x);
418:   VecRestoreArray(F,&f);
419:   return 0;
420: }
421: /* --------------------  Evaluate Jacobian F'(x) -------------------- */

425: int FormJacobian2(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *pptr)
426: {
427:   AppCtx       *user = (AppCtx*)pptr;
428:   int          i, j, row, mx, my, col, ierr;
429:   PetscScalar  two = 2.0, one = 1.0, zero = 0.0, pb, pbb,pbr, pl,pll,p,pr,prr;
430:   PetscScalar  ptl,pt,ptt,dpdy,dpdx,pblap,ptlap,rey,pbl,ptr,pllap,plap,prlap;
431:   PetscScalar  val,four = 4.0, three = 3.0,*x;
432:   PetscReal    hx, hy,hx2, hy2, hxhy2;

434:   mx         = user->mx;
435:   my         = user->my;
436:   hx     = 1.0 / (PetscReal)(mx-1);
437:   hy     = 1.0 / (PetscReal)(my-1);
438:   hx2    = hx*hx;
439:   hy2    = hy*hy;
440:   hxhy2  = hx2*hy2;
441:   rey    = user->param;

443:   MatZeroEntries(*J);
444:   VecGetArray(X,&x);
445:   for (j=0; j<my; j++) {
446:     for (i=0; i<mx; i++) {
447:       row = i + j*mx;
448:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
449:         MatSetValues(*J,1,&row,1,&row,&one,ADD_VALUES);
450:         continue;
451:       }
452:       if (i == 1 || j == 1) {
453:            pbl = zero;
454:       }
455:       else {
456:            pbl = x[row-mx-1];
457:       }
458:       if (j == 1) {
459:            pb = zero;
460:            pbb = x[row];
461:       }
462:       else if (j == 2) {
463:            pb = x[row-mx];
464:            pbb = zero;
465:       }
466:       else {
467:            pb = x[row-mx];
468:            pbb = x[row-2*mx];
469:       }
470:       if (j == 1 || i == mx-2) {
471:            pbr = zero;
472:       }
473:       else {
474:            pbr = x[row-mx+1];
475:       }
476:       if (i == 1) {
477:            pl = zero;
478:            pll = x[row];
479:       }
480:       else if (i == 2) {
481:            pl = x[row-1];
482:            pll = zero;
483:       }
484:       else {
485:            pl = x[row-1];
486:            pll = x[row-2];
487:       }
488:       p = x[row];
489:       if (i == mx-3) {
490:            pr = x[row+1];
491:            prr = zero;
492:       }
493:       else if (i == mx-2) {
494:            pr = zero;
495:            prr = x[row];
496:       }
497:       else {
498:            pr = x[row+1];
499:            prr = x[row+2];
500:       }
501:       if (j == my-2 || i == 1) {
502:            ptl = zero;
503:       }
504:       else {
505:            ptl = x[row+mx-1];
506:       }
507:       if (j == my-3) {
508:            pt = x[row+mx];
509:            ptt = zero;
510:       }
511:       else if (j == my-2) {
512:            pt = zero;
513:            ptt = x[row] + two*hy;
514:       }
515:       else {
516:            pt = x[row+mx];
517:            ptt = x[row+2*mx];
518:       }
519:       if (j == my-2 || i == mx-2) {
520:            ptr = zero;
521:       }
522:       else {
523:            ptr = x[row+mx+1];
524:       }

526:       dpdy = (pt - pb)/(two*hy);
527:       dpdx = (pr - pl)/(two*hx);

529:       /*  Laplacians at each point in the 5 point stencil */
530:       pblap = (pbr - two*pb + pbl)/hx2 + (p   - two*pb + pbb)/hy2;
531:       pllap = (p   - two*pl + pll)/hx2 + (ptl - two*pl + pbl)/hy2;
532:       plap =  (pr  - two*p  + pl)/hx2 + (pt  - two*p  + pb)/hy2;
533:       prlap = (prr - two*pr + p)/hx2 + (ptr - two*pr + pbr)/hy2;
534:       ptlap = (ptr - two*pt + ptl)/hx2 + (ptt - two*pt + p)/hy2;

536:       if (j > 2) {
537:         val = hxhy2*(one/hy2/hy2 - rey*dpdx/hy2/(two*hy));
538:         col = row - 2*mx;
539:         MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
540:       }
541:       if (i > 2) {
542:         val = hxhy2*(one/hx2/hx2 + rey*dpdy/hx2/(two*hx));
543:         col = row - 2;
544:         MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
545:       }
546:       if (i < mx-3) {
547:         val = hxhy2*(one/hx2/hx2 - rey*dpdy/hx2/(two*hx));
548:         col = row + 2;
549:         MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
550:       }
551:       if (j < my-3) {
552:         val = hxhy2*(one/hy2/hy2 + rey*dpdx/hy2/(two*hy));
553:         col = row + 2*mx;
554:         MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
555:       }
556:       if (i != 1 && j != 1) {
557:         val = hxhy2*(two/hy2/hx2 + rey*(dpdy/hy2/(two*hx) - dpdx/hx2/(two*hy)));
558:         col = row - mx - 1;
559:         MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
560:       }
561:       if (j != 1 && i != mx-2) {
562:         val = hxhy2*(two/hy2/hx2 - rey*(dpdy/hy2/(two*hx) + dpdx/hx2/(two*hy)));
563:         col = row - mx + 1;
564:         MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
565:       }
566:       if (j != my-2 && i != 1) {
567:         val = hxhy2*(two/hy2/hx2 + rey*(dpdy/hy2/(two*hx) + dpdx/hx2/(two*hy)));
568:         col = row + mx - 1;
569:         MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
570:       }
571:       if (j != my-2 && i != mx-2) {
572:         val = hxhy2*(two/hy2/hx2 - rey*(dpdy/hy2/(two*hx) - dpdx/hx2/(two*hy)));
573:         col = row + mx + 1;
574:         MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
575:       }
576:       if (j != 1) {
577:         val = hxhy2*(-four*(one/hy2/hx2 + one/hy2/hy2)
578:                      + rey*((prlap - pllap)/(two*hx)/(two*hy)
579:                      + dpdx*(one/hx2 + one/hy2)/hy));
580:         col = row - mx;
581:         MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
582:       }
583:       if (i != 1) {
584:         val = hxhy2*(-four*(one/hy2/hx2 + one/hx2/hx2)
585:                      - rey*((ptlap - pblap)/(two*hx)/(two*hy)
586:                      + dpdy*(one/hx2 + one/hy2)/hx));
587:         col = row - 1;
588:         MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
589:       }
590:       if (i != mx-2) {
591:         val = hxhy2*(-four*(one/hy2/hx2 + one/hx2/hx2)
592:                      + rey*((ptlap - pblap)/(two*hx)/(two*hy)
593:                      + dpdy*(one/hx2 + one/hy2)/hx));
594:         col = row + 1;
595:         MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
596:       }
597:       if (j != my-2) {
598:         val = hxhy2*(-four*(one/hy2/hx2 + one/hy2/hy2)
599:                      - rey*((prlap - pllap)/(two*hx)/(two*hy)
600:                      + dpdx*(one/hx2 + one/hy2)/hy));
601:         col = row + mx;
602:         MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
603:       }
604:       val = hxhy2*(two*(four/hx2/hy2 + three/hx2/hx2 + three/hy2/hy2));
605:       MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
606:       if (j == 1) {
607:         val = hxhy2*(one/hy2/hy2 - rey*(dpdx/hy2/(two*hy)));
608:         MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
609:       }
610:       if (i == 1) {
611:         val = hxhy2*(one/hx2/hx2 + rey*(dpdy/hx2/(two*hx)));
612:         MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
613:       }
614:       if (i == mx-2) {
615:         val = hxhy2*(one/hx2/hx2 - rey*(dpdy/hx2/(two*hx)));
616:         MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
617:       }
618:       if (j == my-2) {
619:         val = hxhy2*(one/hy2/hy2 + rey*(dpdx/hy2/(two*hy)));
620:         MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
621:       }
622:     }
623:   }
624:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
625:   VecRestoreArray(X,&x);
626:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
627:   *flag = SAME_NONZERO_PATTERN;
628:   return 0;
629: }