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);if (ierr) return ierr;
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;
138: if (user.p != 3) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Picard iteration is only supported for p == 3");
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 ierr;
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;
278: PetscReal temp1,temp,hx,hy;
279: PetscScalar **x;
282: 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);
284: hx = 1.0/(PetscReal)(Mx-1);
285: hy = 1.0/(PetscReal)(My-1);
286: temp1 = user->lambda / (user->lambda + 1.);
288: /*
289: Get a pointer to vector data.
290: - For default PETSc vectors, VecGetArray() returns a pointer to
291: the data array. Otherwise, the routine is implementation dependent.
292: - You MUST call VecRestoreArray() when you no longer need access to
293: the array.
294: */
295: DMDAVecGetArray(da,X,&x);
297: /*
298: Get local grid boundaries (for 2-dimensional DA):
299: xs, ys - starting grid indices (no ghost points)
300: xm, ym - widths of local grid (no ghost points)
302: */
303: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
305: /*
306: Compute initial guess over the locally owned part of the grid
307: */
308: for (j=ys; j<ys+ym; j++) {
309: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
310: for (i=xs; i<xs+xm; i++) {
311: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
312: /* boundary conditions are all zero Dirichlet */
313: x[j][i] = 0.0;
314: } else {
315: if (user->initial == -1) {
316: if (user->lambda != 0) {
317: x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
318: } else {
319: /* The solution above is an exact solution for lambda=0, this avoids "accidentally" starting
320: * with an exact solution. */
321: const PetscReal
322: xx = 2*(PetscReal)i/(Mx-1) - 1,
323: yy = 2*(PetscReal)j/(My-1) - 1;
324: x[j][i] = (1 - xx*xx) * (1-yy*yy) * xx * yy;
325: }
326: } else if (user->initial == 0) {
327: x[j][i] = 0.;
328: } else if (user->initial == 1) {
329: const PetscReal
330: xx = 2*(PetscReal)i/(Mx-1) - 1,
331: yy = 2*(PetscReal)j/(My-1) - 1;
332: x[j][i] = (1 - xx*xx) * (1-yy*yy) * xx * yy;
333: } else {
334: if (user->lambda != 0) {
335: x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
336: } else {
337: x[j][i] = 0.5*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
338: }
339: }
340: }
341: }
342: }
343: /*
344: Restore vector
345: */
346: DMDAVecRestoreArray(da,X,&x);
347: return(0);
348: }
350: /*
351: FormRHS - Forms constant RHS for the problem.
353: Input Parameters:
354: user - user-defined application context
355: B - RHS vector
357: Output Parameter:
358: B - vector
359: */
360: static PetscErrorCode FormRHS(AppCtx *user,DM da,Vec B)
361: {
362: PetscInt i,j,Mx,My,xs,ys,xm,ym;
364: PetscReal hx,hy;
365: PetscScalar **b;
368: 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);
370: hx = 1.0/(PetscReal)(Mx-1);
371: hy = 1.0/(PetscReal)(My-1);
372: DMDAVecGetArray(da,B,&b);
373: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
374: for (j=ys; j<ys+ym; j++) {
375: for (i=xs; i<xs+xm; i++) {
376: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
377: b[j][i] = 0.0;
378: } else {
379: b[j][i] = hx*hy*user->source;
380: }
381: }
382: }
383: DMDAVecRestoreArray(da,B,&b);
384: return(0);
385: }
387: PETSC_STATIC_INLINE PetscReal kappa(const AppCtx *ctx,PetscReal x,PetscReal y)
388: {
389: return (((PetscInt)(x*ctx->blocks[0])) + ((PetscInt)(y*ctx->blocks[1]))) % 2 ? ctx->kappa : 1.0;
390: }
391: /* p-Laplacian diffusivity */
392: PETSC_STATIC_INLINE PetscScalar eta(const AppCtx *ctx,PetscReal x,PetscReal y,PetscScalar ux,PetscScalar uy)
393: {
394: return kappa(ctx,x,y) * PetscPowScalar(PetscSqr(ctx->epsilon)+0.5*(ux*ux + uy*uy),0.5*(ctx->p-2.));
395: }
396: PETSC_STATIC_INLINE PetscScalar deta(const AppCtx *ctx,PetscReal x,PetscReal y,PetscScalar ux,PetscScalar uy)
397: {
398: return (ctx->p == 2)
399: ? 0
400: : kappa(ctx,x,y)*PetscPowScalar(PetscSqr(ctx->epsilon)+0.5*(ux*ux + uy*uy),0.5*(ctx->p-4)) * 0.5 * (ctx->p-2.);
401: }
403: /* ------------------------------------------------------------------- */
404: /*
405: FormFunctionLocal - Evaluates nonlinear function, F(x).
406: */
407: static PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
408: {
409: PetscReal hx,hy,dhx,dhy,sc;
410: PetscInt i,j;
411: PetscScalar eu;
415: hx = 1.0/(PetscReal)(info->mx-1);
416: hy = 1.0/(PetscReal)(info->my-1);
417: sc = hx*hy*user->lambda;
418: dhx = 1/hx;
419: dhy = 1/hy;
420: /*
421: Compute function over the locally owned part of the grid
422: */
423: for (j=info->ys; j<info->ys+info->ym; j++) {
424: for (i=info->xs; i<info->xs+info->xm; i++) {
425: PetscReal xx = i*hx,yy = j*hy;
426: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
427: f[j][i] = x[j][i];
428: } else {
429: const PetscScalar
430: u = x[j][i],
431: ux_E = dhx*(x[j][i+1]-x[j][i]),
432: uy_E = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]),
433: ux_W = dhx*(x[j][i]-x[j][i-1]),
434: uy_W = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]),
435: ux_N = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]),
436: uy_N = dhy*(x[j+1][i]-x[j][i]),
437: ux_S = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]),
438: uy_S = dhy*(x[j][i]-x[j-1][i]),
439: e_E = eta(user,xx,yy,ux_E,uy_E),
440: e_W = eta(user,xx,yy,ux_W,uy_W),
441: e_N = eta(user,xx,yy,ux_N,uy_N),
442: e_S = eta(user,xx,yy,ux_S,uy_S),
443: uxx = -hy * (e_E*ux_E - e_W*ux_W),
444: uyy = -hx * (e_N*uy_N - e_S*uy_S);
445: if (sc) eu = PetscExpScalar(u);
446: else eu = 0.;
447: /* For p=2, these terms decay to:
448: uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx
449: uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy
450: */
451: f[j][i] = uxx + uyy - sc*eu;
452: }
453: }
454: }
455: PetscLogFlops(info->xm*info->ym*(72.0));
456: return(0);
457: }
459: /*
460: This is the opposite sign of the part of FormFunctionLocal that excludes the A(x) x part of the operation,
461: 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)
463: */
464: static PetscErrorCode FormFunctionPicardLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
465: {
466: PetscReal hx,hy,sc;
467: PetscInt i,j;
471: hx = 1.0/(PetscReal)(info->mx-1);
472: hy = 1.0/(PetscReal)(info->my-1);
473: sc = hx*hy*user->lambda;
474: /*
475: Compute function over the locally owned part of the grid
476: */
477: for (j=info->ys; j<info->ys+info->ym; j++) {
478: for (i=info->xs; i<info->xs+info->xm; i++) {
479: if (!(i == 0 || j == 0 || i == info->mx-1 || j == info->my-1)) {
480: const PetscScalar u = x[j][i];
481: f[j][i] = sc*PetscExpScalar(u);
482: } else {
483: f[j][i] = 0.0; /* this is zero because the A(x) x term forces the x to be zero on the boundary */
484: }
485: }
486: }
487: PetscLogFlops(info->xm*info->ym*2.0);
488: return(0);
489: }
491: /*
492: FormJacobianLocal - Evaluates Jacobian matrix.
493: */
494: static PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat J,Mat B,AppCtx *user)
495: {
497: PetscInt i,j;
498: MatStencil col[9],row;
499: PetscScalar v[9];
500: PetscReal hx,hy,hxdhy,hydhx,dhx,dhy,sc;
503: MatZeroEntries(B);
504: hx = 1.0/(PetscReal)(info->mx-1);
505: hy = 1.0/(PetscReal)(info->my-1);
506: sc = hx*hy*user->lambda;
507: dhx = 1/hx;
508: dhy = 1/hy;
509: hxdhy = hx/hy;
510: hydhx = hy/hx;
512: /*
513: Compute entries for the locally owned part of the Jacobian.
514: - PETSc parallel matrix formats are partitioned by
515: contiguous chunks of rows across the processors.
516: - Each processor needs to insert only elements that it owns
517: locally (but any non-local elements will be sent to the
518: appropriate processor during matrix assembly).
519: - Here, we set all entries for a particular row at once.
520: */
521: for (j=info->ys; j<info->ys+info->ym; j++) {
522: for (i=info->xs; i<info->xs+info->xm; i++) {
523: PetscReal xx = i*hx,yy = j*hy;
524: row.j = j; row.i = i;
525: /* boundary points */
526: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
527: v[0] = 1.0;
528: MatSetValuesStencil(B,1,&row,1,&row,v,INSERT_VALUES);
529: } else {
530: /* interior grid points */
531: const PetscScalar
532: ux_E = dhx*(x[j][i+1]-x[j][i]),
533: uy_E = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]),
534: ux_W = dhx*(x[j][i]-x[j][i-1]),
535: uy_W = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]),
536: ux_N = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]),
537: uy_N = dhy*(x[j+1][i]-x[j][i]),
538: ux_S = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]),
539: uy_S = dhy*(x[j][i]-x[j-1][i]),
540: u = x[j][i],
541: e_E = eta(user,xx,yy,ux_E,uy_E),
542: e_W = eta(user,xx,yy,ux_W,uy_W),
543: e_N = eta(user,xx,yy,ux_N,uy_N),
544: e_S = eta(user,xx,yy,ux_S,uy_S),
545: de_E = deta(user,xx,yy,ux_E,uy_E),
546: de_W = deta(user,xx,yy,ux_W,uy_W),
547: de_N = deta(user,xx,yy,ux_N,uy_N),
548: de_S = deta(user,xx,yy,ux_S,uy_S),
549: skew_E = de_E*ux_E*uy_E,
550: skew_W = de_W*ux_W*uy_W,
551: skew_N = de_N*ux_N*uy_N,
552: skew_S = de_S*ux_S*uy_S,
553: cross_EW = 0.25*(skew_E - skew_W),
554: cross_NS = 0.25*(skew_N - skew_S),
555: newt_E = e_E+de_E*PetscSqr(ux_E),
556: newt_W = e_W+de_W*PetscSqr(ux_W),
557: newt_N = e_N+de_N*PetscSqr(uy_N),
558: newt_S = e_S+de_S*PetscSqr(uy_S);
559: /* interior grid points */
560: switch (user->jtype) {
561: case JAC_BRATU:
562: /* Jacobian from p=2 */
563: v[0] = -hxdhy; col[0].j = j-1; col[0].i = i;
564: v[1] = -hydhx; col[1].j = j; col[1].i = i-1;
565: v[2] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(u); col[2].j = row.j; col[2].i = row.i;
566: v[3] = -hydhx; col[3].j = j; col[3].i = i+1;
567: v[4] = -hxdhy; col[4].j = j+1; col[4].i = i;
568: MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);
569: break;
570: case JAC_PICARD:
571: /* Jacobian arising from Picard linearization */
572: v[0] = -hxdhy*e_S; col[0].j = j-1; col[0].i = i;
573: v[1] = -hydhx*e_W; col[1].j = j; col[1].i = i-1;
574: v[2] = (e_W+e_E)*hydhx + (e_S+e_N)*hxdhy; col[2].j = row.j; col[2].i = row.i;
575: v[3] = -hydhx*e_E; col[3].j = j; col[3].i = i+1;
576: v[4] = -hxdhy*e_N; col[4].j = j+1; col[4].i = i;
577: MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);
578: break;
579: case JAC_STAR:
580: /* Full Jacobian, but only a star stencil */
581: col[0].j = j-1; col[0].i = i;
582: col[1].j = j; col[1].i = i-1;
583: col[2].j = j; col[2].i = i;
584: col[3].j = j; col[3].i = i+1;
585: col[4].j = j+1; col[4].i = i;
586: v[0] = -hxdhy*newt_S + cross_EW;
587: v[1] = -hydhx*newt_W + cross_NS;
588: v[2] = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*PetscExpScalar(u);
589: v[3] = -hydhx*newt_E - cross_NS;
590: v[4] = -hxdhy*newt_N - cross_EW;
591: MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);
592: break;
593: case JAC_NEWTON:
594: /* The Jacobian is
596: -div [ eta (grad u) + deta (grad u0 . grad u) grad u0 ] - (eE u0) u
598: */
599: col[0].j = j-1; col[0].i = i-1;
600: col[1].j = j-1; col[1].i = i;
601: col[2].j = j-1; col[2].i = i+1;
602: col[3].j = j; col[3].i = i-1;
603: col[4].j = j; col[4].i = i;
604: col[5].j = j; col[5].i = i+1;
605: col[6].j = j+1; col[6].i = i-1;
606: col[7].j = j+1; col[7].i = i;
607: col[8].j = j+1; col[8].i = i+1;
608: v[0] = -0.25*(skew_S + skew_W);
609: v[1] = -hxdhy*newt_S + cross_EW;
610: v[2] = 0.25*(skew_S + skew_E);
611: v[3] = -hydhx*newt_W + cross_NS;
612: v[4] = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*PetscExpScalar(u);
613: v[5] = -hydhx*newt_E - cross_NS;
614: v[6] = 0.25*(skew_N + skew_W);
615: v[7] = -hxdhy*newt_N - cross_EW;
616: v[8] = -0.25*(skew_N + skew_E);
617: MatSetValuesStencil(B,1,&row,9,col,v,INSERT_VALUES);
618: break;
619: default:
620: SETERRQ1(PetscObjectComm((PetscObject)info->da),PETSC_ERR_SUP,"Jacobian type %d not implemented",user->jtype);
621: }
622: }
623: }
624: }
626: /*
627: Assemble matrix, using the 2-step process:
628: MatAssemblyBegin(), MatAssemblyEnd().
629: */
630: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
631: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
633: if (J != B) {
634: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
635: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
636: }
638: /*
639: Tell the matrix we will never add a new nonzero location to the
640: matrix. If we do, it will generate an error.
641: */
642: if (user->jtype == JAC_NEWTON) {
643: PetscLogFlops(info->xm*info->ym*(131.0));
644: }
645: return(0);
646: }
648: /***********************************************************
649: * PreCheck implementation
650: ***********************************************************/
651: PetscErrorCode PreCheckSetFromOptions(PreCheck precheck)
652: {
654: PetscBool flg;
657: PetscOptionsBegin(precheck->comm,NULL,"PreCheck Options","none");
658: PetscOptionsReal("-precheck_angle","Angle in degrees between successive search directions necessary to activate step correction","",precheck->angle,&precheck->angle,NULL);
659: flg = PETSC_FALSE;
660: PetscOptionsBool("-precheck_monitor","Monitor choices made by precheck routine","",flg,&flg,NULL);
661: if (flg) {
662: PetscViewerASCIIOpen(precheck->comm,"stdout",&precheck->monitor);
663: }
664: PetscOptionsEnd();
665: return(0);
666: }
668: /*
669: Compare the direction of the current and previous step, modify the current step accordingly
670: */
671: PetscErrorCode PreCheckFunction(SNESLineSearch linesearch,Vec X,Vec Y,PetscBool *changed, void *ctx)
672: {
674: PreCheck precheck;
675: Vec Ylast;
676: PetscScalar dot;
677: PetscInt iter;
678: PetscReal ynorm,ylastnorm,theta,angle_radians;
679: SNES snes;
682: SNESLineSearchGetSNES(linesearch, &snes);
683: precheck = (PreCheck)ctx;
684: if (!precheck->Ylast) {VecDuplicate(Y,&precheck->Ylast);}
685: Ylast = precheck->Ylast;
686: SNESGetIterationNumber(snes,&iter);
687: if (iter < 1) {
688: VecCopy(Y,Ylast);
689: *changed = PETSC_FALSE;
690: return(0);
691: }
693: VecDot(Y,Ylast,&dot);
694: VecNorm(Y,NORM_2,&ynorm);
695: VecNorm(Ylast,NORM_2,&ylastnorm);
696: /* Compute the angle between the vectors Y and Ylast, clip to keep inside the domain of acos() */
697: theta = PetscAcosReal((PetscReal)PetscClipInterval(PetscAbsScalar(dot) / (ynorm * ylastnorm),-1.0,1.0));
698: angle_radians = precheck->angle * PETSC_PI / 180.;
699: if (PetscAbsReal(theta) < angle_radians || PetscAbsReal(theta - PETSC_PI) < angle_radians) {
700: /* Modify the step Y */
701: PetscReal alpha,ydiffnorm;
702: VecAXPY(Ylast,-1.0,Y);
703: VecNorm(Ylast,NORM_2,&ydiffnorm);
704: alpha = ylastnorm / ydiffnorm;
705: VecCopy(Y,Ylast);
706: VecScale(Y,alpha);
707: if (precheck->monitor) {
708: 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);
709: }
710: } else {
711: VecCopy(Y,Ylast);
712: *changed = PETSC_FALSE;
713: if (precheck->monitor) {
714: PetscViewerASCIIPrintf(precheck->monitor,"Angle %E degrees exceeds threshold %g, no correction applied\n",(double)(theta*180./PETSC_PI),(double)precheck->angle);
715: }
716: }
717: return(0);
718: }
720: PetscErrorCode PreCheckDestroy(PreCheck *precheck)
721: {
725: if (!*precheck) return(0);
726: VecDestroy(&(*precheck)->Ylast);
727: PetscViewerDestroy(&(*precheck)->monitor);
728: PetscFree(*precheck);
729: return(0);
730: }
732: PetscErrorCode PreCheckCreate(MPI_Comm comm,PreCheck *precheck)
733: {
737: PetscNew(precheck);
739: (*precheck)->comm = comm;
740: (*precheck)->angle = 10.; /* only active if angle is less than 10 degrees */
741: return(0);
742: }
744: /*
745: Applies some sweeps on nonlinear Gauss-Seidel on each process
747: */
748: PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx)
749: {
750: PetscInt i,j,k,xs,ys,xm,ym,its,tot_its,sweeps,l,m;
752: PetscReal hx,hy,hxdhy,hydhx,dhx,dhy,sc;
753: PetscScalar **x,**b,bij,F,F0=0,J,y,u,eu;
754: PetscReal atol,rtol,stol;
755: DM da;
756: AppCtx *user = (AppCtx*)ctx;
757: Vec localX,localB;
758: DMDALocalInfo info;
761: SNESGetDM(snes,&da);
762: DMDAGetLocalInfo(da,&info);
764: hx = 1.0/(PetscReal)(info.mx-1);
765: hy = 1.0/(PetscReal)(info.my-1);
766: sc = hx*hy*user->lambda;
767: dhx = 1/hx;
768: dhy = 1/hy;
769: hxdhy = hx/hy;
770: hydhx = hy/hx;
772: tot_its = 0;
773: SNESNGSGetSweeps(snes,&sweeps);
774: SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);
775: DMGetLocalVector(da,&localX);
776: if (B) {
777: DMGetLocalVector(da,&localB);
778: }
779: if (B) {
780: DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);
781: DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);
782: }
783: if (B) {DMDAVecGetArrayRead(da,localB,&b);}
784: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
785: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
786: DMDAVecGetArray(da,localX,&x);
787: for (l=0; l<sweeps; l++) {
788: /*
789: Get local grid boundaries (for 2-dimensional DMDA):
790: xs, ys - starting grid indices (no ghost points)
791: xm, ym - widths of local grid (no ghost points)
792: */
793: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
794: for (m=0; m<2; m++) {
795: for (j=ys; j<ys+ym; j++) {
796: for (i=xs+(m+j)%2; i<xs+xm; i+=2) {
797: PetscReal xx = i*hx,yy = j*hy;
798: if (B) bij = b[j][i];
799: else bij = 0.;
801: if (i == 0 || j == 0 || i == info.mx-1 || j == info.my-1) {
802: /* boundary conditions are all zero Dirichlet */
803: x[j][i] = 0.0 + bij;
804: } else {
805: const PetscScalar
806: u_E = x[j][i+1],
807: u_W = x[j][i-1],
808: u_N = x[j+1][i],
809: u_S = x[j-1][i];
810: const PetscScalar
811: uy_E = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]),
812: uy_W = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]),
813: ux_N = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]),
814: ux_S = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]);
815: u = x[j][i];
816: for (k=0; k<its; k++) {
817: const PetscScalar
818: ux_E = dhx*(u_E-u),
819: ux_W = dhx*(u-u_W),
820: uy_N = dhy*(u_N-u),
821: uy_S = dhy*(u-u_S),
822: e_E = eta(user,xx,yy,ux_E,uy_E),
823: e_W = eta(user,xx,yy,ux_W,uy_W),
824: e_N = eta(user,xx,yy,ux_N,uy_N),
825: e_S = eta(user,xx,yy,ux_S,uy_S),
826: de_E = deta(user,xx,yy,ux_E,uy_E),
827: de_W = deta(user,xx,yy,ux_W,uy_W),
828: de_N = deta(user,xx,yy,ux_N,uy_N),
829: de_S = deta(user,xx,yy,ux_S,uy_S),
830: newt_E = e_E+de_E*PetscSqr(ux_E),
831: newt_W = e_W+de_W*PetscSqr(ux_W),
832: newt_N = e_N+de_N*PetscSqr(uy_N),
833: newt_S = e_S+de_S*PetscSqr(uy_S),
834: uxx = -hy * (e_E*ux_E - e_W*ux_W),
835: uyy = -hx * (e_N*uy_N - e_S*uy_S);
837: if (sc) eu = PetscExpScalar(u);
838: else eu = 0;
840: F = uxx + uyy - sc*eu - bij;
841: if (k == 0) F0 = F;
842: J = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*eu;
843: y = F/J;
844: u -= y;
845: tot_its++;
846: if (atol > PetscAbsReal(PetscRealPart(F)) ||
847: rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) ||
848: stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) {
849: break;
850: }
851: }
852: x[j][i] = u;
853: }
854: }
855: }
856: }
857: /*
858: x Restore vector
859: */
860: }
861: DMDAVecRestoreArray(da,localX,&x);
862: DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
863: DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
864: PetscLogFlops(tot_its*(118.0));
865: DMRestoreLocalVector(da,&localX);
866: if (B) {
867: DMDAVecRestoreArrayRead(da,localB,&b);
868: DMRestoreLocalVector(da,&localB);
869: }
870: return(0);
871: }
873: /*TEST
875: test:
876: nsize: 2
877: args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype NEWTON
878: requires: !single
880: test:
881: suffix: 2
882: nsize: 2
883: args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -precheck 1
884: requires: !single
886: test:
887: suffix: 3
888: nsize: 2
889: args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -picard -precheck 1 -p 3
890: requires: !single
892: test:
893: suffix: 3_star
894: nsize: 2
895: 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
896: output_file: output/ex15_3.out
897: requires: !single
899: test:
900: suffix: 4
901: 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
902: requires: !single
904: test:
905: suffix: lag_jac
906: nsize: 4
907: 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
908: requires: !single
910: test:
911: suffix: lag_pc
912: nsize: 4
913: 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
914: requires: !single
916: test:
917: suffix: nleqerr
918: 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
919: requires: !single
921: test:
922: suffix: mf
923: args: -snes_monitor_short -pc_type lu -da_refine 4 -p 3 -ksp_rtol 1.e-12 -snes_mf_operator
924: requires: !single
926: test:
927: suffix: mf_picard
928: args: -snes_monitor_short -pc_type lu -da_refine 4 -p 3 -ksp_rtol 1.e-12 -snes_mf_operator -picard
929: requires: !single
930: output_file: output/ex15_mf.out
932: test:
933: suffix: fd_picard
934: args: -snes_monitor_short -pc_type lu -da_refine 2 -p 3 -ksp_rtol 1.e-12 -snes_fd -picard
935: requires: !single
937: test:
938: suffix: fd_color_picard
939: args: -snes_monitor_short -pc_type lu -da_refine 4 -p 3 -ksp_rtol 1.e-12 -snes_fd_color -picard
940: requires: !single
941: output_file: output/ex15_mf.out
943: TEST*/