Actual source code: ex5.c

petsc-3.6.0 2015-06-09
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";

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

 16: /* ------------------------------------------------------------------------

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

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

 23:     with boundary conditions

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

 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.


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

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

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

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

 44:   ------------------------------------------------------------------------- */

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

 55: /*
 56:    User-defined application context - contains data needed by the
 57:    application-provided call-back routines, FormJacobianLocal() and
 58:    FormFunctionLocal().
 59: */
 60: typedef struct {
 61:   PassiveReal param;          /* test problem parameter */
 62: } AppCtx;

 64: /*
 65:    User-defined routines
 66: */
 67: extern PetscErrorCode FormInitialGuess(DM,AppCtx*,Vec);
 68: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
 69: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,Mat,AppCtx*);
 70: extern PetscErrorCode FormObjectiveLocal(DMDALocalInfo*,PetscScalar**,PetscReal*,AppCtx*);
 71: #if defined(PETSC_HAVE_MATLAB_ENGINE)
 72: extern PetscErrorCode FormFunctionMatlab(SNES,Vec,Vec,void*);
 73: #endif
 74: extern PetscErrorCode NonlinearGS(SNES,Vec,Vec,void*);

 78: int main(int argc,char **argv)
 79: {
 80:   SNES           snes;                         /* nonlinear solver */
 81:   Vec            x;                            /* solution vector */
 82:   AppCtx         user;                         /* user-defined work context */
 83:   PetscInt       its;                          /* iterations for convergence */
 85:   PetscReal      bratu_lambda_max = 6.81;
 86:   PetscReal      bratu_lambda_min = 0.;
 87:   PetscBool      flg              = PETSC_FALSE;
 88:   DM             da;
 89: #if defined(PETSC_HAVE_MATLAB_ENGINE)
 90:   Vec            r               = NULL;
 91:   PetscBool      matlab_function = PETSC_FALSE;
 92: #endif
 93:   KSP            ksp;
 94:   PetscInt       lits,slits;

 96:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 97:      Initialize program
 98:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

102:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103:      Initialize problem parameters
104:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105:   user.param = 6.0;
106:   PetscOptionsGetReal(NULL,"-par",&user.param,NULL);
107:   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);

109:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110:      Create nonlinear solver context
111:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112:   SNESCreate(PETSC_COMM_WORLD,&snes);
113:   SNESSetCountersReset(snes,PETSC_FALSE);
114:   SNESSetNGS(snes, NonlinearGS, NULL);

116:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117:      Create distributed array (DMDA) to manage parallel grid and vectors
118:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
120:   DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
121:   DMSetApplicationContext(da,&user);
122:   SNESSetDM(snes,da);
123:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124:      Extract global vectors from DMDA; then duplicate for remaining
125:      vectors that are the same types
126:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
127:   DMCreateGlobalVector(da,&x);

129:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130:      Set local function evaluation routine
131:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132:   DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocal,&user);
133:   PetscOptionsGetBool(NULL,"-fd",&flg,NULL);
134:   if (!flg) {
135:     DMDASNESSetJacobianLocal(da,(DMDASNESJacobian)FormJacobianLocal,&user);
136:   }

138:   PetscOptionsGetBool(NULL,"-obj",&flg,NULL);
139:   if (flg) {
140:     DMDASNESSetObjectiveLocal(da,(DMDASNESObjective)FormObjectiveLocal,&user);
141:   }

143: #if defined(PETSC_HAVE_MATLAB_ENGINE)
144:   PetscOptionsGetBool(NULL,"-matlab_function",&matlab_function,0);
145:   if (matlab_function) {
146:     VecDuplicate(x,&r);
147:     SNESSetFunction(snes,r,FormFunctionMatlab,&user);
148:   }
149: #endif

151:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152:      Customize nonlinear solver; set runtime options
153:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154:   SNESSetFromOptions(snes);

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

158:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159:      Solve nonlinear system
160:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161:   SNESSolve(snes,NULL,x);
162:   SNESGetIterationNumber(snes,&its);

164:   SNESGetLinearSolveIterations(snes,&slits);
165:   SNESGetKSP(snes,&ksp);
166:   KSPGetTotalIterations(ksp,&lits);
167:   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);
168:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

171:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172:      Free work space.  All PETSc objects should be destroyed when they
173:      are no longer needed.
174:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175: #if defined(PETSC_HAVE_MATLAB_ENGINE)
176:   VecDestroy(&r);
177: #endif
178:   VecDestroy(&x);
179:   SNESDestroy(&snes);
180:   DMDestroy(&da);
181:   PetscFinalize();
182:   return 0;
183: }
184: /* ------------------------------------------------------------------- */
187: /*
188:    FormInitialGuess - Forms initial approximation.

190:    Input Parameters:
191:    user - user-defined application context
192:    X - vector

194:    Output Parameter:
195:    X - vector
196:  */
197: PetscErrorCode FormInitialGuess(DM da,AppCtx *user,Vec X)
198: {
199:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
201:   PetscReal      lambda,temp1,temp,hx,hy;
202:   PetscScalar    **x;

205:   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);

207:   lambda = user->param;
208:   hx     = 1.0/(PetscReal)(Mx-1);
209:   hy     = 1.0/(PetscReal)(My-1);
210:   temp1  = lambda/(lambda + 1.0);

212:   /*
213:      Get a pointer to vector data.
214:        - For default PETSc vectors, VecGetArray() returns a pointer to
215:          the data array.  Otherwise, the routine is implementation dependent.
216:        - You MUST call VecRestoreArray() when you no longer need access to
217:          the array.
218:   */
219:   DMDAVecGetArray(da,X,&x);

221:   /*
222:      Get local grid boundaries (for 2-dimensional DMDA):
223:        xs, ys   - starting grid indices (no ghost points)
224:        xm, ym   - widths of local grid (no ghost points)

226:   */
227:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

229:   /*
230:      Compute initial guess over the locally owned part of the grid
231:   */
232:   for (j=ys; j<ys+ym; j++) {
233:     temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
234:     for (i=xs; i<xs+xm; i++) {
235:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
236:         /* boundary conditions are all zero Dirichlet */
237:         x[j][i] = 0.0;
238:       } else {
239:         x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
240:       }
241:     }
242:   }

244:   /*
245:      Restore vector
246:   */
247:   DMDAVecRestoreArray(da,X,&x);
248:   return(0);
249: }
250: /* ------------------------------------------------------------------- */
253: /*
254:    FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch


257:  */
258: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
259: {
261:   PetscInt       i,j;
262:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc;
263:   PetscScalar    u,ue,uw,un,us,uxx,uyy;

266:   lambda = user->param;
267:   hx     = 1.0/(PetscReal)(info->mx-1);
268:   hy     = 1.0/(PetscReal)(info->my-1);
269:   sc     = hx*hy*lambda;
270:   hxdhy  = hx/hy;
271:   hydhx  = hy/hx;
272:   /*
273:      Compute function over the locally owned part of the grid
274:   */
275:   for (j=info->ys; j<info->ys+info->ym; j++) {
276:     for (i=info->xs; i<info->xs+info->xm; i++) {
277:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
278:         f[j][i] = 2.0*(hydhx+hxdhy)*x[j][i];
279:       } else {
280:         u  = x[j][i];
281:         uw = x[j][i-1];
282:         ue = x[j][i+1];
283:         un = x[j-1][i];
284:         us = x[j+1][i];

286:         if (i-1 == 0) uw = 0.;
287:         if (i+1 == info->mx-1) ue = 0.;
288:         if (j-1 == 0) un = 0.;
289:         if (j+1 == info->my-1) us = 0.;

291:         uxx     = (2.0*u - uw - ue)*hydhx;
292:         uyy     = (2.0*u - un - us)*hxdhy;
293:         f[j][i] = uxx + uyy - sc*PetscExpScalar(u);
294:       }
295:     }
296:   }
297:   PetscLogFlops(11.0*info->ym*info->xm);
298:   return(0);
299: }


304: /*
305:    FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch


308:  */
309: PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info,PetscScalar **x,PetscReal *obj,AppCtx *user)
310: {
312:   PetscInt       i,j;
313:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc,lobj=0;
314:   PetscScalar    u,ue,uw,un,us,uxux,uyuy;
315:   MPI_Comm       comm;

318:   *obj   = 0;
319:   PetscObjectGetComm((PetscObject)info,&comm);
320:   lambda = user->param;
321:   hx     = 1.0/(PetscReal)(info->mx-1);
322:   hy     = 1.0/(PetscReal)(info->my-1);
323:   sc     = hx*hy*lambda;
324:   hxdhy  = hx/hy;
325:   hydhx  = hy/hx;
326:   /*
327:      Compute function over the locally owned part of the grid
328:   */
329:   for (j=info->ys; j<info->ys+info->ym; j++) {
330:     for (i=info->xs; i<info->xs+info->xm; i++) {
331:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
332:         lobj += PetscRealPart((hydhx + hxdhy)*x[j][i]*x[j][i]);
333:       } else {
334:         u  = x[j][i];
335:         uw = x[j][i-1];
336:         ue = x[j][i+1];
337:         un = x[j-1][i];
338:         us = x[j+1][i];

340:         if (i-1 == 0) uw = 0.;
341:         if (i+1 == info->mx-1) ue = 0.;
342:         if (j-1 == 0) un = 0.;
343:         if (j+1 == info->my-1) us = 0.;

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

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

350:         lobj += PetscRealPart(0.5*(uxux + uyuy) - sc*PetscExpScalar(u));
351:       }
352:     }
353:   }
354:   PetscLogFlops(12.0*info->ym*info->xm);
355:   MPI_Allreduce(&lobj,obj,1,MPIU_REAL,MPIU_SUM,comm);
356:   return(0);
357: }

