Actual source code: ex5.c

petsc-3.7.6 2017-04-24
Report Typos and Errors
  2: static char help[] = "Bratu nonlinear PDE in 2d.\n\
3: We solve the  Bratu (SFI - solid fuel ignition) problem in a 2D rectangular\n\
4: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\
5: The command line options include:\n\
6:   -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
7:      problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)\n\n\
8:   -m_par/n_par <parameter>, where <parameter> indicates an integer\n \
9:       that MMS3 will be evaluated with 2^m_par, 2^n_par";

11: /*T
12:    Concepts: SNES^parallel Bratu example
13:    Concepts: DMDA^using distributed arrays;
14:    Concepts: IS coloirng types;
15:    Processors: n
16: T*/

18: /* ------------------------------------------------------------------------

20:     Solid Fuel Ignition (SFI) problem.  This problem is modeled by
21:     the partial differential equation

23:             -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,

25:     with boundary conditions

27:              u = 0  for  x = 0, x = 1, y = 0, y = 1.

29:     A finite difference approximation with the usual 5-point stencil
30:     is used to discretize the boundary value problem to obtain a nonlinear
31:     system of equations.

34:       This example shows how geometric multigrid can be run transparently with a nonlinear solver so long
35:       as SNESSetDM() is provided. Example usage

37:       ./ex5 -pc_type mg -ksp_monitor  -snes_view -pc_mg_levels 3 -pc_mg_galerkin -da_grid_x 17 -da_grid_y 17
38:              -mg_levels_ksp_monitor -snes_monitor -mg_levels_pc_type sor -pc_mg_type full

40:       or to run with grid sequencing on the nonlinear problem (note that you do not need to provide the number of
41:          multigrid levels, it will be determined automatically based on the number of refinements done)

43:       ./ex5 -pc_type mg -ksp_monitor  -snes_view -pc_mg_galerkin -snes_grid_sequence 3
44:              -mg_levels_ksp_monitor -snes_monitor -mg_levels_pc_type sor -pc_mg_type full

46:   ------------------------------------------------------------------------- */

48: /*
49:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
50:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
51: */
52: #include <petscdm.h>
53: #include <petscdmda.h>
54: #include <petscsnes.h>
55: #include <petscmatlab.h>

57: /*
58:    User-defined application context - contains data needed by the
59:    application-provided call-back routines, FormJacobianLocal() and
60:    FormFunctionLocal().
61: */
62: typedef struct {
63:   PetscReal param;          /* test problem parameter */
64:   PetscInt  mPar;           /* MMS3 m parameter */
65:   PetscInt  nPar;           /* MMS3 n parameter */
66: } AppCtx;

68: /*
69:    User-defined routines
70: */
71: extern PetscErrorCode FormInitialGuess(DM,AppCtx*,Vec);
72: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
73: extern PetscErrorCode FormExactSolution1(DM,AppCtx*,Vec);
74: extern PetscErrorCode FormFunctionLocalMMS1(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
75: extern PetscErrorCode FormExactSolution2(DM,AppCtx*,Vec);
76: extern PetscErrorCode FormFunctionLocalMMS2(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
77: extern PetscErrorCode FormExactSolution3(DM,AppCtx*,Vec);
78: extern PetscErrorCode FormFunctionLocalMMS3(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
79: extern PetscErrorCode FormExactSolution4(DM,AppCtx*,Vec);
80: extern PetscErrorCode FormFunctionLocalMMS4(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
81: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,Mat,AppCtx*);
82: extern PetscErrorCode FormObjectiveLocal(DMDALocalInfo*,PetscScalar**,PetscReal*,AppCtx*);
83: #if defined(PETSC_HAVE_MATLAB_ENGINE)
84: extern PetscErrorCode FormFunctionMatlab(SNES,Vec,Vec,void*);
85: #endif
86: extern PetscErrorCode NonlinearGS(SNES,Vec,Vec,void*);

90: int main(int argc,char **argv)
91: {
92:   SNES           snes;                         /* nonlinear solver */
93:   Vec            x;                            /* solution vector */
94:   AppCtx         user;                         /* user-defined work context */
95:   PetscInt       its;                          /* iterations for convergence */
97:   PetscReal      bratu_lambda_max = 6.81;
98:   PetscReal      bratu_lambda_min = 0.;
99:   PetscInt       MMS              = 0;
100:   PetscBool      flg              = PETSC_FALSE;
101:   DM             da;
102: #if defined(PETSC_HAVE_MATLAB_ENGINE)
103:   Vec            r               = NULL;
104:   PetscBool      matlab_function = PETSC_FALSE;
105: #endif
106:   KSP            ksp;
107:   PetscInt       lits,slits;

109:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110:      Initialize program
111:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

113:   PetscInitialize(&argc,&argv,(char*)0,help);

115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116:      Initialize problem parameters
117:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118:   user.param = 6.0;
119:   PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
120:   if (user.param > bratu_lambda_max || user.param < bratu_lambda_min) SETERRQ3(PETSC_COMM_SELF,1,"Lambda, %g, is out of range, [%g, %g]", user.param, bratu_lambda_min, bratu_lambda_max);
121:   PetscOptionsGetInt(NULL,NULL,"-mms",&MMS,NULL);
122:   if (MMS == 3) {
123:     user.mPar = 2;
124:     user.nPar = 1;
125:     PetscOptionsGetInt(NULL,NULL,"-m_par",&user.mPar,NULL);
126:     PetscOptionsGetInt(NULL,NULL,"-n_par",&user.nPar,NULL);
127:   }

129:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130:      Create nonlinear solver context
131:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132:   SNESCreate(PETSC_COMM_WORLD,&snes);
133:   SNESSetCountersReset(snes,PETSC_FALSE);
134:   SNESSetNGS(snes, NonlinearGS, NULL);

136:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137:      Create distributed array (DMDA) to manage parallel grid and vectors
138:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
140:   DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
141:   DMSetApplicationContext(da,&user);
142:   SNESSetDM(snes,da);
143:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144:      Extract global vectors from DMDA; then duplicate for remaining
145:      vectors that are the same types
146:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147:   DMCreateGlobalVector(da,&x);

149:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150:      Set local function evaluation routine
151:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152:   switch (MMS) {
153:   case 1:  DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocalMMS1,&user);break;
154:   case 2:  DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocalMMS2,&user);break;
155:   case 3:  DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocalMMS3,&user);break;
156:   case 4:  DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocalMMS4,&user);break;
157:   default: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocal,&user);
158:   }
159:   PetscOptionsGetBool(NULL,NULL,"-fd",&flg,NULL);
160:   if (!flg) {
161:     DMDASNESSetJacobianLocal(da,(DMDASNESJacobian)FormJacobianLocal,&user);
162:   }

