Actual source code: ex5.c

petsc-3.4.4 2014-03-13
  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 <petscdmda.h>
 51: #include <petscsnes.h>
 52: #include <petscmatlab.h>

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

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

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

 93:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 94:      Initialize program
 95:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

106:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107:      Create nonlinear solver context
108:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109:   SNESCreate(PETSC_COMM_WORLD,&snes);
110:   SNESSetGS(snes, NonlinearGS, NULL);

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

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

134:   PetscOptionsGetBool(NULL,"-obj",&flg,NULL);
135:   if (flg) {
136:     DMDASNESSetObjectiveLocal(da,(DMDASNESObjective)FormObjectiveLocal,&user);
137:   }

139: #if defined(PETSC_HAVE_MATLAB_ENGINE)
140:   PetscOptionsGetBool(NULL,"-matlab_function",&matlab_function,0);
141:   if (matlab_function) {
142:     VecDuplicate(x,&r);
143:     SNESSetFunction(snes,r,FormFunctionMatlab,&user);
144:   }
145: #endif

147:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148:      Customize nonlinear solver; set runtime options
149:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150:   SNESSetFromOptions(snes);

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

154:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155:      Solve nonlinear system
156:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157:   SNESSolve(snes,NULL,x);
158:   SNESGetIterationNumber(snes,&its);

160:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

163:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164:      Free work space.  All PETSc objects should be destroyed when they
165:      are no longer needed.
166:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167: #if defined(PETSC_HAVE_MATLAB_ENGINE)
168:   VecDestroy(&r);
169: #endif
170:   VecDestroy(&x);
171:   SNESDestroy(&snes);
172:   DMDestroy(&da);
173:   PetscFinalize();
174:   return 0;
175: }
176: /* ------------------------------------------------------------------- */
179: /*
180:    FormInitialGuess - Forms initial approximation.

182:    Input Parameters:
183:    user - user-defined application context
184:    X - vector

186:    Output Parameter:
187:    X - vector
188:  */
189: PetscErrorCode FormInitialGuess(DM da,AppCtx *user,Vec X)
190: {
191:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
193:   PetscReal      lambda,temp1,temp,hx,hy;
194:   PetscScalar    **x;

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

199:   lambda = user->param;
200:   hx     = 1.0/(PetscReal)(Mx-1);
201:   hy     = 1.0/(PetscReal)(My-1);
202:   temp1  = lambda/(lambda + 1.0);

204:   /*
205:      Get a pointer to vector data.
206:        - For default PETSc vectors, VecGetArray() returns a pointer to
207:          the data array.  Otherwise, the routine is implementation dependent.
208:        - You MUST call VecRestoreArray() when you no longer need access to
209:          the array.
210:   */
211:   DMDAVecGetArray(da,X,&x);

213:   /*
214:      Get local grid boundaries (for 2-dimensional DMDA):
215:        xs, ys   - starting grid indices (no ghost points)
216:        xm, ym   - widths of local grid (no ghost points)

218:   */
219:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

221:   /*
222:      Compute initial guess over the locally owned part of the grid
223:   */
224:   for (j=ys; j<ys+ym; j++) {
225:     temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
226:     for (i=xs; i<xs+xm; i++) {
227:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
228:         /* boundary conditions are all zero Dirichlet */
229:         x[j][i] = 0.0;
230:       } else {
231:         x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
232:       }
233:     }
234:   }

236:   /*
237:      Restore vector
238:   */
239:   DMDAVecRestoreArray(da,X,&x);
240:   return(0);
241: }
242: /* ------------------------------------------------------------------- */
245: /*
246:    FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch


249:  */
250: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
251: {
253:   PetscInt       i,j;
254:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc;
255:   PetscScalar    u,ue,uw,un,us,uxx,uyy;

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

278:         if (i-1 == 0) uw = 0.;
279:         if (i+1 == info->mx-1) ue = 0.;
280:         if (j-1 == 0) un = 0.;
281:         if (j+1 == info->my-1) us = 0.;

283:         uxx     = (2.0*u - uw - ue)*hydhx;
284:         uyy     = (2.0*u - un - us)*hxdhy;
285:         f[j][i] = uxx + uyy - sc*PetscExpScalar(u);
286:       }
287:     }
288:   }
289:   PetscLogFlops(11.0*info->ym*info->xm);
290:   return(0);
291: }


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


300:  */
301: PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info,PetscScalar **x,PetscReal *obj,AppCtx *user)
302: {
304:   PetscInt       i,j;
305:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc,lobj=0;
306:   PetscScalar    u,ue,uw,un,us,uxux,uyuy;
307:   MPI_Comm       comm;

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

332:         if (i-1 == 0) uw = 0.;
333:         if (i+1 == info->mx-1) ue = 0.;
334:         if (j-1 == 0) un = 0.;
335:         if (j+1 == info->my-1) us = 0.;

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

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

342:         lobj += PetscRealPart(0.5*(uxux + uyuy) - sc*PetscExpScalar(u));
343:       }
344:     }
345:   }
346:   PetscLogFlops(12.0*info->ym*info->xm);
347:   MPI_Allreduce(&lobj,obj,1,MPIU_REAL,MPIU_SUM,comm);
348:   return(0);
349: }

