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