Actual source code: ex5.c

petsc-master 2017-06-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*);

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

107:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108:      Initialize program
109:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

111:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

113:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114:      Initialize problem parameters
115:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116:   user.param = 6.0;
117:   PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
118:   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);
119:   PetscOptionsGetInt(NULL,NULL,"-mms",&MMS,NULL);
120:   if (MMS == 3) {
121:     user.mPar = 2;
122:     user.nPar = 1;
123:     PetscOptionsGetInt(NULL,NULL,"-m_par",&user.mPar,NULL);
124:     PetscOptionsGetInt(NULL,NULL,"-n_par",&user.nPar,NULL);
125:   }

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

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

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

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

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

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

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

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

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

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

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

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

238:    Input Parameters:
239:    da - The DM
240:    user - user-defined application context

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

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

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

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

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

274:   */
275:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

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

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

299: /*
300:   FormExactSolution1 - Forms initial approximation.

302:   Input Parameters:
303:   da - The DM
304:   user - user-defined application context

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

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

337: /*
338:   FormExactSolution2 - Forms initial approximation.

340:   Input Parameters:
341:   da - The DM
342:   user - user-defined application context

344:   Output Parameter:
345:   X - vector
346:  */
347: PetscErrorCode FormExactSolution2(DM da, AppCtx *user, Vec U)
348: {
349:   DM             coordDA;
350:   Vec            coordinates;
351:   DMDACoor2d   **coords;
352:   PetscScalar  **u;
353:   PetscReal      x, y;
354:   PetscInt       xs, ys, xm, ym, i, j;

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

375: /*
376:   FormExactSolution3 - Forms initial approximation.

378:   Input Parameters:
379:   da - The DM
380:   user - user-defined application context

382:   Output Parameter:
383:   X - vector
384:  */
385: PetscErrorCode FormExactSolution3(DM da, AppCtx *user, Vec U)
386: {
387:   DM             coordDA;
388:   Vec            coordinates;
389:   DMDACoor2d   **coords;
390:   PetscScalar  **u;
391:   PetscReal      x, y;
392:   PetscInt       xs, ys, xm, ym, i, j, m, n;

395:   m = PetscPowReal(2,user->mPar);
396:   n = PetscPowReal(2,user->nPar);

399:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
400:   DMGetCoordinateDM(da, &coordDA);
401:   DMGetCoordinates(da, &coordinates);
402:   DMDAVecGetArray(coordDA, coordinates, &coords);
403:   DMDAVecGetArray(da, U, &u);

405:   for (j = ys; j < ys+ym; ++j) {
406:     for (i = xs; i < xs+xm; ++i) {
407:       x = PetscRealPart(coords[j][i].x);
408:       y = PetscRealPart(coords[j][i].y);

410:       u[j][i] = PetscSinReal(m*PETSC_PI*x*(1-y))*PetscSinReal(n*PETSC_PI*y*(1-x));
411:     }
412:   }
413:   DMDAVecRestoreArray(da, U, &u);
414:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
415:   return(0);
416: }
417: /* ------------------------------------------------------------------- */

419: /*
420:   FormExactSolution4 - Forms initial approximation.

422:   Input Parameters:
423:   da - The DM
424:   user - user-defined application context

426:   Output Parameter:
427:   X - vector
428:  */
429: PetscErrorCode FormExactSolution4(DM da, AppCtx *user, Vec U)
430: {
431:   DM             coordDA;
432:   Vec            coordinates;
433:   DMDACoor2d   **coords;
434:   PetscScalar  **u;
435:   PetscReal      x, y, Lx, Ly;
436:   PetscInt       xs, ys, xm, ym, i, j;

440:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
441:   DMGetCoordinateDM(da, &coordDA);
442:   DMGetCoordinates(da, &coordinates);
443:   DMDAVecGetArray(coordDA, coordinates, &coords);
444:   DMDAVecGetArray(da, U, &u);

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

449:   for (j = ys; j < ys+ym; ++j) {
450:     for (i = xs; i < xs+xm; ++i) {
451:       x = PetscRealPart(coords[j][i].x);
452:       y = PetscRealPart(coords[j][i].y);
453:       u[j][i] = (PetscPowReal(x,4)-PetscSqr(Lx)*PetscSqr(x))*(PetscPowReal(y,4)-PetscSqr(Ly)*PetscSqr(y));
454:     }
455:   }
456:   DMDAVecRestoreArray(da, U, &u);
457:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
458:   return(0);
459: }
460: /* ------------------------------------------------------------------- */
461: /*
462:    FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch

465:  */
466: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
467: {
469:   PetscInt       i,j;
470:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc;
471:   PetscScalar    u,ue,uw,un,us,uxx,uyy;

474:   lambda = user->param;
475:   hx     = 1.0/(PetscReal)(info->mx-1);
476:   hy     = 1.0/(PetscReal)(info->my-1);
477:   sc     = hx*hy*lambda;
478:   hxdhy  = hx/hy;
479:   hydhx  = hy/hx;
480:   /*
481:      Compute function over the locally owned part of the grid
482:   */
483:   for (j=info->ys; j<info->ys+info->ym; j++) {
484:     for (i=info->xs; i<info->xs+info->xm; i++) {
485:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
486:         f[j][i] = 2.0*(hydhx+hxdhy)*x[j][i];
487:       } else {
488:         u  = x[j][i];
489:         uw = x[j][i-1];
490:         ue = x[j][i+1];
491:         un = x[j-1][i];
492:         us = x[j+1][i];

494:         if (i-1 == 0) uw = 0.;
495:         if (i+1 == info->mx-1) ue = 0.;
496:         if (j-1 == 0) un = 0.;
497:         if (j+1 == info->my-1) us = 0.;

499:         uxx     = (2.0*u - uw - ue)*hydhx;
500:         uyy     = (2.0*u - un - us)*hxdhy;
501:         f[j][i] = uxx + uyy - sc*PetscExpScalar(u);
502:       }
503:     }
504:   }
505:   PetscLogFlops(11.0*info->ym*info->xm);
506:   return(0);
507: }

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

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

