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