Actual source code: ex15.c

petsc-master 2018-04-19
Report Typos and Errors
  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.
 35: /*
 36:       mpiexec -n 2 ./ex15 -snes_monitor -ksp_monitor log_summary
 37: */

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

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

 54: /* These functions _should_ be internal, but currently have a reverse dependency so cannot be set with
 55:  * DMDASNESSetPicardLocal.  This hack needs to be fixed in PETSc. */
 56: PETSC_EXTERN PetscErrorCode SNESPicardComputeFunction(SNES,Vec,Vec,void*);
 57: PETSC_EXTERN PetscErrorCode SNESPicardComputeJacobian(SNES,Vec,Mat,Mat,void*);

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

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

 79: /*
 80:    User-defined routines
 81: */
 82: static PetscErrorCode FormRHS(AppCtx*,DM,Vec);
 83: static PetscErrorCode FormInitialGuess(AppCtx*,DM,Vec);
 84: static PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
 85: static PetscErrorCode FormFunctionPicardLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
 86: static PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,Mat,AppCtx*);
 87: static PetscErrorCode NonlinearGS(SNES,Vec,Vec,void*);

 89: typedef struct _n_PreCheck *PreCheck;
 90: struct _n_PreCheck {
 91:   MPI_Comm    comm;
 92:   PetscReal   angle;
 93:   Vec         Ylast;
 94:   PetscViewer monitor;
 95: };
 96: PetscErrorCode PreCheckCreate(MPI_Comm,PreCheck*);
 97: PetscErrorCode PreCheckDestroy(PreCheck*);
 98: PetscErrorCode PreCheckFunction(SNESLineSearch,Vec,Vec,PetscBool*,void*);
 99: PetscErrorCode PreCheckSetFromOptions(PreCheck);

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

120:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121:      Initialize program
122:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

164:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165:      Create distributed array (DMDA) to manage parallel grid and vectors
166:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,4,4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
168:   DMSetFromOptions(da);
169:   DMSetUp(da);
170:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&dastar);
171:   DMSetFromOptions(dastar);
172:   DMSetUp(dastar);

174:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175:      Extract global vectors from DM; then duplicate for remaining
176:      vectors that are the same types
177:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
178:   DMCreateGlobalVector(da,&x);
179:   VecDuplicate(x,&r);
180:   VecDuplicate(x,&b);

182:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183:      Create matrix data structure; set Jacobian evaluation routine

185:      Set Jacobian matrix data structure and default Jacobian evaluation
186:      routine. User can override with:
187:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
188:                 (unless user explicitly sets preconditioner)
189:      -snes_mf_operator : form preconditioning matrix as set by the user,
190:                          but use matrix-free approx for Jacobian-vector
191:                          products within Newton-Krylov method

193:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
194:   /* B can be type of MATAIJ,MATBAIJ or MATSBAIJ */
195:   DMCreateMatrix(alloc_star ? dastar : da,&B);
196:   A    = B;

198:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199:      Set local function evaluation routine
200:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
201:   DMSetApplicationContext(da, &user);
202:   SNESSetDM(snes,da);
203:   if (user.picard) {
204:     /*
205:         This is not really right requiring the user to call SNESSetFunction/Jacobian but the DMDASNESSetPicardLocal() cannot access
206:         the SNES to set it
207:     */
208:     DMDASNESSetPicardLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionPicardLocal,
209:                                   (PetscErrorCode (*)(DMDALocalInfo*,void*,Mat,Mat,void*))FormJacobianLocal,&user);
210:     SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user);
211:     SNESSetJacobian(snes,NULL,NULL,SNESPicardComputeJacobian,&user);
212:   } else {
213:     DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);
214:     DMDASNESSetJacobianLocal(da,(PetscErrorCode (*)(DMDALocalInfo*,void*,Mat,Mat,void*))FormJacobianLocal,&user);
215:   }


218:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
219:      Customize nonlinear solver; set runtime options
220:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
221:   SNESSetFromOptions(snes);
222:   SNESSetNGS(snes,NonlinearGS,&user);
223:   SNESGetLineSearch(snes, &linesearch);
224:   /* Set up the precheck context if requested */
225:   if (use_precheck == 1) {      /* Use the precheck routines in this file */
226:     PreCheckCreate(PETSC_COMM_WORLD,&precheck);
227:     PreCheckSetFromOptions(precheck);
228:     SNESLineSearchSetPreCheck(linesearch,PreCheckFunction,precheck);
229:   } else if (use_precheck == 2) { /* Use the version provided by the library */
230:     SNESLineSearchSetPreCheck(linesearch,SNESLineSearchPreCheckPicard,&precheck_angle);
231:   }

233:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
234:      Evaluate initial guess
235:      Note: The user should initialize the vector, x, with the initial guess
236:      for the nonlinear solver prior to calling SNESSolve().  In particular,
237:      to employ an initial guess of zero, the user should explicitly set
238:      this vector to zero by calling VecSet().
239:   */

241:   FormInitialGuess(&user,da,x);
242:   FormRHS(&user,da,b);

244:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
245:      Solve nonlinear system
246:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
247:   SNESSolve(snes,b,x);
248:   SNESGetIterationNumber(snes,&its);
249:   SNESGetConvergedReason(snes,&reason);

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

253:   if (write_output) {
254:     PetscViewer viewer;
255:     PetscViewerVTKOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer);
256:     VecView(x,viewer);
257:     PetscViewerDestroy(&viewer);
258:   }

260:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
261:      Free work space.  All PETSc objects should be destroyed when they
262:      are no longer needed.
263:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

265:   if (A != B) {
266:     MatDestroy(&A);
267:   }
268:   MatDestroy(&B);
269:   VecDestroy(&x);
270:   VecDestroy(&r);
271:   VecDestroy(&b);
272:   SNESDestroy(&snes);
273:   DMDestroy(&da);
274:   DMDestroy(&dastar);
275:   PreCheckDestroy(&precheck);
276:   PetscFinalize();
277:   return ierr;
278: }

280: /* ------------------------------------------------------------------- */
281: /*
282:    FormInitialGuess - Forms initial approximation.

284:    Input Parameters:
285:    user - user-defined application context
286:    X - vector

288:    Output Parameter:
289:    X - vector
290:  */
291: static PetscErrorCode FormInitialGuess(AppCtx *user,DM da,Vec X)
292: {
293:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
295:   PetscReal      temp1,temp,hx,hy;
296:   PetscScalar    **x;

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

301:   hx    = 1.0/(PetscReal)(Mx-1);
302:   hy    = 1.0/(PetscReal)(My-1);
303:   temp1 = user->lambda / (user->lambda + 1.);

305:   /*
306:      Get a pointer to vector data.
307:        - For default PETSc vectors, VecGetArray() returns a pointer to
308:          the data array.  Otherwise, the routine is implementation dependent.
309:        - You MUST call VecRestoreArray() when you no longer need access to
310:          the array.
311:   */
312:   DMDAVecGetArray(da,X,&x);

314:   /*
315:      Get local grid boundaries (for 2-dimensional DA):
316:        xs, ys   - starting grid indices (no ghost points)
317:        xm, ym   - widths of local grid (no ghost points)

319:   */
320:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

322:   /*
323:      Compute initial guess over the locally owned part of the grid
324:   */
325:   for (j=ys; j<ys+ym; j++) {
326:     temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
327:     for (i=xs; i<xs+xm; i++) {
328:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
329:         /* boundary conditions are all zero Dirichlet */
330:         x[j][i] = 0.0;
331:       } else {
332:         if (user->initial == -1) {
333:           if (user->lambda != 0) {
334:             x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
335:           } else {
336:             /* The solution above is an exact solution for lambda=0, this avoids "accidentally" starting
337:              * with an exact solution. */
338:             const PetscReal
339:               xx = 2*(PetscReal)i/(Mx-1) - 1,
340:               yy = 2*(PetscReal)j/(My-1) - 1;
341:             x[j][i] = (1 - xx*xx) * (1-yy*yy) * xx * yy;
342:           }
343:         } else if (user->initial == 0) {
344:           x[j][i] = 0.;
345:         } else if (user->initial == 1) {
346:           const PetscReal
347:             xx = 2*(PetscReal)i/(Mx-1) - 1,
348:             yy = 2*(PetscReal)j/(My-1) - 1;
349:           x[j][i] = (1 - xx*xx) * (1-yy*yy) * xx * yy;
350:         } else {
351:           if (user->lambda != 0) {
352:             x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
353:           } else {
354:             x[j][i] = 0.5*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
355:           }
356:         }
357:       }
358:     }
359:   }
360:   /*
361:      Restore vector
362:   */
363:   DMDAVecRestoreArray(da,X,&x);
364:   return(0);
365: }