514:  -laplacian u* - lambda exp(u*) = 2x(1-x) + 2y(1-y) - lambda exp(x(1-x)y(1-y))
515:
516:  Remark: the above is subtracted from the residual
517: -----------------------------------------------------------------------------------*/
518: PetscErrorCode FormFunctionLocalMMS1(DMDALocalInfo *info,PetscScalar **vx,PetscScalar **f,AppCtx *user)
519: {
521:   PetscInt       i,j;
522:   PetscReal      lambda,hx,hy,hxdhy,hydhx;
523:   PetscScalar    u,ue,uw,un,us,uxx,uyy;
524:   PetscReal      x,y;
525:   DM             coordDA;
526:   Vec            coordinates;
527:   DMDACoor2d   **coords;
528:   Vec            bcv = NULL;
529:   PetscScalar  **bcx = NULL;

532:   lambda = user->param;
533:   /* Extract coordinates */
534:   DMGetCoordinateDM(info->da, &coordDA);
535:   DMGetCoordinates(info->da, &coordinates);
536:   DMDAVecGetArray(coordDA, coordinates, &coords);
537:   hx     = info->xm > 1 ? PetscRealPart(coords[info->ys][info->xs+1].x) - PetscRealPart(coords[info->ys][info->xs].x) : 1.0;
538:   hy     = info->ym > 1 ? PetscRealPart(coords[info->ys+1][info->xs].y) - PetscRealPart(coords[info->ys][info->xs].y) : 1.0;
539:   hxdhy  = hx/hy;
540:   hydhx  = hy/hx;
541:   DMGetNamedLocalVector(info->da, "_petsc_boundary_conditions_", &bcv);
542:   DMDAVecGetArray(info->da, bcv, &bcx);
543:   /* Compute function over the locally owned part of the grid */
544:   for (j=info->ys; j<info->ys+info->ym; j++) {
545:     for (i=info->xs; i<info->xs+info->xm; i++) {
546:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
547:         f[j][i] = 2.0*(hydhx+hxdhy)*(vx[j][i] - bcx[j][i]);
548:       } else {
549:         x  = PetscRealPart(coords[j][i].x);
550:         y  = PetscRealPart(coords[j][i].y);
551:         u  = vx[j][i];
552:         uw = vx[j][i-1];
553:         ue = vx[j][i+1];
554:         un = vx[j-1][i];
555:         us = vx[j+1][i];

557:         if (i-1 == 0)          uw = bcx[j][i-1];
558:         if (i+1 == info->mx-1) ue = bcx[j][i+1];
559:         if (j-1 == 0)          un = bcx[j-1][i];
560:         if (j+1 == info->my-1) us = bcx[j+1][i];

562:         uxx     = (2.0*u - uw - ue)*hydhx;
563:         uyy     = (2.0*u - un - us)*hxdhy;
564:         f[j][i] = uxx + uyy - hx*hy*(lambda*PetscExpScalar(u) + 2*x*(1 - x) + 2*y*(1 - y) - lambda*PetscExpReal(x*(1 - x)*y*(1 - y)));
565:       }
566:     }
567:   }
568:   DMDAVecRestoreArray(info->da, bcv, &bcx);
569:   DMRestoreNamedLocalVector(info->da, "_petsc_boundary_conditions_", &bcv);
570:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
571:   PetscLogFlops(11.0*info->ym*info->xm);
572:   return(0);
573: }

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

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

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

