Actual source code: ex15.c

  1: static const char help[] = "p-Bratu nonlinear PDE in 2d.\n\
2: We solve the  p-Laplacian (nonlinear diffusion) combined with\n\
3: the Bratu (solid fuel ignition) nonlinearity in a 2D rectangular\n\
4: domain, using distributed arrays (DAs) to partition the parallel grid.\n\
5: The command line options include:\n\
6:   -p <2>: p' in p-Laplacian term\n\
7:   -epsilon <1e-05>: Strain-regularization in p-Laplacian\n\
8:   -lambda <6>: Bratu parameter\n\
9:   -blocks <bx,by>: number of coefficient interfaces in x and y direction\n\
10:   -kappa <1e-3>: diffusivity in odd regions\n\
11: \n";


The $p$-Bratu problem is a combination of the $p$-Laplacian (nonlinear diffusion) and the Brutu solid fuel ignition problem.
This problem is modeled by the partial differential equation

\begin{equation*}
-\nabla\cdot (\eta \nabla u) - \lambda \exp(u) = 0
\end{equation*}

on $\Omega = (-1,1)^2$ with closure

\begin{align*}
\eta(\gamma) &= (\epsilon^2 + \gamma)^{(p-2)/2} & \gamma &= \frac 1 2 |\nabla u|^2
\end{align*}

and boundary conditions $u = 0$ for $(x,y) \in \partial \Omega$

A 9-point finite difference stencil is used to discretize
the boundary value problem to obtain a nonlinear system of equations.
This would be a 5-point stencil if not for the $p$-Laplacian's nonlinearity.
 34: /*
35:       mpiexec -n 2 ./ex15 -snes_monitor -ksp_monitor log_summary
36: */

38: /*
39:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
40:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
41:    file automatically includes:
42:      petsc.h       - base PETSc routines   petscvec.h - vectors
43:      petscsys.h    - system routines       petscmat.h - matrices
44:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
45:      petscviewer.h - viewers               petscpc.h  - preconditioners
46:      petscksp.h   - linear solvers
47: */

49: #include <petscdm.h>
50: #include <petscdmda.h>
51: #include <petscsnes.h>

53: typedef enum {JAC_BRATU,JAC_PICARD,JAC_STAR,JAC_NEWTON} JacType;
54: static const char *const JacTypes[] = {"BRATU","PICARD","STAR","NEWTON","JacType","JAC_",0};

56: /*
57:    User-defined application context - contains data needed by the
58:    application-provided call-back routines, FormJacobianLocal() and
59:    FormFunctionLocal().
60: */
61: typedef struct {
62:   PetscReal   lambda;         /* Bratu parameter */
63:   PetscReal   p;              /* Exponent in p-Laplacian */
64:   PetscReal   epsilon;        /* Regularization */
65:   PetscReal   source;         /* Source term */
66:   JacType     jtype;          /* What type of Jacobian to assemble */
67:   PetscBool   picard;
68:   PetscInt    blocks[2];
69:   PetscReal   kappa;
70:   PetscInt    initial;        /* initial conditions type */
71: } AppCtx;

73: /*
74:    User-defined routines
75: */
76: static PetscErrorCode FormRHS(AppCtx*,DM,Vec);
77: static PetscErrorCode FormInitialGuess(AppCtx*,DM,Vec);
78: static PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
79: static PetscErrorCode FormFunctionPicardLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
80: static PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,Mat,AppCtx*);
81: static PetscErrorCode NonlinearGS(SNES,Vec,Vec,void*);

83: typedef struct _n_PreCheck *PreCheck;
84: struct _n_PreCheck {
85:   MPI_Comm    comm;
86:   PetscReal   angle;
87:   Vec         Ylast;
88:   PetscViewer monitor;
89: };
90: PetscErrorCode PreCheckCreate(MPI_Comm,PreCheck*);
91: PetscErrorCode PreCheckDestroy(PreCheck*);
92: PetscErrorCode PreCheckFunction(SNESLineSearch,Vec,Vec,PetscBool*,void*);
93: PetscErrorCode PreCheckSetFromOptions(PreCheck);