367: /*
368:    FormRHS - Forms constant RHS for the problem.

370:    Input Parameters:
371:    user - user-defined application context
372:    B - RHS vector

374:    Output Parameter:
375:    B - vector
376:  */
377: static PetscErrorCode FormRHS(AppCtx *user,DM da,Vec B)
378: {
379:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
381:   PetscReal      hx,hy;
382:   PetscScalar    **b;

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

387:   hx    = 1.0/(PetscReal)(Mx-1);
388:   hy    = 1.0/(PetscReal)(My-1);
389:   DMDAVecGetArray(da,B,&b);
390:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
391:   for (j=ys; j<ys+ym; j++) {
392:     for (i=xs; i<xs+xm; i++) {
393:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
394:         b[j][i] = 0.0;
395:       } else {
396:         b[j][i] = hx*hy*user->source;
397:       }
398:     }
399:   }
400:   DMDAVecRestoreArray(da,B,&b);
401:   return(0);
402: }

404: PETSC_STATIC_INLINE PetscReal kappa(const AppCtx *ctx,PetscReal x,PetscReal y)
405: {
406:   return (((PetscInt)(x*ctx->blocks[0])) + ((PetscInt)(y*ctx->blocks[1]))) % 2 ? ctx->kappa : 1.0;
407: }
408: /* p-Laplacian diffusivity */
409: PETSC_STATIC_INLINE PetscScalar eta(const AppCtx *ctx,PetscReal x,PetscReal y,PetscScalar ux,PetscScalar uy)
410: {
411:   return kappa(ctx,x,y) * PetscPowScalar(PetscSqr(ctx->epsilon)+0.5*(ux*ux + uy*uy),0.5*(ctx->p-2.));
412: }
413: PETSC_STATIC_INLINE PetscScalar deta(const AppCtx *ctx,PetscReal x,PetscReal y,PetscScalar ux,PetscScalar uy)
414: {
415:   return (ctx->p == 2)
416:          ? 0
417:          : kappa(ctx,x,y)*PetscPowScalar(PetscSqr(ctx->epsilon)+0.5*(ux*ux + uy*uy),0.5*(ctx->p-4)) * 0.5 * (ctx->p-2.);
418: }