582:  Remark: the above is subtracted from the residual
583: -----------------------------------------------------------------------------------*/
584: PetscErrorCode FormFunctionLocalMMS2(DMDALocalInfo *info,PetscScalar **vx,PetscScalar **f,AppCtx *user)
585: {
587:   PetscInt       i,j;
588:   PetscReal      lambda,hx,hy,hxdhy,hydhx;
589:   PetscScalar    u,ue,uw,un,us,uxx,uyy;
590:   PetscReal      x,y;
591:   DM             coordDA;
592:   Vec            coordinates;
593:   DMDACoor2d   **coords;
594:   Vec            bcv = NULL;
595:   PetscScalar  **bcx = NULL;

598:   lambda = user->param;
599:   /* Extract coordinates */
600:   DMGetCoordinateDM(info->da, &coordDA);
601:   DMGetCoordinates(info->da, &coordinates);
602:   DMDAVecGetArray(coordDA, coordinates, &coords);
603:   hx     = info->xm > 1 ? PetscRealPart(coords[info->ys][info->xs+1].x) - PetscRealPart(coords[info->ys][info->xs].x) : 1.0;
604:   hy     = info->ym > 1 ? PetscRealPart(coords[info->ys+1][info->xs].y) - PetscRealPart(coords[info->ys][info->xs].y) : 1.0;
605:   hxdhy  = hx/hy;
606:   hydhx  = hy/hx;
607:   DMGetNamedLocalVector(info->da, "_petsc_boundary_conditions_", &bcv);
608:   DMDAVecGetArray(info->da, bcv, &bcx);
609:   /* Compute function over the locally owned part of the grid */
610:   for (j=info->ys; j<info->ys+info->ym; j++) {
611:     for (i=info->xs; i<info->xs+info->xm; i++) {
612:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
613:         f[j][i] = 2.0*(hydhx+hxdhy)*(vx[j][i] - bcx[j][i]);
614:       } else {
615:         x  = PetscRealPart(coords[j][i].x);
616:         y  = PetscRealPart(coords[j][i].y);
617:         u  = vx[j][i];
618:         uw = vx[j][i-1];
619:         ue = vx[j][i+1];
620:         un = vx[j-1][i];
621:         us = vx[j+1][i];

623:         if (i-1 == 0)          uw = bcx[j][i-1];
624:         if (i+1 == info->mx-1) ue = bcx[j][i+1];
625:         if (j-1 == 0)          un = bcx[j-1][i];
626:         if (j+1 == info->my-1) us = bcx[j+1][i];

628:         uxx     = (2.0*u - uw - ue)*hydhx;
629:         uyy     = (2.0*u - un - us)*hxdhy;
630:         f[j][i] = uxx + uyy - hx*hy*(lambda*PetscExpScalar(u) + 2*PetscSqr(PETSC_PI)*PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y) - lambda*PetscExpReal(PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y)));
631:       }
632:     }
633:   }
634:   DMDAVecRestoreArray(info->da, bcv, &bcx);
635:   DMRestoreNamedLocalVector(info->da, "_petsc_boundary_conditions_", &bcv);
636:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
637:   PetscLogFlops(11.0*info->ym*info->xm);
638:   return(0);
639: }

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

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