95: int main(int argc,char **argv)
96: {
97:   SNES                snes;                    /* nonlinear solver */
98:   Vec                 x,r,b;                   /* solution, residual, rhs vectors */
99:   AppCtx              user;                    /* user-defined work context */
100:   PetscInt            its;                     /* iterations for convergence */
101:   SNESConvergedReason reason;                  /* Check convergence */
102:   PetscBool           alloc_star;              /* Only allocate for the STAR stencil  */
103:   PetscBool           write_output;
104:   char                filename[PETSC_MAX_PATH_LEN] = "ex15.vts";
105:   PetscReal           bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
106:   DM                  da;                      /* distributed array data structure */
107:   PreCheck            precheck = NULL;         /* precheck context for version in this file */
108:   PetscInt            use_precheck;            /* 0=none, 1=version in this file, 2=SNES-provided version */
109:   PetscReal           precheck_angle;          /* When manually setting the SNES-provided precheck function */
110:   PetscErrorCode      ierr;
111:   SNESLineSearch      linesearch;

113:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114:      Initialize program
115:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

119:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120:      Initialize problem parameters
121:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122:   user.lambda    = 0.0; user.p = 2.0; user.epsilon = 1e-5; user.source = 0.1; user.jtype = JAC_NEWTON;user.initial=-1;
123:   user.blocks[0] = 1; user.blocks[1] = 1; user.kappa = 1e-3;
124:   alloc_star     = PETSC_FALSE;
125:   use_precheck   = 0; precheck_angle = 10.;
126:   user.picard    = PETSC_FALSE;
127:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"p-Bratu options",__FILE__);
128:   {
129:     PetscInt two=2;
130:     PetscOptionsReal("-lambda","Bratu parameter","",user.lambda,&user.lambda,NULL);
131:     PetscOptionsReal("-p","Exponent p' in p-Laplacian","",user.p,&user.p,NULL);
132:     PetscOptionsReal("-epsilon","Strain-regularization in p-Laplacian","",user.epsilon,&user.epsilon,NULL);
133:     PetscOptionsReal("-source","Constant source term","",user.source,&user.source,NULL);
134:     PetscOptionsEnum("-jtype","Jacobian approximation to assemble","",JacTypes,(PetscEnum)user.jtype,(PetscEnum*)&user.jtype,NULL);
135:     PetscOptionsName("-picard","Solve with defect-correction Picard iteration","",&user.picard);
136:     if (user.picard) {
137:       user.jtype = JAC_PICARD;
139:       /* the Picard linearization only requires a 5 point stencil, while the Newton linearization requires a 9 point stencil */
140:       /* hence allocating the 5 point stencil gives the same convergence as the 9 point stencil since the extra stencil points are not used */
141:       PetscOptionsBool("-alloc_star","Allocate for STAR stencil (5-point)","",alloc_star,&alloc_star,NULL);
142:     }
143:     PetscOptionsInt("-precheck","Use a pre-check correction intended for use with Picard iteration 1=this version, 2=library","",use_precheck,&use_precheck,NULL);
144:     PetscOptionsInt("-initial","Initial conditions type (-1: default, 0: zero-valued, 1: peaked guess)","",user.initial,&user.initial,NULL);
145:     if (use_precheck == 2) {    /* Using library version, get the angle */
146:       PetscOptionsReal("-precheck_angle","Angle in degrees between successive search directions necessary to activate step correction","",precheck_angle,&precheck_angle,NULL);
147:     }
148:     PetscOptionsIntArray("-blocks","number of coefficient interfaces in x and y direction","",user.blocks,&two,NULL);
149:     PetscOptionsReal("-kappa","diffusivity in odd regions","",user.kappa,&user.kappa,NULL);
150:     PetscOptionsString("-o","Output solution in vts format","",filename,filename,sizeof(filename),&write_output);
151:   }
152:   PetscOptionsEnd();
153:   if (user.lambda > bratu_lambda_max || user.lambda < bratu_lambda_min) {
154:     PetscPrintf(PETSC_COMM_WORLD,"WARNING: lambda %g out of range for p=2\n",(double)user.lambda);
155:   }

157:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158:      Create nonlinear solver context
159:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160:   SNESCreate(PETSC_COMM_WORLD,&snes);

162:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163:      Create distributed array (DMDA) to manage parallel grid and vectors
164:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,alloc_star ? DMDA_STENCIL_STAR : DMDA_STENCIL_BOX,4,4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
166:   DMSetFromOptions(da);
167:   DMSetUp(da);

169:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
170:      Extract global vectors from DM; then duplicate for remaining
171:      vectors that are the same types
172:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173:   DMCreateGlobalVector(da,&x);
174:   VecDuplicate(x,&r);
175:   VecDuplicate(x,&b);

177:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178:      User can override with:
179:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
180:                 (unless user explicitly sets preconditioner)
181:      -snes_mf_operator : form preconditioning matrix as set by the user,
182:                          but use matrix-free approx for Jacobian-vector
183:                          products within Newton-Krylov method

185:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

187:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188:      Set local function evaluation routine
189:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
190:   DMSetApplicationContext(da, &user);
191:   SNESSetDM(snes,da);
192:   if (user.picard) {
193:     /*
194:         This is not really right requiring the user to call SNESSetFunction/Jacobian but the DMDASNESSetPicardLocal() cannot access
195:         the SNES to set it
196:     */
197:     DMDASNESSetPicardLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionPicardLocal,
198:                                   (PetscErrorCode (*)(DMDALocalInfo*,void*,Mat,Mat,void*))FormJacobianLocal,&user);
199:     SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user);
200:     SNESSetJacobian(snes,NULL,NULL,SNESPicardComputeJacobian,&user);
201:   } else {
202:     DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);
203:     DMDASNESSetJacobianLocal(da,(PetscErrorCode (*)(DMDALocalInfo*,void*,Mat,Mat,void*))FormJacobianLocal,&user);
204:   }

206:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
207:      Customize nonlinear solver; set runtime options
208:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
209:   SNESSetFromOptions(snes);
210:   SNESSetNGS(snes,NonlinearGS,&user);
211:   SNESGetLineSearch(snes, &linesearch);
212:   /* Set up the precheck context if requested */
213:   if (use_precheck == 1) {      /* Use the precheck routines in this file */
214:     PreCheckCreate(PETSC_COMM_WORLD,&precheck);
215:     PreCheckSetFromOptions(precheck);
216:     SNESLineSearchSetPreCheck(linesearch,PreCheckFunction,precheck);
217:   } else if (use_precheck == 2) { /* Use the version provided by the library */
218:     SNESLineSearchSetPreCheck(linesearch,SNESLineSearchPreCheckPicard,&precheck_angle);
219:   }

221:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
222:      Evaluate initial guess
223:      Note: The user should initialize the vector, x, with the initial guess
224:      for the nonlinear solver prior to calling SNESSolve().  In particular,
225:      to employ an initial guess of zero, the user should explicitly set
226:      this vector to zero by calling VecSet().
227:   */

229:   FormInitialGuess(&user,da,x);
230:   FormRHS(&user,da,b);

232:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
233:      Solve nonlinear system
234:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
235:   SNESSolve(snes,b,x);
236:   SNESGetIterationNumber(snes,&its);
237:   SNESGetConvergedReason(snes,&reason);

239:   PetscPrintf(PETSC_COMM_WORLD,"%s Number of nonlinear iterations = %D\n",SNESConvergedReasons[reason],its);

241:   if (write_output) {
242:     PetscViewer viewer;
243:     PetscViewerVTKOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer);
244:     VecView(x,viewer);
245:     PetscViewerDestroy(&viewer);
246:   }

248:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
249:      Free work space.  All PETSc objects should be destroyed when they
250:      are no longer needed.
251:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

253:   VecDestroy(&x);
254:   VecDestroy(&r);
255:   VecDestroy(&b);
256:   SNESDestroy(&snes);
257:   DMDestroy(&da);
258:   PreCheckDestroy(&precheck);
259:   PetscFinalize();
260:   return 0;
261: }

263: /* ------------------------------------------------------------------- */
264: /*
265:    FormInitialGuess - Forms initial approximation.

267:    Input Parameters:
268:    user - user-defined application context
269:    X - vector

271:    Output Parameter:
272:    X - vector
273:  */
274: static PetscErrorCode FormInitialGuess(AppCtx *user,DM da,Vec X)
275: {
276:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
277:   PetscReal      temp1,temp,hx,hy;
278:   PetscScalar    **x;

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

283:   hx    = 1.0/(PetscReal)(Mx-1);
284:   hy    = 1.0/(PetscReal)(My-1);
285:   temp1 = user->lambda / (user->lambda + 1.);

287:   /*
288:      Get a pointer to vector data.
289:        - For default PETSc vectors, VecGetArray() returns a pointer to
290:          the data array.  Otherwise, the routine is implementation dependent.
291:        - You MUST call VecRestoreArray() when you no longer need access to
292:          the array.
293:   */
294:   DMDAVecGetArray(da,X,&x);

296:   /*
297:      Get local grid boundaries (for 2-dimensional DA):
298:        xs, ys   - starting grid indices (no ghost points)
299:        xm, ym   - widths of local grid (no ghost points)

301:   */
302:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

304:   /*
305:      Compute initial guess over the locally owned part of the grid
306:   */
307:   for (j=ys; j<ys+ym; j++) {
308:     temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
309:     for (i=xs; i<xs+xm; i++) {
310:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
311:         /* boundary conditions are all zero Dirichlet */
312:         x[j][i] = 0.0;
313:       } else {
314:         if (user->initial == -1) {
315:           if (user->lambda != 0) {
316:             x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
317:           } else {
318:             /* The solution above is an exact solution for lambda=0, this avoids "accidentally" starting
319:              * with an exact solution. */
320:             const PetscReal
321:               xx = 2*(PetscReal)i/(Mx-1) - 1,
322:               yy = 2*(PetscReal)j/(My-1) - 1;
323:             x[j][i] = (1 - xx*xx) * (1-yy*yy) * xx * yy;
324:           }
325:         } else if (user->initial == 0) {
326:           x[j][i] = 0.;
327:         } else if (user->initial == 1) {
328:           const PetscReal
329:             xx = 2*(PetscReal)i/(Mx-1) - 1,
330:             yy = 2*(PetscReal)j/(My-1) - 1;
331:           x[j][i] = (1 - xx*xx) * (1-yy*yy) * xx * yy;
332:         } else {
333:           if (user->lambda != 0) {
334:             x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
335:           } else {
336:             x[j][i] = 0.5*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
337:           }
338:         }
339:       }
340:     }
341:   }
342:   /*
343:      Restore vector
344:   */
345:   DMDAVecRestoreArray(da,X,&x);
346:   return 0;
347: }

