Actual source code: ex5.c
petsc-3.4.5 2014-06-29
2: static char help[] = "Solves the modified Bratu problem in a 2D rectangular domain,\n\
3: using distributed arrays (DMDAs) to partition the parallel grid.\n\
4: The command line options include:\n\
5: -lambda <parameter>, where <parameter> indicates the problem's nonlinearity\n\
6: -kappa <parameter>, where <parameter> indicates the problem's nonlinearity\n\
7: -mx <xg>, where <xg> = number of grid points in the x-direction\n\
8: -my <yg>, where <yg> = number of grid points in the y-direction\n\
9: -Nx <npx>, where <npx> = number of processors in the x-direction\n\
10: -Ny <npy>, where <npy> = number of processors in the y-direction\n\n";
12: /*T
13: Concepts: SNES^solving a system of nonlinear equations (parallel Bratu example);
14: Concepts: DM^using distributed arrays;
15: Processors: n
16: T*/
18: /* ------------------------------------------------------------------------
20: Modified Solid Fuel Ignition problem. This problem is modeled by
21: the partial differential equation
23: -Laplacian u - kappa*\PartDer{u}{x} - lambda*exp(u) = 0,
25: where
27: 0 < x,y < 1,
29: with boundary conditions
31: u = 0 for x = 0, x = 1, y = 0, y = 1.
33: A finite difference approximation with the usual 5-point stencil
34: is used to discretize the boundary value problem to obtain a nonlinear
35: system of equations.
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: petscsys.h - base PETSc routines petscvec.h - vectors
44: 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 <petscdmda.h>
50: #include <petscsnes.h>
52: /*
53: User-defined application context - contains data needed by the
54: application-provided call-back routines, FormJacobian() and
55: FormFunction().
56: */
57: typedef struct {
58: PetscReal param; /* test problem parameter */
59: PetscReal param2; /* test problem parameter */
60: PetscInt mx,my; /* discretization in x, y directions */
61: Vec localX,localF; /* ghosted local vector */
62: DM da; /* distributed array data structure */
63: PetscMPIInt rank; /* processor rank */
64: } AppCtx;
66: /*
67: User-defined routines
68: */
69: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*),FormInitialGuess(AppCtx*,Vec);
70: extern PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
74: int main(int argc,char **argv)
75: {
76: SNES snes; /* nonlinear solver */
77: Vec x,r; /* solution, residual vectors */
78: Mat J; /* Jacobian matrix */
79: AppCtx user; /* user-defined work context */
80: PetscInt its; /* iterations for convergence */
81: PetscInt Nx,Ny; /* number of preocessors in x- and y- directions */
82: PetscBool matrix_free = PETSC_FALSE; /* flag - 1 indicates matrix-free version */
83: PetscMPIInt size; /* number of processors */
84: PetscInt m,N;
86: PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
87: PetscReal bratu_kappa_max = 10000,bratu_kappa_min = 0.;
89: PetscInitialize(&argc,&argv,(char*)0,help);
90: MPI_Comm_rank(PETSC_COMM_WORLD,&user.rank);
92: /*
93: Initialize problem parameters
94: */
95: user.mx = 4; user.my = 4; user.param = 6.0; user.param2 = 0.0;
96: PetscOptionsGetInt(NULL,"-mx",&user.mx,NULL);
97: PetscOptionsGetInt(NULL,"-my",&user.my,NULL);
98: PetscOptionsGetReal(NULL,"-lambda",&user.param,NULL);
99: if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,1,"Lambda is out of range");
100: PetscOptionsGetReal(NULL,"-kappa",&user.param2,NULL);
101: if (user.param2 >= bratu_kappa_max || user.param2 < bratu_kappa_min) SETERRQ(PETSC_COMM_SELF,1,"Kappa is out of range");
102: PetscPrintf(PETSC_COMM_WORLD,"Solving the Bratu problem with lambda=%G, kappa=%G\n",user.param,user.param2);
104: N = user.mx*user.my;
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Create nonlinear solver context
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: SNESCreate(PETSC_COMM_WORLD,&snes);
112: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: Create vector data structures; set function evaluation routine
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116: /*
117: Create distributed array (DMDA) to manage parallel grid and vectors
118: */
119: MPI_Comm_size(PETSC_COMM_WORLD,&size);
120: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
121: PetscOptionsGetInt(NULL,"-Nx",&Nx,NULL);
122: PetscOptionsGetInt(NULL,"-Ny",&Ny,NULL);
123: if (Nx*Ny != size && (Nx != PETSC_DECIDE || Ny != PETSC_DECIDE)) SETERRQ(PETSC_COMM_SELF,1,"Incompatible number of processors: Nx * Ny != size");
124: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.mx,user.my,Nx,Ny,1,1,NULL,NULL,&user.da);
125: SNESSetDM(snes,user.da);
126: /*
127: Visualize the distribution of the array across the processors
128: */
129: /* DMView(user.da,PETSC_VIEWER_DRAW_WORLD); */
132: /*
133: Extract global and local vectors from DMDA; then duplicate for remaining
134: vectors that are the same types
135: */
136: DMCreateGlobalVector(user.da,&x);
137: DMCreateLocalVector(user.da,&user.localX);
138: VecDuplicate(x,&r);
139: VecDuplicate(user.localX,&user.localF);
141: /*
142: Set function evaluation routine and vector
143: */
144: SNESSetFunction(snes,r,FormFunction,(void*)&user);
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Create matrix data structure; set Jacobian evaluation routine
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150: /*
151: Set Jacobian matrix data structure and default Jacobian evaluation
152: routine. User can override with:
153: -snes_fd : default finite differencing approximation of Jacobian
154: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
155: (unless user explicitly sets preconditioner)
156: -snes_mf_operator : form preconditioning matrix as set by the user,
157: but use matrix-free approx for Jacobian-vector
158: products within Newton-Krylov method
160: Note: For the parallel case, vectors and matrices MUST be partitioned
161: accordingly. When using distributed arrays (DMDAs) to create vectors,
162: the DMDAs determine the problem partitioning. We must explicitly
163: specify the local matrix dimensions upon its creation for compatibility
164: with the vector distribution. Thus, the generic MatCreate() routine
165: is NOT sufficient when working with distributed arrays.
167: Note: Here we only approximately preallocate storage space for the
168: Jacobian. See the users manual for a discussion of better techniques
169: for preallocating matrix memory.
170: */
171: PetscOptionsGetBool(NULL,"-snes_mf",&matrix_free,NULL);
172: if (!matrix_free) {
173: PetscBool matrix_free_operator = PETSC_FALSE;
174: PetscOptionsGetBool(NULL,"-snes_mf_operator",&matrix_free_operator,NULL);
175: if (matrix_free_operator) matrix_free = PETSC_FALSE;
176: }
177: if (!matrix_free) {
178: if (size == 1) {
179: MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,5,NULL,&J);
180: } else {
181: VecGetLocalSize(x,&m);
182: MatCreateAIJ(PETSC_COMM_WORLD,m,m,N,N,5,NULL,3,NULL,&J);
183: }
184: SNESSetJacobian(snes,J,J,FormJacobian,&user);
185: }
187: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188: Customize nonlinear solver; set runtime options
189: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
191: /*
192: Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
193: */
194: SNESSetFromOptions(snes);
196: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
197: Evaluate initial guess; then solve nonlinear system
198: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
199: /*
200: Note: The user should initialize the vector, x, with the initial guess
201: for the nonlinear solver prior to calling SNESSolve(). In particular,
202: to employ an initial guess of zero, the user should explicitly set
203: this vector to zero by calling VecSet().
204: */
205: FormInitialGuess(&user,x);
206: SNESSolve(snes,NULL,x);
207: SNESGetIterationNumber(snes,&its);
208: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
210: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
211: Free work space. All PETSc objects should be destroyed when they
212: are no longer needed.
213: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
215: if (!matrix_free) {
216: MatDestroy(&J);
217: }
218: VecDestroy(&user.localX); VecDestroy(&x);
219: VecDestroy(&user.localF); VecDestroy(&r);
220: SNESDestroy(&snes); DMDestroy(&user.da);
221: PetscFinalize();
223: return 0;
224: }
225: /* ------------------------------------------------------------------- */
228: /*
229: FormInitialGuess - Forms initial approximation.
231: Input Parameters:
232: user - user-defined application context
233: X - vector
235: Output Parameter:
236: X - vector
237: */
238: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
239: {
240: PetscInt i,j,row,mx,my,xs,ys,xm,ym,gxm,gym,gxs,gys;
242: PetscReal one = 1.0,lambda,temp1,temp,hx,hy,hxdhy,hydhx,sc;
243: PetscScalar *x;
244: Vec localX = user->localX;
246: mx = user->mx; my = user->my; lambda = user->param;
247: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
248: sc = hx*hy*lambda; hxdhy = hx/hy; hydhx = hy/hx;
249: temp1 = lambda/(lambda + one);
251: /*
252: Get a pointer to vector data.
253: - For default PETSc vectors,VecGetArray() returns a pointer to
254: the data array. Otherwise, the routine is implementation dependent.
255: - You MUST call VecRestoreArray() when you no longer need access to
256: the array.
257: */
258: VecGetArray(localX,&x);
260: /*
261: Get local grid boundaries (for 2-dimensional DMDA):
262: xs, ys - starting grid indices (no ghost points)
263: xm, ym - widths of local grid (no ghost points)
264: gxs, gys - starting grid indices (including ghost points)
265: gxm, gym - widths of local grid (including ghost points)
266: */
267: DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);
268: DMDAGetGhostCorners(user->da,&gxs,&gys,NULL,&gxm,&gym,NULL);
270: /*
271: Compute initial guess over the locally owned part of the grid
272: */
273: for (j=ys; j<ys+ym; j++) {
274: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
275: for (i=xs; i<xs+xm; i++) {
276: row = i - gxs + (j - gys)*gxm;
277: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
278: x[row] = 0.0;
279: continue;
280: }
281: x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
282: }
283: }
285: /*
286: Restore vector
287: */
288: VecRestoreArray(localX,&x);
290: /*
291: Insert values into global vector
292: */
293: DMLocalToGlobalBegin(user->da,localX,INSERT_VALUES,X);
294: DMLocalToGlobalEnd(user->da,localX,INSERT_VALUES,X);
295: return 0;
296: }
297: /* ------------------------------------------------------------------- */
300: /*
301: FormFunction - Evaluates nonlinear function, F(x).
303: Input Parameters:
304: . snes - the SNES context
305: . X - input vector
306: . ptr - optional user-defined context, as set by SNESSetFunction()
308: Output Parameter:
309: . F - function vector
310: */
311: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
312: {
313: AppCtx *user = (AppCtx*)ptr;
315: PetscInt i,j,row,mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym;
316: PetscReal two = 2.0,one = 1.0,half = 0.5;
317: PetscReal lambda,hx,hy,hxdhy,hydhx,sc;
318: PetscScalar u,ux,uxx,uyy,*x,*f,kappa;
319: Vec localX = user->localX,localF = user->localF;
321: mx = user->mx; my = user->my; lambda = user->param;
322: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
323: sc = hx*hy*lambda; hxdhy = hx/hy; hydhx = hy/hx;
324: kappa = user->param2;
326: /*
327: Scatter ghost points to local vector, using the 2-step process
328: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
329: By placing code between these two statements, computations can be
330: done while messages are in transition.
331: */
332: DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
333: DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
335: /*
336: Get pointers to vector data
337: */
338: VecGetArray(localX,&x);
339: VecGetArray(localF,&f);
341: /*
342: Get local grid boundaries
343: */
344: DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);
345: DMDAGetGhostCorners(user->da,&gxs,&gys,NULL,&gxm,&gym,NULL);
347: /*
348: Compute function over the locally owned part of the grid
349: */
350: for (j=ys; j<ys+ym; j++) {
351: row = (j - gys)*gxm + xs - gxs - 1;
352: for (i=xs; i<xs+xm; i++) {
353: row++;
354: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
355: f[row] = x[row];
356: continue;
357: }
358: u = x[row];
359: ux = (x[row+1] - x[row-1])*half*hy;
360: uxx = (two*u - x[row-1] - x[row+1])*hydhx;
361: uyy = (two*u - x[row-gxm] - x[row+gxm])*hxdhy;
362: f[row] = uxx + uyy - kappa*ux - sc*exp(u);
363: }
364: }
366: /*
367: Restore vectors
368: */
369: VecRestoreArray(localX,&x);
370: VecRestoreArray(localF,&f);
372: /*
373: Insert values into global vector
374: */
375: DMLocalToGlobalBegin(user->da,localF,INSERT_VALUES,F);
376: DMLocalToGlobalEnd(user->da,localF,INSERT_VALUES,F);
377: PetscLogFlops(11.0*ym*xm);
378: return 0;
379: }
380: /* ------------------------------------------------------------------- */
383: /*
384: FormJacobian - Evaluates Jacobian matrix.
386: Input Parameters:
387: . snes - the SNES context
388: . x - input vector
389: . ptr - optional user-defined context, as set by SNESSetJacobian()
391: Output Parameters:
392: . A - Jacobian matrix
393: . B - optionally different preconditioning matrix
394: . flag - flag indicating matrix structure
396: Notes:
397: Due to grid point reordering with DMDAs, we must always work
398: with the local grid points, and then transform them to the new
399: global numbering with the "ltog" mapping (via DMDAGetGlobalIndices()).
400: We cannot work directly with the global numbers for the original
401: uniprocessor grid!
402: */
403: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
404: {
405: AppCtx *user = (AppCtx*)ptr; /* user-defined application context */
406: Mat jac = *B; /* Jacobian matrix */
407: Vec localX = user->localX; /* local vector */
409: PetscInt *ltog; /* local-to-global mapping */
410: PetscInt i,j,row,mx,my,col[5];
411: PetscInt nloc,xs,ys,xm,ym,gxs,gys,gxm,gym,grow;
412: PetscScalar two = 2.0,one = 1.0,lambda,v[5],hx,hy,hxdhy,hydhx,sc,*x;
414: mx = user->mx; my = user->my; lambda = user->param;
415: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
416: sc = hx*hy; hxdhy = hx/hy; hydhx = hy/hx;
418: /*
419: Scatter ghost points to local vector,using the 2-step process
420: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
421: By placing code between these two statements, computations can be
422: done while messages are in transition.
423: */
424: DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
425: DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
427: /*
428: Get pointer to vector data
429: */
430: VecGetArray(localX,&x);
432: /*
433: Get local grid boundaries
434: */
435: DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);
436: DMDAGetGhostCorners(user->da,&gxs,&gys,NULL,&gxm,&gym,NULL);
438: /*
439: Get the global node numbers for all local nodes, including ghost points
440: */
441: DMDAGetGlobalIndices(user->da,&nloc,<og);
443: /*
444: Compute entries for the locally owned part of the Jacobian.
445: - Currently, all PETSc parallel matrix formats are partitioned by
446: contiguous chunks of rows across the processors. The "grow"
447: parameter computed below specifies the global row number
448: corresponding to each local grid point.
449: - Each processor needs to insert only elements that it owns
450: locally (but any non-local elements will be sent to the
451: appropriate processor during matrix assembly).
452: - Always specify global row and columns of matrix entries.
453: - Here, we set all entries for a particular row at once.
454: */
455: for (j=ys; j<ys+ym; j++) {
456: row = (j - gys)*gxm + xs - gxs - 1;
457: for (i=xs; i<xs+xm; i++) {
458: row++;
459: grow = ltog[row];
460: /* boundary points */
461: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
462: MatSetValues(jac,1,&grow,1,&grow,&one,INSERT_VALUES);
463: continue;
464: }
465: /* interior grid points */
466: v[0] = -hxdhy; col[0] = ltog[row - gxm];
467: v[1] = -hydhx; col[1] = ltog[row - 1];
468: v[2] = two*(hydhx + hxdhy) - sc*lambda*exp(x[row]); col[2] = grow;
469: v[3] = -hydhx; col[3] = ltog[row + 1];
470: v[4] = -hxdhy; col[4] = ltog[row + gxm];
471: MatSetValues(jac,1,&grow,5,col,v,INSERT_VALUES);
472: }
473: }
475: /*
476: Assemble matrix, using the 2-step process:
477: MatAssemblyBegin(), MatAssemblyEnd().
478: By placing code between these two statements, computations can be
479: done while messages are in transition.
480: */
481: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
482: VecRestoreArray(localX,&x);
483: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
485: /*
486: Set flag to indicate that the Jacobian matrix retains an identical
487: nonzero structure throughout all nonlinear iterations (although the
488: values of the entries change). Thus, we can save some work in setting
489: up the preconditioner (e.g., no need to redo symbolic factorization for
490: ILU/ICC preconditioners).
491: - If the nonzero structure of the matrix is different during
492: successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
493: must be used instead. If you are unsure whether the matrix
494: structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
495: - Caution: If you specify SAME_NONZERO_PATTERN, PETSc
496: believes your assertion and does not check the structure
497: of the matrix. If you erroneously claim that the structure
498: is the same when it actually is not, the new preconditioner
499: will not function correctly. Thus, use this optimization
500: feature with caution!
501: */
502: *flag = SAME_NONZERO_PATTERN;
503: /*
504: Tell the matrix we will never add a new nonzero location to the
505: matrix. If we do it will generate an error.
506: */
507: /* MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE); */
508: return 0;
509: }