164:   PetscOptionsGetBool(NULL,NULL,"-obj",&flg,NULL);
165:   if (flg) {
166:     DMDASNESSetObjectiveLocal(da,(DMDASNESObjective)FormObjectiveLocal,&user);
167:   }

169: #if defined(PETSC_HAVE_MATLAB_ENGINE)
170:   PetscOptionsGetBool(NULL,NULL,"-matlab_function",&matlab_function,0);
171:   if (matlab_function) {
172:     VecDuplicate(x,&r);
173:     SNESSetFunction(snes,r,FormFunctionMatlab,&user);
174:   }
175: #endif

177:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178:      Customize nonlinear solver; set runtime options
179:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180:   SNESSetFromOptions(snes);

182:   FormInitialGuess(da,&user,x);

184:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185:      Solve nonlinear system
186:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
187:   SNESSolve(snes,NULL,x);
188:   SNESGetIterationNumber(snes,&its);

190:   SNESGetLinearSolveIterations(snes,&slits);
191:   SNESGetKSP(snes,&ksp);
192:   KSPGetTotalIterations(ksp,&lits);
193:   if (lits != slits) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Number of total linear iterations reported by SNES %D does not match reported by KSP %D",slits,lits);
194:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

197:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198:      If using MMS, check the l_2 error
199:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
200:   if (MMS) {
201:     Vec       e;
202:     PetscReal errorl2, errorinf;
203:     PetscInt  N;

205:     VecDuplicate(x, &e);
206:     PetscObjectViewFromOptions((PetscObject) x, NULL, "-sol_view");
207:     if (MMS == 1)      {FormExactSolution1(da, &user, e);}
208:     else if (MMS == 2) {FormExactSolution2(da, &user, e);}
209:     else if (MMS == 3) {FormExactSolution3(da, &user, e);}
210:     else               {FormExactSolution4(da, &user, e);}
211:     PetscObjectViewFromOptions((PetscObject) e, NULL, "-exact_view");
212:     VecAXPY(e, -1.0, x);
213:     PetscObjectViewFromOptions((PetscObject) e, NULL, "-error_view");
214:     VecNorm(e, NORM_2, &errorl2);
215:     VecNorm(e, NORM_INFINITY, &errorinf);
216:     VecGetSize(e, &N);
217:     PetscPrintf(PETSC_COMM_WORLD, "N: %D error l2 %g inf %g\n", N, (double) errorl2/N, (double) errorinf);
218:     VecDestroy(&e);
219:   }

221:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
222:      Free work space.  All PETSc objects should be destroyed when they
223:      are no longer needed.
224:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
225: #if defined(PETSC_HAVE_MATLAB_ENGINE)
226:   VecDestroy(&r);
227: #endif
228:   VecDestroy(&x);
229:   SNESDestroy(&snes);
230:   DMDestroy(&da);
231:   PetscFinalize();
232:   return 0;
233: }
234: /* ------------------------------------------------------------------- */
237: /*
238:    FormInitialGuess - Forms initial approximation.

240:    Input Parameters:
241:    da - The DM
242:    user - user-defined application context

244:    Output Parameter:
245:    X - vector
246:  */
247: PetscErrorCode FormInitialGuess(DM da,AppCtx *user,Vec X)
248: {
249:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
251:   PetscReal      lambda,temp1,temp,hx,hy;
252:   PetscScalar    **x;

255:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

257:   lambda = user->param;
258:   hx     = 1.0/(PetscReal)(Mx-1);
259:   hy     = 1.0/(PetscReal)(My-1);
260:   temp1  = lambda/(lambda + 1.0);

262:   /*
263:      Get a pointer to vector data.
264:        - For default PETSc vectors, VecGetArray() returns a pointer to
265:          the data array.  Otherwise, the routine is implementation dependent.
266:        - You MUST call VecRestoreArray() when you no longer need access to
267:          the array.
268:   */
269:   DMDAVecGetArray(da,X,&x);

271:   /*
272:      Get local grid boundaries (for 2-dimensional DMDA):
273:        xs, ys   - starting grid indices (no ghost points)
274:        xm, ym   - widths of local grid (no ghost points)

276:   */
277:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

279:   /*
280:      Compute initial guess over the locally owned part of the grid
281:   */
282:   for (j=ys; j<ys+ym; j++) {
283:     temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
284:     for (i=xs; i<xs+xm; i++) {
285:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
286:         /* boundary conditions are all zero Dirichlet */
287:         x[j][i] = 0.0;
288:       } else {
289:         x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
290:       }
291:     }
292:   }

