Actual source code: ex5.c

petsc-master 2015-07-31
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 FormExactSolution1(DM,AppCtx*,Vec);
 70: extern PetscErrorCode FormFunctionLocalMMS1(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
 71: extern PetscErrorCode FormExactSolution2(DM,AppCtx*,Vec);
 72: extern PetscErrorCode FormFunctionLocalMMS2(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
 73: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,Mat,AppCtx*);
 74: extern PetscErrorCode FormObjectiveLocal(DMDALocalInfo*,PetscScalar**,PetscReal*,AppCtx*);
 75: #if defined(PETSC_HAVE_MATLAB_ENGINE)
 76: extern PetscErrorCode FormFunctionMatlab(SNES,Vec,Vec,void*);
 77: #endif
 78: extern PetscErrorCode NonlinearGS(SNES,Vec,Vec,void*);

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

101:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102:      Initialize program
103:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

107:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108:      Initialize problem parameters
109:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110:   user.param = 6.0;
111:   PetscOptionsGetReal(NULL,"-par",&user.param,NULL);
112:   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);
113:   PetscOptionsGetInt(NULL,"-mms",&MMS,NULL);

115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116:      Create nonlinear solver context
117:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118:   SNESCreate(PETSC_COMM_WORLD,&snes);
119:   SNESSetCountersReset(snes,PETSC_FALSE);
120:   SNESSetNGS(snes, NonlinearGS, NULL);

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

135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136:      Set local function evaluation routine
137:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138:   switch (MMS) {
139:   case 2:  DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocalMMS2,&user);break;
140:   case 1:  DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocalMMS1,&user);break;
141:   default: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocal,&user);
142:   }
143:   PetscOptionsGetBool(NULL,"-fd",&flg,NULL);
144:   if (!flg) {
145:     DMDASNESSetJacobianLocal(da,(DMDASNESJacobian)FormJacobianLocal,&user);
146:   }

148:   PetscOptionsGetBool(NULL,"-obj",&flg,NULL);
149:   if (flg) {
150:     DMDASNESSetObjectiveLocal(da,(DMDASNESObjective)FormObjectiveLocal,&user);
151:   }

153: #if defined(PETSC_HAVE_MATLAB_ENGINE)
154:   PetscOptionsGetBool(NULL,"-matlab_function",&matlab_function,0);
155:   if (matlab_function) {
156:     VecDuplicate(x,&r);
157:     SNESSetFunction(snes,r,FormFunctionMatlab,&user);
158:   }
159: #endif

161:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162:      Customize nonlinear solver; set runtime options
163:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
164:   SNESSetFromOptions(snes);

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

168:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169:      Solve nonlinear system
170:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
171:   SNESSolve(snes,NULL,x);
172:   SNESGetIterationNumber(snes,&its);

174:   SNESGetLinearSolveIterations(snes,&slits);
175:   SNESGetKSP(snes,&ksp);
176:   KSPGetTotalIterations(ksp,&lits);
177:   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);
178:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

181:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182:      If using MMS, check the l_2 error
183:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
184:   if (MMS) {
185:     Vec       e;
186:     PetscReal errorl2, errorinf;
187:     PetscInt  N;

189:     VecDuplicate(x, &e);
190:     PetscObjectViewFromOptions((PetscObject) x, NULL, "-sol_view");
191:     if (MMS == 1) {FormExactSolution1(da, &user, e);}
192:     else          {FormExactSolution2(da, &user, e);}
193:     PetscObjectViewFromOptions((PetscObject) e, NULL, "-exact_view");
194:     VecAXPY(e, -1.0, x);
195:     PetscObjectViewFromOptions((PetscObject) e, NULL, "-error_view");
196:     VecNorm(e, NORM_2, &errorl2);
197:     VecNorm(e, NORM_INFINITY, &errorinf);
198:     VecGetSize(e, &N);
199:     PetscPrintf(PETSC_COMM_WORLD, "N: %D error l2 %g inf %g\n", N, (double) errorl2/N, (double) errorinf);
200:     VecDestroy(&e);
201:   }

203:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204:      Free work space.  All PETSc objects should be destroyed when they
205:      are no longer needed.
206:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
207: #if defined(PETSC_HAVE_MATLAB_ENGINE)
208:   VecDestroy(&r);
209: #endif
210:   VecDestroy(&x);
211:   SNESDestroy(&snes);
212:   DMDestroy(&da);
213:   PetscFinalize();
214:   return 0;
215: }
216: /* ------------------------------------------------------------------- */
219: /*
220:    FormInitialGuess - Forms initial approximation.

222:    Input Parameters:
223:    da - The DM
224:    user - user-defined application context

226:    Output Parameter:
227:    X - vector
228:  */
229: PetscErrorCode FormInitialGuess(DM da,AppCtx *user,Vec X)
230: {
231:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
233:   PetscReal      lambda,temp1,temp,hx,hy;
234:   PetscScalar    **x;

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

239:   lambda = user->param;
240:   hx     = 1.0/(PetscReal)(Mx-1);
241:   hy     = 1.0/(PetscReal)(My-1);
242:   temp1  = lambda/(lambda + 1.0);

244:   /*
245:      Get a pointer to vector data.
246:        - For default PETSc vectors, VecGetArray() returns a pointer to
247:          the data array.  Otherwise, the routine is implementation dependent.
248:        - You MUST call VecRestoreArray() when you no longer need access to
249:          the array.
250:   */
251:   DMDAVecGetArray(da,X,&x);

253:   /*
254:      Get local grid boundaries (for 2-dimensional DMDA):
255:        xs, ys   - starting grid indices (no ghost points)
256:        xm, ym   - widths of local grid (no ghost points)

258:   */
259:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

261:   /*
262:      Compute initial guess over the locally owned part of the grid
263:   */
264:   for (j=ys; j<ys+ym; j++) {
265:     temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
266:     for (i=xs; i<xs+xm; i++) {
267:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
268:         /* boundary conditions are all zero Dirichlet */
269:         x[j][i] = 0.0;
270:       } else {
271:         x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
272:       }
273:     }
274:   }

276:   /*
277:      Restore vector
278:   */
279:   DMDAVecRestoreArray(da,X,&x);
280:   return(0);
281: }

285: /*
286:   FormExactSolution1 - Forms initial approximation.

288:   Input Parameters:
289:   da - The DM
290:   user - user-defined application context

292:   Output Parameter:
293:   X - vector
294:  */
295: PetscErrorCode FormExactSolution1(DM da, AppCtx *user, Vec U)
296: {
297:   DM             coordDA;
298:   Vec            coordinates;
299:   DMDACoor2d   **coords;
300:   PetscScalar  **u;
301:   PetscReal      x, y;
302:   PetscInt       xs, ys, xm, ym, i, j;

306:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
307:   DMGetCoordinateDM(da, &coordDA);
308:   DMGetCoordinates(da, &coordinates);
309:   DMDAVecGetArray(coordDA, coordinates, &coords);
310:   DMDAVecGetArray(da, U, &u);
311:   for (j = ys; j < ys+ym; ++j) {
312:     for (i = xs; i < xs+xm; ++i) {
313:       x = PetscRealPart(coords[j][i].x);
314:       y = PetscRealPart(coords[j][i].y);
315:       u[j][i] = x*(1 - x)*y*(1 - y);
316:     }
317:   }
318:   DMDAVecRestoreArray(da, U, &u);
319:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
320:   return(0);
321: }

