Actual source code: ex15.c
petsc-3.11.1 2019-04-12
1: static const char help[] = " p-Bratu nonlinear PDE in 2d.\n\
2: We solve the p-Laplacian ( nonlinear diffusion) combined with\n\
3: the Bratu ( solid fuel ignition) nonlinearity in a 2D rectangular\n\
4: domain, using distributed arrays ( DAs) to partition the parallel grid.\n\
5: The command line options include:\n\
6: -p <2>: `p' in p-Laplacian term\n\
7: -epsilon <1e-05>: Strain-regularization in p-Laplacian\n\
8: -lambda <6>: Bratu parameter\n\
9: -blocks <bx,by>: number of coefficient interfaces in x and y direction\n\
10: -kappa <1e-3>: diffusivity in odd regions\n\
11 : \n" ;
The $p$-Bratu problem is a combination of the $p$-Laplacian ( nonlinear diffusion) and the Brutu solid fuel ignition problem.
This problem is modeled by the partial differential equation
\begin{equation*}
-\nabla\cdot ( \eta \nabla u) - \lambda \exp( u) = 0
\end{equation*}
on $\Omega = ( -1,1)^2$ with closure
\begin{align*}
\eta( \gamma) &= ( \epsilon^2 + \gamma)^{( p-2)/2} & \gamma &= \frac 1 2 |\nabla u|^2
\end{align*}
and boundary conditions $u = 0$ for $( x,y) \in \partial \Omega$
A 9-point finite difference stencil is used to discretize
the boundary value problem to obtain a nonlinear system of equations.
This would be a 5-point stencil if not for the $p$-Laplacian's nonlinearity.
35: /*
36: mpiexec -n 2 ./ex15 -snes_monitor -ksp_monitor log_summary
37: */
39: /*
40: Include " petscdmda.h" so that we can use distributed arrays ( DMDAs).
41: Include " petscsnes.h" so that we can use SNES solvers. Note that this
42: file automatically includes:
43: petsc.h - base PETSc routines petscvec.h - vectors
44: petscsys.h - system routines petscmat.h - matrices
45: petscis.h - index sets petscksp.h - Krylov subspace methods
46: petscviewer.h - viewers petscpc.h - preconditioners
47: petscksp.h - linear solvers
48: */
50: #include <petscdm.h>
51: #include <petscdmda.h>
52: #include <petscsnes.h>
54: typedef enum {JAC_BRATU,JAC_PICARD,JAC_STAR,JAC_NEWTON} JacType;
55: static const char *const JacTypes[] = {" BRATU" ," PICARD" ," STAR" ," NEWTON" ," JacType" ," JAC_" ,0};
57: /*
58: User-defined application context - contains data needed by the
59: application-provided call-back routines, FormJacobianLocal( ) and
60: FormFunctionLocal( ).
61: */
62: typedef struct {
63: PetscReal lambda; /* Bratu parameter */
64: PetscReal p; /* Exponent in p-Laplacian */
65: PetscReal epsilon; /* Regularization */
66: PetscReal source; /* Source term */
67: JacType jtype; /* What type of Jacobian to assemble */
68: PetscBool picard;
69: PetscInt blocks[2];
70: PetscReal kappa;
71: PetscInt initial; /* initial conditions type */
72: } AppCtx;
74: /*
75: User-defined routines
76: */
77: static PetscErrorCode FormRHS( AppCtx*,DM ,Vec ) ;
78: static PetscErrorCode FormInitialGuess( AppCtx*,DM ,Vec ) ;
79: static PetscErrorCode FormFunctionLocal( DMDALocalInfo *,PetscScalar **,PetscScalar **,AppCtx*) ;
80: static PetscErrorCode FormFunctionPicardLocal( DMDALocalInfo *,PetscScalar **,PetscScalar **,AppCtx*) ;
81: static PetscErrorCode FormJacobianLocal( DMDALocalInfo *,PetscScalar **,Mat ,Mat ,AppCtx*) ;
82: static PetscErrorCode NonlinearGS( SNES ,Vec ,Vec ,void*) ;
84: typedef struct _n_PreCheck *PreCheck;
85: struct _n_PreCheck {
86: MPI_Comm comm;
87: PetscReal angle;
88: Vec Ylast;
89: PetscViewer monitor;
90: };
91: PetscErrorCode PreCheckCreate( MPI_Comm ,PreCheck*) ;
92: PetscErrorCode PreCheckDestroy( PreCheck*) ;
93: PetscErrorCode PreCheckFunction( SNESLineSearch ,Vec ,Vec ,PetscBool *,void*) ;
94: PetscErrorCode PreCheckSetFromOptions( PreCheck) ;
96: int main( int argc,char **argv)
97: {
98: SNES snes; /* nonlinear solver */
99: Vec x,r,b; /* solution, residual, rhs vectors */
100: Mat A,B; /* Jacobian and preconditioning matrices */
101: AppCtx user; /* user-defined work context */
102: PetscInt its; /* iterations for convergence */
103: SNESConvergedReason reason; /* Check convergence */
104: PetscBool alloc_star; /* Only allocate for the STAR stencil */
105: PetscBool write_output;
106: char filename[PETSC_MAX_PATH_LEN] = " ex15.vts" ;
107: PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
108: DM da,dastar; /* distributed array data structure */
109: PreCheck precheck = NULL; /* precheck context for version in this file */
110: PetscInt use_precheck; /* 0= none, 1= version in this file, 2= SNES -provided version */
111: PetscReal precheck_angle; /* When manually setting the SNES -provided precheck function */
112: PetscErrorCode ierr;
113: SNESLineSearch linesearch;
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Initialize program
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: PetscInitialize ( &argc,&argv,( char*)0,help);if ( ierr) return ierr;
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Initialize problem parameters
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124: user.lambda = 0.0; user.p = 2.0; user.epsilon = 1e-5; user.source = 0.1; user.jtype = JAC_NEWTON;user.initial= -1;
125: user.blocks[0] = 1; user.blocks[1] = 1; user.kappa = 1e-3;
126: alloc_star = PETSC_FALSE ;
127: use_precheck = 0; precheck_angle = 10.;
128: user.picard = PETSC_FALSE ;
129: PetscOptionsBegin ( PETSC_COMM_WORLD ,NULL," p-Bratu options" ,__FILE__);
130: {
131: PetscInt two= 2;
132: PetscOptionsReal ( " -lambda" ," Bratu parameter" ," " ,user.lambda,&user.lambda,NULL);
133: PetscOptionsReal ( " -p" ," Exponent `p' in p-Laplacian" ," " ,user.p,&user.p,NULL);
134: PetscOptionsReal ( " -epsilon" ," Strain-regularization in p-Laplacian" ," " ,user.epsilon,&user.epsilon,NULL);
135: PetscOptionsReal ( " -source" ," Constant source term" ," " ,user.source,&user.source,NULL);
136: PetscOptionsEnum ( " -jtype" ," Jacobian approximation to assemble" ," " ,JacTypes,( PetscEnum )user.jtype,( PetscEnum *)&user.jtype,NULL);
137: PetscOptionsName ( " -picard" ," Solve with defect-correction Picard iteration" ," " ,&user.picard);
138: if ( user.picard) {user.jtype = JAC_PICARD; user.p = 3;}
139: PetscOptionsBool ( " -alloc_star" ," Allocate for STAR stencil ( 5-point)" ," " ,alloc_star,&alloc_star,NULL);
140: PetscOptionsInt ( " -precheck" ," Use a pre-check correction intended for use with Picard iteration 1= this version, 2= library" ," " ,use_precheck,&use_precheck,NULL);
141: PetscOptionsInt ( " -initial" ," Initial conditions type ( -1: default, 0: zero-valued, 1: peaked guess)" ," " ,user.initial,&user.initial,NULL);
142: if ( use_precheck = = 2) { /* Using library version, get the angle */
143: PetscOptionsReal ( " -precheck_angle" ," Angle in degrees between successive search directions necessary to activate step correction" ," " ,precheck_angle,&precheck_angle,NULL);
144: }
145: PetscOptionsIntArray ( " -blocks" ," number of coefficient interfaces in x and y direction" ," " ,user.blocks,&two,NULL);
146: PetscOptionsReal ( " -kappa" ," diffusivity in odd regions" ," " ,user.kappa,&user.kappa,NULL);
147: PetscOptionsString ( " -o" ," Output solution in vts format" ," " ,filename,filename,sizeof ( filename),&write_output);
148: }
149: PetscOptionsEnd ( );
150: if ( user.lambda > bratu_lambda_max || user.lambda < bratu_lambda_min) {
151: PetscPrintf ( PETSC_COMM_WORLD ," WARNING: lambda %g out of range for p= 2\n" ,user.lambda);
152: }
154: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: Create nonlinear solver context
156: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157: SNESCreate ( PETSC_COMM_WORLD ,&snes);
159: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: Create distributed array ( DMDA ) to manage parallel grid and vectors
161: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
162: DMDACreate2d ( PETSC_COMM_WORLD ,DM_BOUNDARY_NONE ,DM_BOUNDARY_NONE ,DMDA_STENCIL_BOX ,4,4,PETSC_DECIDE ,PETSC_DECIDE ,1,1,NULL,NULL,&da);
163: DMSetFromOptions ( da);
164: DMSetUp ( da);
165: DMDACreate2d ( PETSC_COMM_WORLD ,DM_BOUNDARY_NONE ,DM_BOUNDARY_NONE ,DMDA_STENCIL_STAR ,4,4,PETSC_DECIDE ,PETSC_DECIDE ,1,1,NULL,NULL,&dastar);
166: DMSetFromOptions ( dastar);
167: DMSetUp ( dastar);
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: Create matrix data structure; set Jacobian evaluation routine
180: Set Jacobian matrix data structure and default Jacobian evaluation
181: routine. User can override with:
182: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
183: ( unless user explicitly sets preconditioner)
184: -snes_mf_operator : form preconditioning matrix as set by the user,
185: but use matrix-free approx for Jacobian-vector
186: products within Newton-Krylov method
188: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
189: /* B can be type of MATAIJ ,MATBAIJ or MATSBAIJ */
190: DMCreateMatrix ( alloc_star ? dastar : da,&B);
191: A = B;
193: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
194: Set local function evaluation routine
195: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
196: DMSetApplicationContext ( da, &user);
197: SNESSetDM ( snes,da);
198: if ( user.picard) {
199: /*
200: This is not really right requiring the user to call SNESSetFunction /Jacobian but the DMDASNESSetPicardLocal ( ) cannot access
201: the SNES to set it
202: */
203: DMDASNESSetPicardLocal ( da,INSERT_VALUES ,( PetscErrorCode ( *)( DMDALocalInfo *,void*,void*,void*))FormFunctionPicardLocal,
204: ( PetscErrorCode ( *)( DMDALocalInfo *,void*,Mat ,Mat ,void*))FormJacobianLocal,&user);
205: SNESSetFunction ( snes,NULL,SNESPicardComputeFunction,&user);
206: SNESSetJacobian ( snes,NULL,NULL,SNESPicardComputeJacobian,&user);
207: } else {
208: DMDASNESSetFunctionLocal ( da,INSERT_VALUES ,( PetscErrorCode ( *)( DMDALocalInfo *,void*,void*,void*))FormFunctionLocal,&user);
209: DMDASNESSetJacobianLocal ( da,( PetscErrorCode ( *)( DMDALocalInfo *,void*,Mat ,Mat ,void*))FormJacobianLocal,&user);
210: }
213: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
214: Customize nonlinear solver; set runtime options
215: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216: SNESSetFromOptions ( snes);
217: SNESSetNGS ( snes,NonlinearGS,&user);
218: SNESGetLineSearch ( snes, &linesearch);
219: /* Set up the precheck context if requested */
220: if ( use_precheck = = 1) { /* Use the precheck routines in this file */
221: PreCheckCreate( PETSC_COMM_WORLD ,&precheck);
222: PreCheckSetFromOptions( precheck);
223: SNESLineSearchSetPreCheck ( linesearch,PreCheckFunction,precheck);
224: } else if ( use_precheck = = 2) { /* Use the version provided by the library */
225: SNESLineSearchSetPreCheck ( linesearch,SNESLineSearchPreCheckPicard ,&precheck_angle);
226: }
228: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
229: Evaluate initial guess
230: Note: The user should initialize the vector, x, with the initial guess
231: for the nonlinear solver prior to calling SNESSolve ( ). In particular,
232: to employ an initial guess of zero, the user should explicitly set
233: this vector to zero by calling VecSet ( ).
234: */
236: FormInitialGuess( &user,da,x);
237: FormRHS( &user,da,b);
239: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
240: Solve nonlinear system
241: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
242: SNESSolve ( snes,b,x);
243: SNESGetIterationNumber ( snes,&its);
244: SNESGetConvergedReason ( snes,&reason);
246: PetscPrintf ( PETSC_COMM_WORLD ," %s Number of nonlinear iterations = %D\n" ,SNESConvergedReasons[reason],its);
248: if ( write_output) {
249: PetscViewer viewer;
250: PetscViewerVTKOpen ( PETSC_COMM_WORLD ,filename,FILE_MODE_WRITE ,&viewer);
251: VecView ( x,viewer);
252: PetscViewerDestroy ( &viewer);
253: }
255: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
256: Free work space. All PETSc objects should be destroyed when they
257: are no longer needed.
258: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
260: if ( A != B) {
261: MatDestroy ( &A);
262: }
263: MatDestroy ( &B);
264: VecDestroy ( &x);
265: VecDestroy ( &r);
266: VecDestroy ( &b);
267: SNESDestroy ( &snes);
268: DMDestroy ( &da);
269: DMDestroy ( &dastar);
270: PreCheckDestroy( &precheck);
271: PetscFinalize ( );
272: return ierr;
273: }
275: /* ------------------------------------------------------------------- */
276: /*
277: FormInitialGuess - Forms initial approximation.
279: Input Parameters:
280: user - user-defined application context
281: X - vector
283: Output Parameter:
284: X - vector
285: */
286: static PetscErrorCode FormInitialGuess( AppCtx *user,DM da,Vec X)
287: {
288: PetscInt i,j,Mx,My,xs,ys,xm,ym;
290: PetscReal temp1,temp,hx,hy;
291: PetscScalar **x;
294: 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 );
296: hx = 1.0/( PetscReal )( Mx-1);
297: hy = 1.0/( PetscReal )( My-1);
298: temp1 = user->lambda / ( user->lambda + 1.);
300: /*
301: Get a pointer to vector data.
302: - For default PETSc vectors, VecGetArray ( ) returns a pointer to
303: the data array. Otherwise, the routine is implementation dependent.
304: - You MUST call VecRestoreArray ( ) when you no longer need access to
305: the array.
306: */
307: DMDAVecGetArray ( da,X,&x);
309: /*
310: Get local grid boundaries ( for 2-dimensional DA):
311: xs, ys - starting grid indices ( no ghost points)
312: xm, ym - widths of local grid ( no ghost points)
314: */
315: DMDAGetCorners ( da,&xs,&ys,NULL,&xm,&ym,NULL);
317: /*
318: Compute initial guess over the locally owned part of the grid
319: */
320: for ( j= ys; j<ys+ym; j++) {
321: temp = ( PetscReal )( PetscMin ( j,My-j-1))*hy;
322: for ( i= xs; i<xs+xm; i++) {
323: if ( i = = 0 || j = = 0 || i = = Mx-1 || j = = My-1) {
324: /* boundary conditions are all zero Dirichlet */
325: x[j][i] = 0.0;
326: } else {
327: if ( user->initial = = -1) {
328: if ( user->lambda != 0) {
329: x[j][i] = temp1*PetscSqrtReal( PetscMin ( ( PetscReal )( PetscMin ( i,Mx-i-1))*hx,temp));
330: } else {
331: /* The solution above is an exact solution for lambda= 0, this avoids " accidentally" starting
332: * with an exact solution. */
333: const PetscReal
334: xx = 2*( PetscReal )i/( Mx-1) - 1,
335: yy = 2*( PetscReal )j/( My-1) - 1;
336: x[j][i] = ( 1 - xx*xx) * ( 1-yy*yy) * xx * yy;
337: }
338: } else if ( user->initial = = 0) {
339: x[j][i] = 0.;
340: } else if ( user->initial = = 1) {
341: const PetscReal
342: xx = 2*( PetscReal )i/( Mx-1) - 1,
343: yy = 2*( PetscReal )j/( My-1) - 1;
344: x[j][i] = ( 1 - xx*xx) * ( 1-yy*yy) * xx * yy;
345: } else {
346: if ( user->lambda != 0) {
347: x[j][i] = temp1*PetscSqrtReal( PetscMin ( ( PetscReal )( PetscMin ( i,Mx-i-1))*hx,temp));
348: } else {
349: x[j][i] = 0.5*PetscSqrtReal( PetscMin ( ( PetscReal )( PetscMin ( i,Mx-i-1))*hx,temp));
350: }
351: }
352: }
353: }
354: }
355: /*
356: Restore vector
357: */
358: DMDAVecRestoreArray ( da,X,&x);
359: return ( 0);
360: }
362: /*
363: FormRHS - Forms constant RHS for the problem.
365: Input Parameters:
366: user - user-defined application context
367: B - RHS vector
369: Output Parameter:
370: B - vector
371: */
372: static PetscErrorCode FormRHS( AppCtx *user,DM da,Vec B)
373: {
374: PetscInt i,j,Mx,My,xs,ys,xm,ym;
376: PetscReal hx,hy;
377: PetscScalar **b;
380: 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 );
382: hx = 1.0/( PetscReal )( Mx-1);
383: hy = 1.0/( PetscReal )( My-1);
384: DMDAVecGetArray ( da,B,&b);
385: DMDAGetCorners ( da,&xs,&ys,NULL,&xm,&ym,NULL);
386: for ( j= ys; j<ys+ym; j++) {
387: for ( i= xs; i<xs+xm; i++) {
388: if ( i = = 0 || j = = 0 || i = = Mx-1 || j = = My-1) {
389: b[j][i] = 0.0;
390: } else {
391: b[j][i] = hx*hy*user->source;
392: }
393: }
394: }
395: DMDAVecRestoreArray ( da,B,&b);
396: return ( 0);
397: }
399: PETSC_STATIC_INLINE PetscReal kappa( const AppCtx *ctx,PetscReal x,PetscReal y)
400: {
401: return ( ( ( PetscInt )( x*ctx->blocks[0])) + ( ( PetscInt )( y*ctx->blocks[1]))) % 2 ? ctx->kappa : 1.0;
402: }
403: /* p-Laplacian diffusivity */
404: PETSC_STATIC_INLINE PetscScalar eta( const AppCtx *ctx,PetscReal x,PetscReal y,PetscScalar ux,PetscScalar uy)
405: {
406: return kappa( ctx,x,y) * PetscPowScalar( PetscSqr ( ctx->epsilon)+0.5*( ux*ux + uy*uy),0.5*( ctx->p-2.));
407: }
408: PETSC_STATIC_INLINE PetscScalar deta( const AppCtx *ctx,PetscReal x,PetscReal y,PetscScalar ux,PetscScalar uy)
409: {
410: return ( ctx->p = = 2)
411: ? 0
412: : kappa( ctx,x,y)*PetscPowScalar( PetscSqr ( ctx->epsilon)+0.5*( ux*ux + uy*uy),0.5*( ctx->p-4)) * 0.5 * ( ctx->p-2.);
413: }
416: /* ------------------------------------------------------------------- */
417: /*
418: FormFunctionLocal - Evaluates nonlinear function, F( x).
419: */
420: static PetscErrorCode FormFunctionLocal( DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
421: {
422: PetscReal hx,hy,dhx,dhy,sc;
423: PetscInt i,j;
424: PetscScalar eu;
429: hx = 1.0/( PetscReal )( info->mx-1);
430: hy = 1.0/( PetscReal )( info->my-1);
431: sc = hx*hy*user->lambda;
432: dhx = 1/hx;
433: dhy = 1/hy;
434: /*
435: Compute function over the locally owned part of the grid
436: */
437: for ( j= info->ys; j<info->ys+info->ym; j++) {
438: for ( i= info->xs; i<info->xs+info->xm; i++) {
439: PetscReal xx = i*hx,yy = j*hy;
440: if ( i = = 0 || j = = 0 || i = = info->mx-1 || j = = info->my-1) {
441: f[j][i] = x[j][i];
442: } else {
443: const PetscScalar
444: u = x[j][i],
445: ux_E = dhx*( x[j][i+1]-x[j][i]),
446: uy_E = 0.25*dhy*( x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]),
447: ux_W = dhx*( x[j][i]-x[j][i-1]),
448: uy_W = 0.25*dhy*( x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]),
449: ux_N = 0.25*dhx*( x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]),
450: uy_N = dhy*( x[j+1][i]-x[j][i]),
451: ux_S = 0.25*dhx*( x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]),
452: uy_S = dhy*( x[j][i]-x[j-1][i]),
453: e_E = eta( user,xx,yy,ux_E,uy_E),
454: e_W = eta( user,xx,yy,ux_W,uy_W),
455: e_N = eta( user,xx,yy,ux_N,uy_N),
456: e_S = eta( user,xx,yy,ux_S,uy_S),
457: uxx = -hy * ( e_E*ux_E - e_W*ux_W),
458: uyy = -hx * ( e_N*uy_N - e_S*uy_S);
459: if ( sc) eu = PetscExpScalar( u);
460: else eu = 0.;
461: /** For p= 2, these terms decay to:
462: * uxx = ( 2.0*u - x[j][i-1] - x[j][i+1])*hydhx
463: * uyy = ( 2.0*u - x[j-1][i] - x[j+1][i])*hxdhy
464: **/
465: f[j][i] = uxx + uyy - sc*eu;
466: }
467: }
468: }
469: PetscLogFlops ( info->xm*info->ym*( 72.0));
470: return ( 0);
471: }
473: /*
474: This is the opposite sign of the part of FormFunctionLocal that excludes the A( x) x part of the operation,
475: 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)
477: */
478: static PetscErrorCode FormFunctionPicardLocal( DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
479: {
480: PetscReal hx,hy,sc;
481: PetscInt i,j;
485: hx = 1.0/( PetscReal )( info->mx-1);
486: hy = 1.0/( PetscReal )( info->my-1);
487: sc = hx*hy*user->lambda;
488: /*
489: Compute function over the locally owned part of the grid
490: */
491: for ( j= info->ys; j<info->ys+info->ym; j++) {
492: for ( i= info->xs; i<info->xs+info->xm; i++) {
493: if ( !( i = = 0 || j = = 0 || i = = info->mx-1 || j = = info->my-1)) {
494: const PetscScalar u = x[j][i];
495: f[j][i] = sc*PetscExpScalar( u);
496: }
497: }
498: }
499: PetscLogFlops ( info->xm*info->ym*2.0);
500: return ( 0);
501: }
503: /*
504: FormJacobianLocal - Evaluates Jacobian matrix.
505: */
506: static PetscErrorCode FormJacobianLocal( DMDALocalInfo *info,PetscScalar **x,Mat J,Mat B,AppCtx *user)
507: {
509: PetscInt i,j;
510: MatStencil col[9],row;
511: PetscScalar v[9];
512: PetscReal hx,hy,hxdhy,hydhx,dhx,dhy,sc;
515: hx = 1.0/( PetscReal )( info->mx-1);
516: hy = 1.0/( PetscReal )( info->my-1);
517: sc = hx*hy*user->lambda;
518: dhx = 1/hx;
519: dhy = 1/hy;
520: hxdhy = hx/hy;
521: hydhx = hy/hx;
523: /*
524: Compute entries for the locally owned part of the Jacobian.
525: - PETSc parallel matrix formats are partitioned by
526: contiguous chunks of rows across the processors.
527: - Each processor needs to insert only elements that it owns
528: locally ( but any non-local elements will be sent to the
529: appropriate processor during matrix assembly).
530: - Here, we set all entries for a particular row at once.
531: */
532: for ( j= info->ys; j<info->ys+info->ym; j++) {
533: for ( i= info->xs; i<info->xs+info->xm; i++) {
534: PetscReal xx = i*hx,yy = j*hy;
535: row.j = j; row.i = i;
536: /* boundary points */
537: if ( i = = 0 || j = = 0 || i = = info->mx-1 || j = = info->my-1) {
538: v[0] = 1.0;
539: MatSetValuesStencil ( B,1,&row,1,&row,v,INSERT_VALUES );
540: } else {
541: /* interior grid points */
542: const PetscScalar
543: ux_E = dhx*( x[j][i+1]-x[j][i]),
544: uy_E = 0.25*dhy*( x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]),
545: ux_W = dhx*( x[j][i]-x[j][i-1]),
546: uy_W = 0.25*dhy*( x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]),
547: ux_N = 0.25*dhx*( x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]),
548: uy_N = dhy*( x[j+1][i]-x[j][i]),
549: ux_S = 0.25*dhx*( x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]),
550: uy_S = dhy*( x[j][i]-x[j-1][i]),
551: u = x[j][i],
552: e_E = eta( user,xx,yy,ux_E,uy_E),
553: e_W = eta( user,xx,yy,ux_W,uy_W),
554: e_N = eta( user,xx,yy,ux_N,uy_N),
555: e_S = eta( user,xx,yy,ux_S,uy_S),
556: de_E = deta( user,xx,yy,ux_E,uy_E),
557: de_W = deta( user,xx,yy,ux_W,uy_W),
558: de_N = deta( user,xx,yy,ux_N,uy_N),
559: de_S = deta( user,xx,yy,ux_S,uy_S),
560: skew_E = de_E*ux_E*uy_E,
561: skew_W = de_W*ux_W*uy_W,
562: skew_N = de_N*ux_N*uy_N,
563: skew_S = de_S*ux_S*uy_S,
564: cross_EW = 0.25*( skew_E - skew_W),
565: cross_NS = 0.25*( skew_N - skew_S),
566: newt_E = e_E+de_E*PetscSqr ( ux_E),
567: newt_W = e_W+de_W*PetscSqr ( ux_W),
568: newt_N = e_N+de_N*PetscSqr ( uy_N),
569: newt_S = e_S+de_S*PetscSqr ( uy_S);
570: /* interior grid points */
571: switch ( user->jtype) {
572: case JAC_BRATU:
573: /* Jacobian from p= 2 */
574: v[0] = -hxdhy; col[0].j = j-1; col[0].i = i;
575: v[1] = -hydhx; col[1].j = j; col[1].i = i-1;
576: v[2] = 2.0*( hydhx + hxdhy) - sc*PetscExpScalar( u); col[2].j = row.j; col[2].i = row.i;
577: v[3] = -hydhx; col[3].j = j; col[3].i = i+1;
578: v[4] = -hxdhy; col[4].j = j+1; col[4].i = i;
579: MatSetValuesStencil ( B,1,&row,5,col,v,INSERT_VALUES );
580: break ;
581: case JAC_PICARD:
582: /* Jacobian arising from Picard linearization */
583: v[0] = -hxdhy*e_S; col[0].j = j-1; col[0].i = i;
584: v[1] = -hydhx*e_W; col[1].j = j; col[1].i = i-1;
585: v[2] = ( e_W+e_E)*hydhx + ( e_S+e_N)*hxdhy; col[2].j = row.j; col[2].i = row.i;
586: v[3] = -hydhx*e_E; col[3].j = j; col[3].i = i+1;
587: v[4] = -hxdhy*e_N; col[4].j = j+1; col[4].i = i;
588: MatSetValuesStencil ( B,1,&row,5,col,v,INSERT_VALUES );
589: break ;
590: case JAC_STAR:
591: /* Full Jacobian, but only a star stencil */
592: col[0].j = j-1; col[0].i = i;
593: col[1].j = j; col[1].i = i-1;
594: col[2].j = j; col[2].i = i;
595: col[3].j = j; col[3].i = i+1;
596: col[4].j = j+1; col[4].i = i;
597: v[0] = -hxdhy*newt_S + cross_EW;
598: v[1] = -hydhx*newt_W + cross_NS;
599: v[2] = hxdhy*( newt_N + newt_S) + hydhx*( newt_E + newt_W) - sc*PetscExpScalar( u);
600: v[3] = -hydhx*newt_E - cross_NS;
601: v[4] = -hxdhy*newt_N - cross_EW;
602: MatSetValuesStencil ( B,1,&row,5,col,v,INSERT_VALUES );
603: break ;
604: case JAC_NEWTON:
605: /** The Jacobian is
606: *
607: * -div [ eta ( grad u) + deta ( grad u0 . grad u) grad u0 ] - ( eE u0) u
608: *
609: **/
610: col[0].j = j-1; col[0].i = i-1;
611: col[1].j = j-1; col[1].i = i;
612: col[2].j = j-1; col[2].i = i+1;
613: col[3].j = j; col[3].i = i-1;
614: col[4].j = j; col[4].i = i;
615: col[5].j = j; col[5].i = i+1;
616: col[6].j = j+1; col[6].i = i-1;
617: col[7].j = j+1; col[7].i = i;
618: col[8].j = j+1; col[8].i = i+1;
619: v[0] = -0.25*( skew_S + skew_W);
620: v[1] = -hxdhy*newt_S + cross_EW;
621: v[2] = 0.25*( skew_S + skew_E);
622: v[3] = -hydhx*newt_W + cross_NS;
623: v[4] = hxdhy*( newt_N + newt_S) + hydhx*( newt_E + newt_W) - sc*PetscExpScalar( u);
624: v[5] = -hydhx*newt_E - cross_NS;
625: v[6] = 0.25*( skew_N + skew_W);
626: v[7] = -hxdhy*newt_N - cross_EW;
627: v[8] = -0.25*( skew_N + skew_E);
628: MatSetValuesStencil ( B,1,&row,9,col,v,INSERT_VALUES );
629: break ;
630: default:
631: SETERRQ1 ( PetscObjectComm ( ( PetscObject )info->da),PETSC_ERR_SUP," Jacobian type %d not implemented" ,user->jtype);
632: }
633: }
634: }
635: }
637: /*
638: Assemble matrix, using the 2-step process:
639: MatAssemblyBegin ( ), MatAssemblyEnd ( ).
640: */
641: MatAssemblyBegin ( B,MAT_FINAL_ASSEMBLY );
642: MatAssemblyEnd ( B,MAT_FINAL_ASSEMBLY );
644: if ( J != B) {
645: MatAssemblyBegin ( J,MAT_FINAL_ASSEMBLY );
646: MatAssemblyEnd ( J,MAT_FINAL_ASSEMBLY );
647: }
649: /*
650: Tell the matrix we will never add a new nonzero location to the
651: matrix. If we do, it will generate an error.
652: */
653: if ( user->jtype = = JAC_NEWTON) {
654: PetscLogFlops ( info->xm*info->ym*( 131.0));
655: }
656: MatSetOption ( B,MAT_NEW_NONZERO_LOCATION_ERR ,PETSC_TRUE );
657: return ( 0);
658: }
660: /***********************************************************
661: * PreCheck implementation
662: ***********************************************************/
663: PetscErrorCode PreCheckSetFromOptions( PreCheck precheck)
664: {
666: PetscBool flg;
669: PetscOptionsBegin ( precheck->comm,NULL," PreCheck Options" ," none" );
670: PetscOptionsReal ( " -precheck_angle" ," Angle in degrees between successive search directions necessary to activate step correction" ," " ,precheck->angle,&precheck->angle,NULL);
671: flg = PETSC_FALSE ;
672: PetscOptionsBool ( " -precheck_monitor" ," Monitor choices made by precheck routine" ," " ,flg,&flg,NULL);
673: if ( flg) {
674: PetscViewerASCIIOpen ( precheck->comm," stdout" ,&precheck->monitor);
675: }
676: PetscOptionsEnd ( );
677: return ( 0);
678: }
680: /*
681: Compare the direction of the current and previous step, modify the current step accordingly
682: */
683: PetscErrorCode PreCheckFunction( SNESLineSearch linesearch,Vec X,Vec Y,PetscBool *changed, void *ctx)
684: {
686: PreCheck precheck;
687: Vec Ylast;
688: PetscScalar dot;
689: PetscInt iter;
690: PetscReal ynorm,ylastnorm,theta,angle_radians;
691: SNES snes;
694: SNESLineSearchGetSNES ( linesearch, &snes);
695: precheck = ( PreCheck)ctx;
696: if ( !precheck->Ylast) {VecDuplicate ( Y,&precheck->Ylast);}
697: Ylast = precheck->Ylast;
698: SNESGetIterationNumber ( snes,&iter);
699: if ( iter < 1) {
700: VecCopy ( Y,Ylast);
701: *changed = PETSC_FALSE ;
702: return ( 0);
703: }
705: VecDot ( Y,Ylast,&dot);
706: VecNorm ( Y,NORM_2 ,&ynorm);
707: VecNorm ( Ylast,NORM_2 ,&ylastnorm);
708: /* Compute the angle between the vectors Y and Ylast, clip to keep inside the domain of acos( ) */
709: theta = PetscAcosReal( ( PetscReal )PetscClipInterval ( PetscAbsScalar( dot) / ( ynorm * ylastnorm),-1.0,1.0));
710: angle_radians = precheck->angle * PETSC_PI / 180.;
711: if ( PetscAbsReal ( theta) < angle_radians || PetscAbsReal ( theta - PETSC_PI) < angle_radians) {
712: /* Modify the step Y */
713: PetscReal alpha,ydiffnorm;
714: VecAXPY ( Ylast,-1.0,Y);
715: VecNorm ( Ylast,NORM_2 ,&ydiffnorm);
716: alpha = ylastnorm / ydiffnorm;
717: VecCopy ( Y,Ylast);
718: VecScale ( Y,alpha);
719: if ( precheck->monitor) {
720: 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);
721: }
722: } else {
723: VecCopy ( Y,Ylast);
724: *changed = PETSC_FALSE ;
725: if ( precheck->monitor) {
726: PetscViewerASCIIPrintf ( precheck->monitor," Angle %E degrees exceeds threshold %g, no correction applied\n" ,( double)( theta*180./PETSC_PI),( double)precheck->angle);
727: }
728: }
729: return ( 0);
730: }
732: PetscErrorCode PreCheckDestroy( PreCheck *precheck)
733: {
737: if ( !*precheck) return ( 0);
738: VecDestroy ( &( *precheck)->Ylast);
739: PetscViewerDestroy ( &( *precheck)->monitor);
740: PetscFree ( *precheck);
741: return ( 0);
742: }
744: PetscErrorCode PreCheckCreate( MPI_Comm comm,PreCheck *precheck)
745: {
749: PetscNew ( precheck);
751: ( *precheck)->comm = comm;
752: ( *precheck)->angle = 10.; /* only active if angle is less than 10 degrees */
753: return ( 0);
754: }
756: /*
757: Applies some sweeps on nonlinear Gauss-Seidel on each process
759: */
760: PetscErrorCode NonlinearGS( SNES snes,Vec X, Vec B, void *ctx)
761: {
762: PetscInt i,j,k,xs,ys,xm,ym,its,tot_its,sweeps,l,m;
764: PetscReal hx,hy,hxdhy,hydhx,dhx,dhy,sc;
765: PetscScalar **x,**b,bij,F,F0= 0,J,y,u,eu;
766: PetscReal atol,rtol,stol;
767: DM da;
768: AppCtx *user = ( AppCtx*)ctx;
769: Vec localX,localB;
770: DMDALocalInfo info;
773: SNESGetDM ( snes,&da);
774: DMDAGetLocalInfo ( da,&info);
776: hx = 1.0/( PetscReal )( info.mx-1);
777: hy = 1.0/( PetscReal )( info.my-1);
778: sc = hx*hy*user->lambda;
779: dhx = 1/hx;
780: dhy = 1/hy;
781: hxdhy = hx/hy;
782: hydhx = hy/hx;
784: tot_its = 0;
785: SNESNGSGetSweeps ( snes,&sweeps);
786: SNESNGSGetTolerances ( snes,&atol,&rtol,&stol,&its);
787: DMGetLocalVector ( da,&localX);
788: if ( B) {
789: DMGetLocalVector ( da,&localB);
790: }
791: if ( B) {
792: DMGlobalToLocalBegin ( da,B,INSERT_VALUES ,localB);
793: DMGlobalToLocalEnd ( da,B,INSERT_VALUES ,localB);
794: }
795: if ( B) DMDAVecGetArrayRead ( da,localB,&b);
796: DMGlobalToLocalBegin ( da,X,INSERT_VALUES ,localX);
797: DMGlobalToLocalEnd ( da,X,INSERT_VALUES ,localX);
798: DMDAVecGetArray ( da,localX,&x);
799: for ( l= 0; l<sweeps; l++) {
800: /*
801: Get local grid boundaries ( for 2-dimensional DMDA ):
802: xs, ys - starting grid indices ( no ghost points)
803: xm, ym - widths of local grid ( no ghost points)
804: */
805: DMDAGetCorners ( da,&xs,&ys,NULL,&xm,&ym,NULL);
806: for ( m= 0; m<2; m++) {
807: for ( j= ys; j<ys+ym; j++) {
808: for ( i= xs+( m+j)%2; i<xs+xm; i+= 2) {
809: PetscReal xx = i*hx,yy = j*hy;
810: if ( B) bij = b[j][i];
811: else bij = 0.;
813: if ( i = = 0 || j = = 0 || i = = info.mx-1 || j = = info.my-1) {
814: /* boundary conditions are all zero Dirichlet */
815: x[j][i] = 0.0 + bij;
816: } else {
817: const PetscScalar
818: u_E = x[j][i+1],
819: u_W = x[j][i-1],
820: u_N = x[j+1][i],
821: u_S = x[j-1][i];
822: const PetscScalar
823: uy_E = 0.25*dhy*( x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]),
824: uy_W = 0.25*dhy*( x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]),
825: ux_N = 0.25*dhx*( x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]),
826: ux_S = 0.25*dhx*( x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]);
827: u = x[j][i];
828: for ( k= 0; k<its; k++) {
829: const PetscScalar
830: ux_E = dhx*( u_E-u),
831: ux_W = dhx*( u-u_W),
832: uy_N = dhy*( u_N-u),
833: uy_S = dhy*( u-u_S),
834: e_E = eta( user,xx,yy,ux_E,uy_E),
835: e_W = eta( user,xx,yy,ux_W,uy_W),
836: e_N = eta( user,xx,yy,ux_N,uy_N),
837: e_S = eta( user,xx,yy,ux_S,uy_S),
838: de_E = deta( user,xx,yy,ux_E,uy_E),
839: de_W = deta( user,xx,yy,ux_W,uy_W),
840: de_N = deta( user,xx,yy,ux_N,uy_N),
841: de_S = deta( user,xx,yy,ux_S,uy_S),
842: newt_E = e_E+de_E*PetscSqr ( ux_E),
843: newt_W = e_W+de_W*PetscSqr ( ux_W),
844: newt_N = e_N+de_N*PetscSqr ( uy_N),
845: newt_S = e_S+de_S*PetscSqr ( uy_S),
846: uxx = -hy * ( e_E*ux_E - e_W*ux_W),
847: uyy = -hx * ( e_N*uy_N - e_S*uy_S);
849: if ( sc) eu = PetscExpScalar( u);
850: else eu = 0;
852: F = uxx + uyy - sc*eu - bij;
853: if ( k = = 0) F0 = F;
854: J = hxdhy*( newt_N + newt_S) + hydhx*( newt_E + newt_W) - sc*eu;
855: y = F/J;
856: u -= y;
857: tot_its++;
858: if ( atol > PetscAbsReal ( PetscRealPart ( F)) ||
859: rtol*PetscAbsReal ( PetscRealPart ( F0)) > PetscAbsReal ( PetscRealPart ( F)) ||
860: stol*PetscAbsReal ( PetscRealPart ( u)) > PetscAbsReal ( PetscRealPart ( y))) {
861: break ;
862: }
863: }
864: x[j][i] = u;
865: }
866: }
867: }
868: }
869: /*
870: x Restore vector
871: */
872: }
873: DMDAVecRestoreArray ( da,localX,&x);
874: DMLocalToGlobalBegin ( da,localX,INSERT_VALUES ,X);
875: DMLocalToGlobalEnd ( da,localX,INSERT_VALUES ,X);
876: PetscLogFlops ( tot_its*( 118.0));
877: DMRestoreLocalVector ( da,&localX);
878: if ( B) {
879: DMDAVecRestoreArrayRead ( da,localB,&b);
880: DMRestoreLocalVector ( da,&localB);
881: }
882: return ( 0);
883: }
886: /*TEST
888: test:
889: nsize: 2
890: args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype NEWTON
891: requires: !single
893: test:
894: suffix: 2
895: nsize: 2
896: args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -precheck 1
897: requires: !single
899: test:
900: suffix: 3
901: nsize: 2
902: args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -picard -precheck 1
903: requires: !single
905: test:
906: suffix: 4
907: 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
908: requires: !single
910: test:
911: suffix: lag_jac
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_jacobian 5 -npc_pc_type asm -npc_ksp_converged_reason -npc_snes_lag_jacobian_persists
914: requires: !single
916: test:
917: suffix: lag_pc
918: nsize: 4
919: 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
920: requires: !single
922: test:
923: suffix: nleqerr
924: 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
925: requires: !single
927: TEST*/