646:  -laplacian u* - lambda exp(u*) = -(exp(Sin[m*Pi*x*(1 - y)]*Sin[n*Pi*(1 - x)*y])*lambda) +
647:                                   Pi^2*(-2*m*n*((-1 + x)*x + (-1 + y)*y)*Cos[m*Pi*x*(-1 + y)]*
648:                                   Cos[n*Pi*(-1 + x)*y] + (m^2*(x^2 + (-1 + y)^2) + n^2*((-1 + x)^2 + y^2))*
649:                                   Sin[m*Pi*x*(-1 + y)]*Sin[n*Pi*(-1 + x)*y])

651:   Remark: the above is subtracted from the residual
652: -----------------------------------------------------------------------------------*/
653: PetscErrorCode FormFunctionLocalMMS3(DMDALocalInfo *info,PetscScalar **vx,PetscScalar **f,AppCtx *user)
654: {
656:   PetscInt       i,j;
657:   PetscReal      lambda,hx,hy,hxdhy,hydhx,m,n;
658:   PetscScalar    u,ue,uw,un,us,uxx,uyy;
659:   PetscReal      x,y;
660:   DM             coordDA;
661:   Vec            coordinates;
662:   DMDACoor2d   **coords;

664:   m = PetscPowReal(2,user->mPar);
665:   n = PetscPowReal(2,user->nPar);

668:   lambda = user->param;
669:   hx     = 1.0/(PetscReal)(info->mx-1);
670:   hy     = 1.0/(PetscReal)(info->my-1);
671:   hxdhy  = hx/hy;
672:   hydhx  = hy/hx;
673:   /* Extract coordinates */
674:   DMGetCoordinateDM(info->da, &coordDA);
675:   DMGetCoordinates(info->da, &coordinates);
676:   DMDAVecGetArray(coordDA, coordinates, &coords);
677:   /* Compute function over the locally owned part of the grid */
678:   for (j=info->ys; j<info->ys+info->ym; j++) {
679:     for (i=info->xs; i<info->xs+info->xm; i++) {
680:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
681:         f[j][i] = 2.0*(hydhx+hxdhy)*vx[j][i];
682:       } else {
683:         x  = PetscRealPart(coords[j][i].x);
684:         y  = PetscRealPart(coords[j][i].y);
685:         u  = vx[j][i];
686:         uw = vx[j][i-1];
687:         ue = vx[j][i+1];
688:         un = vx[j-1][i];
689:         us = vx[j+1][i];

691:         if (i-1 == 0) uw = 0.;
692:         if (i+1 == info->mx-1) ue = 0.;
693:         if (j-1 == 0) un = 0.;
694:         if (j+1 == info->my-1) us = 0.;

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

699:         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)));
700:       }
701:     }
702:   }
703:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
704:   PetscLogFlops(11.0*info->ym*info->xm);
705:   return(0);

707: }

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

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

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