421: /* ------------------------------------------------------------------- */
422: /*
423:    FormFunctionLocal - Evaluates nonlinear function, F(x).
424:  */
425: static PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
426: {
427:   PetscReal      hx,hy,dhx,dhy,sc;
428:   PetscInt       i,j;
429:   PetscScalar    eu;


434:   hx     = 1.0/(PetscReal)(info->mx-1);
435:   hy     = 1.0/(PetscReal)(info->my-1);
436:   sc     = hx*hy*user->lambda;
437:   dhx    = 1/hx;
438:   dhy    = 1/hy;
439:   /*
440:      Compute function over the locally owned part of the grid
441:   */
442:   for (j=info->ys; j<info->ys+info->ym; j++) {
443:     for (i=info->xs; i<info->xs+info->xm; i++) {
444:       PetscReal xx = i*hx,yy = j*hy;
445:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
446:         f[j][i] = x[j][i];
447:       } else {
448:         const PetscScalar
449:           u    = x[j][i],
450:           ux_E = dhx*(x[j][i+1]-x[j][i]),
451:           uy_E = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]),
452:           ux_W = dhx*(x[j][i]-x[j][i-1]),
453:           uy_W = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]),
454:           ux_N = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]),
455:           uy_N = dhy*(x[j+1][i]-x[j][i]),
456:           ux_S = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]),
457:           uy_S = dhy*(x[j][i]-x[j-1][i]),
458:           e_E  = eta(user,xx,yy,ux_E,uy_E),
459:           e_W  = eta(user,xx,yy,ux_W,uy_W),
460:           e_N  = eta(user,xx,yy,ux_N,uy_N),
461:           e_S  = eta(user,xx,yy,ux_S,uy_S),
462:           uxx  = -hy * (e_E*ux_E - e_W*ux_W),
463:           uyy  = -hx * (e_N*uy_N - e_S*uy_S);
464:         if (sc) eu = PetscExpScalar(u);
465:         else    eu = 0.;
466:         /** For p=2, these terms decay to:
467:         * uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx
468:         * uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy
469:         **/
470:         f[j][i] = uxx + uyy - sc*eu;
471:       }
472:     }
473:   }
474:   PetscLogFlops(info->xm*info->ym*(72.0));
475:   return(0);
476: }