325: /*
326:   FormExactSolution2 - Forms initial approximation.

328:   Input Parameters:
329:   da - The DM
330:   user - user-defined application context

332:   Output Parameter:
333:   X - vector
334:  */
335: PetscErrorCode FormExactSolution2(DM da, AppCtx *user, Vec U)
336: {
337:   DM             coordDA;
338:   Vec            coordinates;
339:   DMDACoor2d   **coords;
340:   PetscScalar  **u;
341:   PetscReal      x, y;
342:   PetscInt       xs, ys, xm, ym, i, j;

346:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
347:   DMGetCoordinateDM(da, &coordDA);
348:   DMGetCoordinates(da, &coordinates);
349:   DMDAVecGetArray(coordDA, coordinates, &coords);
350:   DMDAVecGetArray(da, U, &u);
351:   for (j = ys; j < ys+ym; ++j) {
352:     for (i = xs; i < xs+xm; ++i) {
353:       x = PetscRealPart(coords[j][i].x);
354:       y = PetscRealPart(coords[j][i].y);
355:       u[j][i] = PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y);
356:     }
357:   }
358:   DMDAVecRestoreArray(da, U, &u);
359:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
360:   return(0);
361: }
362: /* ------------------------------------------------------------------- */
365: /*
366:    FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch


369:  */
370: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
371: {
373:   PetscInt       i,j;
374:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc;
375:   PetscScalar    u,ue,uw,un,us,uxx,uyy;

378:   lambda = user->param;
379:   hx     = 1.0/(PetscReal)(info->mx-1);
380:   hy     = 1.0/(PetscReal)(info->my-1);
381:   sc     = hx*hy*lambda;
382:   hxdhy  = hx/hy;
383:   hydhx  = hy/hx;
384:   /*
385:      Compute function over the locally owned part of the grid
386:   */
387:   for (j=info->ys; j<info->ys+info->ym; j++) {
388:     for (i=info->xs; i<info->xs+info->xm; i++) {
389:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
390:         f[j][i] = 2.0*(hydhx+hxdhy)*x[j][i];
391:       } else {
392:         u  = x[j][i];
393:         uw = x[j][i-1];
394:         ue = x[j][i+1];
395:         un = x[j-1][i];
396:         us = x[j+1][i];

398:         if (i-1 == 0) uw = 0.;
399:         if (i+1 == info->mx-1) ue = 0.;
400:         if (j-1 == 0) un = 0.;
401:         if (j+1 == info->my-1) us = 0.;

403:         uxx     = (2.0*u - uw - ue)*hydhx;
404:         uyy     = (2.0*u - un - us)*hxdhy;
405:         f[j][i] = uxx + uyy - sc*PetscExpScalar(u);
406:       }
407:     }
408:   }
409:   PetscLogFlops(11.0*info->ym*info->xm);
410:   return(0);
411: }

415: /* FormFunctionLocalMMS1 - Evaluates nonlinear function, F(x) on local process patch */
416: PetscErrorCode FormFunctionLocalMMS1(DMDALocalInfo *info,PetscScalar **vx,PetscScalar **f,AppCtx *user)
417: {
419:   PetscInt       i,j;
420:   PetscReal      lambda,hx,hy,hxdhy,hydhx;
421:   PetscScalar    u,ue,uw,un,us,uxx,uyy;
422:   PetscReal      x,y;
423:   DM             coordDA;
424:   Vec            coordinates;
425:   DMDACoor2d   **coords;

428:   lambda = user->param;
429:   hx     = 1.0/(PetscReal)(info->mx-1);
430:   hy     = 1.0/(PetscReal)(info->my-1);
431:   hxdhy  = hx/hy;
432:   hydhx  = hy/hx;
433:   /* Extract coordinates */
434:   DMGetCoordinateDM(info->da, &coordDA);
435:   DMGetCoordinates(info->da, &coordinates);
436:   DMDAVecGetArray(coordDA, coordinates, &coords);
437:   /* Compute function over the locally owned part of the grid */
438:   for (j=info->ys; j<info->ys+info->ym; j++) {
439:     for (i=info->xs; i<info->xs+info->xm; i++) {
440:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
441:         f[j][i] = 2.0*(hydhx+hxdhy)*vx[j][i];
442:       } else {
443:         x  = PetscRealPart(coords[j][i].x);
444:         y  = PetscRealPart(coords[j][i].y);
445:         u  = vx[j][i];
446:         uw = vx[j][i-1];
447:         ue = vx[j][i+1];
448:         un = vx[j-1][i];
449:         us = vx[j+1][i];

451:         if (i-1 == 0) uw = 0.;
452:         if (i+1 == info->mx-1) ue = 0.;
453:         if (j-1 == 0) un = 0.;
454:         if (j+1 == info->my-1) us = 0.;

456:         uxx     = (2.0*u - uw - ue)*hydhx;
457:         uyy     = (2.0*u - un - us)*hxdhy;
458:         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)));
459:       }
460:     }
461:   }
462:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
463:   PetscLogFlops(11.0*info->ym*info->xm);
464:   return(0);
465: }

