Actual source code: ex15.c

petsc-3.3-p7 2013-05-11
  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: \n";

12: /*F

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.
 31: F*/

33: /*
34:     Program usage:  mpiexec -n <procs> ./pbratu [-help] [all PETSc options]
35:      e.g.,
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: */
49: #include <petscdmda.h>
50: #include <petscsnes.h>

52: /*
53:    User-defined application context - contains data needed by the
54:    application-provided call-back routines, FormJacobianLocal() and
55:    FormFunctionLocal().
56: */
57: typedef struct {
58:   PassiveReal lambda;         /* Bratu parameter */
59:   PassiveReal p;              /* Exponent in p-Laplacian */
60:   PassiveReal epsilon;        /* Regularization */
61:   PassiveReal source;         /* Source term */
62:   PetscInt    jtype;          /* What type of Jacobian to assemble */
63:   PetscBool   picard;
64: } AppCtx;

66: /*
67:    User-defined routines
68: */
69: static PetscErrorCode FormInitialGuess(AppCtx*,DM,Vec);
70: static PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
71: static PetscErrorCode FormFunctionPicardLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
72: static PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,AppCtx*);

74: typedef struct _n_PreCheck *PreCheck;
75: struct _n_PreCheck {
76:   MPI_Comm comm;
77:   PetscReal angle;
78:   Vec Ylast;
79:   PetscViewer monitor;
80: };
81: PetscErrorCode PreCheckCreate(MPI_Comm,PreCheck*);
82: PetscErrorCode PreCheckDestroy(PreCheck*);
83: PetscErrorCode PreCheckFunction(SNESLineSearch,Vec,Vec,PetscBool*,void*);
84: PetscErrorCode PreCheckSetFromOptions(PreCheck);

88: int main(int argc,char **argv)
89: {
90:   SNES                   snes;                 /* nonlinear solver */
91:   Vec                    x,r;                  /* solution, residual vectors */
92:   Mat                    A,B;                  /* Jacobian and preconditioning matrices */
93:   AppCtx                 user;                 /* user-defined work context */
94:   PetscInt               its;                  /* iterations for convergence */
95:   SNESConvergedReason    reason;               /* Check convergence */
96:   PetscBool              alloc_star;           /* Only allocate for the STAR stencil  */
97:   PetscReal              bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
98:   DM                     da,dastar;            /* distributed array data structure */
99:   PreCheck               precheck = PETSC_NULL; /* precheck context for version in this file */
100:   PetscInt               use_precheck;   /* 0=none, 1=version in this file, 2=SNES-provided version */
101:   PetscReal              precheck_angle; /* When manually setting the SNES-provided precheck function */
102:   PetscErrorCode         ierr;
103:   SNESLineSearch        linesearch;

105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106:      Initialize program
107:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

111:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112:      Initialize problem parameters
113:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114:   user.lambda = 6.0; user.p = 2.0; user.epsilon = 1e-5; user.source = 0.1; user.jtype = 4; alloc_star = PETSC_FALSE;
115:   use_precheck = 0; precheck_angle = 10.;
116:   user.picard = PETSC_FALSE;
117:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"p-Bratu options",__FILE__);
118:   {
119:     PetscOptionsReal("-lambda","Bratu parameter","",user.lambda,&user.lambda,NULL);
120:     PetscOptionsReal("-p","Exponent p' in p-Laplacian","",user.p,&user.p,NULL);
121:     PetscOptionsReal("-epsilon","Strain-regularization in p-Laplacian","",user.epsilon,&user.epsilon,NULL);
122:     PetscOptionsReal("-source","Constant source term","",user.source,&user.source,NULL);
123:     PetscOptionsInt("-jtype","Jacobian type, 1=plain, 2=first term 3=star 4=full","",user.jtype,&user.jtype,NULL);
124:     PetscOptionsName("-picard","Solve with defect-correction Picard iteration","",&user.picard);
125:     if (user.picard) {user.jtype = 2; user.p = 3;}
126:     PetscOptionsBool("-alloc_star","Allocate for STAR stencil (5-point)","",alloc_star,&alloc_star,NULL);
127:     PetscOptionsInt("-precheck","Use a pre-check correction intended for use with Picard iteration 1=this version, 2=library","",use_precheck,&use_precheck,NULL);
128:     if (use_precheck == 2) {    /* Using library version, get the angle */
129:       PetscOptionsReal("-precheck_angle","Angle in degrees between successive search directions necessary to activate step correction","",precheck_angle,&precheck_angle,PETSC_NULL);
130:     }
131:   }
132:   PetscOptionsEnd();
133:   if (user.lambda > bratu_lambda_max || user.lambda < bratu_lambda_min) {
134:     PetscPrintf(PETSC_COMM_WORLD,"WARNING: lambda %g out of range for p=2\n",user.lambda);
135:   }

137:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138:      Create nonlinear solver context
139:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140:   SNESCreate(PETSC_COMM_WORLD,&snes);

142:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143:      Create distributed array (DMDA) to manage parallel grid and vectors
144:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE,
146:                       1,1,PETSC_NULL,PETSC_NULL,&da);
147:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,
148:                       1,1,PETSC_NULL,PETSC_NULL,&dastar);

151:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152:      Extract global vectors from DM; then duplicate for remaining
153:      vectors that are the same types
154:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
155:   DMCreateGlobalVector(da,&x);
156:   VecDuplicate(x,&r);

158:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159:      Create matrix data structure; set Jacobian evaluation routine

161:      Set Jacobian matrix data structure and default Jacobian evaluation
162:      routine. User can override with:
163:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
164:                 (unless user explicitly sets preconditioner)
165:      -snes_mf_operator : form preconditioning matrix as set by the user,
166:                          but use matrix-free approx for Jacobian-vector
167:                          products within Newton-Krylov method

169:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170:   /* B can be type of MATAIJ,MATBAIJ or MATSBAIJ */
171:   DMCreateMatrix(alloc_star ? dastar : da,MATAIJ,&B);
172:   A    = B;

174:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175:      Set local function evaluation routine
176:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177:   DMSetApplicationContext(da, &user);
178:   SNESSetDM(snes,da);
179:   if (user.picard) {
180:     DMDASetLocalFunction(da,(DMDALocalFunction1)FormFunctionPicardLocal);
181:     SNESSetPicard(snes,r,PETSC_NULL,A,B,PETSC_NULL,&user);
182:   } else {
183:     DMDASetLocalFunction(da,(DMDALocalFunction1)FormFunctionLocal);
184:   }
185:  DMDASetLocalJacobian(da,(DMDALocalFunction1)FormJacobianLocal);

187:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188:      Customize nonlinear solver; set runtime options
189:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
190:   SNESSetFromOptions(snes);
191:   SNESGetSNESLineSearch(snes, &linesearch);
192:   /* Set up the precheck context if requested */
193:   if (use_precheck == 1) {      /* Use the precheck routines in this file */
194:     PreCheckCreate(PETSC_COMM_WORLD,&precheck);
195:     PreCheckSetFromOptions(precheck);
196:     SNESLineSearchSetPreCheck(linesearch,PreCheckFunction,precheck);
197:   } else if (use_precheck == 2) { /* Use the version provided by the library */
198:     SNESLineSearchSetPreCheck(linesearch,SNESLineSearchPreCheckPicard,&precheck_angle);
199:   }

201:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202:      Evaluate initial guess
203:      Note: The user should initialize the vector, x, with the initial guess
204:      for the nonlinear solver prior to calling SNESSolve().  In particular,
205:      to employ an initial guess of zero, the user should explicitly set
206:      this vector to zero by calling VecSet().
207:   */

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

211:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
212:      Solve nonlinear system
213:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
214:   SNESSolve(snes,PETSC_NULL,x);
215:   SNESGetIterationNumber(snes,&its);
216:   SNESGetConvergedReason(snes,&reason);

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

220:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
221:      Free work space.  All PETSc objects should be destroyed when they
222:      are no longer needed.
223:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

225:   if (A != B) {
226:     MatDestroy(&A);
227:   }
228:   MatDestroy(&B);
229:   VecDestroy(&x);
230:   VecDestroy(&r);
231:   SNESDestroy(&snes);
232:   DMDestroy(&da);
233:   DMDestroy(&dastar);
234:   PreCheckDestroy(&precheck);
235:   PetscFinalize();
236:   return 0;
237: }