478: /*
479:     This is the opposite sign of the part of FormFunctionLocal that excludes the A(x) x part of the operation,
480:     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)

482: */
483: static PetscErrorCode FormFunctionPicardLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
484: {
485:   PetscReal hx,hy,sc;
486:   PetscInt  i,j;

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

508: /*
509:    FormJacobianLocal - Evaluates Jacobian matrix.
510: */
511: static PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat J,Mat B,AppCtx *user)
512: {
514:   PetscInt       i,j;
515:   MatStencil     col[9],row;
516:   PetscScalar    v[9];
517:   PetscReal      hx,hy,hxdhy,hydhx,dhx,dhy,sc;

520:   hx    = 1.0/(PetscReal)(info->mx-1);
521:   hy    = 1.0/(PetscReal)(info->my-1);
522:   sc    = hx*hy*user->lambda;
523:   dhx   = 1/hx;
524:   dhy   = 1/hy;
525:   hxdhy = hx/hy;
526:   hydhx = hy/hx;

528:   /*
529:      Compute entries for the locally owned part of the Jacobian.
530:       - PETSc parallel matrix formats are partitioned by
531:         contiguous chunks of rows across the processors.
532:       - Each processor needs to insert only elements that it owns
533:         locally (but any non-local elements will be sent to the
534:         appropriate processor during matrix assembly).
535:       - Here, we set all entries for a particular row at once.
536:   */
537:   for (j=info->ys; j<info->ys+info->ym; j++) {
538:     for (i=info->xs; i<info->xs+info->xm; i++) {
539:       PetscReal xx = i*hx,yy = j*hy;
540:       row.j = j; row.i = i;
541:       /* boundary points */
542:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
543:         v[0] = 1.0;
544:         MatSetValuesStencil(B,1,&row,1,&row,v,INSERT_VALUES);
545:       } else {
546:         /* interior grid points */
547:         const PetscScalar
548:           ux_E     = dhx*(x[j][i+1]-x[j][i]),
549:           uy_E     = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]),
550:           ux_W     = dhx*(x[j][i]-x[j][i-1]),
551:           uy_W     = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]),
552:           ux_N     = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]),
553:           uy_N     = dhy*(x[j+1][i]-x[j][i]),
554:           ux_S     = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]),
555:           uy_S     = dhy*(x[j][i]-x[j-1][i]),
556:           u        = x[j][i],
557:           e_E      = eta(user,xx,yy,ux_E,uy_E),
558:           e_W      = eta(user,xx,yy,ux_W,uy_W),
559:           e_N      = eta(user,xx,yy,ux_N,uy_N),
560:           e_S      = eta(user,xx,yy,ux_S,uy_S),
561:           de_E     = deta(user,xx,yy,ux_E,uy_E),
562:           de_W     = deta(user,xx,yy,ux_W,uy_W),
563:           de_N     = deta(user,xx,yy,ux_N,uy_N),
564:           de_S     = deta(user,xx,yy,ux_S,uy_S),
565:           skew_E   = de_E*ux_E*uy_E,
566:           skew_W   = de_W*ux_W*uy_W,
567:           skew_N   = de_N*ux_N*uy_N,
568:           skew_S   = de_S*ux_S*uy_S,
569:           cross_EW = 0.25*(skew_E - skew_W),
570:           cross_NS = 0.25*(skew_N - skew_S),
571:           newt_E   = e_E+de_E*PetscSqr(ux_E),
572:           newt_W   = e_W+de_W*PetscSqr(ux_W),
573:           newt_N   = e_N+de_N*PetscSqr(uy_N),
574:           newt_S   = e_S+de_S*PetscSqr(uy_S);
575:         /* interior grid points */
576:         switch (user->jtype) {
577:         case JAC_BRATU:
578:           /* Jacobian from p=2 */
579:           v[0] = -hxdhy;                                           col[0].j = j-1;   col[0].i = i;
580:           v[1] = -hydhx;                                           col[1].j = j;     col[1].i = i-1;
581:           v[2] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(u);       col[2].j = row.j; col[2].i = row.i;
582:           v[3] = -hydhx;                                           col[3].j = j;     col[3].i = i+1;
583:           v[4] = -hxdhy;                                           col[4].j = j+1;   col[4].i = i;
584:           MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);
585:           break;
586:         case JAC_PICARD:
587:           /* Jacobian arising from Picard linearization */
588:           v[0] = -hxdhy*e_S;                                           col[0].j = j-1;   col[0].i = i;
589:           v[1] = -hydhx*e_W;                                           col[1].j = j;     col[1].i = i-1;
590:           v[2] = (e_W+e_E)*hydhx + (e_S+e_N)*hxdhy;                    col[2].j = row.j; col[2].i = row.i;
591:           v[3] = -hydhx*e_E;                                           col[3].j = j;     col[3].i = i+1;
592:           v[4] = -hxdhy*e_N;                                           col[4].j = j+1;   col[4].i = i;
593:           MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);
594:           break;
595:         case JAC_STAR:
596:           /* Full Jacobian, but only a star stencil */
597:           col[0].j = j-1; col[0].i = i;
598:           col[1].j = j;   col[1].i = i-1;
599:           col[2].j = j;   col[2].i = i;
600:           col[3].j = j;   col[3].i = i+1;
601:           col[4].j = j+1; col[4].i = i;
602:           v[0]     = -hxdhy*newt_S + cross_EW;
603:           v[1]     = -hydhx*newt_W + cross_NS;
604:           v[2]     = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*PetscExpScalar(u);
605:           v[3]     = -hydhx*newt_E - cross_NS;
606:           v[4]     = -hxdhy*newt_N - cross_EW;
607:           MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);
608:           break;
609:         case JAC_NEWTON:
610:           /** The Jacobian is
611:           *
612:           * -div [ eta (grad u) + deta (grad u0 . grad u) grad u0 ] - (eE u0) u
613:           *
614:           **/
615:           col[0].j = j-1; col[0].i = i-1;
616:           col[1].j = j-1; col[1].i = i;
617:           col[2].j = j-1; col[2].i = i+1;
618:           col[3].j = j;   col[3].i = i-1;
619:           col[4].j = j;   col[4].i = i;
620:           col[5].j = j;   col[5].i = i+1;
621:           col[6].j = j+1; col[6].i = i-1;
622:           col[7].j = j+1; col[7].i = i;
623:           col[8].j = j+1; col[8].i = i+1;
624:           v[0]     = -0.25*(skew_S + skew_W);
625:           v[1]     = -hxdhy*newt_S + cross_EW;
626:           v[2]     =  0.25*(skew_S + skew_E);
627:           v[3]     = -hydhx*newt_W + cross_NS;
628:           v[4]     = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*PetscExpScalar(u);
629:           v[5]     = -hydhx*newt_E - cross_NS;
630:           v[6]     =  0.25*(skew_N + skew_W);
631:           v[7]     = -hxdhy*newt_N - cross_EW;
632:           v[8]     = -0.25*(skew_N + skew_E);
633:           MatSetValuesStencil(B,1,&row,9,col,v,INSERT_VALUES);
634:           break;
635:         default:
636:           SETERRQ1(PetscObjectComm((PetscObject)info->da),PETSC_ERR_SUP,"Jacobian type %d not implemented",user->jtype);
637:         }
638:       }
639:     }
640:   }