361: /*
362:    FormJacobianLocal - Evaluates Jacobian matrix on local process patch
363: */
364: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,Mat jacpre,AppCtx *user)
365: {
367:   PetscInt       i,j,k;
368:   MatStencil     col[5],row;
369:   PetscScalar    lambda,v[5],hx,hy,hxdhy,hydhx,sc;

372:   lambda = user->param;
373:   hx     = 1.0/(PetscReal)(info->mx-1);
374:   hy     = 1.0/(PetscReal)(info->my-1);
375:   sc     = hx*hy*lambda;
376:   hxdhy  = hx/hy;
377:   hydhx  = hy/hx;


380:   /*
381:      Compute entries for the locally owned part of the Jacobian.
382:       - Currently, all PETSc parallel matrix formats are partitioned by
383:         contiguous chunks of rows across the processors.
384:       - Each processor needs to insert only elements that it owns
385:         locally (but any non-local elements will be sent to the
386:         appropriate processor during matrix assembly).
387:       - Here, we set all entries for a particular row at once.
388:       - We can set matrix entries either using either
389:         MatSetValuesLocal() or MatSetValues(), as discussed above.
390:   */
391:   for (j=info->ys; j<info->ys+info->ym; j++) {
392:     for (i=info->xs; i<info->xs+info->xm; i++) {
393:       row.j = j; row.i = i;
394:       /* boundary points */
395:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
396:         v[0] =  2.0*(hydhx + hxdhy);
397:         MatSetValuesStencil(jacpre,1,&row,1,&row,v,INSERT_VALUES);
398:       } else {
399:         k = 0;
400:         /* interior grid points */
401:         if (j-1 != 0) {
402:           v[k]     = -hxdhy;
403:           col[k].j = j - 1; col[k].i = i;
404:           k++;
405:         }
406:         if (i-1 != 0) {
407:           v[k]     = -hydhx;
408:           col[k].j = j;     col[k].i = i-1;
409:           k++;
410:         }

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

414:         if (i+1 != info->mx-1) {
415:           v[k]     = -hydhx;
416:           col[k].j = j;     col[k].i = i+1;
417:           k++;
418:         }
419:         if (j+1 != info->mx-1) {
420:           v[k]     = -hxdhy;
421:           col[k].j = j + 1; col[k].i = i;
422:           k++;
423:         }
424:         MatSetValuesStencil(jacpre,1,&row,k,col,v,INSERT_VALUES);
425:       }
426:     }
427:   }

429:   /*
430:      Assemble matrix, using the 2-step process:
431:        MatAssemblyBegin(), MatAssemblyEnd().
432:   */
433:   MatAssemblyBegin(jacpre,MAT_FINAL_ASSEMBLY);
434:   MatAssemblyEnd(jacpre,MAT_FINAL_ASSEMBLY);
435:   /*
436:      Tell the matrix we will never add a new nonzero location to the
437:      matrix. If we do, it will generate an error.
438:   */
439:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
440:   return(0);
441: }

443: #if defined(PETSC_HAVE_MATLAB_ENGINE)
446: PetscErrorCode FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
447: {
448:   AppCtx         *user = (AppCtx*)ptr;
450:   PetscInt       Mx,My;
451:   PetscReal      lambda,hx,hy;
452:   Vec            localX,localF;
453:   MPI_Comm       comm;
454:   DM             da;

457:   SNESGetDM(snes,&da);
458:   DMGetLocalVector(da,&localX);
459:   DMGetLocalVector(da,&localF);
460:   PetscObjectSetName((PetscObject)localX,"localX");
461:   PetscObjectSetName((PetscObject)localF,"localF");
462:   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);