294:   /*
295:      Restore vector
296:   */
297:   DMDAVecRestoreArray(da,X,&x);
298:   return(0);
299: }

303: /*
304:   FormExactSolution1 - Forms initial approximation.

306:   Input Parameters:
307:   da - The DM
308:   user - user-defined application context

310:   Output Parameter:
311:   X - vector
312:  */
313: PetscErrorCode FormExactSolution1(DM da, AppCtx *user, Vec U)
314: {
315:   DM             coordDA;
316:   Vec            coordinates;
317:   DMDACoor2d   **coords;
318:   PetscScalar  **u;
319:   PetscReal      x, y;
320:   PetscInt       xs, ys, xm, ym, i, j;

324:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
325:   DMGetCoordinateDM(da, &coordDA);
326:   DMGetCoordinates(da, &coordinates);
327:   DMDAVecGetArray(coordDA, coordinates, &coords);
328:   DMDAVecGetArray(da, U, &u);
329:   for (j = ys; j < ys+ym; ++j) {
330:     for (i = xs; i < xs+xm; ++i) {
331:       x = PetscRealPart(coords[j][i].x);
332:       y = PetscRealPart(coords[j][i].y);
333:       u[j][i] = x*(1 - x)*y*(1 - y);
334:     }
335:   }
336:   DMDAVecRestoreArray(da, U, &u);
337:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
338:   return(0);
339: }

343: /*
344:   FormExactSolution2 - Forms initial approximation.

346:   Input Parameters:
347:   da - The DM
348:   user - user-defined application context

350:   Output Parameter:
351:   X - vector
352:  */
353: PetscErrorCode FormExactSolution2(DM da, AppCtx *user, Vec U)
354: {
355:   DM             coordDA;
356:   Vec            coordinates;
357:   DMDACoor2d   **coords;
358:   PetscScalar  **u;
359:   PetscReal      x, y;
360:   PetscInt       xs, ys, xm, ym, i, j;

364:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
365:   DMGetCoordinateDM(da, &coordDA);
366:   DMGetCoordinates(da, &coordinates);
367:   DMDAVecGetArray(coordDA, coordinates, &coords);
368:   DMDAVecGetArray(da, U, &u);
369:   for (j = ys; j < ys+ym; ++j) {
370:     for (i = xs; i < xs+xm; ++i) {
371:       x = PetscRealPart(coords[j][i].x);
372:       y = PetscRealPart(coords[j][i].y);
373:       u[j][i] = PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y);
374:     }
375:   }
376:   DMDAVecRestoreArray(da, U, &u);
377:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
378:   return(0);
379: }

383: /*
384:   FormExactSolution3 - Forms initial approximation.

386:   Input Parameters:
387:   da - The DM
388:   user - user-defined application context

390:   Output Parameter:
391:   X - vector
392:  */
393: PetscErrorCode FormExactSolution3(DM da, AppCtx *user, Vec U)
394: {
395:   DM             coordDA;
396:   Vec            coordinates;
397:   DMDACoor2d   **coords;
398:   PetscScalar  **u;
399:   PetscReal      x, y;
400:   PetscInt       xs, ys, xm, ym, i, j, m, n;

403:   m = PetscPowReal(2,user->mPar);
404:   n = PetscPowReal(2,user->nPar);

407:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
408:   DMGetCoordinateDM(da, &coordDA);
409:   DMGetCoordinates(da, &coordinates);
410:   DMDAVecGetArray(coordDA, coordinates, &coords);
411:   DMDAVecGetArray(da, U, &u);

413:   for (j = ys; j < ys+ym; ++j) {
414:     for (i = xs; i < xs+xm; ++i) {
415:       x = PetscRealPart(coords[j][i].x);
416:       y = PetscRealPart(coords[j][i].y);

418:       u[j][i] = PetscSinReal(m*PETSC_PI*x*(1-y))*PetscSinReal(n*PETSC_PI*y*(1-x));
419:     }
420:   }
421:   DMDAVecRestoreArray(da, U, &u);
422:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
423:   return(0);
424: }
425: /* ------------------------------------------------------------------- */

