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,&its);
108:   SNESGetNumberUnsuccessfulSteps(snes,&nfails);

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

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

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

128: /* --------------------  Form initial approximation ----------------- */

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

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

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

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

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

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

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

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


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

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

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

257: /* --------------------  Form initial approximation ----------------- */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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