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