239: /* ------------------------------------------------------------------- */
242: /*
243:    FormInitialGuess - Forms initial approximation.

245:    Input Parameters:
246:    user - user-defined application context
247:    X - vector

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

260:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
261:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

263:   hx     = 1.0/(PetscReal)(Mx-1);
264:   hy     = 1.0/(PetscReal)(My-1);
265:   temp1  = user->lambda / (user->lambda + 1.);

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

276:   /*
277:      Get local grid boundaries (for 2-dimensional DA):
278:        xs, ys   - starting grid indices (no ghost points)
279:        xm, ym   - widths of local grid (no ghost points)

281:   */
282:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

284:   /*
285:      Compute initial guess over the locally owned part of the grid
286:   */
287:   for (j=ys; j<ys+ym; j++) {
288:     temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
289:     for (i=xs; i<xs+xm; i++) {
290:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
291:         /* boundary conditions are all zero Dirichlet */
292:         x[j][i] = 0.0;
293:       } else {
294:         if (user->lambda != 0) {
295:           x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
296:         } else {
297:           /* The solution above is an exact solution for lambda=0, this avoids "accidentally" starting
298:            * with an exact solution. */
299:           const PetscReal
300:             xx = 2*(PetscReal)i/(Mx-1) - 1,
301:             yy = 2*(PetscReal)j/(My-1) - 1;
302:           x[j][i] = (1 - xx*xx) * (1-yy*yy) * xx * yy;
303:         }
304:       }
305:     }
306:   }

308:   /*
309:      Restore vector
310:   */
311:   DMDAVecRestoreArray(da,X,&x);

313:   return(0);
314: }

316: /* p-Laplacian diffusivity */
317: PETSC_STATIC_INLINE PetscScalar eta(const AppCtx *ctx,PetscScalar ux,PetscScalar uy)
318: {return pow(PetscSqr(ctx->epsilon)+0.5*(ux*ux + uy*uy),0.5*(ctx->p-2.));}
319: PETSC_STATIC_INLINE PetscScalar deta(const AppCtx *ctx,PetscScalar ux,PetscScalar uy)
320: {
321:   return (ctx->p == 2)
322:     ? 0
323:     : pow(PetscSqr(ctx->epsilon)+0.5*(ux*ux + uy*uy),0.5*(ctx->p-4)) * 0.5 * (ctx->p-2.);
324: }

327: /* ------------------------------------------------------------------- */
330: /*
331:    FormFunctionLocal - Evaluates nonlinear function, F(x).
332:  */
333: static PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
334: {
335:   PetscReal      hx,hy,dhx,dhy,sc,source;
336:   PetscInt       i,j;

339:   hx     = 1.0/(PetscReal)(info->mx-1);
340:   hy     = 1.0/(PetscReal)(info->my-1);
341:   sc     = hx*hy*user->lambda;
342:   source = hx*hy*user->source;
343:   dhx    = 1/hx;
344:   dhy    = 1/hy;
345:   /*
346:      Compute function over the locally owned part of the grid
347:   */
348:   for (j=info->ys; j<info->ys+info->ym; j++) {
349:     for (i=info->xs; i<info->xs+info->xm; i++) {
350:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
351:         f[j][i] = x[j][i];
352:       } else {
353:         const PetscScalar
354:           u = x[j][i],
355:           ux_E = dhx*(x[j][i+1]-x[j][i]),
356:           uy_E = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]),
357:           ux_W = dhx*(x[j][i]-x[j][i-1]),
358:           uy_W = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]),
359:           ux_N = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]),
360:           uy_N = dhy*(x[j+1][i]-x[j][i]),
361:           ux_S = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]),
362:           uy_S = dhy*(x[j][i]-x[j-1][i]),
363:           e_E = eta(user,ux_E,uy_E),
364:           e_W = eta(user,ux_W,uy_W),
365:           e_N = eta(user,ux_N,uy_N),
366:           e_S = eta(user,ux_S,uy_S),
367:           uxx = -hy * (e_E*ux_E - e_W*ux_W),
368:           uyy = -hx * (e_N*uy_N - e_S*uy_S);
369:         /** For p=2, these terms decay to:
370:         * uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx
371:         * uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy
372:         **/
373:         f[j][i] = uxx + uyy - sc*PetscExpScalar(u) - source;
374:       }
375:     }
376:   }
377:   return(0);
378: }

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