349: /*
350:    FormRHS - Forms constant RHS for the problem.

352:    Input Parameters:
353:    user - user-defined application context
354:    B - RHS vector

356:    Output Parameter:
357:    B - vector
358:  */
359: static PetscErrorCode FormRHS(AppCtx *user,DM da,Vec B)
360: {
361:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
362:   PetscReal      hx,hy;
363:   PetscScalar    **b;

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

368:   hx    = 1.0/(PetscReal)(Mx-1);
369:   hy    = 1.0/(PetscReal)(My-1);
370:   DMDAVecGetArray(da,B,&b);
371:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
372:   for (j=ys; j<ys+ym; j++) {
373:     for (i=xs; i<xs+xm; i++) {
374:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
375:         b[j][i] = 0.0;
376:       } else {
377:         b[j][i] = hx*hy*user->source;
378:       }
379:     }
380:   }
381:   DMDAVecRestoreArray(da,B,&b);
382:   return 0;
383: }

385: static inline PetscReal kappa(const AppCtx *ctx,PetscReal x,PetscReal y)
386: {
387:   return (((PetscInt)(x*ctx->blocks[0])) + ((PetscInt)(y*ctx->blocks[1]))) % 2 ? ctx->kappa : 1.0;
388: }
389: /* p-Laplacian diffusivity */
390: static inline PetscScalar eta(const AppCtx *ctx,PetscReal x,PetscReal y,PetscScalar ux,PetscScalar uy)
391: {
392:   return kappa(ctx,x,y) * PetscPowScalar(PetscSqr(ctx->epsilon)+0.5*(ux*ux + uy*uy),0.5*(ctx->p-2.));
393: }
394: static inline PetscScalar deta(const AppCtx *ctx,PetscReal x,PetscReal y,PetscScalar ux,PetscScalar uy)
395: {
396:   return (ctx->p == 2)
397:          ? 0
398:          : kappa(ctx,x,y)*PetscPowScalar(PetscSqr(ctx->epsilon)+0.5*(ux*ux + uy*uy),0.5*(ctx->p-4)) * 0.5 * (ctx->p-2.);
399: }

401: /* ------------------------------------------------------------------- */
402: /*
403:    FormFunctionLocal - Evaluates nonlinear function, F(x).
404:  */
405: static PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
406: {
407:   PetscReal      hx,hy,dhx,dhy,sc;
408:   PetscInt       i,j;
409:   PetscScalar    eu;

412:   hx     = 1.0/(PetscReal)(info->mx-1);
413:   hy     = 1.0/(PetscReal)(info->my-1);
414:   sc     = hx*hy*user->lambda;
415:   dhx    = 1/hx;
416:   dhy    = 1/hy;
417:   /*
418:      Compute function over the locally owned part of the grid
419:   */
420:   for (j=info->ys; j<info->ys+info->ym; j++) {
421:     for (i=info->xs; i<info->xs+info->xm; i++) {
422:       PetscReal xx = i*hx,yy = j*hy;
423:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
424:         f[j][i] = x[j][i];
425:       } else {
426:         const PetscScalar
427:           u    = x[j][i],
428:           ux_E = dhx*(x[j][i+1]-x[j][i]),
429:           uy_E = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]),
430:           ux_W = dhx*(x[j][i]-x[j][i-1]),
431:           uy_W = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]),
432:           ux_N = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]),
433:           uy_N = dhy*(x[j+1][i]-x[j][i]),
434:           ux_S = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]),
435:           uy_S = dhy*(x[j][i]-x[j-1][i]),
436:           e_E  = eta(user,xx,yy,ux_E,uy_E),
437:           e_W  = eta(user,xx,yy,ux_W,uy_W),
438:           e_N  = eta(user,xx,yy,ux_N,uy_N),
439:           e_S  = eta(user,xx,yy,ux_S,uy_S),
440:           uxx  = -hy * (e_E*ux_E - e_W*ux_W),
441:           uyy  = -hx * (e_N*uy_N - e_S*uy_S);
442:         if (sc) eu = PetscExpScalar(u);
443:         else    eu = 0.;
444:         /* For p=2, these terms decay to:
445:          uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx
446:          uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy
447:         */
448:         f[j][i] = uxx + uyy - sc*eu;
449:       }
450:     }
451:   }
452:   PetscLogFlops(info->xm*info->ym*(72.0));
453:   return 0;
454: }