429: /*
430:   FormExactSolution4 - Forms initial approximation.

432:   Input Parameters:
433:   da - The DM
434:   user - user-defined application context

436:   Output Parameter:
437:   X - vector
438:  */
439: PetscErrorCode FormExactSolution4(DM da, AppCtx *user, Vec U)
440: {
441:   DM             coordDA;
442:   Vec            coordinates;
443:   DMDACoor2d   **coords;
444:   PetscScalar  **u;
445:   PetscReal      x, y, Lx, Ly;
446:   PetscInt       xs, ys, xm, ym, i, j;

450:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
451:   DMGetCoordinateDM(da, &coordDA);
452:   DMGetCoordinates(da, &coordinates);
453:   DMDAVecGetArray(coordDA, coordinates, &coords);
454:   DMDAVecGetArray(da, U, &u);

456:   Lx = PetscRealPart(coords[ys][xs+xm-1].x - coords[ys][xs].x);
457:   Ly = PetscRealPart(coords[ys+ym-1][xs].y - coords[ys][xs].y);

459:   for (j = ys; j < ys+ym; ++j) {
460:     for (i = xs; i < xs+xm; ++i) {
461:       x = PetscRealPart(coords[j][i].x);
462:       y = PetscRealPart(coords[j][i].y);
463:       u[j][i] = (PetscPowReal(x,4)-PetscSqr(Lx)*PetscSqr(x))*(PetscPowReal(y,4)-PetscSqr(Ly)*PetscSqr(y));
464:     }
465:   }
466:   DMDAVecRestoreArray(da, U, &u);
467:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
468:   return(0);
469: }
470: /* ------------------------------------------------------------------- */
473: /*
474:    FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch

477:  */
478: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
479: {
481:   PetscInt       i,j;
482:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc;
483:   PetscScalar    u,ue,uw,un,us,uxx,uyy;

486:   lambda = user->param;
487:   hx     = 1.0/(PetscReal)(info->mx-1);
488:   hy     = 1.0/(PetscReal)(info->my-1);
489:   sc     = hx*hy*lambda;
490:   hxdhy  = hx/hy;
491:   hydhx  = hy/hx;
492:   /*
493:      Compute function over the locally owned part of the grid
494:   */
495:   for (j=info->ys; j<info->ys+info->ym; j++) {
496:     for (i=info->xs; i<info->xs+info->xm; i++) {
497:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
498:         f[j][i] = 2.0*(hydhx+hxdhy)*x[j][i];
499:       } else {
500:         u  = x[j][i];
501:         uw = x[j][i-1];
502:         ue = x[j][i+1];
503:         un = x[j-1][i];
504:         us = x[j+1][i];

506:         if (i-1 == 0) uw = 0.;
507:         if (i+1 == info->mx-1) ue = 0.;
508:         if (j-1 == 0) un = 0.;
509:         if (j+1 == info->my-1) us = 0.;

511:         uxx     = (2.0*u - uw - ue)*hydhx;
512:         uyy     = (2.0*u - un - us)*hxdhy;
513:         f[j][i] = uxx + uyy - sc*PetscExpScalar(u);
514:       }
515:     }
516:   }
517:   PetscLogFlops(11.0*info->ym*info->xm);
518:   return(0);
519: }

523: /* ---------------------------------------------------------------------------------
524:  FormFunctionLocalMMS1 - Evaluates nonlinear function, F(x) on local process patch

526:  u(x,y) = x(1-x)y(1-y)

528:  -laplacian u* - lambda exp(u*) = 2x(1-x) + 2y(1-y) - lambda exp(x(1-x)y(1-y))
529:
530:  Remark: the above is subtracted from the residual
531: -----------------------------------------------------------------------------------*/
532: PetscErrorCode FormFunctionLocalMMS1(DMDALocalInfo *info,PetscScalar **vx,PetscScalar **f,AppCtx *user)
533: {
535:   PetscInt       i,j;
536:   PetscReal      lambda,hx,hy,hxdhy,hydhx;
537:   PetscScalar    u,ue,uw,un,us,uxx,uyy;
538:   PetscReal      x,y;
539:   DM             coordDA;
540:   Vec            coordinates;
541:   DMDACoor2d   **coords;
542:   Vec            bcv = NULL;
543:   PetscScalar  **bcx = NULL;

546:   lambda = user->param;
547:   /* Extract coordinates */
548:   DMGetCoordinateDM(info->da, &coordDA);
549:   DMGetCoordinates(info->da, &coordinates);
550:   DMDAVecGetArray(coordDA, coordinates, &coords);
551:   hx     = info->xm > 1 ? PetscRealPart(coords[info->ys][info->xs+1].x) - PetscRealPart(coords[info->ys][info->xs].x) : 1.0;
552:   hy     = info->ym > 1 ? PetscRealPart(coords[info->ys+1][info->xs].y) - PetscRealPart(coords[info->ys][info->xs].y) : 1.0;
553:   hxdhy  = hx/hy;
554:   hydhx  = hy/hx;
555:   DMGetNamedLocalVector(info->da, "_petsc_boundary_conditions_", &bcv);
556:   DMDAVecGetArray(info->da, bcv, &bcx);
557:   /* Compute function over the locally owned part of the grid */
558:   for (j=info->ys; j<info->ys+info->ym; j++) {
559:     for (i=info->xs; i<info->xs+info->xm; i++) {
560:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
561:         f[j][i] = 2.0*(hydhx+hxdhy)*(vx[j][i] - bcx[j][i]);
562:       } else {
563:         x  = PetscRealPart(coords[j][i].x);
564:         y  = PetscRealPart(coords[j][i].y);
565:         u  = vx[j][i];
566:         uw = vx[j][i-1];
567:         ue = vx[j][i+1];
568:         un = vx[j-1][i];
569:         us = vx[j+1][i];

571:         if (i-1 == 0)          uw = bcx[j][i-1];
572:         if (i+1 == info->mx-1) ue = bcx[j][i+1];
573:         if (j-1 == 0)          un = bcx[j-1][i];
574:         if (j+1 == info->my-1) us = bcx[j+1][i];

576:         uxx     = (2.0*u - uw - ue)*hydhx;
577:         uyy     = (2.0*u - un - us)*hxdhy;
578:         f[j][i] = uxx + uyy - hx*hy*(lambda*PetscExpScalar(u) + 2*x*(1 - x) + 2*y*(1 - y) - lambda*exp(x*(1 - x)*y*(1 - y)));
579:       }
580:     }
581:   }
582:   DMDAVecRestoreArray(info->da, bcv, &bcx);
583:   DMRestoreNamedLocalVector(info->da, "_petsc_boundary_conditions_", &bcv);
584:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
585:   PetscLogFlops(11.0*info->ym*info->xm);
586:   return(0);
587: }