386: */
387: static PetscErrorCode FormFunctionPicardLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
388: {
389:   PetscReal      hx,hy,sc,source;
390:   PetscInt       i,j;

393:   hx     = 1.0/(PetscReal)(info->mx-1);
394:   hy     = 1.0/(PetscReal)(info->my-1);
395:   sc     = hx*hy*user->lambda;
396:   source = hx*hy*user->source;
397:   /*
398:      Compute function over the locally owned part of the grid
399:   */
400:   for (j=info->ys; j<info->ys+info->ym; j++) {
401:     for (i=info->xs; i<info->xs+info->xm; i++) {
402:       if (!(i == 0 || j == 0 || i == info->mx-1 || j == info->my-1)) {
403:         const PetscScalar u = x[j][i];
404:         f[j][i] = sc*PetscExpScalar(u) + source;
405:       }
406:     }
407:   }
408:   return(0);
409: }

413: /*
414:    FormJacobianLocal - Evaluates Jacobian matrix.
415: */
416: static PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat B,AppCtx *user)
417: {
419:   PetscInt       i,j;
420:   MatStencil     col[9],row;
421:   PetscScalar    v[9];
422:   PetscReal      hx,hy,hxdhy,hydhx,dhx,dhy,sc;

425:   hx     = 1.0/(PetscReal)(info->mx-1);
426:   hy     = 1.0/(PetscReal)(info->my-1);
427:   sc     = hx*hy*user->lambda;
428:   dhx    = 1/hx;
429:   dhy    = 1/hy;
430:   hxdhy  = hx/hy;
431:   hydhx  = hy/hx;

433:   /*
434:      Compute entries for the locally owned part of the Jacobian.
435:       - PETSc parallel matrix formats are partitioned by
436:         contiguous chunks of rows across the processors.
437:       - Each processor needs to insert only elements that it owns
438:         locally (but any non-local elements will be sent to the
439:         appropriate processor during matrix assembly).
440:       - Here, we set all entries for a particular row at once.
441:   */
442:   for (j=info->ys; j<info->ys+info->ym; j++) {
443:     for (i=info->xs; i<info->xs+info->xm; i++) {
444:       row.j = j; row.i = i;
445:       /* boundary points */
446:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
447:         v[0] = 1.0;
448:         MatSetValuesStencil(B,1,&row,1,&row,v,INSERT_VALUES);
449:       } else {
450:       /* interior grid points */
451:         const PetscScalar
452:           ux_E = dhx*(x[j][i+1]-x[j][i]),
453:           uy_E = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]),
454:           ux_W = dhx*(x[j][i]-x[j][i-1]),
455:           uy_W = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]),
456:           ux_N = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]),
457:           uy_N = dhy*(x[j+1][i]-x[j][i]),
458:           ux_S = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]),
459:           uy_S = dhy*(x[j][i]-x[j-1][i]),
460:           u = x[j][i],
461:           e_E = eta(user,ux_E,uy_E),
462:           e_W = eta(user,ux_W,uy_W),
463:           e_N = eta(user,ux_N,uy_N),
464:           e_S = eta(user,ux_S,uy_S),
465:           de_E = deta(user,ux_E,uy_E),
466:           de_W = deta(user,ux_W,uy_W),
467:           de_N = deta(user,ux_N,uy_N),
468:           de_S = deta(user,ux_S,uy_S),
469:           skew_E = de_E*ux_E*uy_E,
470:           skew_W = de_W*ux_W*uy_W,
471:           skew_N = de_N*ux_N*uy_N,
472:           skew_S = de_S*ux_S*uy_S,
473:           cross_EW = 0.25*(skew_E - skew_W),
474:           cross_NS = 0.25*(skew_N - skew_S),
475:           newt_E = e_E+de_E*PetscSqr(ux_E),
476:           newt_W = e_W+de_W*PetscSqr(ux_W),
477:           newt_N = e_N+de_N*PetscSqr(uy_N),
478:           newt_S = e_S+de_S*PetscSqr(uy_S);
479:       /* interior grid points */
480:         switch (user->jtype) {
481:           case 1:
482:             /* Jacobian from p=2 */
483:             v[0] = -hxdhy;                                           col[0].j = j-1;   col[0].i = i;
484:             v[1] = -hydhx;                                           col[1].j = j;     col[1].i = i-1;
485:             v[2] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(u);       col[2].j = row.j; col[2].i = row.i;
486:             v[3] = -hydhx;                                           col[3].j = j;     col[3].i = i+1;
487:             v[4] = -hxdhy;                                           col[4].j = j+1;   col[4].i = i;
488:             MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);
489:             break;
490:           case 2:
491:             /* Jacobian arising from Picard linearization */
492:             v[0] = -hxdhy*e_S;                                           col[0].j = j-1;   col[0].i = i;
493:             v[1] = -hydhx*e_W;                                           col[1].j = j;     col[1].i = i-1;
494:             v[2] = (e_W+e_E)*hydhx + (e_S+e_N)*hxdhy;                    col[2].j = row.j; col[2].i = row.i;
495:             v[3] = -hydhx*e_E;                                           col[3].j = j;     col[3].i = i+1;
496:             v[4] = -hxdhy*e_N;                                           col[4].j = j+1;   col[4].i = i;
497:             MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);
498:             break;
499:           case 3:
500:             /* Full Jacobian, but only a star stencil */
501:             col[0].j = j-1; col[0].i = i;
502:             col[1].j = j;   col[1].i = i-1;
503:             col[2].j = j;   col[2].i = i;
504:             col[3].j = j;   col[3].i = i+1;
505:             col[4].j = j+1; col[4].i = i;
506:             v[0] = -hxdhy*newt_S + cross_EW;
507:             v[1] = -hydhx*newt_W + cross_NS;
508:             v[2] = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*PetscExpScalar(u);
509:             v[3] = -hydhx*newt_E - cross_NS;
510:             v[4] = -hxdhy*newt_N - cross_EW;
511:             MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);
512:             break;
513:           case 4:
514:             /** The Jacobian is
515:             *
517:             *
518:             **/
519:             col[0].j = j-1; col[0].i = i-1;
520:             col[1].j = j-1; col[1].i = i;
521:             col[2].j = j-1; col[2].i = i+1;
522:             col[3].j = j;   col[3].i = i-1;
523:             col[4].j = j;   col[4].i = i;
524:             col[5].j = j;   col[5].i = i+1;
525:             col[6].j = j+1; col[6].i = i-1;
526:             col[7].j = j+1; col[7].i = i;
527:             col[8].j = j+1; col[8].i = i+1;
528:             v[0] = -0.25*(skew_S + skew_W);
529:             v[1] = -hxdhy*newt_S + cross_EW;
530:             v[2] =  0.25*(skew_S + skew_E);
531:             v[3] = -hydhx*newt_W + cross_NS;
532:             v[4] = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*PetscExpScalar(u);
533:             v[5] = -hydhx*newt_E - cross_NS;
534:             v[6] =  0.25*(skew_N + skew_W);
535:             v[7] = -hxdhy*newt_N - cross_EW;
536:             v[8] = -0.25*(skew_N + skew_E);
537:             MatSetValuesStencil(B,1,&row,9,col,v,INSERT_VALUES);
538:             break;
539:           default:
540:             SETERRQ1(((PetscObject)info->da)->comm,PETSC_ERR_SUP,"Jacobian type %d not implemented",user->jtype);
541:         }
542:       }
543:     }
544:   }