456: /*
457:     This is the opposite sign of the part of FormFunctionLocal that excludes the A(x) x part of the operation,
458:     that is FormFunction applies A(x) x - b(x) while this applies b(x) because for Picard we think of it as solving A(x) x = b(x)

460: */
461: static PetscErrorCode FormFunctionPicardLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
462: {
463:   PetscReal hx,hy,sc;
464:   PetscInt  i,j;

467:   hx     = 1.0/(PetscReal)(info->mx-1);
468:   hy     = 1.0/(PetscReal)(info->my-1);
469:   sc     = hx*hy*user->lambda;
470:   /*
471:      Compute function over the locally owned part of the grid
472:   */
473:   for (j=info->ys; j<info->ys+info->ym; j++) {
474:     for (i=info->xs; i<info->xs+info->xm; i++) {
475:       if (!(i == 0 || j == 0 || i == info->mx-1 || j == info->my-1)) {
476:         const PetscScalar u = x[j][i];
477:         f[j][i] = sc*PetscExpScalar(u);
478:       } else {
479:         f[j][i] = 0.0; /* this is zero because the A(x) x term forces the x to be zero on the boundary */
480:       }
481:     }
482:   }
483:   PetscLogFlops(info->xm*info->ym*2.0);
484:   return 0;
485: }

487: /*
488:    FormJacobianLocal - Evaluates Jacobian matrix.
489: */
490: static PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat J,Mat B,AppCtx *user)
491: {
492:   PetscInt       i,j;
493:   MatStencil     col[9],row;
494:   PetscScalar    v[9];
495:   PetscReal      hx,hy,hxdhy,hydhx,dhx,dhy,sc;

498:   MatZeroEntries(B);
499:   hx    = 1.0/(PetscReal)(info->mx-1);
500:   hy    = 1.0/(PetscReal)(info->my-1);
501:   sc    = hx*hy*user->lambda;
502:   dhx   = 1/hx;
503:   dhy   = 1/hy;
504:   hxdhy = hx/hy;
505:   hydhx = hy/hx;

507:   /*
508:      Compute entries for the locally owned part of the Jacobian.
509:       - PETSc parallel matrix formats are partitioned by
510:         contiguous chunks of rows across the processors.
511:       - Each processor needs to insert only elements that it owns
512:         locally (but any non-local elements will be sent to the
513:         appropriate processor during matrix assembly).
514:       - Here, we set all entries for a particular row at once.
515:   */
516:   for (j=info->ys; j<info->ys+info->ym; j++) {
517:     for (i=info->xs; i<info->xs+info->xm; i++) {
518:       PetscReal xx = i*hx,yy = j*hy;
519:       row.j = j; row.i = i;
520:       /* boundary points */
521:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
522:         v[0] = 1.0;
523:         MatSetValuesStencil(B,1,&row,1,&row,v,INSERT_VALUES);
524:       } else {
525:         /* interior grid points */
526:         const PetscScalar
527:           ux_E     = dhx*(x[j][i+1]-x[j][i]),
528:           uy_E     = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]),
529:           ux_W     = dhx*(x[j][i]-x[j][i-1]),
530:           uy_W     = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]),
531:           ux_N     = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]),
532:           uy_N     = dhy*(x[j+1][i]-x[j][i]),
533:           ux_S     = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]),
534:           uy_S     = dhy*(x[j][i]-x[j-1][i]),
535:           u        = x[j][i],
536:           e_E      = eta(user,xx,yy,ux_E,uy_E),
537:           e_W      = eta(user,xx,yy,ux_W,uy_W),
538:           e_N      = eta(user,xx,yy,ux_N,uy_N),
539:           e_S      = eta(user,xx,yy,ux_S,uy_S),
540:           de_E     = deta(user,xx,yy,ux_E,uy_E),
541:           de_W     = deta(user,xx,yy,ux_W,uy_W),
542:           de_N     = deta(user,xx,yy,ux_N,uy_N),
543:           de_S     = deta(user,xx,yy,ux_S,uy_S),
544:           skew_E   = de_E*ux_E*uy_E,
545:           skew_W   = de_W*ux_W*uy_W,
546:           skew_N   = de_N*ux_N*uy_N,
547:           skew_S   = de_S*ux_S*uy_S,
548:           cross_EW = 0.25*(skew_E - skew_W),
549:           cross_NS = 0.25*(skew_N - skew_S),
550:           newt_E   = e_E+de_E*PetscSqr(ux_E),
551:           newt_W   = e_W+de_W*PetscSqr(ux_W),
552:           newt_N   = e_N+de_N*PetscSqr(uy_N),
553:           newt_S   = e_S+de_S*PetscSqr(uy_S);
554:         /* interior grid points */
555:         switch (user->jtype) {
556:         case JAC_BRATU:
557:           /* Jacobian from p=2 */
558:           v[0] = -hxdhy;                                           col[0].j = j-1;   col[0].i = i;
559:           v[1] = -hydhx;                                           col[1].j = j;     col[1].i = i-1;
560:           v[2] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(u);       col[2].j = row.j; col[2].i = row.i;
561:           v[3] = -hydhx;                                           col[3].j = j;     col[3].i = i+1;
562:           v[4] = -hxdhy;                                           col[4].j = j+1;   col[4].i = i;
563:           MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);
564:           break;
565:         case JAC_PICARD:
566:           /* Jacobian arising from Picard linearization */
567:           v[0] = -hxdhy*e_S;                                           col[0].j = j-1;   col[0].i = i;
568:           v[1] = -hydhx*e_W;                                           col[1].j = j;     col[1].i = i-1;
569:           v[2] = (e_W+e_E)*hydhx + (e_S+e_N)*hxdhy;                    col[2].j = row.j; col[2].i = row.i;
570:           v[3] = -hydhx*e_E;                                           col[3].j = j;     col[3].i = i+1;
571:           v[4] = -hxdhy*e_N;                                           col[4].j = j+1;   col[4].i = i;
572:           MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);
573:           break;
574:         case JAC_STAR:
575:           /* Full Jacobian, but only a star stencil */
576:           col[0].j = j-1; col[0].i = i;
577:           col[1].j = j;   col[1].i = i-1;
578:           col[2].j = j;   col[2].i = i;
579:           col[3].j = j;   col[3].i = i+1;
580:           col[4].j = j+1; col[4].i = i;
581:           v[0]     = -hxdhy*newt_S + cross_EW;
582:           v[1]     = -hydhx*newt_W + cross_NS;
583:           v[2]     = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*PetscExpScalar(u);
584:           v[3]     = -hydhx*newt_E - cross_NS;
585:           v[4]     = -hxdhy*newt_N - cross_EW;
586:           MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);
587:           break;
588:         case JAC_NEWTON:
589:           /* The Jacobian is

593:           */
594:           col[0].j = j-1; col[0].i = i-1;
595:           col[1].j = j-1; col[1].i = i;
596:           col[2].j = j-1; col[2].i = i+1;
597:           col[3].j = j;   col[3].i = i-1;
598:           col[4].j = j;   col[4].i = i;
599:           col[5].j = j;   col[5].i = i+1;
600:           col[6].j = j+1; col[6].i = i-1;
601:           col[7].j = j+1; col[7].i = i;
602:           col[8].j = j+1; col[8].i = i+1;
603:           v[0]     = -0.25*(skew_S + skew_W);
604:           v[1]     = -hxdhy*newt_S + cross_EW;
605:           v[2]     =  0.25*(skew_S + skew_E);
606:           v[3]     = -hydhx*newt_W + cross_NS;
607:           v[4]     = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*PetscExpScalar(u);
608:           v[5]     = -hydhx*newt_E - cross_NS;
609:           v[6]     =  0.25*(skew_N + skew_W);
610:           v[7]     = -hxdhy*newt_N - cross_EW;
611:           v[8]     = -0.25*(skew_N + skew_E);
612:           MatSetValuesStencil(B,1,&row,9,col,v,INSERT_VALUES);
613:           break;
614:         default:
615:           SETERRQ(PetscObjectComm((PetscObject)info->da),PETSC_ERR_SUP,"Jacobian type %d not implemented",user->jtype);
616:         }
617:       }
618:     }
619:   }