591: /* ---------------------------------------------------------------------------------
592:  FormFunctionLocalMMS2 - Evaluates nonlinear function, F(x) on local process patch

594:  u(x,y) = sin(pi x)sin(pi y)

596:  -laplacian u* - lambda exp(u*) = 2 pi^2 sin(pi x) sin(pi y) - lambda exp(sin(pi x)sin(pi y))

598:  Remark: the above is subtracted from the residual
599: -----------------------------------------------------------------------------------*/
600: PetscErrorCode FormFunctionLocalMMS2(DMDALocalInfo *info,PetscScalar **vx,PetscScalar **f,AppCtx *user)
601: {
603:   PetscInt       i,j;
604:   PetscReal      lambda,hx,hy,hxdhy,hydhx;
605:   PetscScalar    u,ue,uw,un,us,uxx,uyy;
606:   PetscReal      x,y;
607:   DM             coordDA;
608:   Vec            coordinates;
609:   DMDACoor2d   **coords;
610:   Vec            bcv = NULL;
611:   PetscScalar  **bcx = NULL;

614:   lambda = user->param;
615:   /* Extract coordinates */
616:   DMGetCoordinateDM(info->da, &coordDA);
617:   DMGetCoordinates(info->da, &coordinates);
618:   DMDAVecGetArray(coordDA, coordinates, &coords);
619:   hx     = info->xm > 1 ? PetscRealPart(coords[info->ys][info->xs+1].x) - PetscRealPart(coords[info->ys][info->xs].x) : 1.0;
620:   hy     = info->ym > 1 ? PetscRealPart(coords[info->ys+1][info->xs].y) - PetscRealPart(coords[info->ys][info->xs].y) : 1.0;
621:   hxdhy  = hx/hy;
622:   hydhx  = hy/hx;
623:   DMGetNamedLocalVector(info->da, "_petsc_boundary_conditions_", &bcv);
624:   DMDAVecGetArray(info->da, bcv, &bcx);
625:   /* Compute function over the locally owned part of the grid */
626:   for (j=info->ys; j<info->ys+info->ym; j++) {
627:     for (i=info->xs; i<info->xs+info->xm; i++) {
628:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
629:         f[j][i] = 2.0*(hydhx+hxdhy)*(vx[j][i] - bcx[j][i]);
630:       } else {
631:         x  = PetscRealPart(coords[j][i].x);
632:         y  = PetscRealPart(coords[j][i].y);
633:         u  = vx[j][i];
634:         uw = vx[j][i-1];
635:         ue = vx[j][i+1];
636:         un = vx[j-1][i];
637:         us = vx[j+1][i];

639:         if (i-1 == 0)          uw = bcx[j][i-1];
640:         if (i+1 == info->mx-1) ue = bcx[j][i+1];
641:         if (j-1 == 0)          un = bcx[j-1][i];
642:         if (j+1 == info->my-1) us = bcx[j+1][i];

644:         uxx     = (2.0*u - uw - ue)*hydhx;
645:         uyy     = (2.0*u - un - us)*hxdhy;
646:         f[j][i] = uxx + uyy - hx*hy*(lambda*PetscExpScalar(u) + 2*PetscSqr(PETSC_PI)*PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y) - lambda*exp(PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y)));
647:       }
648:     }
649:   }
650:   DMDAVecRestoreArray(info->da, bcv, &bcx);
651:   DMRestoreNamedLocalVector(info->da, "_petsc_boundary_conditions_", &bcv);
652:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
653:   PetscLogFlops(11.0*info->ym*info->xm);
654:   return(0);
655: }

659: /* ---------------------------------------------------------------------------------
660:  FormFunctionLocalMMS3 - Evaluates nonlinear function, F(x) on local process patch

662:  u(x,y) = sin(m pi x(1-y))sin(n pi y(1-x))

664:  -laplacian u* - lambda exp(u*) = -(exp(Sin[m*Pi*x*(1 - y)]*Sin[n*Pi*(1 - x)*y])*lambda) +
665:                                   Pi^2*(-2*m*n*((-1 + x)*x + (-1 + y)*y)*Cos[m*Pi*x*(-1 + y)]*
666:                                   Cos[n*Pi*(-1 + x)*y] + (m^2*(x^2 + (-1 + y)^2) + n^2*((-1 + x)^2 + y^2))*
667:                                   Sin[m*Pi*x*(-1 + y)]*Sin[n*Pi*(-1 + x)*y])

669:   Remark: the above is subtracted from the residual
670: -----------------------------------------------------------------------------------*/
671: PetscErrorCode FormFunctionLocalMMS3(DMDALocalInfo *info,PetscScalar **vx,PetscScalar **f,AppCtx *user)
672: {
674:   PetscInt       i,j;
675:   PetscReal      lambda,hx,hy,hxdhy,hydhx,m,n;
676:   PetscScalar    u,ue,uw,un,us,uxx,uyy;
677:   PetscReal      x,y;
678:   DM             coordDA;
679:   Vec            coordinates;
680:   DMDACoor2d   **coords;

682:   m = PetscPowReal(2,user->mPar);
683:   n = PetscPowReal(2,user->nPar);

686:   lambda = user->param;
687:   hx     = 1.0/(PetscReal)(info->mx-1);
688:   hy     = 1.0/(PetscReal)(info->my-1);
689:   hxdhy  = hx/hy;
690:   hydhx  = hy/hx;
691:   /* Extract coordinates */
692:   DMGetCoordinateDM(info->da, &coordDA);
693:   DMGetCoordinates(info->da, &coordinates);
694:   DMDAVecGetArray(coordDA, coordinates, &coords);
695:   /* Compute function over the locally owned part of the grid */
696:   for (j=info->ys; j<info->ys+info->ym; j++) {
697:     for (i=info->xs; i<info->xs+info->xm; i++) {
698:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
699:         f[j][i] = 2.0*(hydhx+hxdhy)*vx[j][i];
700:       } else {
701:         x  = PetscRealPart(coords[j][i].x);
702:         y  = PetscRealPart(coords[j][i].y);
703:         u  = vx[j][i];
704:         uw = vx[j][i-1];
705:         ue = vx[j][i+1];
706:         un = vx[j-1][i];
707:         us = vx[j+1][i];

709:         if (i-1 == 0) uw = 0.;
710:         if (i+1 == info->mx-1) ue = 0.;
711:         if (j-1 == 0) un = 0.;
712:         if (j+1 == info->my-1) us = 0.;

714:         uxx     = (2.0*u - uw - ue)*hydhx;
715:         uyy     = (2.0*u - un - us)*hxdhy;

717:         f[j][i] = uxx + uyy - hx*hy*(lambda*PetscExpScalar(u)-(PetscExpReal(PetscSinReal(m*PETSC_PI*x*(1 - y))*PetscSinReal(n*PETSC_PI*(1 - x)*y))*lambda) + PetscSqr(PETSC_PI)*(-2*m*n*((-1 + x)*x + (-1 + y)*y)*PetscCosReal(m*PETSC_PI*x*(-1 + y))*PetscCosReal(n*PETSC_PI*(-1 + x)*y) + (PetscSqr(m)*(PetscSqr(x) + PetscSqr(-1 + y)) + PetscSqr(n)*(PetscSqr(-1 + x) + PetscSqr(y)))*PetscSinReal(m*PETSC_PI*x*(-1 + y))*PetscSinReal(n*PETSC_PI*(-1 + x)*y)));
718:       }
719:     }
720:   }
721:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
722:   PetscLogFlops(11.0*info->ym*info->xm);
723:   return(0);