642:   /*
643:      Assemble matrix, using the 2-step process:
644:        MatAssemblyBegin(), MatAssemblyEnd().
645:   */
646:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
647:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

649:   if (J != B) {
650:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
651:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
652:   }

654:   /*
655:      Tell the matrix we will never add a new nonzero location to the
656:      matrix. If we do, it will generate an error.
657:   */
658:   if (user->jtype == JAC_NEWTON) {
659:     PetscLogFlops(info->xm*info->ym*(131.0));
660:   }
661:   MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
662:   return(0);
663: }

665: /***********************************************************
666:  * PreCheck implementation
667:  ***********************************************************/
668: PetscErrorCode PreCheckSetFromOptions(PreCheck precheck)
669: {
671:   PetscBool      flg;

674:   PetscOptionsBegin(precheck->comm,NULL,"PreCheck Options","none");
675:   PetscOptionsReal("-precheck_angle","Angle in degrees between successive search directions necessary to activate step correction","",precheck->angle,&precheck->angle,NULL);
676:   flg  = PETSC_FALSE;
677:   PetscOptionsBool("-precheck_monitor","Monitor choices made by precheck routine","",flg,&flg,NULL);
678:   if (flg) {
679:     PetscViewerASCIIOpen(precheck->comm,"stdout",&precheck->monitor);
680:   }
681:   PetscOptionsEnd();
682:   return(0);
683: }

685: /*
686:   Compare the direction of the current and previous step, modify the current step accordingly
687: */
688: PetscErrorCode PreCheckFunction(SNESLineSearch linesearch,Vec X,Vec Y,PetscBool *changed, void *ctx)
689: {
691:   PreCheck       precheck;
692:   Vec            Ylast;
693:   PetscScalar    dot;
694:   PetscInt       iter;
695:   PetscReal      ynorm,ylastnorm,theta,angle_radians;
696:   SNES           snes;

699:   SNESLineSearchGetSNES(linesearch, &snes);
700:   precheck = (PreCheck)ctx;
701:   if (!precheck->Ylast) {VecDuplicate(Y,&precheck->Ylast);}
702:   Ylast = precheck->Ylast;
703:   SNESGetIterationNumber(snes,&iter);
704:   if (iter < 1) {
705:     VecCopy(Y,Ylast);
706:     *changed = PETSC_FALSE;
707:     return(0);
708:   }

710:   VecDot(Y,Ylast,&dot);
711:   VecNorm(Y,NORM_2,&ynorm);
712:   VecNorm(Ylast,NORM_2,&ylastnorm);
713:   /* Compute the angle between the vectors Y and Ylast, clip to keep inside the domain of acos() */
714:   theta         = PetscAcosReal((PetscReal)PetscClipInterval(PetscAbsScalar(dot) / (ynorm * ylastnorm),-1.0,1.0));
715:   angle_radians = precheck->angle * PETSC_PI / 180.;
716:   if (PetscAbsReal(theta) < angle_radians || PetscAbsReal(theta - PETSC_PI) < angle_radians) {
717:     /* Modify the step Y */
718:     PetscReal alpha,ydiffnorm;
719:     VecAXPY(Ylast,-1.0,Y);
720:     VecNorm(Ylast,NORM_2,&ydiffnorm);
721:     alpha = ylastnorm / ydiffnorm;
722:     VecCopy(Y,Ylast);
723:     VecScale(Y,alpha);
724:     if (precheck->monitor) {
725:       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);
726:     }
727:   } else {
728:     VecCopy(Y,Ylast);
729:     *changed = PETSC_FALSE;
730:     if (precheck->monitor) {
731:       PetscViewerASCIIPrintf(precheck->monitor,"Angle %E degrees exceeds threshold %g, no correction applied\n",(double)(theta*180./PETSC_PI),(double)precheck->angle);
732:     }
733:   }
734:   return(0);
735: }