621:   /*
622:      Assemble matrix, using the 2-step process:
623:        MatAssemblyBegin(), MatAssemblyEnd().
624:   */
625:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
626:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

628:   if (J != B) {
629:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
630:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
631:   }

633:   /*
634:      Tell the matrix we will never add a new nonzero location to the
635:      matrix. If we do, it will generate an error.
636:   */
637:   if (user->jtype == JAC_NEWTON) {
638:     PetscLogFlops(info->xm*info->ym*(131.0));
639:   }
640:   return 0;
641: }

643: /***********************************************************
644:  * PreCheck implementation
645:  ***********************************************************/
646: PetscErrorCode PreCheckSetFromOptions(PreCheck precheck)
647: {
649:   PetscBool      flg;

652:   PetscOptionsBegin(precheck->comm,NULL,"PreCheck Options","none");
653:   PetscOptionsReal("-precheck_angle","Angle in degrees between successive search directions necessary to activate step correction","",precheck->angle,&precheck->angle,NULL);
654:   flg  = PETSC_FALSE;
655:   PetscOptionsBool("-precheck_monitor","Monitor choices made by precheck routine","",flg,&flg,NULL);
656:   if (flg) {
657:     PetscViewerASCIIOpen(precheck->comm,"stdout",&precheck->monitor);
658:   }
659:   PetscOptionsEnd();
660:   return 0;
661: }

