Actual source code: ex5.c

petsc-master 2016-09-28
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 pmat -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 pmat -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);if (ierr) return ierr;

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:   DMSetFromOptions(da);
141:   DMSetUp(da);
142:   DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
143:   DMSetApplicationContext(da,&user);
144:   SNESSetDM(snes,da);
145:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146:      Extract global vectors from DMDA; then duplicate for remaining
147:      vectors that are the same types
148:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149:   DMCreateGlobalVector(da,&x);

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

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

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

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

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

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

192:   SNESGetLinearSolveIterations(snes,&slits);
193:   SNESGetKSP(snes,&ksp);
194:   KSPGetTotalIterations(ksp,&lits);
195:   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);
196:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
197:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

242:    Input Parameters:
243:    da - The DM
244:    user - user-defined application context

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

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

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

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

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

278:   */
279:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

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

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

305: /*
306:   FormExactSolution1 - Forms initial approximation.

308:   Input Parameters:
309:   da - The DM
310:   user - user-defined application context

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

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

345: /*
346:   FormExactSolution2 - Forms initial approximation.

348:   Input Parameters:
349:   da - The DM
350:   user - user-defined application context

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

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

385: /*
386:   FormExactSolution3 - Forms initial approximation.

388:   Input Parameters:
389:   da - The DM
390:   user - user-defined application context

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

405:   m = PetscPowReal(2,user->mPar);
406:   n = PetscPowReal(2,user->nPar);

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

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

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

431: /*
432:   FormExactSolution4 - Forms initial approximation.

434:   Input Parameters:
435:   da - The DM
436:   user - user-defined application context

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

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

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

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


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

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

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

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

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

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

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

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

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

578:         uxx     = (2.0*u - uw - ue)*hydhx;
579:         uyy     = (2.0*u - un - us)*hxdhy;
580:         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)));
581:       }
582:     }
583:   }
584:   DMDAVecRestoreArray(info->da, bcv, &bcx);
585:   DMRestoreNamedLocalVector(info->da, "_petsc_boundary_conditions_", &bcv);
586:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
587:   PetscLogFlops(11.0*info->ym*info->xm);
588:   return(0);
589: }

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

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

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

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

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

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

646:         uxx     = (2.0*u - uw - ue)*hydhx;
647:         uyy     = (2.0*u - un - us)*hxdhy;
648:         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)));
649:       }
650:     }
651:   }
652:   DMDAVecRestoreArray(info->da, bcv, &bcx);
653:   DMRestoreNamedLocalVector(info->da, "_petsc_boundary_conditions_", &bcv);
654:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
655:   PetscLogFlops(11.0*info->ym*info->xm);
656:   return(0);
657: }

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

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

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

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

684:   m = PetscPowReal(2,user->mPar);
685:   n = PetscPowReal(2,user->nPar);

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

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

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

719:         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)));
720:       }
721:     }
722:   }
723:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
724:   PetscLogFlops(11.0*info->ym*info->xm);
725:   return(0);

727: }

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

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

736:  -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)
737:                                   -lambda exp((x^4 - Lx^2 x^2)(y^2-Ly^2 y^2))

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

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

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

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

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

785:         uxx     = (2.0*u - uw - ue)*hydhx;
786:         uyy     = (2.0*u - un - us)*hxdhy;
787:         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))));
788:       }
789:     }
790:   }
791:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
792:   PetscLogFlops(11.0*info->ym*info->xm);
793:   return(0);

795: }

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

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

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

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

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

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

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

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


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

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

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

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

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

956:   SNESGetDM(snes,&da);
957:   DMGetLocalVector(da,&localX);
958:   DMGetLocalVector(da,&localF);
959:   PetscObjectSetName((PetscObject)localX,"localX");
960:   PetscObjectSetName((PetscObject)localF,"localF");
961:   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);

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

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

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

991: /* ------------------------------------------------------------------- */
994: /*
995:       Applies some sweeps on nonlinear Gauss-Seidel on each process

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

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

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

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


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

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

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

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

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

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