Actual source code: ex4.c

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

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

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

 36: #include <petscsnes.h>

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

 45:                       FormFunction1(SNES,Vec,Vec,void*),
 46:                       FormInitialGuess1(AppCtx*,Vec);
 48:                       FormFunction2(SNES,Vec,Vec,void*),
 49:                       FormInitialGuess2(AppCtx*,Vec);
 50: 
 53: int main(int argc, char **argv)
 54: {
 55:   SNES           snes;                 /* SNES context */
 56:   const SNESType method = SNESLS;      /* default nonlinear solution method */
 57:   Vec            x, r;                 /* solution, residual vectors */
 58:   Mat            J;                    /* Jacobian matrix */
 59:   AppCtx         user;                 /* user-defined application context */
 60:   PetscDraw      draw;                 /* drawing context */
 61:   PetscInt       its, N, nfails;
 63:   PetscBool      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)) SETERRQ(PETSC_COMM_SELF,1,"Lambda is out of range");
 81:   N = user.mx*user.my;
 82: 
 83:   /* Set up data structures */
 84:   VecCreateSeq(PETSC_COMM_SELF,N,&x);
 85:   VecDuplicate(x,&r);
 86:   MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,5,PETSC_NULL,&J);

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

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

103:   /* Set options and solve nonlinear system */
104:   SNESSetFromOptions(snes);
105:   SNESSolve(snes,PETSC_NULL,x);
106:   SNESGetIterationNumber(snes,&its);
107:   SNESGetNonlinearStepFailures(snes,&nfails);

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

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

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

127: /* --------------------  Form initial approximation ----------------- */

131: PetscErrorCode  FormInitialGuess1(AppCtx *user,Vec X)
132: {
133:   PetscInt       i, j, row, mx, my;
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*PetscSqrtReal(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: PetscErrorCode  FormFunction1(SNES snes,Vec X,Vec F,void *ptr)
169: {
170:   AppCtx         *user = (AppCtx*)ptr;
171:   PetscInt       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: PetscErrorCode  FormJacobian1(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
214: {
215:   AppCtx         *user = (AppCtx*)ptr;
216:   Mat            jac = *J;
217:   PetscInt       i, j, row, mx, my, col[5];
219:   PetscScalar    two = 2.0, one = 1.0, lambda, v[5],sc, *x;
220:   PetscReal      hx, hy, hxdhy, hydhx;


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

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

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

259: /* --------------------  Form initial approximation ----------------- */

263: PetscErrorCode  FormInitialGuess2(AppCtx *user,Vec X)
264: {
266:   PetscInt       i, j, row, mx, my;
267:   PetscScalar    xx,yy,*x;
268:   PetscReal      hx, hy;

270:   mx         = user->mx;
271:   my         = user->my;

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

276:   hx    = 1.0 / (PetscReal)(mx-1);
277:   hy    = 1.0 / (PetscReal)(my-1);

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

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

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

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

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

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

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

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

439:   mx         = user->mx;
440:   my         = user->my;
441:   hx     = 1.0 / (PetscReal)(mx-1);
442:   hy     = 1.0 / (PetscReal)(my-1);
443:   hx2    = hx*hx;
444:   hy2    = hy*hy;
445:   hxhy2  = hx2*hy2;
446:   rey    = user->param;

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

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

537:       /*  Laplacians at each point in the 5 point stencil */
538:       pblap = (pbr - two*pb + pbl)/hx2 + (p   - two*pb + pbb)/hy2;
539:       pllap = (p   - two*pl + pll)/hx2 + (ptl - two*pl + pbl)/hy2;
540:       plap =  (pr  - two*p  + pl)/hx2 + (pt  - two*p  + pb)/hy2;
541:       prlap = (prr - two*pr + p)/hx2 + (ptr - two*pr + pbr)/hy2;
542:       ptlap = (ptr - two*pt + ptl)/hx2 + (ptt - two*pt + p)/hy2;

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