737: PetscErrorCode PreCheckDestroy(PreCheck *precheck)
738: {

742:   if (!*precheck) return(0);
743:   VecDestroy(&(*precheck)->Ylast);
744:   PetscViewerDestroy(&(*precheck)->monitor);
745:   PetscFree(*precheck);
746:   return(0);
747: }

749: PetscErrorCode PreCheckCreate(MPI_Comm comm,PreCheck *precheck)
750: {

754:   PetscNew(precheck);

756:   (*precheck)->comm  = comm;
757:   (*precheck)->angle = 10.;     /* only active if angle is less than 10 degrees */
758:   return(0);
759: }

761: /*
762:       Applies some sweeps on nonlinear Gauss-Seidel on each process

764:  */
765: PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx)
766: {
767:   PetscInt       i,j,k,xs,ys,xm,ym,its,tot_its,sweeps,l,m;
769:   PetscReal      hx,hy,hxdhy,hydhx,dhx,dhy,sc;
770:   PetscScalar    **x,**b,bij,F,F0=0,J,y,u,eu;
771:   PetscReal      atol,rtol,stol;
772:   DM             da;
773:   AppCtx         *user = (AppCtx*)ctx;
774:   Vec            localX,localB;
775:   DMDALocalInfo  info;

778:   SNESGetDM(snes,&da);
779:   DMDAGetLocalInfo(da,&info);

781:   hx     = 1.0/(PetscReal)(info.mx-1);
782:   hy     = 1.0/(PetscReal)(info.my-1);
783:   sc     = hx*hy*user->lambda;
784:   dhx    = 1/hx;
785:   dhy    = 1/hy;
786:   hxdhy  = hx/hy;
787:   hydhx  = hy/hx;

789:   tot_its = 0;
790:   SNESNGSGetSweeps(snes,&sweeps);
791:   SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);
792:   DMGetLocalVector(da,&localX);
793:   if (B) {
794:     DMGetLocalVector(da,&localB);
795:   }
796:   if (B) {
797:     DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);
798:     DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);
799:   }
800:   if (B) DMDAVecGetArrayRead(da,localB,&b);
801:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
802:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
803:   DMDAVecGetArray(da,localX,&x);
804:   for (l=0; l<sweeps; l++) {
805:     /*
806:      Get local grid boundaries (for 2-dimensional DMDA):
807:      xs, ys   - starting grid indices (no ghost points)
808:      xm, ym   - widths of local grid (no ghost points)
809:      */
810:     DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
811:     for (m=0; m<2; m++) {
812:       for (j=ys; j<ys+ym; j++) {
813:         for (i=xs+(m+j)%2; i<xs+xm; i+=2) {
814:           PetscReal xx = i*hx,yy = j*hy;
815:           if (B) bij = b[j][i];
816:           else   bij = 0.;

818:           if (i == 0 || j == 0 || i == info.mx-1 || j == info.my-1) {
819:             /* boundary conditions are all zero Dirichlet */
820:             x[j][i] = 0.0 + bij;
821:           } else {
822:             const PetscScalar
823:               u_E = x[j][i+1],
824:               u_W = x[j][i-1],
825:               u_N = x[j+1][i],
826:               u_S = x[j-1][i];
827:             const PetscScalar
828:               uy_E   = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]),
829:               uy_W   = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]),
830:               ux_N   = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]),
831:               ux_S   = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]);
832:             u = x[j][i];
833:             for (k=0; k<its; k++) {
834:               const PetscScalar
835:                 ux_E   = dhx*(u_E-u),
836:                 ux_W   = dhx*(u-u_W),
837:                 uy_N   = dhy*(u_N-u),
838:                 uy_S   = dhy*(u-u_S),
839:                 e_E    = eta(user,xx,yy,ux_E,uy_E),
840:                 e_W    = eta(user,xx,yy,ux_W,uy_W),
841:                 e_N    = eta(user,xx,yy,ux_N,uy_N),
842:                 e_S    = eta(user,xx,yy,ux_S,uy_S),
843:                 de_E   = deta(user,xx,yy,ux_E,uy_E),
844:                 de_W   = deta(user,xx,yy,ux_W,uy_W),
845:                 de_N   = deta(user,xx,yy,ux_N,uy_N),
846:                 de_S   = deta(user,xx,yy,ux_S,uy_S),
847:                 newt_E = e_E+de_E*PetscSqr(ux_E),
848:                 newt_W = e_W+de_W*PetscSqr(ux_W),
849:                 newt_N = e_N+de_N*PetscSqr(uy_N),
850:                 newt_S = e_S+de_S*PetscSqr(uy_S),
851:                 uxx    = -hy * (e_E*ux_E - e_W*ux_W),
852:                 uyy    = -hx * (e_N*uy_N - e_S*uy_S);

854:               if (sc) eu = PetscExpScalar(u);
855:               else    eu = 0;

857:               F = uxx + uyy - sc*eu - bij;
858:               if (k == 0) F0 = F;
859:               J  = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*eu;
860:               y  = F/J;
861:               u -= y;
862:               tot_its++;
863:               if (atol > PetscAbsReal(PetscRealPart(F)) ||
864:                   rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) ||
865:                   stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) {
866:                 break;
867:               }
868:             }
869:             x[j][i] = u;
870:           }
871:         }
872:       }
873:     }
874:     /*
875: x     Restore vector
876:      */
877:   }
878:   DMDAVecRestoreArray(da,localX,&x);
879:   DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
880:   DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
881:   PetscLogFlops(tot_its*(118.0));
882:   DMRestoreLocalVector(da,&localX);
883:   if (B) {
884:     DMDAVecRestoreArrayRead(da,localB,&b);
885:     DMRestoreLocalVector(da,&localB);
886:   }
887:   return(0);
888: }


891: /*TEST

893:    test:
894:       nsize: 2
895:       args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype NEWTON
896:       requires: !single

898:    test:
899:       suffix: 2
900:       nsize: 2
901:       args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -precheck 1
902:       requires: !single

904:    test:
905:       suffix: 3
906:       nsize: 2
907:       args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -picard -precheck 1
908:       requires: !single

910:    test:
911:       suffix: 4
912:       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
913:       requires: !single

915:    test:
916:       suffix: lag_jac
917:       nsize: 4
918:       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
919:       requires: !single

921:    test:
922:       suffix: lag_pc
923:       nsize: 4
924:       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
925:       requires: !single

927:    test:
928:       suffix: nleqerr
929:       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
930:       requires: !single

932: TEST*/