464:   lambda = user->param;
465:   hx     = 1.0/(PetscReal)(Mx-1);
466:   hy     = 1.0/(PetscReal)(My-1);

468:   PetscObjectGetComm((PetscObject)snes,&comm);
469:   /*
470:      Scatter ghost points to local vector,using the 2-step process
471:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
472:      By placing code between these two statements, computations can be
473:      done while messages are in transition.
474:   */
475:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
476:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
477:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
478:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
479:   PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);

481:   /*
482:      Insert values into global vector
483:   */
484:   DMLocalToGlobalBegin(da,localF,INSERT_VALUES,F);
485:   DMLocalToGlobalEnd(da,localF,INSERT_VALUES,F);
486:   DMRestoreLocalVector(da,&localX);
487:   DMRestoreLocalVector(da,&localF);
488:   return(0);
489: }
490: #endif

492: /* ------------------------------------------------------------------- */
495: /*
496:       Applies some sweeps on nonlinear Gauss-Seidel on each process

498:  */
499: PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx)
500: {
501:   PetscInt       i,j,k,Mx,My,xs,ys,xm,ym,its,tot_its,sweeps,l;
503:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc;
504:   PetscScalar    **x,**b,bij,F,F0=0,J,u,un,us,ue,eu,uw,uxx,uyy,y;
505:   PetscReal      atol,rtol,stol;
506:   DM             da;
507:   AppCtx         *user;
508:   Vec            localX,localB;

511:   tot_its = 0;
512:   SNESNGSGetSweeps(snes,&sweeps);
513:   SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);
514:   SNESGetDM(snes,&da);
515:   DMGetApplicationContext(da,(void**)&user);

517:   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);

519:   lambda = user->param;
520:   hx     = 1.0/(PetscReal)(Mx-1);
521:   hy     = 1.0/(PetscReal)(My-1);
522:   sc     = hx*hy*lambda;
523:   hxdhy  = hx/hy;
524:   hydhx  = hy/hx;


527:   DMGetLocalVector(da,&localX);
528:   if (B) {
529:     DMGetLocalVector(da,&localB);
530:   }
531:   for (l=0; l<sweeps; l++) {

533:     DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
534:     DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
535:     if (B) {
536:       DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);
537:       DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);
538:     }
539:     /*
540:      Get a pointer to vector data.
541:      - For default PETSc vectors, VecGetArray() returns a pointer to
542:      the data array.  Otherwise, the routine is implementation dependent.
543:      - You MUST call VecRestoreArray() when you no longer need access to
544:      the array.
545:      */
546:     DMDAVecGetArray(da,localX,&x);
547:     if (B) DMDAVecGetArray(da,localB,&b);
548:     /*
549:      Get local grid boundaries (for 2-dimensional DMDA):
550:      xs, ys   - starting grid indices (no ghost points)
551:      xm, ym   - widths of local grid (no ghost points)
552:      */
553:     DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

555:     for (j=ys; j<ys+ym; j++) {
556:       for (i=xs; i<xs+xm; i++) {
557:         if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
558:           /* boundary conditions are all zero Dirichlet */
559:           x[j][i] = 0.0;
560:         } else {
561:           if (B) bij = b[j][i];
562:           else   bij = 0.;

564:           u  = x[j][i];
565:           un = x[j-1][i];
566:           us = x[j+1][i];
567:           ue = x[j][i-1];
568:           uw = x[j][i+1];

570:           for (k=0; k<its; k++) {
571:             eu  = PetscExpScalar(u);
572:             uxx = (2.0*u - ue - uw)*hydhx;
573:             uyy = (2.0*u - un - us)*hxdhy;
574:             F   = uxx + uyy - sc*eu - bij;
575:             if (k == 0) F0 = F;
576:             J  = 2.0*(hydhx + hxdhy) - sc*eu;
577:             y  = F/J;
578:             u -= y;
579:             tot_its++;

581:             if (atol > PetscAbsReal(PetscRealPart(F)) ||
582:                 rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) ||
583:                 stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) {
584:               break;
585:             }
586:           }
587:           x[j][i] = u;
588:         }
589:       }
590:     }
591:     /*
592:      Restore vector
593:      */
594:     DMDAVecRestoreArray(da,localX,&x);
595:     DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
596:     DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
597:   }
598:   PetscLogFlops(tot_its*(21.0));
599:   DMRestoreLocalVector(da,&localX);
600:   if (B) {
601:     DMDAVecRestoreArray(da,localB,&b);
602:     DMRestoreLocalVector(da,&localB);
603:   }
604:   return(0);
605: }