Actual source code: ex5.c

petsc-master 2014-12-20
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

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

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

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

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

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

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

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

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

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

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

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

161:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

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

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

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

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

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


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

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

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

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


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


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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

485: /* ------------------------------------------------------------------- */
488: /*
489:       Applies some sweeps on nonlinear Gauss-Seidel on each process

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

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

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

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


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

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

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

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

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

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