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
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 PetscReal300: 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 PetscScalar354: 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 PetscScalar452: 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: *
516: * -div [ eta (grad u) + deta (grad u0 . grad u) grad u0 ] - (eE u0) u
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;
595: PetscReal ynorm,ylastnorm,theta,angle_radians;
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: }