469: /* FormFunctionLocalMMS2 - Evaluates nonlinear function, F(x) on local process patch */
470: PetscErrorCode FormFunctionLocalMMS2(DMDALocalInfo *info,PetscScalar **vx,PetscScalar **f,AppCtx *user)
471: {
473:   PetscInt       i,j;
474:   PetscReal      lambda,hx,hy,hxdhy,hydhx;
475:   PetscScalar    u,ue,uw,un,us,uxx,uyy;
476:   PetscReal      x,y;
477:   DM             coordDA;
478:   Vec            coordinates;
479:   DMDACoor2d   **coords;

482:   lambda = user->param;
483:   hx     = 1.0/(PetscReal)(info->mx-1);
484:   hy     = 1.0/(PetscReal)(info->my-1);
485:   hxdhy  = hx/hy;
486:   hydhx  = hy/hx;
487:   /* Extract coordinates */
488:   DMGetCoordinateDM(info->da, &coordDA);
489:   DMGetCoordinates(info->da, &coordinates);
490:   DMDAVecGetArray(coordDA, coordinates, &coords);
491:   /* Compute function over the locally owned part of the grid */
492:   for (j=info->ys; j<info->ys+info->ym; j++) {
493:     for (i=info->xs; i<info->xs+info->xm; i++) {
494:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
495:         f[j][i] = 2.0*(hydhx+hxdhy)*vx[j][i];
496:       } else {
497:         x  = PetscRealPart(coords[j][i].x);
498:         y  = PetscRealPart(coords[j][i].y);
499:         u  = vx[j][i];
500:         uw = vx[j][i-1];
501:         ue = vx[j][i+1];
502:         un = vx[j-1][i];
503:         us = vx[j+1][i];

505:         if (i-1 == 0) uw = 0.;
506:         if (i+1 == info->mx-1) ue = 0.;
507:         if (j-1 == 0) un = 0.;
508:         if (j+1 == info->my-1) us = 0.;

510:         uxx     = (2.0*u - uw - ue)*hydhx;
511:         uyy     = (2.0*u - un - us)*hxdhy;
512:         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)));
513:       }
514:     }
515:   }
516:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
517:   PetscLogFlops(11.0*info->ym*info->xm);
518:   return(0);
519: }

523: /* FormObjectiveLocal - Evaluates nonlinear function, F(x) on local process patch */
524: PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info,PetscScalar **x,PetscReal *obj,AppCtx *user)
525: {
527:   PetscInt       i,j;
528:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc,lobj=0;
529:   PetscScalar    u,ue,uw,un,us,uxux,uyuy;
530:   MPI_Comm       comm;

533:   *obj   = 0;
534:   PetscObjectGetComm((PetscObject)info,&comm);
535:   lambda = user->param;
536:   hx     = 1.0/(PetscReal)(info->mx-1);
537:   hy     = 1.0/(PetscReal)(info->my-1);
538:   sc     = hx*hy*lambda;
539:   hxdhy  = hx/hy;
540:   hydhx  = hy/hx;
541:   /*
542:      Compute function over the locally owned part of the grid
543:   */
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:         lobj += PetscRealPart((hydhx + hxdhy)*x[j][i]*x[j][i]);
548:       } else {
549:         u  = x[j][i];
550:         uw = x[j][i-1];
551:         ue = x[j][i+1];
552:         un = x[j-1][i];
553:         us = x[j+1][i];

555:         if (i-1 == 0) uw = 0.;
556:         if (i+1 == info->mx-1) ue = 0.;
557:         if (j-1 == 0) un = 0.;
558:         if (j+1 == info->my-1) us = 0.;

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

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

565:         lobj += PetscRealPart(0.5*(uxux + uyuy) - sc*PetscExpScalar(u));
566:       }
567:     }
568:   }
569:   PetscLogFlops(12.0*info->ym*info->xm);
570:   MPI_Allreduce(&lobj,obj,1,MPIU_REAL,MPIU_SUM,comm);
571:   return(0);
572: }