663: /*
664:   Compare the direction of the current and previous step, modify the current step accordingly
665: */
666: PetscErrorCode PreCheckFunction(SNESLineSearch linesearch,Vec X,Vec Y,PetscBool *changed, void *ctx)
667: {
668:   PreCheck       precheck;
669:   Vec            Ylast;
670:   PetscScalar    dot;
671:   PetscInt       iter;
673:   SNES           snes;

676:   SNESLineSearchGetSNES(linesearch, &snes);
677:   precheck = (PreCheck)ctx;
678:   if (!precheck->Ylast) VecDuplicate(Y,&precheck->Ylast);
679:   Ylast = precheck->Ylast;
680:   SNESGetIterationNumber(snes,&iter);
681:   if (iter < 1) {
682:     VecCopy(Y,Ylast);
683:     *changed = PETSC_FALSE;
684:     return 0;
685:   }

687:   VecDot(Y,Ylast,&dot);
688:   VecNorm(Y,NORM_2,&ynorm);
689:   VecNorm(Ylast,NORM_2,&ylastnorm);
690:   /* Compute the angle between the vectors Y and Ylast, clip to keep inside the domain of acos() */
691:   theta         = PetscAcosReal((PetscReal)PetscClipInterval(PetscAbsScalar(dot) / (ynorm * ylastnorm),-1.0,1.0));
692:   angle_radians = precheck->angle * PETSC_PI / 180.;
693:   if (PetscAbsReal(theta) < angle_radians || PetscAbsReal(theta - PETSC_PI) < angle_radians) {
694:     /* Modify the step Y */
695:     PetscReal alpha,ydiffnorm;
696:     VecAXPY(Ylast,-1.0,Y);
697:     VecNorm(Ylast,NORM_2,&ydiffnorm);
698:     alpha = ylastnorm / ydiffnorm;
699:     VecCopy(Y,Ylast);
700:     VecScale(Y,alpha);
701:     if (precheck->monitor) {
702:       PetscViewerASCIIPrintf(precheck->monitor,"Angle %E degrees less than threshold %g, corrected step by alpha=%g\n",(double)(theta*180./PETSC_PI),(double)precheck->angle,(double)alpha);
703:     }
704:   } else {
705:     VecCopy(Y,Ylast);
706:     *changed = PETSC_FALSE;
707:     if (precheck->monitor) {
708:       PetscViewerASCIIPrintf(precheck->monitor,"Angle %E degrees exceeds threshold %g, no correction applied\n",(double)(theta*180./PETSC_PI),(double)precheck->angle);
709:     }
710:   }
711:   return 0;
712: }

714: PetscErrorCode PreCheckDestroy(PreCheck *precheck)
715: {
717:   if (!*precheck) return 0;
718:   VecDestroy(&(*precheck)->Ylast);
719:   PetscViewerDestroy(&(*precheck)->monitor);
720:   PetscFree(*precheck);
721:   return 0;
722: }

724: PetscErrorCode PreCheckCreate(MPI_Comm comm,PreCheck *precheck)
725: {
727:   PetscNew(precheck);

729:   (*precheck)->comm  = comm;
730:   (*precheck)->angle = 10.;     /* only active if angle is less than 10 degrees */
731:   return 0;
732: }

734: /*
735:       Applies some sweeps on nonlinear Gauss-Seidel on each process

737:  */
738: PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx)
739: {
740:   PetscInt       i,j,k,xs,ys,xm,ym,its,tot_its,sweeps,l,m;
741:   PetscReal      hx,hy,hxdhy,hydhx,dhx,dhy,sc;
742:   PetscScalar    **x,**b,bij,F,F0=0,J,y,u,eu;
743:   PetscReal      atol,rtol,stol;
744:   DM             da;
745:   AppCtx         *user = (AppCtx*)ctx;
746:   Vec            localX,localB;
747:   DMDALocalInfo  info;

750:   SNESGetDM(snes,&da);
751:   DMDAGetLocalInfo(da,&info);

753:   hx     = 1.0/(PetscReal)(info.mx-1);
754:   hy     = 1.0/(PetscReal)(info.my-1);
755:   sc     = hx*hy*user->lambda;
756:   dhx    = 1/hx;
757:   dhy    = 1/hy;
758:   hxdhy  = hx/hy;
759:   hydhx  = hy/hx;

761:   tot_its = 0;
762:   SNESNGSGetSweeps(snes,&sweeps);
763:   SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);
764:   DMGetLocalVector(da,&localX);
765:   if (B) {
766:     DMGetLocalVector(da,&localB);
767:   }
768:   if (B) {
769:     DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);
770:     DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);
771:   }
773:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
774:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
775:   DMDAVecGetArray(da,localX,&x);
776:   for (l=0; l<sweeps; l++) {
777:     /*
778:      Get local grid boundaries (for 2-dimensional DMDA):
779:      xs, ys   - starting grid indices (no ghost points)
780:      xm, ym   - widths of local grid (no ghost points)
781:      */
782:     DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
783:     for (m=0; m<2; m++) {
784:       for (j=ys; j<ys+ym; j++) {
785:         for (i=xs+(m+j)%2; i<xs+xm; i+=2) {
786:           PetscReal xx = i*hx,yy = j*hy;
787:           if (B) bij = b[j][i];
788:           else   bij = 0.;

790:           if (i == 0 || j == 0 || i == info.mx-1 || j == info.my-1) {
791:             /* boundary conditions are all zero Dirichlet */
792:             x[j][i] = 0.0 + bij;
793:           } else {
794:             const PetscScalar
795:               u_E = x[j][i+1],
796:               u_W = x[j][i-1],
797:               u_N = x[j+1][i],
798:               u_S = x[j-1][i];
799:             const PetscScalar
800:               uy_E   = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]),
801:               uy_W   = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]),
802:               ux_N   = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]),
803:               ux_S   = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]);
804:             u = x[j][i];
805:             for (k=0; k<its; k++) {
806:               const PetscScalar
807:                 ux_E   = dhx*(u_E-u),
808:                 ux_W   = dhx*(u-u_W),
809:                 uy_N   = dhy*(u_N-u),
810:                 uy_S   = dhy*(u-u_S),
811:                 e_E    = eta(user,xx,yy,ux_E,uy_E),
812:                 e_W    = eta(user,xx,yy,ux_W,uy_W),
813:                 e_N    = eta(user,xx,yy,ux_N,uy_N),
814:                 e_S    = eta(user,xx,yy,ux_S,uy_S),
815:                 de_E   = deta(user,xx,yy,ux_E,uy_E),
816:                 de_W   = deta(user,xx,yy,ux_W,uy_W),
817:                 de_N   = deta(user,xx,yy,ux_N,uy_N),
818:                 de_S   = deta(user,xx,yy,ux_S,uy_S),
819:                 newt_E = e_E+de_E*PetscSqr(ux_E),
820:                 newt_W = e_W+de_W*PetscSqr(ux_W),
821:                 newt_N = e_N+de_N*PetscSqr(uy_N),
822:                 newt_S = e_S+de_S*PetscSqr(uy_S),
823:                 uxx    = -hy * (e_E*ux_E - e_W*ux_W),
824:                 uyy    = -hx * (e_N*uy_N - e_S*uy_S);

826:               if (sc) eu = PetscExpScalar(u);
827:               else    eu = 0;

829:               F = uxx + uyy - sc*eu - bij;
830:               if (k == 0) F0 = F;
831:               J  = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*eu;
832:               y  = F/J;
833:               u -= y;
834:               tot_its++;
835:               if (atol > PetscAbsReal(PetscRealPart(F)) ||
836:                   rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) ||
837:                   stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) {
838:                 break;
839:               }
840:             }
841:             x[j][i] = u;
842:           }
843:         }
844:       }
845:     }
846:     /*
847: x     Restore vector
848:      */
849:   }
850:   DMDAVecRestoreArray(da,localX,&x);
851:   DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
852:   DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
853:   PetscLogFlops(tot_its*(118.0));
854:   DMRestoreLocalVector(da,&localX);
855:   if (B) {
857:     DMRestoreLocalVector(da,&localB);
858:   }
859:   return 0;
860: }