725: }

729: /* ---------------------------------------------------------------------------------
730:  FormFunctionLocalMMS4 - Evaluates nonlinear function, F(x) on local process patch

732:  u(x,y) = (x^4 - Lx^2 x^2)(y^2-Ly^2 y^2)

734:  -laplacian u* - lambda exp(u*) = 2x^2(x^2-Lx^2)(Ly^2-6y^2)+2y^2(Lx^2-6x^2)(y^2-Ly^2)
735:                                   -lambda exp((x^4 - Lx^2 x^2)(y^2-Ly^2 y^2))

737:   Remark: the above is subtracted from the residual
738: -----------------------------------------------------------------------------------*/
739: PetscErrorCode FormFunctionLocalMMS4(DMDALocalInfo *info,PetscScalar **vx,PetscScalar **f,AppCtx *user)
740: {
742:   PetscInt       i,j;
743:   PetscReal      lambda,hx,hy,hxdhy,hydhx,Lx,Ly;
744:   PetscScalar    u,ue,uw,un,us,uxx,uyy;
745:   PetscReal      x,y;
746:   DM             coordDA;
747:   Vec            coordinates;
748:   DMDACoor2d   **coords;

751:   lambda = user->param;
752:   hx     = 1.0/(PetscReal)(info->mx-1);
753:   hy     = 1.0/(PetscReal)(info->my-1);
754:   hxdhy  = hx/hy;
755:   hydhx  = hy/hx;
756:   /* Extract coordinates */
757:   DMGetCoordinateDM(info->da, &coordDA);
758:   DMGetCoordinates(info->da, &coordinates);
759:   DMDAVecGetArray(coordDA, coordinates, &coords);

761:   Lx = PetscRealPart(coords[info->ys][info->xs+info->xm-1].x - coords[info->ys][info->xs].x);
762:   Ly = PetscRealPart(coords[info->ys+info->ym-1][info->xs].y - coords[info->ys][info->xs].y);

764:   /* Compute function over the locally owned part of the grid */
765:   for (j=info->ys; j<info->ys+info->ym; j++) {
766:     for (i=info->xs; i<info->xs+info->xm; i++) {
767:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
768:         f[j][i] = 2.0*(hydhx+hxdhy)*vx[j][i];
769:       } else {
770:         x  = PetscRealPart(coords[j][i].x);
771:         y  = PetscRealPart(coords[j][i].y);
772:         u  = vx[j][i];
773:         uw = vx[j][i-1];
774:         ue = vx[j][i+1];
775:         un = vx[j-1][i];
776:         us = vx[j+1][i];

778:         if (i-1 == 0) uw = 0.;
779:         if (i+1 == info->mx-1) ue = 0.;
780:         if (j-1 == 0) un = 0.;
781:         if (j+1 == info->my-1) us = 0.;

783:         uxx     = (2.0*u - uw - ue)*hydhx;
784:         uyy     = (2.0*u - un - us)*hxdhy;
785:         f[j][i] = uxx + uyy - hx*hy*(lambda*PetscExpScalar(u)+ 2*PetscSqr(x)*(PetscSqr(x)-PetscSqr(Lx))*(PetscSqr(Ly)-6*PetscSqr(y))+2*PetscSqr(y)*(PetscSqr(Lx)-6*PetscSqr(x))*(PetscSqr(y)-PetscSqr(Ly))-lambda*PetscExpReal((PetscPowReal(x,4)-PetscSqr(Lx)*PetscSqr(x))*(PetscPowReal(y,4)-PetscSqr(Ly)*PetscSqr(y))));
786:       }
787:     }
788:   }
789:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
790:   PetscLogFlops(11.0*info->ym*info->xm);
791:   return(0);

793: }

