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