Actual source code: ex5.c
petsc-3.10.5 2019-03-28
2: static char help[] = "Bratu nonlinear PDE in 2d.\n\
3: We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular\n\
4: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\
5: The command line options include:\n\
6: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
7: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\n\
8: -m_par/n_par <parameter>, where <parameter> indicates an integer\n \
9: that MMS3 will be evaluated with 2^m_par, 2^n_par";
11: /*T
12: Concepts: SNES^parallel Bratu example
13: Concepts: DMDA^using distributed arrays;
14: Concepts: IS coloirng types;
15: Processors: n
16: T*/
20: /* ------------------------------------------------------------------------
22: Solid Fuel Ignition (SFI) problem. This problem is modeled by
23: the partial differential equation
25: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
27: with boundary conditions
29: u = 0 for x = 0, x = 1, y = 0, y = 1.
31: A finite difference approximation with the usual 5-point stencil
32: is used to discretize the boundary value problem to obtain a nonlinear
33: system of equations.
36: This example shows how geometric multigrid can be run transparently with a nonlinear solver so long
37: as SNESSetDM() is provided. Example usage
39: ./ex5 -pc_type mg -ksp_monitor -snes_view -pc_mg_levels 3 -pc_mg_galerkin pmat -da_grid_x 17 -da_grid_y 17
40: -mg_levels_ksp_monitor -snes_monitor -mg_levels_pc_type sor -pc_mg_type full
42: or to run with grid sequencing on the nonlinear problem (note that you do not need to provide the number of
43: multigrid levels, it will be determined automatically based on the number of refinements done)
45: ./ex5 -pc_type mg -ksp_monitor -snes_view -pc_mg_galerkin pmat -snes_grid_sequence 3
46: -mg_levels_ksp_monitor -snes_monitor -mg_levels_pc_type sor -pc_mg_type full
48: ------------------------------------------------------------------------- */
50: /*
51: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
52: Include "petscsnes.h" so that we can use SNES solvers. Note that this
53: */
54: #include <petscdm.h>
55: #include <petscdmda.h>
56: #include <petscsnes.h>
57: #include <petscmatlab.h>
59: /*
60: User-defined application context - contains data needed by the
61: application-provided call-back routines, FormJacobianLocal() and
62: FormFunctionLocal().
63: */
64: typedef struct AppCtx AppCtx;
65: struct AppCtx {
66: PetscReal param; /* test problem parameter */
67: PetscInt m,n; /* MMS3 parameters */
68: PetscErrorCode (*mms_solution)(AppCtx*,const DMDACoor2d*,PetscScalar*);
69: PetscErrorCode (*mms_forcing)(AppCtx*,const DMDACoor2d*,PetscScalar*);
70: };
72: /*
73: User-defined routines
74: */
75: extern PetscErrorCode FormInitialGuess(DM,AppCtx*,Vec);
76: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
77: extern PetscErrorCode FormExactSolution(DM,AppCtx*,Vec);
78: extern PetscErrorCode ZeroBCSolution(AppCtx*,const DMDACoor2d*,PetscScalar*);
79: extern PetscErrorCode MMSSolution1(AppCtx*,const DMDACoor2d*,PetscScalar*);
80: extern PetscErrorCode MMSForcing1(AppCtx*,const DMDACoor2d*,PetscScalar*);
81: extern PetscErrorCode MMSSolution2(AppCtx*,const DMDACoor2d*,PetscScalar*);
82: extern PetscErrorCode MMSForcing2(AppCtx*,const DMDACoor2d*,PetscScalar*);
83: extern PetscErrorCode MMSSolution3(AppCtx*,const DMDACoor2d*,PetscScalar*);
84: extern PetscErrorCode MMSForcing3(AppCtx*,const DMDACoor2d*,PetscScalar*);
85: extern PetscErrorCode MMSSolution4(AppCtx*,const DMDACoor2d*,PetscScalar*);
86: extern PetscErrorCode MMSForcing4(AppCtx*,const DMDACoor2d*,PetscScalar*);
87: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,Mat,AppCtx*);
88: extern PetscErrorCode FormObjectiveLocal(DMDALocalInfo*,PetscScalar**,PetscReal*,AppCtx*);
89: #if defined(PETSC_HAVE_MATLAB_ENGINE)
90: extern PetscErrorCode FormFunctionMatlab(SNES,Vec,Vec,void*);
91: #endif
92: extern PetscErrorCode NonlinearGS(SNES,Vec,Vec,void*);
94: int main(int argc,char **argv)
95: {
96: SNES snes; /* nonlinear solver */
97: Vec x; /* solution vector */
98: AppCtx user; /* user-defined work context */
99: PetscInt its; /* iterations for convergence */
101: PetscReal bratu_lambda_max = 6.81;
102: PetscReal bratu_lambda_min = 0.;
103: PetscInt MMS = 0;
104: PetscBool flg = PETSC_FALSE;
105: DM da;
106: #if defined(PETSC_HAVE_MATLAB_ENGINE)
107: Vec r = NULL;
108: PetscBool matlab_function = PETSC_FALSE;
109: #endif
110: KSP ksp;
111: PetscInt lits,slits;
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Initialize program
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
119: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: Initialize problem parameters
121: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122: user.param = 6.0;
123: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
124: if (user.param > bratu_lambda_max || user.param < bratu_lambda_min) SETERRQ3(PETSC_COMM_SELF,1,"Lambda, %g, is out of range, [%g, %g]", user.param, bratu_lambda_min, bratu_lambda_max);
125: PetscOptionsGetInt(NULL,NULL,"-mms",&MMS,NULL);
126: if (MMS == 3) {
127: PetscInt mPar = 2, nPar = 1;
128: PetscOptionsGetInt(NULL,NULL,"-m_par",&mPar,NULL);
129: PetscOptionsGetInt(NULL,NULL,"-n_par",&nPar,NULL);
130: user.m = PetscPowInt(2,mPar);
131: user.n = PetscPowInt(2,nPar);
132: }
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Create nonlinear solver context
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137: SNESCreate(PETSC_COMM_WORLD,&snes);
138: SNESSetCountersReset(snes,PETSC_FALSE);
139: SNESSetNGS(snes, NonlinearGS, NULL);
141: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: Create distributed array (DMDA) to manage parallel grid and vectors
143: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
145: DMSetFromOptions(da);
146: DMSetUp(da);
147: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
148: DMSetApplicationContext(da,&user);
149: SNESSetDM(snes,da);
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: Extract global vectors from DMDA; then duplicate for remaining
152: vectors that are the same types
153: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154: DMCreateGlobalVector(da,&x);
156: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: Set local function evaluation routine
158: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
159: user.mms_solution = ZeroBCSolution;
160: switch (MMS) {
161: case 0: user.mms_solution = NULL; user.mms_forcing = NULL;
162: case 1: user.mms_solution = MMSSolution1; user.mms_forcing = MMSForcing1; break;
163: case 2: user.mms_solution = MMSSolution2; user.mms_forcing = MMSForcing2; break;
164: case 3: user.mms_solution = MMSSolution3; user.mms_forcing = MMSForcing3; break;
165: case 4: user.mms_solution = MMSSolution4; user.mms_forcing = MMSForcing4; break;
166: default: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_USER,"Unknown MMS type %d",MMS);
167: }
168: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocal,&user);
169: PetscOptionsGetBool(NULL,NULL,"-fd",&flg,NULL);
170: if (!flg) {
171: DMDASNESSetJacobianLocal(da,(DMDASNESJacobian)FormJacobianLocal,&user);
172: }
174: PetscOptionsGetBool(NULL,NULL,"-obj",&flg,NULL);
175: if (flg) {
176: DMDASNESSetObjectiveLocal(da,(DMDASNESObjective)FormObjectiveLocal,&user);
177: }
179: #if defined(PETSC_HAVE_MATLAB_ENGINE)
180: PetscOptionsGetBool(NULL,NULL,"-matlab_function",&matlab_function,0);
181: if (matlab_function) {
182: VecDuplicate(x,&r);
183: SNESSetFunction(snes,r,FormFunctionMatlab,&user);
184: }
185: #endif
187: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188: Customize nonlinear solver; set runtime options
189: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
190: SNESSetFromOptions(snes);
192: FormInitialGuess(da,&user,x);
194: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195: Solve nonlinear system
196: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
197: SNESSolve(snes,NULL,x);
198: SNESGetIterationNumber(snes,&its);
200: SNESGetLinearSolveIterations(snes,&slits);
201: SNESGetKSP(snes,&ksp);
202: KSPGetTotalIterations(ksp,&lits);
203: if (lits != slits) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Number of total linear iterations reported by SNES %D does not match reported by KSP %D",slits,lits);
204: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
207: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
208: If using MMS, check the l_2 error
209: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
210: if (MMS) {
211: Vec e;
212: PetscReal errorl2, errorinf;
213: PetscInt N;
215: VecDuplicate(x, &e);
216: PetscObjectViewFromOptions((PetscObject) x, NULL, "-sol_view");
217: FormExactSolution(da, &user, e);
218: PetscObjectViewFromOptions((PetscObject) e, NULL, "-exact_view");
219: VecAXPY(e, -1.0, x);
220: PetscObjectViewFromOptions((PetscObject) e, NULL, "-error_view");
221: VecNorm(e, NORM_2, &errorl2);
222: VecNorm(e, NORM_INFINITY, &errorinf);
223: VecGetSize(e, &N);
224: PetscPrintf(PETSC_COMM_WORLD, "N: %D error L2 %g inf %g\n", N, (double) errorl2/PetscSqrtReal(N), (double) errorinf);
225: VecDestroy(&e);
226: }
228: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
229: Free work space. All PETSc objects should be destroyed when they
230: are no longer needed.
231: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
232: #if defined(PETSC_HAVE_MATLAB_ENGINE)
233: VecDestroy(&r);
234: #endif
235: VecDestroy(&x);
236: SNESDestroy(&snes);
237: DMDestroy(&da);
238: PetscFinalize();
239: return ierr;
240: }
241: /* ------------------------------------------------------------------- */
242: /*
243: FormInitialGuess - Forms initial approximation.
245: Input Parameters:
246: da - The DM
247: user - user-defined application context
249: Output Parameter:
250: X - vector
251: */
252: PetscErrorCode FormInitialGuess(DM da,AppCtx *user,Vec X)
253: {
254: PetscInt i,j,Mx,My,xs,ys,xm,ym;
256: PetscReal lambda,temp1,temp,hx,hy;
257: PetscScalar **x;
260: 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);
262: lambda = user->param;
263: hx = 1.0/(PetscReal)(Mx-1);
264: hy = 1.0/(PetscReal)(My-1);
265: temp1 = lambda/(lambda + 1.0);
267: /*
268: Get a pointer to vector data.
269: - For default PETSc vectors, VecGetArray() returns a pointer to
270: the data array. Otherwise, the routine is implementation dependent.
271: - You MUST call VecRestoreArray() when you no longer need access to
272: the array.
273: */
274: DMDAVecGetArray(da,X,&x);
276: /*
277: Get local grid boundaries (for 2-dimensional DMDA):
278: xs, ys - starting grid indices (no ghost points)
279: xm, ym - widths of local grid (no ghost points)
281: */
282: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
284: /*
285: Compute initial guess over the locally owned part of the grid
286: */
287: for (j=ys; j<ys+ym; j++) {
288: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
289: for (i=xs; i<xs+xm; i++) {
290: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
291: /* boundary conditions are all zero Dirichlet */
292: x[j][i] = 0.0;
293: } else {
294: x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
295: }
296: }
297: }
299: /*
300: Restore vector
301: */
302: DMDAVecRestoreArray(da,X,&x);
303: return(0);
304: }
306: /*
307: FormExactSolution - Forms MMS solution
309: Input Parameters:
310: da - The DM
311: user - user-defined application context
313: Output Parameter:
314: X - vector
315: */
316: PetscErrorCode FormExactSolution(DM da, AppCtx *user, Vec U)
317: {
318: DM coordDA;
319: Vec coordinates;
320: DMDACoor2d **coords;
321: PetscScalar **u;
322: PetscInt xs, ys, xm, ym, i, j;
326: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
327: DMGetCoordinateDM(da, &coordDA);
328: DMGetCoordinates(da, &coordinates);
329: DMDAVecGetArray(coordDA, coordinates, &coords);
330: DMDAVecGetArray(da, U, &u);
331: for (j = ys; j < ys+ym; ++j) {
332: for (i = xs; i < xs+xm; ++i) {
333: user->mms_solution(user,&coords[j][i],&u[j][i]);
334: }
335: }
336: DMDAVecRestoreArray(da, U, &u);
337: DMDAVecRestoreArray(coordDA, coordinates, &coords);
338: return(0);
339: }
341: PetscErrorCode ZeroBCSolution(AppCtx *user,const DMDACoor2d *c,PetscScalar *u)
342: {
343: u[0] = 0.;
344: return 0;
345: }
347: /* The functions below evaluate the MMS solution u(x,y) and associated forcing
349: f(x,y) = -u_xx - u_yy - lambda exp(u)
351: such that u(x,y) is an exact solution with f(x,y) as the right hand side forcing term.
352: */
353: PetscErrorCode MMSSolution1(AppCtx *user,const DMDACoor2d *c,PetscScalar *u)
354: {
355: PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
356: u[0] = x*(1 - x)*y*(1 - y);
357: PetscLogFlops(5);
358: return 0;
359: }
360: PetscErrorCode MMSForcing1(AppCtx *user,const DMDACoor2d *c,PetscScalar *f)
361: {
362: PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
363: f[0] = 2*x*(1 - x) + 2*y*(1 - y) - user->param*PetscExpReal(x*(1 - x)*y*(1 - y));
364: return 0;
365: }
367: PetscErrorCode MMSSolution2(AppCtx *user,const DMDACoor2d *c,PetscScalar *u)
368: {
369: PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
370: u[0] = PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y);
371: PetscLogFlops(5);
372: return 0;
373: }
374: PetscErrorCode MMSForcing2(AppCtx *user,const DMDACoor2d *c,PetscScalar *f)
375: {
376: PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
377: f[0] = 2*PetscSqr(PETSC_PI)*PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y) - user->param*PetscExpReal(PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y));
378: return 0;
379: }
381: PetscErrorCode MMSSolution3(AppCtx *user,const DMDACoor2d *c,PetscScalar *u)
382: {
383: PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
384: u[0] = PetscSinReal(user->m*PETSC_PI*x*(1-y))*PetscSinReal(user->n*PETSC_PI*y*(1-x));
385: PetscLogFlops(5);
386: return 0;
387: }
388: PetscErrorCode MMSForcing3(AppCtx *user,const DMDACoor2d *c,PetscScalar *f)
389: {
390: PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
391: PetscReal m = user->m, n = user->n, lambda = user->param;
392: f[0] = (-(PetscExpReal(PetscSinReal(m*PETSC_PI*x*(1 - y))*PetscSinReal(n*PETSC_PI*(1 - x)*y))*lambda)
393: + PetscSqr(PETSC_PI)*(-2*m*n*((-1 + x)*x + (-1 + y)*y)*PetscCosReal(m*PETSC_PI*x*(-1 + y))*PetscCosReal(n*PETSC_PI*(-1 + x)*y)
394: + (PetscSqr(m)*(PetscSqr(x) + PetscSqr(-1 + y)) + PetscSqr(n)*(PetscSqr(-1 + x) + PetscSqr(y)))
395: *PetscSinReal(m*PETSC_PI*x*(-1 + y))*PetscSinReal(n*PETSC_PI*(-1 + x)*y)));
396: return 0;
397: }
399: PetscErrorCode MMSSolution4(AppCtx *user,const DMDACoor2d *c,PetscScalar *u)
400: {
401: const PetscReal Lx = 1.,Ly = 1.;
402: PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
403: u[0] = (PetscPowReal(x,4)-PetscSqr(Lx)*PetscSqr(x))*(PetscPowReal(y,4)-PetscSqr(Ly)*PetscSqr(y));
404: PetscLogFlops(9);
405: return 0;
406: }
407: PetscErrorCode MMSForcing4(AppCtx *user,const DMDACoor2d *c,PetscScalar *f)
408: {
409: const PetscReal Lx = 1.,Ly = 1.;
410: PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
411: f[0] = (2*PetscSqr(x)*(PetscSqr(x)-PetscSqr(Lx))*(PetscSqr(Ly)-6*PetscSqr(y))
412: + 2*PetscSqr(y)*(PetscSqr(Lx)-6*PetscSqr(x))*(PetscSqr(y)-PetscSqr(Ly))
413: - user->param*PetscExpReal((PetscPowReal(x,4)-PetscSqr(Lx)*PetscSqr(x))*(PetscPowReal(y,4)-PetscSqr(Ly)*PetscSqr(y))));
414: return 0;
415: }
417: /* ------------------------------------------------------------------- */
418: /*
419: FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch
422: */
423: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
424: {
426: PetscInt i,j;
427: PetscReal lambda,hx,hy,hxdhy,hydhx;
428: PetscScalar u,ue,uw,un,us,uxx,uyy,mms_solution,mms_forcing;
429: DMDACoor2d c;
432: lambda = user->param;
433: hx = 1.0/(PetscReal)(info->mx-1);
434: hy = 1.0/(PetscReal)(info->my-1);
435: hxdhy = hx/hy;
436: hydhx = hy/hx;
437: /*
438: Compute function over the locally owned part of the grid
439: */
440: for (j=info->ys; j<info->ys+info->ym; j++) {
441: for (i=info->xs; i<info->xs+info->xm; i++) {
442: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
443: c.x = i*hx; c.y = j*hy;
444: user->mms_solution(user,&c,&mms_solution);
445: f[j][i] = 2.0*(hydhx+hxdhy)*(x[j][i] - mms_solution);
446: } else {
447: u = x[j][i];
448: uw = x[j][i-1];
449: ue = x[j][i+1];
450: un = x[j-1][i];
451: us = x[j+1][i];
453: /* Enforce boundary conditions at neighboring points -- setting these values causes the Jacobian to be symmetric. */
454: if (i-1 == 0) {c.x = (i-1)*hx; c.y = j*hy; user->mms_solution(user,&c,&uw);}
455: if (i+1 == info->mx-1) {c.x = (i+1)*hx; c.y = j*hy; user->mms_solution(user,&c,&ue);}
456: if (j-1 == 0) {c.x = i*hx; c.y = (j-1)*hy; user->mms_solution(user,&c,&un);}
457: if (j+1 == info->my-1) {c.x = i*hx; c.y = (j+1)*hy; user->mms_solution(user,&c,&us);}
459: uxx = (2.0*u - uw - ue)*hydhx;
460: uyy = (2.0*u - un - us)*hxdhy;
461: mms_forcing = 0;
462: c.x = i*hx; c.y = j*hy;
463: if (user->mms_forcing) {user->mms_forcing(user,&c,&mms_forcing);}
464: f[j][i] = uxx + uyy - hx*hy*(lambda*PetscExpScalar(u) + mms_forcing);
465: }
466: }
467: }
468: PetscLogFlops(11.0*info->ym*info->xm);
469: return(0);
470: }
472: /* FormObjectiveLocal - Evaluates nonlinear function, F(x) on local process patch */
473: PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info,PetscScalar **x,PetscReal *obj,AppCtx *user)
474: {
476: PetscInt i,j;
477: PetscReal lambda,hx,hy,hxdhy,hydhx,sc,lobj=0;
478: PetscScalar u,ue,uw,un,us,uxux,uyuy;
479: MPI_Comm comm;
482: *obj = 0;
483: PetscObjectGetComm((PetscObject)info->da,&comm);
484: lambda = user->param;
485: hx = 1.0/(PetscReal)(info->mx-1);
486: hy = 1.0/(PetscReal)(info->my-1);
487: sc = hx*hy*lambda;
488: hxdhy = hx/hy;
489: hydhx = hy/hx;
490: /*
491: Compute function over the locally owned part of the grid
492: */
493: for (j=info->ys; j<info->ys+info->ym; j++) {
494: for (i=info->xs; i<info->xs+info->xm; i++) {
495: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
496: lobj += PetscRealPart((hydhx + hxdhy)*x[j][i]*x[j][i]);
497: } else {
498: u = x[j][i];
499: uw = x[j][i-1];
500: ue = x[j][i+1];
501: un = x[j-1][i];
502: us = x[j+1][i];
504: if (i-1 == 0) uw = 0.;
505: if (i+1 == info->mx-1) ue = 0.;
506: if (j-1 == 0) un = 0.;
507: if (j+1 == info->my-1) us = 0.;
509: /* F[u] = 1/2\int_{\omega}\nabla^2u(x)*u(x)*dx */
511: uxux = u*(2.*u - ue - uw)*hydhx;
512: uyuy = u*(2.*u - un - us)*hxdhy;
514: lobj += PetscRealPart(0.5*(uxux + uyuy) - sc*PetscExpScalar(u));
515: }
516: }
517: }
518: PetscLogFlops(12.0*info->ym*info->xm);
519: MPI_Allreduce(&lobj,obj,1,MPIU_REAL,MPIU_SUM,comm);
520: return(0);
521: }
523: /*
524: FormJacobianLocal - Evaluates Jacobian matrix on local process patch
525: */
526: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,Mat jacpre,AppCtx *user)
527: {
529: PetscInt i,j,k;
530: MatStencil col[5],row;
531: PetscScalar lambda,v[5],hx,hy,hxdhy,hydhx,sc;
532: DM coordDA;
533: Vec coordinates;
534: DMDACoor2d **coords;
537: lambda = user->param;
538: /* Extract coordinates */
539: DMGetCoordinateDM(info->da, &coordDA);
540: DMGetCoordinates(info->da, &coordinates);
541: DMDAVecGetArray(coordDA, coordinates, &coords);
542: hx = info->xm > 1 ? PetscRealPart(coords[info->ys][info->xs+1].x) - PetscRealPart(coords[info->ys][info->xs].x) : 1.0;
543: hy = info->ym > 1 ? PetscRealPart(coords[info->ys+1][info->xs].y) - PetscRealPart(coords[info->ys][info->xs].y) : 1.0;
544: DMDAVecRestoreArray(coordDA, coordinates, &coords);
545: hxdhy = hx/hy;
546: hydhx = hy/hx;
547: sc = hx*hy*lambda;
550: /*
551: Compute entries for the locally owned part of the Jacobian.
552: - Currently, all PETSc parallel matrix formats are partitioned by
553: contiguous chunks of rows across the processors.
554: - Each processor needs to insert only elements that it owns
555: locally (but any non-local elements will be sent to the
556: appropriate processor during matrix assembly).
557: - Here, we set all entries for a particular row at once.
558: - We can set matrix entries either using either
559: MatSetValuesLocal() or MatSetValues(), as discussed above.
560: */
561: for (j=info->ys; j<info->ys+info->ym; j++) {
562: for (i=info->xs; i<info->xs+info->xm; i++) {
563: row.j = j; row.i = i;
564: /* boundary points */
565: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
566: v[0] = 2.0*(hydhx + hxdhy);
567: MatSetValuesStencil(jacpre,1,&row,1,&row,v,INSERT_VALUES);
568: } else {
569: k = 0;
570: /* interior grid points */
571: if (j-1 != 0) {
572: v[k] = -hxdhy;
573: col[k].j = j - 1; col[k].i = i;
574: k++;
575: }
576: if (i-1 != 0) {
577: v[k] = -hydhx;
578: col[k].j = j; col[k].i = i-1;
579: k++;
580: }
582: v[k] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(x[j][i]); col[k].j = row.j; col[k].i = row.i; k++;
584: if (i+1 != info->mx-1) {
585: v[k] = -hydhx;
586: col[k].j = j; col[k].i = i+1;
587: k++;
588: }
589: if (j+1 != info->mx-1) {
590: v[k] = -hxdhy;
591: col[k].j = j + 1; col[k].i = i;
592: k++;
593: }
594: MatSetValuesStencil(jacpre,1,&row,k,col,v,INSERT_VALUES);
595: }
596: }
597: }
599: /*
600: Assemble matrix, using the 2-step process:
601: MatAssemblyBegin(), MatAssemblyEnd().
602: */
603: MatAssemblyBegin(jacpre,MAT_FINAL_ASSEMBLY);
604: MatAssemblyEnd(jacpre,MAT_FINAL_ASSEMBLY);
605: /*
606: Tell the matrix we will never add a new nonzero location to the
607: matrix. If we do, it will generate an error.
608: */
609: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
610: return(0);
611: }
613: #if defined(PETSC_HAVE_MATLAB_ENGINE)
614: PetscErrorCode FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
615: {
616: AppCtx *user = (AppCtx*)ptr;
618: PetscInt Mx,My;
619: PetscReal lambda,hx,hy;
620: Vec localX,localF;
621: MPI_Comm comm;
622: DM da;
625: SNESGetDM(snes,&da);
626: DMGetLocalVector(da,&localX);
627: DMGetLocalVector(da,&localF);
628: PetscObjectSetName((PetscObject)localX,"localX");
629: PetscObjectSetName((PetscObject)localF,"localF");
630: 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);
632: lambda = user->param;
633: hx = 1.0/(PetscReal)(Mx-1);
634: hy = 1.0/(PetscReal)(My-1);
636: PetscObjectGetComm((PetscObject)snes,&comm);
637: /*
638: Scatter ghost points to local vector,using the 2-step process
639: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
640: By placing code between these two statements, computations can be
641: done while messages are in transition.
642: */
643: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
644: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
645: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
646: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
647: PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);
649: /*
650: Insert values into global vector
651: */
652: DMLocalToGlobalBegin(da,localF,INSERT_VALUES,F);
653: DMLocalToGlobalEnd(da,localF,INSERT_VALUES,F);
654: DMRestoreLocalVector(da,&localX);
655: DMRestoreLocalVector(da,&localF);
656: return(0);
657: }
658: #endif
660: /* ------------------------------------------------------------------- */
661: /*
662: Applies some sweeps on nonlinear Gauss-Seidel on each process
664: */
665: PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx)
666: {
667: PetscInt i,j,k,Mx,My,xs,ys,xm,ym,its,tot_its,sweeps,l;
669: PetscReal lambda,hx,hy,hxdhy,hydhx,sc;
670: PetscScalar **x,**b,bij,F,F0=0,J,u,un,us,ue,eu,uw,uxx,uyy,y;
671: PetscReal atol,rtol,stol;
672: DM da;
673: AppCtx *user;
674: Vec localX,localB;
677: tot_its = 0;
678: SNESNGSGetSweeps(snes,&sweeps);
679: SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);
680: SNESGetDM(snes,&da);
681: DMGetApplicationContext(da,(void**)&user);
683: 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);
685: lambda = user->param;
686: hx = 1.0/(PetscReal)(Mx-1);
687: hy = 1.0/(PetscReal)(My-1);
688: sc = hx*hy*lambda;
689: hxdhy = hx/hy;
690: hydhx = hy/hx;
693: DMGetLocalVector(da,&localX);
694: if (B) {
695: DMGetLocalVector(da,&localB);
696: }
697: for (l=0; l<sweeps; l++) {
699: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
700: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
701: if (B) {
702: DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);
703: DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);
704: }
705: /*
706: Get a pointer to vector data.
707: - For default PETSc vectors, VecGetArray() returns a pointer to
708: the data array. Otherwise, the routine is implementation dependent.
709: - You MUST call VecRestoreArray() when you no longer need access to
710: the array.
711: */
712: DMDAVecGetArray(da,localX,&x);
713: if (B) DMDAVecGetArray(da,localB,&b);
714: /*
715: Get local grid boundaries (for 2-dimensional DMDA):
716: xs, ys - starting grid indices (no ghost points)
717: xm, ym - widths of local grid (no ghost points)
718: */
719: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
721: for (j=ys; j<ys+ym; j++) {
722: for (i=xs; i<xs+xm; i++) {
723: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
724: /* boundary conditions are all zero Dirichlet */
725: x[j][i] = 0.0;
726: } else {
727: if (B) bij = b[j][i];
728: else bij = 0.;
730: u = x[j][i];
731: un = x[j-1][i];
732: us = x[j+1][i];
733: ue = x[j][i-1];
734: uw = x[j][i+1];
736: for (k=0; k<its; k++) {
737: eu = PetscExpScalar(u);
738: uxx = (2.0*u - ue - uw)*hydhx;
739: uyy = (2.0*u - un - us)*hxdhy;
740: F = uxx + uyy - sc*eu - bij;
741: if (k == 0) F0 = F;
742: J = 2.0*(hydhx + hxdhy) - sc*eu;
743: y = F/J;
744: u -= y;
745: tot_its++;
747: if (atol > PetscAbsReal(PetscRealPart(F)) ||
748: rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) ||
749: stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) {
750: break;
751: }
752: }
753: x[j][i] = u;
754: }
755: }
756: }
757: /*
758: Restore vector
759: */
760: DMDAVecRestoreArray(da,localX,&x);
761: DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
762: DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
763: }
764: PetscLogFlops(tot_its*(21.0));
765: DMRestoreLocalVector(da,&localX);
766: if (B) {
767: DMDAVecRestoreArray(da,localB,&b);
768: DMRestoreLocalVector(da,&localB);
769: }
770: return(0);
771: }
773: /*TEST
775: test:
776: suffix: asm_0
777: requires: !single
778: args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu
780: test:
781: suffix: msm_0
782: requires: !single
783: args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu
785: test:
786: suffix: asm_1
787: requires: !single
788: args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu -da_grid_x 8
790: test:
791: suffix: msm_1
792: requires: !single
793: args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu -da_grid_x 8
795: test:
796: suffix: asm_2
797: requires: !single
798: args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 3 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu -da_grid_x 8
800: test:
801: suffix: msm_2
802: requires: !single
803: args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 3 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu -da_grid_x 8
805: test:
806: suffix: asm_3
807: requires: !single
808: args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 4 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu -da_grid_x 8
810: test:
811: suffix: msm_3
812: requires: !single
813: args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 4 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu -da_grid_x 8
815: test:
816: suffix: asm_4
817: requires: !single
818: nsize: 2
819: args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu -da_grid_x 8
821: test:
822: suffix: msm_4
823: requires: !single
824: nsize: 2
825: args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu -da_grid_x 8
827: test:
828: suffix: asm_5
829: requires: !single
830: nsize: 2
831: args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 4 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu -da_grid_x 8
833: test:
834: suffix: msm_5
835: requires: !single
836: nsize: 2
837: args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 4 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu -da_grid_x 8
839: test:
840: args: -snes_rtol 1.e-5 -pc_type mg -ksp_monitor_short -snes_view -pc_mg_levels 3 -pc_mg_galerkin pmat -da_grid_x 17 -da_grid_y 17 -mg_levels_ksp_monitor_short -mg_levels_ksp_norm_type unpreconditioned -snes_monitor_short -mg_levels_ksp_chebyshev_esteig 0.5,1.1 -mg_levels_pc_type sor -pc_mg_type full
841: requires: !single
843: test:
844: suffix: 2
845: args: -pc_type mg -ksp_converged_reason -snes_view -pc_mg_galerkin pmat -snes_grid_sequence 3 -mg_levels_ksp_norm_type unpreconditioned -snes_monitor_short -mg_levels_ksp_chebyshev_esteig 0.5,1.1 -mg_levels_pc_type sor -pc_mg_type full -ksp_atol -1.
847: test:
848: suffix: 3
849: nsize: 2
850: args: -snes_grid_sequence 2 -snes_mf_operator -snes_converged_reason -snes_view -pc_type mg -snes_atol -1 -snes_rtol 1.e-2
851: filter: grep -v "otal number of function evaluations"
853: test:
854: suffix: 4
855: nsize: 2
856: args: -snes_grid_sequence 2 -snes_monitor_short -ksp_converged_reason -snes_converged_reason -snes_view -pc_type mg -snes_atol -1 -ksp_atol -1
858: test:
859: suffix: 5_anderson
860: args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type anderson
862: test:
863: suffix: 5_aspin
864: nsize: 4
865: args: -snes_monitor_short -ksp_monitor_short -snes_converged_reason -da_refine 4 -da_overlap 3 -snes_type aspin
867: test:
868: suffix: 5_broyden
869: args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type qn -snes_qn_type broyden -snes_qn_m 10
871: test:
872: suffix: 5_fas
873: args: -fas_coarse_snes_max_it 1 -fas_coarse_pc_type lu -fas_coarse_ksp_type preonly -snes_monitor_short -snes_type fas -fas_coarse_ksp_type richardson -da_refine 6
874: requires: !single
876: test:
877: suffix: 5_fas_additive
878: args: -fas_coarse_snes_max_it 1 -fas_coarse_pc_type lu -fas_coarse_ksp_type preonly -snes_monitor_short -snes_type fas -fas_coarse_ksp_type richardson -da_refine 6 -snes_fas_type additive -snes_max_it 50
880: test:
881: suffix: 5_fas_monitor
882: args: -da_refine 1 -snes_type fas -snes_fas_monitor
883: requires: !single
885: test:
886: suffix: 5_ls
887: args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type newtonls
889: test:
890: suffix: 5_ls_sell_sor
891: args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type newtonls -dm_mat_type sell -pc_type sor
892: output_file: output/ex5_5_ls.out
894: test:
895: suffix: 5_nasm
896: nsize: 4
897: args: -snes_monitor_short -snes_converged_reason -da_refine 4 -da_overlap 3 -snes_type nasm -snes_nasm_type restrict -snes_max_it 10
899: test:
900: suffix: 5_ncg
901: args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type ncg -snes_ncg_type fr
903: test:
904: suffix: 5_newton_asm_dmda
905: nsize: 4
906: args: -snes_monitor_short -ksp_monitor_short -snes_converged_reason -da_refine 4 -da_overlap 3 -snes_type newtonls -pc_type asm -pc_asm_dm_subdomains -malloc_dump
907: requires: !single
909: test:
910: suffix: 5_newton_gasm_dmda
911: nsize: 4
912: args: -snes_monitor_short -ksp_monitor_short -snes_converged_reason -da_refine 4 -da_overlap 3 -snes_type newtonls -pc_type gasm -malloc_dump
913: requires: !single
915: test:
916: suffix: 5_ngmres
917: args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type ngmres -snes_ngmres_m 10
919: test:
920: suffix: 5_ngmres_fas
921: args: -snes_rtol 1.e-4 -snes_type ngmres -npc_fas_coarse_snes_max_it 1 -npc_fas_coarse_snes_type newtonls -npc_fas_coarse_pc_type lu -npc_fas_coarse_ksp_type preonly -snes_ngmres_m 10 -snes_monitor_short -npc_snes_max_it 1 -npc_snes_type fas -npc_fas_coarse_ksp_type richardson -da_refine 6
923: test:
924: suffix: 5_ngmres_ngs
925: args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type ngmres -npc_snes_type ngs -npc_snes_max_it 1
927: test:
928: suffix: 5_ngmres_nrichardson
929: args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type ngmres -snes_ngmres_m 10 -npc_snes_type nrichardson -npc_snes_max_it 3
930: output_file: output/ex5_5_ngmres_richardson.out
932: test:
933: suffix: 5_nrichardson
934: args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type nrichardson
936: test:
937: suffix: 5_qn
938: args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type qn -snes_linesearch_type cp -snes_qn_m 10
940: test:
941: suffix: 6
942: nsize: 4
943: args: -snes_converged_reason -ksp_converged_reason -da_grid_x 129 -da_grid_y 129 -pc_type mg -pc_mg_levels 8 -mg_levels_ksp_type chebyshev -mg_levels_ksp_chebyshev_esteig 0,0.5,0,1.1 -mg_levels_ksp_max_it 2
945: TEST*/