797: /* FormObjectiveLocal - Evaluates nonlinear function, F(x) on local process patch */
798: PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info,PetscScalar **x,PetscReal *obj,AppCtx *user)
799: {
801:   PetscInt       i,j;
802:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc,lobj=0;
803:   PetscScalar    u,ue,uw,un,us,uxux,uyuy;
804:   MPI_Comm       comm;

807:   *obj   = 0;
808:   PetscObjectGetComm((PetscObject)info->da,&comm);
809:   lambda = user->param;
810:   hx     = 1.0/(PetscReal)(info->mx-1);
811:   hy     = 1.0/(PetscReal)(info->my-1);
812:   sc     = hx*hy*lambda;
813:   hxdhy  = hx/hy;
814:   hydhx  = hy/hx;
815:   /*
816:      Compute function over the locally owned part of the grid
817:   */
818:   for (j=info->ys; j<info->ys+info->ym; j++) {
819:     for (i=info->xs; i<info->xs+info->xm; i++) {
820:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
821:         lobj += PetscRealPart((hydhx + hxdhy)*x[j][i]*x[j][i]);
822:       } else {
823:         u  = x[j][i];
824:         uw = x[j][i-1];
825:         ue = x[j][i+1];
826:         un = x[j-1][i];
827:         us = x[j+1][i];

829:         if (i-1 == 0) uw = 0.;
830:         if (i+1 == info->mx-1) ue = 0.;
831:         if (j-1 == 0) un = 0.;
832:         if (j+1 == info->my-1) us = 0.;

834:         /* F[u] = 1/2\int_{\omega}\nabla^2u(x)*u(x)*dx */

836:         uxux = u*(2.*u - ue - uw)*hydhx;
837:         uyuy = u*(2.*u - un - us)*hxdhy;

839:         lobj += PetscRealPart(0.5*(uxux + uyuy) - sc*PetscExpScalar(u));
840:       }
841:     }
842:   }
843:   PetscLogFlops(12.0*info->ym*info->xm);
844:   MPI_Allreduce(&lobj,obj,1,MPIU_REAL,MPIU_SUM,comm);
845:   return(0);
846: }

850: /*
851:    FormJacobianLocal - Evaluates Jacobian matrix on local process patch
852: */
853: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,Mat jacpre,AppCtx *user)
854: {
856:   PetscInt       i,j,k;
857:   MatStencil     col[5],row;
858:   PetscScalar    lambda,v[5],hx,hy,hxdhy,hydhx,sc;
859:   DM             coordDA;
860:   Vec            coordinates;
861:   DMDACoor2d   **coords;

864:   lambda = user->param;
865:   /* Extract coordinates */
866:   DMGetCoordinateDM(info->da, &coordDA);
867:   DMGetCoordinates(info->da, &coordinates);
868:   DMDAVecGetArray(coordDA, coordinates, &coords);
869:   hx     = info->xm > 1 ? PetscRealPart(coords[info->ys][info->xs+1].x) - PetscRealPart(coords[info->ys][info->xs].x) : 1.0;
870:   hy     = info->ym > 1 ? PetscRealPart(coords[info->ys+1][info->xs].y) - PetscRealPart(coords[info->ys][info->xs].y) : 1.0;
871:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
872:   hxdhy  = hx/hy;
873:   hydhx  = hy/hx;
874:   sc     = hx*hy*lambda;

877:   /*
878:      Compute entries for the locally owned part of the Jacobian.
879:       - Currently, all PETSc parallel matrix formats are partitioned by
880:         contiguous chunks of rows across the processors.
881:       - Each processor needs to insert only elements that it owns
882:         locally (but any non-local elements will be sent to the
883:         appropriate processor during matrix assembly).
884:       - Here, we set all entries for a particular row at once.
885:       - We can set matrix entries either using either
886:         MatSetValuesLocal() or MatSetValues(), as discussed above.
887:   */
888:   for (j=info->ys; j<info->ys+info->ym; j++) {
889:     for (i=info->xs; i<info->xs+info->xm; i++) {
890:       row.j = j; row.i = i;
891:       /* boundary points */
892:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
893:         v[0] =  2.0*(hydhx + hxdhy);
894:         MatSetValuesStencil(jacpre,1,&row,1,&row,v,INSERT_VALUES);
895:       } else {
896:         k = 0;
897:         /* interior grid points */
898:         if (j-1 != 0) {
899:           v[k]     = -hxdhy;
900:           col[k].j = j - 1; col[k].i = i;
901:           k++;
902:         }
903:         if (i-1 != 0) {
904:           v[k]     = -hydhx;
905:           col[k].j = j;     col[k].i = i-1;
906:           k++;
907:         }

909:         v[k] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(x[j][i]); col[k].j = row.j; col[k].i = row.i; k++;

911:         if (i+1 != info->mx-1) {
912:           v[k]     = -hydhx;
913:           col[k].j = j;     col[k].i = i+1;
914:           k++;
915:         }
916:         if (j+1 != info->mx-1) {
917:           v[k]     = -hxdhy;
918:           col[k].j = j + 1; col[k].i = i;
919:           k++;
920:         }
921:         MatSetValuesStencil(jacpre,1,&row,k,col,v,INSERT_VALUES);
922:       }
923:     }
924:   }

926:   /*
927:      Assemble matrix, using the 2-step process:
928:        MatAssemblyBegin(), MatAssemblyEnd().
929:   */
930:   MatAssemblyBegin(jacpre,MAT_FINAL_ASSEMBLY);
931:   MatAssemblyEnd(jacpre,MAT_FINAL_ASSEMBLY);
932:   /*
933:      Tell the matrix we will never add a new nonzero location to the
934:      matrix. If we do, it will generate an error.
935:   */
936:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
937:   return(0);
938: }