576: /*
577:    FormJacobianLocal - Evaluates Jacobian matrix on local process patch
578: */
579: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,Mat jacpre,AppCtx *user)
580: {
582:   PetscInt       i,j,k;
583:   MatStencil     col[5],row;
584:   PetscScalar    lambda,v[5],hx,hy,hxdhy,hydhx,sc;

587:   lambda = user->param;
588:   hx     = 1.0/(PetscReal)(info->mx-1);
589:   hy     = 1.0/(PetscReal)(info->my-1);
590:   sc     = hx*hy*lambda;
591:   hxdhy  = hx/hy;
592:   hydhx  = hy/hx;


595:   /*
596:      Compute entries for the locally owned part of the Jacobian.
597:       - Currently, all PETSc parallel matrix formats are partitioned by
598:         contiguous chunks of rows across the processors.
599:       - Each processor needs to insert only elements that it owns
600:         locally (but any non-local elements will be sent to the
601:         appropriate processor during matrix assembly).
602:       - Here, we set all entries for a particular row at once.
603:       - We can set matrix entries either using either
604:         MatSetValuesLocal() or MatSetValues(), as discussed above.
605:   */
606:   for (j=info->ys; j<info->ys+info->ym; j++) {
607:     for (i=info->xs; i<info->xs+info->xm; i++) {
608:       row.j = j; row.i = i;
609:       /* boundary points */
610:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
611:         v[0] =  2.0*(hydhx + hxdhy);
612:         MatSetValuesStencil(jacpre,1,&row,1,&row,v,INSERT_VALUES);
613:       } else {
614:         k = 0;
615:         /* interior grid points */
616:         if (j-1 != 0) {
617:           v[k]     = -hxdhy;
618:           col[k].j = j - 1; col[k].i = i;
619:           k++;
620:         }
621:         if (i-1 != 0) {
622:           v[k]     = -hydhx;
623:           col[k].j = j;     col[k].i = i-1;
624:           k++;
625:         }

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

629:         if (i+1 != info->mx-1) {
630:           v[k]     = -hydhx;
631:           col[k].j = j;     col[k].i = i+1;
632:           k++;
633:         }
634:         if (j+1 != info->mx-1) {
635:           v[k]     = -hxdhy;
636:           col[k].j = j + 1; col[k].i = i;
637:           k++;
638:         }
639:         MatSetValuesStencil(jacpre,1,&row,k,col,v,INSERT_VALUES);
640:       }
641:     }
642:   }

644:   /*
645:      Assemble matrix, using the 2-step process:
646:        MatAssemblyBegin(), MatAssemblyEnd().
647:   */
648:   MatAssemblyBegin(jacpre,MAT_FINAL_ASSEMBLY);
649:   MatAssemblyEnd(jacpre,MAT_FINAL_ASSEMBLY);
650:   /*
651:      Tell the matrix we will never add a new nonzero location to the
652:      matrix. If we do, it will generate an error.
653:   */
654:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
655:   return(0);
656: }

658: #if defined(PETSC_HAVE_MATLAB_ENGINE)
661: PetscErrorCode FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
662: {
663:   AppCtx         *user = (AppCtx*)ptr;
665:   PetscInt       Mx,My;
666:   PetscReal      lambda,hx,hy;
667:   Vec            localX,localF;
668:   MPI_Comm       comm;
669:   DM             da;

672:   SNESGetDM(snes,&da);
673:   DMGetLocalVector(da,&localX);
674:   DMGetLocalVector(da,&localF);
675:   PetscObjectSetName((PetscObject)localX,"localX");
676:   PetscObjectSetName((PetscObject)localF,"localF");
677:   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);

679:   lambda = user->param;
680:   hx     = 1.0/(PetscReal)(Mx-1);
681:   hy     = 1.0/(PetscReal)(My-1);

683:   PetscObjectGetComm((PetscObject)snes,&comm);
684:   /*
685:      Scatter ghost points to local vector,using the 2-step process
686:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
687:      By placing code between these two statements, computations can be
688:      done while messages are in transition.
689:   */
690:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
691:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
692:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
693:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
694:   PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);