717:   Remark: the above is subtracted from the residual
718: -----------------------------------------------------------------------------------*/
719: PetscErrorCode FormFunctionLocalMMS4(DMDALocalInfo *info,PetscScalar **vx,PetscScalar **f,AppCtx *user)
720: {
722:   PetscInt       i,j;
723:   PetscReal      lambda,hx,hy,hxdhy,hydhx,Lx,Ly;
724:   PetscScalar    u,ue,uw,un,us,uxx,uyy;
725:   PetscReal      x,y;
726:   DM             coordDA;
727:   Vec            coordinates;
728:   DMDACoor2d   **coords;

731:   lambda = user->param;
732:   hx     = 1.0/(PetscReal)(info->mx-1);
733:   hy     = 1.0/(PetscReal)(info->my-1);
734:   hxdhy  = hx/hy;
735:   hydhx  = hy/hx;
736:   /* Extract coordinates */
737:   DMGetCoordinateDM(info->da, &coordDA);
738:   DMGetCoordinates(info->da, &coordinates);
739:   DMDAVecGetArray(coordDA, coordinates, &coords);

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

744:   /* Compute function over the locally owned part of the grid */
745:   for (j=info->ys; j<info->ys+info->ym; j++) {
746:     for (i=info->xs; i<info->xs+info->xm; i++) {
747:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
748:         f[j][i] = 2.0*(hydhx+hxdhy)*vx[j][i];
749:       } else {
750:         x  = PetscRealPart(coords[j][i].x);
751:         y  = PetscRealPart(coords[j][i].y);
752:         u  = vx[j][i];
753:         uw = vx[j][i-1];
754:         ue = vx[j][i+1];
755:         un = vx[j-1][i];
756:         us = vx[j+1][i];

758:         if (i-1 == 0) uw = 0.;
759:         if (i+1 == info->mx-1) ue = 0.;
760:         if (j-1 == 0) un = 0.;
761:         if (j+1 == info->my-1) us = 0.;

763:         uxx     = (2.0*u - uw - ue)*hydhx;
764:         uyy     = (2.0*u - un - us)*hxdhy;
765:         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))));
766:       }
767:     }
768:   }
769:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
770:   PetscLogFlops(11.0*info->ym*info->xm);
771:   return(0);

773: }

775: /* FormObjectiveLocal - Evaluates nonlinear function, F(x) on local process patch */
776: PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info,PetscScalar **x,PetscReal *obj,AppCtx *user)
777: {
779:   PetscInt       i,j;
780:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc,lobj=0;
781:   PetscScalar    u,ue,uw,un,us,uxux,uyuy;
782:   MPI_Comm       comm;

785:   *obj   = 0;
786:   PetscObjectGetComm((PetscObject)info->da,&comm);
787:   lambda = user->param;
788:   hx     = 1.0/(PetscReal)(info->mx-1);
789:   hy     = 1.0/(PetscReal)(info->my-1);
790:   sc     = hx*hy*lambda;
791:   hxdhy  = hx/hy;
792:   hydhx  = hy/hx;
793:   /*
794:      Compute function over the locally owned part of the grid
795:   */
796:   for (j=info->ys; j<info->ys+info->ym; j++) {
797:     for (i=info->xs; i<info->xs+info->xm; i++) {
798:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
799:         lobj += PetscRealPart((hydhx + hxdhy)*x[j][i]*x[j][i]);
800:       } else {
801:         u  = x[j][i];
802:         uw = x[j][i-1];
803:         ue = x[j][i+1];
804:         un = x[j-1][i];
805:         us = x[j+1][i];

807:         if (i-1 == 0) uw = 0.;
808:         if (i+1 == info->mx-1) ue = 0.;
809:         if (j-1 == 0) un = 0.;
810:         if (j+1 == info->my-1) us = 0.;

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

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

817:         lobj += PetscRealPart(0.5*(uxux + uyuy) - sc*PetscExpScalar(u));
818:       }
819:     }
820:   }
821:   PetscLogFlops(12.0*info->ym*info->xm);
822:   MPI_Allreduce(&lobj,obj,1,MPIU_REAL,MPIU_SUM,comm);
823:   return(0);
824: }

