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
591: -div [ eta (grad u) + deta (grad u0 . grad u) grad u0 ] - (eE u0) u
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;
672: PetscReal ynorm,ylastnorm,theta,angle_radians;
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: }
772: if (B) DMDAVecGetArrayRead(da,localB,&b);
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) {
856: DMDAVecRestoreArrayRead(da,localB,&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*/