940: #if defined(PETSC_HAVE_MATLAB_ENGINE)
943: PetscErrorCode FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
944: {
945:   AppCtx         *user = (AppCtx*)ptr;
947:   PetscInt       Mx,My;
948:   PetscReal      lambda,hx,hy;
949:   Vec            localX,localF;
950:   MPI_Comm       comm;
951:   DM             da;

954:   SNESGetDM(snes,&da);
955:   DMGetLocalVector(da,&localX);
956:   DMGetLocalVector(da,&localF);
957:   PetscObjectSetName((PetscObject)localX,"localX");
958:   PetscObjectSetName((PetscObject)localF,"localF");
959:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

961:   lambda = user->param;
962:   hx     = 1.0/(PetscReal)(Mx-1);
963:   hy     = 1.0/(PetscReal)(My-1);

965:   PetscObjectGetComm((PetscObject)snes,&comm);
966:   /*
967:      Scatter ghost points to local vector,using the 2-step process
968:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
969:      By placing code between these two statements, computations can be
970:      done while messages are in transition.
971:   */
972:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
973:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
974:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
975:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
976:   PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);

978:   /*
979:      Insert values into global vector
980:   */
981:   DMLocalToGlobalBegin(da,localF,INSERT_VALUES,F);
982:   DMLocalToGlobalEnd(da,localF,INSERT_VALUES,F);
983:   DMRestoreLocalVector(da,&localX);
984:   DMRestoreLocalVector(da,&localF);
985:   return(0);
986: }
987: #endif

989: /* ------------------------------------------------------------------- */
992: /*
993:       Applies some sweeps on nonlinear Gauss-Seidel on each process

995:  */
996: PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx)
997: {
998:   PetscInt       i,j,k,Mx,My,xs,ys,xm,ym,its,tot_its,sweeps,l;
1000:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc;
1001:   PetscScalar    **x,**b,bij,F,F0=0,J,u,un,us,ue,eu,uw,uxx,uyy,y;
1002:   PetscReal      atol,rtol,stol;
1003:   DM             da;
1004:   AppCtx         *user;
1005:   Vec            localX,localB;

1008:   tot_its = 0;
1009:   SNESNGSGetSweeps(snes,&sweeps);
1010:   SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);
1011:   SNESGetDM(snes,&da);
1012:   DMGetApplicationContext(da,(void**)&user);

1014:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

1016:   lambda = user->param;
1017:   hx     = 1.0/(PetscReal)(Mx-1);
1018:   hy     = 1.0/(PetscReal)(My-1);
1019:   sc     = hx*hy*lambda;
1020:   hxdhy  = hx/hy;
1021:   hydhx  = hy/hx;

1024:   DMGetLocalVector(da,&localX);
1025:   if (B) {
1026:     DMGetLocalVector(da,&localB);
1027:   }
1028:   for (l=0; l<sweeps; l++) {

1030:     DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
1031:     DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
1032:     if (B) {
1033:       DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);
1034:       DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);
1035:     }
1036:     /*
1037:      Get a pointer to vector data.
1038:      - For default PETSc vectors, VecGetArray() returns a pointer to
1039:      the data array.  Otherwise, the routine is implementation dependent.
1040:      - You MUST call VecRestoreArray() when you no longer need access to
1041:      the array.
1042:      */
1043:     DMDAVecGetArray(da,localX,&x);
1044:     if (B) DMDAVecGetArray(da,localB,&b);
1045:     /*
1046:      Get local grid boundaries (for 2-dimensional DMDA):
1047:      xs, ys   - starting grid indices (no ghost points)
1048:      xm, ym   - widths of local grid (no ghost points)
1049:      */
1050:     DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

1052:     for (j=ys; j<ys+ym; j++) {
1053:       for (i=xs; i<xs+xm; i++) {
1054:         if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
1055:           /* boundary conditions are all zero Dirichlet */
1056:           x[j][i] = 0.0;
1057:         } else {
1058:           if (B) bij = b[j][i];
1059:           else   bij = 0.;

1061:           u  = x[j][i];
1062:           un = x[j-1][i];
1063:           us = x[j+1][i];
1064:           ue = x[j][i-1];
1065:           uw = x[j][i+1];

1067:           for (k=0; k<its; k++) {
1068:             eu  = PetscExpScalar(u);
1069:             uxx = (2.0*u - ue - uw)*hydhx;
1070:             uyy = (2.0*u - un - us)*hxdhy;
1071:             F   = uxx + uyy - sc*eu - bij;
1072:             if (k == 0) F0 = F;
1073:             J  = 2.0*(hydhx + hxdhy) - sc*eu;
1074:             y  = F/J;
1075:             u -= y;
1076:             tot_its++;

1078:             if (atol > PetscAbsReal(PetscRealPart(F)) ||
1079:                 rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) ||
1080:                 stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) {
1081:               break;
1082:             }
1083:           }
1084:           x[j][i] = u;
1085:         }
1086:       }
1087:     }
1088:     /*
1089:      Restore vector
1090:      */
1091:     DMDAVecRestoreArray(da,localX,&x);
1092:     DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
1093:     DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
1094:   }
1095:   PetscLogFlops(tot_its*(21.0));
1096:   DMRestoreLocalVector(da,&localX);
1097:   if (B) {
1098:     DMDAVecRestoreArray(da,localB,&b);
1099:     DMRestoreLocalVector(da,&localB);
1100:   }
1101:   return(0);
1102: }