826: /*
827:    FormJacobianLocal - Evaluates Jacobian matrix on local process patch
828: */
829: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,Mat jacpre,AppCtx *user)
830: {
832:   PetscInt       i,j,k;
833:   MatStencil     col[5],row;
834:   PetscScalar    lambda,v[5],hx,hy,hxdhy,hydhx,sc;
835:   DM             coordDA;
836:   Vec            coordinates;
837:   DMDACoor2d   **coords;

840:   lambda = user->param;
841:   /* Extract coordinates */
842:   DMGetCoordinateDM(info->da, &coordDA);
843:   DMGetCoordinates(info->da, &coordinates);
844:   DMDAVecGetArray(coordDA, coordinates, &coords);
845:   hx     = info->xm > 1 ? PetscRealPart(coords[info->ys][info->xs+1].x) - PetscRealPart(coords[info->ys][info->xs].x) : 1.0;
846:   hy     = info->ym > 1 ? PetscRealPart(coords[info->ys+1][info->xs].y) - PetscRealPart(coords[info->ys][info->xs].y) : 1.0;
847:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
848:   hxdhy  = hx/hy;
849:   hydhx  = hy/hx;
850:   sc     = hx*hy*lambda;

853:   /*
854:      Compute entries for the locally owned part of the Jacobian.
855:       - Currently, all PETSc parallel matrix formats are partitioned by
856:         contiguous chunks of rows across the processors.
857:       - Each processor needs to insert only elements that it owns
858:         locally (but any non-local elements will be sent to the
859:         appropriate processor during matrix assembly).
860:       - Here, we set all entries for a particular row at once.
861:       - We can set matrix entries either using either
862:         MatSetValuesLocal() or MatSetValues(), as discussed above.
863:   */
864:   for (j=info->ys; j<info->ys+info->ym; j++) {
865:     for (i=info->xs; i<info->xs+info->xm; i++) {
866:       row.j = j; row.i = i;
867:       /* boundary points */
868:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
869:         v[0] =  2.0*(hydhx + hxdhy);
870:         MatSetValuesStencil(jacpre,1,&row,1,&row,v,INSERT_VALUES);
871:       } else {
872:         k = 0;
873:         /* interior grid points */
874:         if (j-1 != 0) {
875:           v[k]     = -hxdhy;
876:           col[k].j = j - 1; col[k].i = i;
877:           k++;
878:         }
879:         if (i-1 != 0) {
880:           v[k]     = -hydhx;
881:           col[k].j = j;     col[k].i = i-1;
882:           k++;
883:         }

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

887:         if (i+1 != info->mx-1) {
888:           v[k]     = -hydhx;
889:           col[k].j = j;     col[k].i = i+1;
890:           k++;
891:         }
892:         if (j+1 != info->mx-1) {
893:           v[k]     = -hxdhy;
894:           col[k].j = j + 1; col[k].i = i;
895:           k++;
896:         }
897:         MatSetValuesStencil(jacpre,1,&row,k,col,v,INSERT_VALUES);
898:       }
899:     }
900:   }

902:   /*
903:      Assemble matrix, using the 2-step process:
904:        MatAssemblyBegin(), MatAssemblyEnd().
905:   */
906:   MatAssemblyBegin(jacpre,MAT_FINAL_ASSEMBLY);
907:   MatAssemblyEnd(jacpre,MAT_FINAL_ASSEMBLY);
908:   /*
909:      Tell the matrix we will never add a new nonzero location to the
910:      matrix. If we do, it will generate an error.
911:   */
912:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
913:   return(0);
914: }