353: /*
354:    FormJacobianLocal - Evaluates Jacobian matrix on local process patch
355: */
356: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jacpre,Mat jac,MatStructure *flg,AppCtx *user)
357: {
359:   PetscInt       i,j,k;
360:   MatStencil     col[5],row;
361:   PetscScalar    lambda,v[5],hx,hy,hxdhy,hydhx,sc;

364:   lambda = user->param;
365:   hx     = 1.0/(PetscReal)(info->mx-1);
366:   hy     = 1.0/(PetscReal)(info->my-1);
367:   sc     = hx*hy*lambda;
368:   hxdhy  = hx/hy;
369:   hydhx  = hy/hx;


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

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

406:         if (i+1 != info->mx-1) {
407:           v[k]     = -hydhx;
408:           col[k].j = j;     col[k].i = i+1;
409:           k++;
410:         }
411:         if (j+1 != info->mx-1) {
412:           v[k]     = -hxdhy;
413:           col[k].j = j + 1; col[k].i = i;
414:           k++;
415:         }
416:         MatSetValuesStencil(jac,1,&row,k,col,v,INSERT_VALUES);
417:       }
418:     }
419:   }

421:   /*
422:      Assemble matrix, using the 2-step process:
423:        MatAssemblyBegin(), MatAssemblyEnd().
424:   */
425:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
426:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
427:   /*
428:      Tell the matrix we will never add a new nonzero location to the
429:      matrix. If we do, it will generate an error.
430:   */
431:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
432:   return(0);
433: }

435: #if defined(PETSC_HAVE_MATLAB_ENGINE)
438: PetscErrorCode FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
439: {
440:   AppCtx         *user = (AppCtx*)ptr;
442:   PetscInt       Mx,My;
443:   PetscReal      lambda,hx,hy;
444:   Vec            localX,localF;
445:   MPI_Comm       comm;
446:   DM             da;

449:   SNESGetDM(snes,&da);
450:   DMGetLocalVector(da,&localX);
451:   DMGetLocalVector(da,&localF);
452:   PetscObjectSetName((PetscObject)localX,"localX");
453:   PetscObjectSetName((PetscObject)localF,"localF");
454:   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);

456:   lambda = user->param;
457:   hx     = 1.0/(PetscReal)(Mx-1);
458:   hy     = 1.0/(PetscReal)(My-1);

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

473:   /*
474:      Insert values into global vector
475:   */
476:   DMLocalToGlobalBegin(da,localF,INSERT_VALUES,F);
477:   DMLocalToGlobalEnd(da,localF,INSERT_VALUES,F);
478:   DMRestoreLocalVector(da,&localX);
479:   DMRestoreLocalVector(da,&localF);
480:   return(0);
481: }
482: #endif

484: /* ------------------------------------------------------------------- */
487: /*
488:       Applies some sweeps on nonlinear Gauss-Seidel on each process

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

503:   tot_its = 0;
504:   SNESGSGetSweeps(snes,&sweeps);
505:   SNESGSGetTolerances(snes,&atol,&rtol,&stol,&its);
506:   SNESGetDM(snes,&da);
507:   DMGetApplicationContext(da,(void**)&user);

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

511:   lambda = user->param;
512:   hx     = 1.0/(PetscReal)(Mx-1);
513:   hy     = 1.0/(PetscReal)(My-1);
514:   sc     = hx*hy*lambda;
515:   hxdhy  = hx/hy;
516:   hydhx  = hy/hx;


519:   DMGetLocalVector(da,&localX);
520:   if (B) {
521:     DMGetLocalVector(da,&localB);
522:   }
523:   for (l=0; l<sweeps; l++) {

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

547:     for (j=ys; j<ys+ym; j++) {
548:       for (i=xs; i<xs+xm; i++) {
549:         if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
550:           /* boundary conditions are all zero Dirichlet */
551:           x[j][i] = 0.0;
552:         } else {
553:           if (B) bij = b[j][i];
554:           else   bij = 0.;

556:           u  = x[j][i];
557:           un = x[j-1][i];
558:           us = x[j+1][i];
559:           ue = x[j][i-1];
560:           uw = x[j][i+1];

562:           for (k=0; k<its; k++) {
563:             eu  = PetscExpScalar(u);
564:             uxx = (2.0*u - ue - uw)*hydhx;
565:             uyy = (2.0*u - un - us)*hxdhy;
566:             F   = uxx + uyy - sc*eu - bij;
567:             if (k == 0) F0 = F;
568:             J  = 2.0*(hydhx + hxdhy) - sc*eu;
569:             y  = F/J;
570:             u -= y;
571:             tot_its++;

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