862: /*TEST

864:    test:
865:       nsize: 2
866:       args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype NEWTON
867:       requires: !single

869:    test:
870:       suffix: 2
871:       nsize: 2
872:       args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -precheck 1
873:       requires: !single

875:    test:
876:       suffix: 3
877:       nsize: 2
878:       args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -picard -precheck 1 -p 3
879:       requires: !single

881:    test:
882:       suffix: 3_star
883:       nsize: 2
884:       args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -picard -precheck 1 -p 3 -alloc_star
885:       output_file: output/ex15_3.out
886:       requires: !single

888:    test:
889:       suffix: 4
890:       args: -snes_monitor_short -snes_type newtonls -npc_snes_type ngs -snes_npc_side left -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -ksp_monitor_short -pc_type none
891:       requires: !single

893:    test:
894:       suffix: lag_jac
895:       nsize: 4
896:       args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 6.0 -lambda 0 -jtype NEWTON -snes_type ngmres -npc_snes_type newtonls -npc_snes_lag_jacobian 5 -npc_pc_type asm -npc_ksp_converged_reason -npc_snes_lag_jacobian_persists
897:       requires: !single

899:    test:
900:       suffix: lag_pc
901:       nsize: 4
902:       args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 6.0 -lambda 0 -jtype NEWTON -snes_type ngmres -npc_snes_type newtonls -npc_snes_lag_preconditioner 5 -npc_pc_type asm -npc_ksp_converged_reason -npc_snes_lag_preconditioner_persists
903:       requires: !single

905:    test:
906:       suffix: nleqerr
907:       args: -snes_monitor_short -snes_type newtonls -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -snes_linesearch_monitor -pc_type lu -snes_linesearch_type nleqerr
908:       requires: !single

910:    test:
911:       suffix: mf
912:       args: -snes_monitor_short -pc_type lu -da_refine 4  -p 3 -ksp_rtol 1.e-12  -snes_mf_operator
913:       requires: !single

915:    test:
916:       suffix: mf_picard
917:       args: -snes_monitor_short -pc_type lu -da_refine 4  -p 3 -ksp_rtol 1.e-12  -snes_mf_operator -picard
918:       requires: !single
919:       output_file: output/ex15_mf.out

921:    test:
922:       suffix: fd_picard
923:       args: -snes_monitor_short -pc_type lu -da_refine 2  -p 3 -ksp_rtol 1.e-12  -snes_fd -picard
924:       requires: !single

926:    test:
927:       suffix: fd_color_picard
928:       args: -snes_monitor_short -pc_type lu -da_refine 4  -p 3 -ksp_rtol 1.e-12  -snes_fd_color -picard
929:       requires: !single
930:       output_file: output/ex15_mf.out

932: TEST*/