916: #if defined(PETSC_HAVE_MATLAB_ENGINE)
917: PetscErrorCode FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
918: {
919:   AppCtx         *user = (AppCtx*)ptr;
921:   PetscInt       Mx,My;
922:   PetscReal      lambda,hx,hy;
923:   Vec            localX,localF;
924:   MPI_Comm       comm;
925:   DM             da;

928:   SNESGetDM(snes,&da);
929:   DMGetLocalVector(da,&localX);
930:   DMGetLocalVector(da,&localF);
931:   PetscObjectSetName((PetscObject)localX,"localX");
932:   PetscObjectSetName((PetscObject)localF,"localF");
933:   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);

935:   lambda = user->param;
936:   hx     = 1.0/(PetscReal)(Mx-1);
937:   hy     = 1.0/(PetscReal)(My-1);

939:   PetscObjectGetComm((PetscObject)snes,&comm);
940:   /*
941:      Scatter ghost points to local vector,using the 2-step process
942:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
943:      By placing code between these two statements, computations can be
944:      done while messages are in transition.
945:   */
946:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
947:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
948:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
949:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
950:   PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);

952:   /*
953:      Insert values into global vector
954:   */
955:   DMLocalToGlobalBegin(da,localF,INSERT_VALUES,F);
956:   DMLocalToGlobalEnd(da,localF,INSERT_VALUES,F);
957:   DMRestoreLocalVector(da,&localX);
958:   DMRestoreLocalVector(da,&localF);
959:   return(0);
960: }
961: #endif

963: /* ------------------------------------------------------------------- */
964: /*
965:       Applies some sweeps on nonlinear Gauss-Seidel on each process

967:  */
968: PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx)
969: {
970:   PetscInt       i,j,k,Mx,My,xs,ys,xm,ym,its,tot_its,sweeps,l;
972:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc;
973:   PetscScalar    **x,**b,bij,F,F0=0,J,u,un,us,ue,eu,uw,uxx,uyy,y;
974:   PetscReal      atol,rtol,stol;
975:   DM             da;
976:   AppCtx         *user;
977:   Vec            localX,localB;

980:   tot_its = 0;
981:   SNESNGSGetSweeps(snes,&sweeps);
982:   SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);
983:   SNESGetDM(snes,&da);
984:   DMGetApplicationContext(da,(void**)&user);

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

988:   lambda = user->param;
989:   hx     = 1.0/(PetscReal)(Mx-1);
990:   hy     = 1.0/(PetscReal)(My-1);
991:   sc     = hx*hy*lambda;
992:   hxdhy  = hx/hy;
993:   hydhx  = hy/hx;

996:   DMGetLocalVector(da,&localX);
997:   if (B) {
998:     DMGetLocalVector(da,&localB);
999:   }
1000:   for (l=0; l<sweeps; l++) {

1002:     DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
1003:     DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
1004:     if (B) {
1005:       DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);
1006:       DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);
1007:     }
1008:     /*
1009:      Get a pointer to vector data.
1010:      - For default PETSc vectors, VecGetArray() returns a pointer to
1011:      the data array.  Otherwise, the routine is implementation dependent.
1012:      - You MUST call VecRestoreArray() when you no longer need access to
1013:      the array.
1014:      */
1015:     DMDAVecGetArray(da,localX,&x);
1016:     if (B) DMDAVecGetArray(da,localB,&b);
1017:     /*
1018:      Get local grid boundaries (for 2-dimensional DMDA):
1019:      xs, ys   - starting grid indices (no ghost points)
1020:      xm, ym   - widths of local grid (no ghost points)
1021:      */
1022:     DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