546:   /*
547:      Assemble matrix, using the 2-step process:
548:        MatAssemblyBegin(), MatAssemblyEnd().
549:   */
550:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
551:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

553:   /*
554:      Tell the matrix we will never add a new nonzero location to the
555:      matrix. If we do, it will generate an error.
556:   */
557:   MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
558:   return(0);
559: }

561: /***********************************************************
562:  * PreCheck implementation
563:  ***********************************************************/
566: PetscErrorCode PreCheckSetFromOptions(PreCheck precheck)
567: {
569:   PetscBool flg;

572:   PetscOptionsBegin(precheck->comm,PETSC_NULL,"PreCheck Options","none");
573:   PetscOptionsReal("-precheck_angle","Angle in degrees between successive search directions necessary to activate step correction","",precheck->angle,&precheck->angle,PETSC_NULL);
574:   flg = PETSC_FALSE;
575:   PetscOptionsBool("-precheck_monitor","Monitor choices made by precheck routine","",flg,&flg,PETSC_NULL);
576:   if (flg) {
577:     PetscViewerASCIIOpen(precheck->comm,"stdout",&precheck->monitor);
578:   }
579:   PetscOptionsEnd();
580:   return(0);
581: }

585: /*
586:   Compare the direction of the current and previous step, modify the current step accordingly
587: */
588: PetscErrorCode PreCheckFunction(SNESLineSearch linesearch,Vec X,Vec Y,PetscBool *changed, void *ctx)
589: {
591:   PreCheck       precheck;
592:   Vec            Ylast;
593:   PetscScalar    dot;
594:   PetscInt       iter;
596:   SNES           snes;

599:   SNESLineSearchGetSNES(linesearch, &snes);
600:   precheck = (PreCheck)ctx;
601:   if (!precheck->Ylast) {VecDuplicate(Y,&precheck->Ylast);}
602:   Ylast = precheck->Ylast;
603:   SNESGetIterationNumber(snes,&iter);
604:   if (iter < 1) {
605:     VecCopy(Y,Ylast);
606:     *changed = PETSC_FALSE;
607:     return(0);
608:   }

610:   VecDot(Y,Ylast,&dot);
611:   VecNorm(Y,NORM_2,&ynorm);
612:   VecNorm(Ylast,NORM_2,&ylastnorm);
613:   /* Compute the angle between the vectors Y and Ylast, clip to keep inside the domain of acos() */
614:   theta = acos((double)PetscClipInterval(PetscAbsScalar(dot) / (ynorm * ylastnorm),-1.0,1.0));
615:   angle_radians = precheck->angle * PETSC_PI / 180.;
616:   if (PetscAbsReal(theta) < angle_radians || PetscAbsReal(theta - PETSC_PI) < angle_radians) {
617:     /* Modify the step Y */
618:     PetscReal alpha,ydiffnorm;
619:     VecAXPY(Ylast,-1.0,Y);
620:     VecNorm(Ylast,NORM_2,&ydiffnorm);
621:     alpha = ylastnorm / ydiffnorm;
622:     VecCopy(Y,Ylast);
623:     VecScale(Y,alpha);
624:     if (precheck->monitor) {
625:       PetscViewerASCIIPrintf(precheck->monitor,"Angle %E degrees less than threshold %G, corrected step by alpha=%G\n",theta*180./PETSC_PI,precheck->angle,alpha);
626:     }
627:   } else {
628:     VecCopy(Y,Ylast);
629:     *changed = PETSC_FALSE;
630:     if (precheck->monitor) {
631:       PetscViewerASCIIPrintf(precheck->monitor,"Angle %E degrees exceeds threshold %G, no correction applied\n",theta*180./PETSC_PI,precheck->angle);
632:     }
633:   }
634:   return(0);
635: }

639: PetscErrorCode PreCheckDestroy(PreCheck *precheck)
640: {

644:   if (!*precheck) return(0);
645:   VecDestroy(&(*precheck)->Ylast);
646:   PetscViewerDestroy(&(*precheck)->monitor);
647:   PetscFree(*precheck);
648:   return(0);
649: }

653: PetscErrorCode PreCheckCreate(MPI_Comm comm,PreCheck *precheck)
654: {

658:   PetscMalloc(sizeof(struct _n_PreCheck),precheck);
659:   PetscMemzero(*precheck,sizeof(struct _n_PreCheck));
660:   (*precheck)->comm = comm;
661:   (*precheck)->angle = 10.;     /* only active if angle is less than 10 degrees */
662:   return(0);
663: }