Actual source code: ex5.c
petsc-3.3-p7 2013-05-11
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";
9: /*T
10: Concepts: SNES^parallel Bratu example
11: Concepts: DMDA^using distributed arrays;
12: Concepts: IS coloirng types;
13: Processors: n
14: T*/
16: /* ------------------------------------------------------------------------
18: Solid Fuel Ignition (SFI) problem. This problem is modeled by
19: the partial differential equation
20:
21: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
22:
23: with boundary conditions
24:
25: u = 0 for x = 0, x = 1, y = 0, y = 1.
26:
27: A finite difference approximation with the usual 5-point stencil
28: is used to discretize the boundary value problem to obtain a nonlinear
29: system of equations.
31: Program usage: mpiexec -n <procs> ex5 [-help] [all PETSc options]
32: e.g.,
33:
34: This example shows how geometric multigrid can be run transparently with a nonlinear solver so long
35: as SNESSetDM() is provided. Example usage
37: ./ex5 -pc_type mg -ksp_monitor -snes_view -pc_mg_levels 3 -pc_mg_galerkin -da_grid_x 17 -da_grid_y 17
38: -mg_levels_ksp_monitor -snes_monitor -mg_levels_pc_type sor -pc_mg_type full
40: or to run with grid sequencing on the nonlinear problem (note that you do not need to provide the number of
41: multigrid levels, it will be determined automatically based on the number of refinements done)
43: ./ex5 -pc_type mg -ksp_monitor -snes_view -pc_mg_galerkin -snes_grid_sequence 3
44: -mg_levels_ksp_monitor -snes_monitor -mg_levels_pc_type sor -pc_mg_type full
46: ------------------------------------------------------------------------- */
48: /*
49: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
50: Include "petscsnes.h" so that we can use SNES solvers. Note that this
51: */
52: #include <petscdmda.h>
53: #include <petscsnes.h>
55: /*
56: User-defined application context - contains data needed by the
57: application-provided call-back routines, FormJacobianLocal() and
58: FormFunctionLocal().
59: */
60: typedef struct {
61: PassiveReal param; /* test problem parameter */
62: } AppCtx;
64: /*
65: User-defined routines
66: */
67: extern PetscErrorCode FormInitialGuess(DM,AppCtx*,Vec);
68: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
69: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,AppCtx*);
70: #if defined(PETSC_HAVE_MATLAB_ENGINE)
71: extern PetscErrorCode FormFunctionMatlab(SNES,Vec,Vec,void *);
72: #endif
73: extern PetscErrorCode NonlinearGS(SNES,Vec,Vec,void*);
77: int main(int argc,char **argv)
78: {
79: SNES snes; /* nonlinear solver */
80: Vec x; /* solution vector */
81: AppCtx user; /* user-defined work context */
82: PetscInt its; /* iterations for convergence */
83: PetscErrorCode ierr;
84: PetscReal bratu_lambda_max = 6.81;
85: PetscReal bratu_lambda_min = 0.;
86: PetscBool flg = PETSC_FALSE;
87: DM da;
88: PetscBool matlab_function = PETSC_FALSE;
90: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: Initialize program
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94: PetscInitialize(&argc,&argv,(char *)0,help);
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Initialize problem parameters
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: user.param = 6.0;
100: PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
101: 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);
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Create nonlinear solver context
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: SNESCreate(PETSC_COMM_WORLD,&snes);
107: SNESSetGS(snes, NonlinearGS, PETSC_NULL);
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Create distributed array (DMDA) to manage parallel grid and vectors
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&da);
113: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
114: DMSetApplicationContext(da,&user);
115: SNESSetDM(snes,da);
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Extract global vectors from DMDA; then duplicate for remaining
118: vectors that are the same types
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: DMCreateGlobalVector(da,&x);
122: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: Set local function evaluation routine
124: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: DMDASetLocalFunction(da,(DMDALocalFunction1)FormFunctionLocal);
126: PetscOptionsGetBool(PETSC_NULL,"-fd",&flg,PETSC_NULL);
127: if (!flg) {
128: DMDASetLocalJacobian(da,(DMDALocalFunction1)FormJacobianLocal);
129: }
131: /* Decide which FormFunction to use */
132: PetscOptionsGetBool(PETSC_NULL,"-matlab_function",&matlab_function,0);
134: #if defined(PETSC_HAVE_MATLAB_ENGINE)
135: Vec r;
136: if (matlab_function) {
137: VecDuplicate(x,&r);
138: SNESSetFunction(snes,r,FormFunctionMatlab,&user);
139: }
140: #endif
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Customize nonlinear solver; set runtime options
144: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145: SNESSetFromOptions(snes);
147: FormInitialGuess(da,&user,x);
149: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: Solve nonlinear system
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: SNESSolve(snes,PETSC_NULL,x);
153: SNESGetIterationNumber(snes,&its);
155: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: Free work space. All PETSc objects should be destroyed when they
160: are no longer needed.
161: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
162: #if defined(PETSC_HAVE_MATLAB_ENGINE)
163: if (r){VecDestroy(&r);}
164: #endif
165: VecDestroy(&x);
166: SNESDestroy(&snes);
167: DMDestroy(&da);
168: PetscFinalize();
169: return 0;
170: }
171: /* ------------------------------------------------------------------- */
174: /*
175: FormInitialGuess - Forms initial approximation.
177: Input Parameters:
178: user - user-defined application context
179: X - vector
181: Output Parameter:
182: X - vector
183: */
184: PetscErrorCode FormInitialGuess(DM da,AppCtx *user,Vec X)
185: {
186: PetscInt i,j,Mx,My,xs,ys,xm,ym;
188: PetscReal lambda,temp1,temp,hx,hy;
189: PetscScalar **x;
192: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
193: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
195: lambda = user->param;
196: hx = 1.0/(PetscReal)(Mx-1);
197: hy = 1.0/(PetscReal)(My-1);
198: temp1 = lambda/(lambda + 1.0);
200: /*
201: Get a pointer to vector data.
202: - For default PETSc vectors, VecGetArray() returns a pointer to
203: the data array. Otherwise, the routine is implementation dependent.
204: - You MUST call VecRestoreArray() when you no longer need access to
205: the array.
206: */
207: DMDAVecGetArray(da,X,&x);
209: /*
210: Get local grid boundaries (for 2-dimensional DMDA):
211: xs, ys - starting grid indices (no ghost points)
212: xm, ym - widths of local grid (no ghost points)
214: */
215: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
217: /*
218: Compute initial guess over the locally owned part of the grid
219: */
220: for (j=ys; j<ys+ym; j++) {
221: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
222: for (i=xs; i<xs+xm; i++) {
223: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
224: /* boundary conditions are all zero Dirichlet */
225: x[j][i] = 0.0;
226: } else {
227: x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
228: }
229: }
230: }
232: /*
233: Restore vector
234: */
235: DMDAVecRestoreArray(da,X,&x);
237: return(0);
238: }
239: /* ------------------------------------------------------------------- */
242: /*
243: FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch
246: */
247: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
248: {
250: PetscInt i,j;
251: PetscReal lambda,hx,hy,hxdhy,hydhx,sc;
252: PetscScalar u,uxx,uyy;
256: lambda = user->param;
257: hx = 1.0/(PetscReal)(info->mx-1);
258: hy = 1.0/(PetscReal)(info->my-1);
259: sc = hx*hy*lambda;
260: hxdhy = hx/hy;
261: hydhx = hy/hx;
262: /*
263: Compute function over the locally owned part of the grid
264: */
265: for (j=info->ys; j<info->ys+info->ym; j++) {
266: for (i=info->xs; i<info->xs+info->xm; i++) {
267: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
268: f[j][i] = 2.0*(hydhx+hxdhy)*x[j][i];
269: } else {
270: u = x[j][i];
271: uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
272: uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
273: f[j][i] = uxx + uyy - sc*PetscExpScalar(u);
274: }
275: }
276: }
277: PetscLogFlops(11.0*info->ym*info->xm);
278: return(0);
279: }
283: /*
284: FormJacobianLocal - Evaluates Jacobian matrix on local process patch
285: */
286: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,AppCtx *user)
287: {
289: PetscInt i,j;
290: MatStencil col[5],row;
291: PetscScalar lambda,v[5],hx,hy,hxdhy,hydhx,sc;
294: lambda = user->param;
295: hx = 1.0/(PetscReal)(info->mx-1);
296: hy = 1.0/(PetscReal)(info->my-1);
297: sc = hx*hy*lambda;
298: hxdhy = hx/hy;
299: hydhx = hy/hx;
302: /*
303: Compute entries for the locally owned part of the Jacobian.
304: - Currently, all PETSc parallel matrix formats are partitioned by
305: contiguous chunks of rows across the processors.
306: - Each processor needs to insert only elements that it owns
307: locally (but any non-local elements will be sent to the
308: appropriate processor during matrix assembly).
309: - Here, we set all entries for a particular row at once.
310: - We can set matrix entries either using either
311: MatSetValuesLocal() or MatSetValues(), as discussed above.
312: */
313: for (j=info->ys; j<info->ys+info->ym; j++) {
314: for (i=info->xs; i<info->xs+info->xm; i++) {
315: row.j = j; row.i = i;
316: /* boundary points */
317: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
318: v[0] = 2.0*(hydhx + hxdhy);
319: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
320: } else {
321: /* interior grid points */
322: v[0] = -hxdhy; col[0].j = j - 1; col[0].i = i;
323: v[1] = -hydhx; col[1].j = j; col[1].i = i-1;
324: v[2] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(x[j][i]); col[2].j = row.j; col[2].i = row.i;
325: v[3] = -hydhx; col[3].j = j; col[3].i = i+1;
326: v[4] = -hxdhy; col[4].j = j + 1; col[4].i = i;
327: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
328: }
329: }
330: }
332: /*
333: Assemble matrix, using the 2-step process:
334: MatAssemblyBegin(), MatAssemblyEnd().
335: */
336: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
337: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
338: /*
339: Tell the matrix we will never add a new nonzero location to the
340: matrix. If we do, it will generate an error.
341: */
342: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
343: return(0);
344: }
346: #if defined(PETSC_HAVE_MATLAB_ENGINE)
349: PetscErrorCode FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
350: {
351: AppCtx *user = (AppCtx*)ptr;
353: PetscInt Mx,My;
354: PetscReal lambda,hx,hy;
355: Vec localX,localF;
356: MPI_Comm comm;
357: DM da;
360: SNESGetDM(snes,&da);
361: DMGetLocalVector(da,&localX);
362: DMGetLocalVector(da,&localF);
363: PetscObjectSetName((PetscObject)localX,"localX");
364: PetscObjectSetName((PetscObject)localF,"localF");
365: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
366: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
368: lambda = user->param;
369: hx = 1.0/(PetscReal)(Mx-1);
370: hy = 1.0/(PetscReal)(My-1);
372: PetscObjectGetComm((PetscObject)snes,&comm);
373: /*
374: Scatter ghost points to local vector,using the 2-step process
375: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
376: By placing code between these two statements, computations can be
377: done while messages are in transition.
378: */
379: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
380: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
381: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
382: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
383: PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);
385: /*
386: Insert values into global vector
387: */
388: DMLocalToGlobalBegin(da,localF,INSERT_VALUES,F);
389: DMLocalToGlobalEnd(da,localF,INSERT_VALUES,F);
390: DMRestoreLocalVector(da,&localX);
391: DMRestoreLocalVector(da,&localF);
392: return(0);
393: }
394: #endif
396: /* ------------------------------------------------------------------- */
399: /*
400: Applies some sweeps on nonlinear Gauss-Seidel on each process
402: */
403: PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void * ctx)
404: {
405: PetscInt i,j,Mx,My,xs,ys,xm,ym,k,its,l;
407: PetscReal lambda,hx,hy,hxdhy,hydhx,sc;
408: PetscScalar **x,**b,bij,F,J,u,uxx,uyy;
409: DM da;
410: AppCtx *user;
411: Vec localX,localB;
414: SNESGetTolerances(snes,PETSC_NULL,PETSC_NULL,PETSC_NULL,&its,PETSC_NULL);
415: SNESGetDM(snes,&da);
416: DMGetApplicationContext(da,(void**)&user);
418: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
419: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
421: lambda = user->param;
422: hx = 1.0/(PetscReal)(Mx-1);
423: hy = 1.0/(PetscReal)(My-1);
424: sc = hx*hy*lambda;
425: hxdhy = hx/hy;
426: hydhx = hy/hx;
429: DMGetLocalVector(da,&localX);
430: if (B) {
431: DMGetLocalVector(da,&localB);
432: }
433: for (l=0; l<1; l++) {
435: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
436: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
437: if (B) {
438: DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);
439: DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);
440: }
441: /*
442: Get a pointer to vector data.
443: - For default PETSc vectors, VecGetArray() returns a pointer to
444: the data array. Otherwise, the routine is implementation dependent.
445: - You MUST call VecRestoreArray() when you no longer need access to
446: the array.
447: */
448: DMDAVecGetArray(da,localX,&x);
449: if (B) DMDAVecGetArray(da,localB,&b);
450: /*
451: Get local grid boundaries (for 2-dimensional DMDA):
452: xs, ys - starting grid indices (no ghost points)
453: xm, ym - widths of local grid (no ghost points)
454: */
455: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
457: for (j=ys; j<ys+ym; j++) {
458: for (i=xs; i<xs+xm; i++) {
459: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
460: /* boundary conditions are all zero Dirichlet */
461: x[j][i] = 0.0;
462: } else {
463: if (B) {
464: bij = b[j][i];
465: } else {
466: bij = 0.;
467: }
468: u = x[j][i];
469: for (k=0; k<1; k++) {
470: uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
471: uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
472: F = uxx + uyy - sc*PetscExpScalar(u) - bij;
473: J = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(u);
474: u = u - F/J;
475: }
476: x[j][i] = u;
477: }
478: }
479: }
480: /*
481: Restore vector
482: */
483: DMDAVecRestoreArray(da,localX,&x);
484: DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
485: DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
486: }
487: PetscLogFlops((11.0 + 5)*ym*xm);
488: DMRestoreLocalVector(da,&localX);
489: if (B) {
490: DMDAVecRestoreArray(da,localB,&b);
491: DMRestoreLocalVector(da,&localB);
492: }
493: return(0);
494: }