696:   /*
697:      Insert values into global vector
698:   */
699:   DMLocalToGlobalBegin(da,localF,INSERT_VALUES,F);
700:   DMLocalToGlobalEnd(da,localF,INSERT_VALUES,F);
701:   DMRestoreLocalVector(da,&localX);
702:   DMRestoreLocalVector(da,&localF);
703:   return(0);
704: }
705: #endif

707: /* ------------------------------------------------------------------- */
710: /*
711:       Applies some sweeps on nonlinear Gauss-Seidel on each process

713:  */
714: PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx)
715: {
716:   PetscInt       i,j,k,Mx,My,xs,ys,xm,ym,its,tot_its,sweeps,l;
718:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc;
719:   PetscScalar    **x,**b,bij,F,F0=0,J,u,un,us,ue,eu,uw,uxx,uyy,y;
720:   PetscReal      atol,rtol,stol;
721:   DM             da;
722:   AppCtx         *user;
723:   Vec            localX,localB;

726:   tot_its = 0;
727:   SNESNGSGetSweeps(snes,&sweeps);
728:   SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);
729:   SNESGetDM(snes,&da);
730:   DMGetApplicationContext(da,(void**)&user);

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

734:   lambda = user->param;
735:   hx     = 1.0/(PetscReal)(Mx-1);
736:   hy     = 1.0/(PetscReal)(My-1);
737:   sc     = hx*hy*lambda;
738:   hxdhy  = hx/hy;
739:   hydhx  = hy/hx;


742:   DMGetLocalVector(da,&localX);
743:   if (B) {
744:     DMGetLocalVector(da,&localB);
745:   }
746:   for (l=0; l<sweeps; l++) {

748:     DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
749:     DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
750:     if (B) {
751:       DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);
752:       DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);
753:     }
754:     /*
755:      Get a pointer to vector data.
756:      - For default PETSc vectors, VecGetArray() returns a pointer to
757:      the data array.  Otherwise, the routine is implementation dependent.
758:      - You MUST call VecRestoreArray() when you no longer need access to
759:      the array.
760:      */
761:     DMDAVecGetArray(da,localX,&x);
762:     if (B) DMDAVecGetArray(da,localB,&b);
763:     /*
764:      Get local grid boundaries (for 2-dimensional DMDA):
765:      xs, ys   - starting grid indices (no ghost points)
766:      xm, ym   - widths of local grid (no ghost points)
767:      */
768:     DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

770:     for (j=ys; j<ys+ym; j++) {
771:       for (i=xs; i<xs+xm; i++) {
772:         if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
773:           /* boundary conditions are all zero Dirichlet */
774:           x[j][i] = 0.0;
775:         } else {
776:           if (B) bij = b[j][i];
777:           else   bij = 0.;

779:           u  = x[j][i];
780:           un = x[j-1][i];
781:           us = x[j+1][i];
782:           ue = x[j][i-1];
783:           uw = x[j][i+1];

785:           for (k=0; k<its; k++) {
786:             eu  = PetscExpScalar(u);
787:             uxx = (2.0*u - ue - uw)*hydhx;
788:             uyy = (2.0*u - un - us)*hxdhy;
789:             F   = uxx + uyy - sc*eu - bij;
790:             if (k == 0) F0 = F;
791:             J  = 2.0*(hydhx + hxdhy) - sc*eu;
792:             y  = F/J;
793:             u -= y;
794:             tot_its++;

796:             if (atol > PetscAbsReal(PetscRealPart(F)) ||
797:                 rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) ||
798:                 stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) {
799:               break;
800:             }
801:           }
802:           x[j][i] = u;
803:         }
804:       }
805:     }
806:     /*
807:      Restore vector
808:      */
809:     DMDAVecRestoreArray(da,localX,&x);
810:     DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
811:     DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
812:   }
813:   PetscLogFlops(tot_its*(21.0));
814:   DMRestoreLocalVector(da,&localX);
815:   if (B) {
816:     DMDAVecRestoreArray(da,localB,&b);
817:     DMRestoreLocalVector(da,&localB);
818:   }
819:   return(0);
820: }