1024:     for (j=ys; j<ys+ym; j++) {
1025:       for (i=xs; i<xs+xm; i++) {
1026:         if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
1027:           /* boundary conditions are all zero Dirichlet */
1028:           x[j][i] = 0.0;
1029:         } else {
1030:           if (B) bij = b[j][i];
1031:           else   bij = 0.;

1033:           u  = x[j][i];
1034:           un = x[j-1][i];
1035:           us = x[j+1][i];
1036:           ue = x[j][i-1];
1037:           uw = x[j][i+1];

1039:           for (k=0; k<its; k++) {
1040:             eu  = PetscExpScalar(u);
1041:             uxx = (2.0*u - ue - uw)*hydhx;
1042:             uyy = (2.0*u - un - us)*hxdhy;
1043:             F   = uxx + uyy - sc*eu - bij;
1044:             if (k == 0) F0 = F;
1045:             J  = 2.0*(hydhx + hxdhy) - sc*eu;
1046:             y  = F/J;
1047:             u -= y;
1048:             tot_its++;

1050:             if (atol > PetscAbsReal(PetscRealPart(F)) ||
1051:                 rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) ||
1052:                 stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) {
1053:               break;
1054:             }
1055:           }
1056:           x[j][i] = u;
1057:         }
1058:       }
1059:     }
1060:     /*
1061:      Restore vector
1062:      */
1063:     DMDAVecRestoreArray(da,localX,&x);
1064:     DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
1065:     DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
1066:   }
1067:   PetscLogFlops(tot_its*(21.0));
1068:   DMRestoreLocalVector(da,&localX);
1069:   if (B) {
1070:     DMDAVecRestoreArray(da,localB,&b);
1071:     DMRestoreLocalVector(da,&localB);
1072:   }
1073:   return(0);
1074: }

1076: /*TEST
1077:   # ASM vs MSM
1078:   test:
1079:     suffix: asm_0
1080:     requires: !single
1081:     args: -mms 1 -par 0.0 -snes_monitor -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu
1082:   test:
1083:     suffix: msm_0
1084:     requires: !single
1085:     args: -mms 1 -par 0.0 -snes_monitor -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu
1086:   test:
1087:     suffix: asm_1
1088:     requires: !single
1089:     args: -mms 1 -par 0.0 -snes_monitor -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu -da_grid_x 8
1090:   test:
1091:     suffix: msm_1
1092:     requires: !single
1093:     args: -mms 1 -par 0.0 -snes_monitor -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu -da_grid_x 8
1094:   test:
1095:     suffix: asm_2
1096:     requires: !single
1097:     args: -mms 1 -par 0.0 -snes_monitor -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor -ksp_type richardson -pc_type asm -pc_asm_blocks 3 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu -da_grid_x 8
1098:   test:
1099:     suffix: msm_2
1100:     requires: !single
1101:     args: -mms 1 -par 0.0 -snes_monitor -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor -ksp_type richardson -pc_type asm -pc_asm_blocks 3 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu -da_grid_x 8
1102:   test:
1103:     suffix: asm_3
1104:     requires: !single
1105:     args: -mms 1 -par 0.0 -snes_monitor -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor -ksp_type richardson -pc_type asm -pc_asm_blocks 4 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu -da_grid_x 8
1106:   test:
1107:     suffix: msm_3
1108:     requires: !single
1109:     args: -mms 1 -par 0.0 -snes_monitor -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor -ksp_type richardson -pc_type asm -pc_asm_blocks 4 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu -da_grid_x 8
1110:   test:
1111:     suffix: asm_4
1112:     requires: !single
1113:     nsize: 2
1114:     args: -mms 1 -par 0.0 -snes_monitor -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu -da_grid_x 8
1115:   test:
1116:     suffix: msm_4
1117:     requires: !single
1118:     nsize: 2
1119:     args: -mms 1 -par 0.0 -snes_monitor -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu -da_grid_x 8
1120:   test:
1121:     suffix: asm_5
1122:     requires: !single
1123:     nsize: 2
1124:     args: -mms 1 -par 0.0 -snes_monitor -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor -ksp_type richardson -pc_type asm -pc_asm_blocks 4 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu -da_grid_x 8
1125:   test:
1126:     suffix: msm_5
1127:     requires: !single
1128:     nsize: 2
1129:     args: -mms 1 -par 0.0 -snes_monitor -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor -ksp_type richardson -pc_type asm -pc_asm_blocks 4 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu -da_grid_x 8

1131: TEST*/