Actual source code: ex14.c
petsc-3.3-p7 2013-05-11
2: /* Program usage: mpiexec -n <procs> ex14 [-help] [all PETSc options] */
4: static char help[] = "Bratu nonlinear PDE in 3d.\n\
5: We solve the Bratu (SFI - solid fuel ignition) problem in a 3D rectangular\n\
6: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\
7: The command line options include:\n\
8: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
9: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\n";
11: /*T
12: Concepts: SNES^parallel Bratu example
13: Concepts: DMDA^using distributed arrays;
14: Processors: n
15: T*/
17: /* ------------------------------------------------------------------------
19: Solid Fuel Ignition (SFI) problem. This problem is modeled by
20: the partial differential equation
21:
22: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
23:
24: with boundary conditions
25:
26: u = 0 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1
27:
28: A finite difference approximation with the usual 7-point stencil
29: is used to discretize the boundary value problem to obtain a nonlinear
30: system of equations.
33: ------------------------------------------------------------------------- */
35: /*
36: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
37: Include "petscsnes.h" so that we can use SNES solvers. Note that this
38: file automatically includes:
39: petscsys.h - base PETSc routines petscvec.h - vectors
40: petscmat.h - matrices
41: petscis.h - index sets petscksp.h - Krylov subspace methods
42: petscviewer.h - viewers petscpc.h - preconditioners
43: petscksp.h - linear solvers
44: */
45: #include <petscdmda.h>
46: #include <petscsnes.h>
49: /*
50: User-defined application context - contains data needed by the
51: application-provided call-back routines, FormJacobian() and
52: FormFunction().
53: */
54: typedef struct {
55: PetscReal param; /* test problem parameter */
56: DM da; /* distributed array data structure */
57: } AppCtx;
59: /*
60: User-defined routines
61: */
62: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*),FormInitialGuess(AppCtx*,Vec);
63: extern PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
67: int main(int argc,char **argv)
68: {
69: SNES snes; /* nonlinear solver */
70: Vec x,r; /* solution, residual vectors */
71: Mat J; /* Jacobian matrix */
72: AppCtx user; /* user-defined work context */
73: PetscInt its; /* iterations for convergence */
74: PetscBool matrix_free = PETSC_FALSE,coloring = PETSC_FALSE;
75: PetscErrorCode ierr;
76: PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.,fnorm;
77: MatFDColoring matfdcoloring;
79: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80: Initialize program
81: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
83: PetscInitialize(&argc,&argv,(char *)0,help);
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Initialize problem parameters
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: user.param = 6.0;
89: PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
90: if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,1,"Lambda is out of range");
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Create nonlinear solver context
94: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
95: SNESCreate(PETSC_COMM_WORLD,&snes);
97: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98: Create distributed array (DMDA) to manage parallel grid and vectors
99: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100: DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,-4,PETSC_DECIDE,PETSC_DECIDE,
101: PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,PETSC_NULL,&user.da);
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Extract global vectors from DMDA; then duplicate for remaining
105: vectors that are the same types
106: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: DMCreateGlobalVector(user.da,&x);
108: VecDuplicate(x,&r);
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Set function evaluation routine and vector
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: SNESSetFunction(snes,r,FormFunction,(void*)&user);
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Create matrix data structure; set Jacobian evaluation routine
118: Set Jacobian matrix data structure and default Jacobian evaluation
119: routine. User can override with:
120: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
121: (unless user explicitly sets preconditioner)
122: -snes_mf_operator : form preconditioning matrix as set by the user,
123: but use matrix-free approx for Jacobian-vector
124: products within Newton-Krylov method
125: -fdcoloring : using finite differences with coloring to compute the Jacobian
127: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128: PetscOptionsGetBool(PETSC_NULL,"-snes_mf",&matrix_free,PETSC_NULL);
129: PetscOptionsGetBool(PETSC_NULL,"-fdcoloring",&coloring,PETSC_NULL);
130: if (!matrix_free) {
131: if (coloring) {
132: ISColoring iscoloring;
134: DMCreateColoring(user.da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
135: DMCreateMatrix(user.da,MATAIJ,&J);
136: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
137: ISColoringDestroy(&iscoloring);
138: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
139: MatFDColoringSetFromOptions(matfdcoloring);
140: SNESSetJacobian(snes,J,J,SNESDefaultComputeJacobianColor,matfdcoloring);
141: } else {
142: DMCreateMatrix(user.da,MATAIJ,&J);
143: SNESSetJacobian(snes,J,J,FormJacobian,&user);
144: }
145: }
147: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: Customize nonlinear solver; set runtime options
149: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150: SNESSetFromOptions(snes);
152: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: Evaluate initial guess
154: Note: The user should initialize the vector, x, with the initial guess
155: for the nonlinear solver prior to calling SNESSolve(). In particular,
156: to employ an initial guess of zero, the user should explicitly set
157: this vector to zero by calling VecSet().
158: */
159: FormInitialGuess(&user,x);
161: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: Solve nonlinear system
163: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
164: SNESSolve(snes,PETSC_NULL,x);
165: SNESGetIterationNumber(snes,&its);
167: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: Explicitly check norm of the residual of the solution
169: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170: FormFunction(snes,x,r,(void*)&user);
171: VecNorm(r,NORM_2,&fnorm);
172: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D fnorm %G\n",its,fnorm);
174: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: Free work space. All PETSc objects should be destroyed when they
176: are no longer needed.
177: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
179: if (!matrix_free) {
180: MatDestroy(&J);
181: }
182: if (coloring) {
183: MatFDColoringDestroy(&matfdcoloring);
184: }
185: VecDestroy(&x);
186: VecDestroy(&r);
187: SNESDestroy(&snes);
188: DMDestroy(&user.da);
189: PetscFinalize();
191: return(0);
192: }
193: /* ------------------------------------------------------------------- */
196: /*
197: FormInitialGuess - Forms initial approximation.
199: Input Parameters:
200: user - user-defined application context
201: X - vector
203: Output Parameter:
204: X - vector
205: */
206: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
207: {
208: PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
210: PetscReal lambda,temp1,hx,hy,hz,tempk,tempj;
211: PetscScalar ***x;
214: DMDAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,
215: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
217: lambda = user->param;
218: hx = 1.0/(PetscReal)(Mx-1);
219: hy = 1.0/(PetscReal)(My-1);
220: hz = 1.0/(PetscReal)(Mz-1);
221: temp1 = lambda/(lambda + 1.0);
223: /*
224: Get a pointer to vector data.
225: - For default PETSc vectors, VecGetArray() returns a pointer to
226: the data array. Otherwise, the routine is implementation dependent.
227: - You MUST call VecRestoreArray() when you no longer need access to
228: the array.
229: */
230: DMDAVecGetArray(user->da,X,&x);
232: /*
233: Get local grid boundaries (for 3-dimensional DMDA):
234: xs, ys, zs - starting grid indices (no ghost points)
235: xm, ym, zm - widths of local grid (no ghost points)
237: */
238: DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);
240: /*
241: Compute initial guess over the locally owned part of the grid
242: */
243: for (k=zs; k<zs+zm; k++) {
244: tempk = (PetscReal)(PetscMin(k,Mz-k-1))*hz;
245: for (j=ys; j<ys+ym; j++) {
246: tempj = PetscMin((PetscReal)(PetscMin(j,My-j-1))*hy,tempk);
247: for (i=xs; i<xs+xm; i++) {
248: if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
249: /* boundary conditions are all zero Dirichlet */
250: x[k][j][i] = 0.0;
251: } else {
252: x[k][j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,tempj));
253: }
254: }
255: }
256: }
258: /*
259: Restore vector
260: */
261: DMDAVecRestoreArray(user->da,X,&x);
262: return(0);
263: }
264: /* ------------------------------------------------------------------- */
267: /*
268: FormFunction - Evaluates nonlinear function, F(x).
270: Input Parameters:
271: . snes - the SNES context
272: . X - input vector
273: . ptr - optional user-defined context, as set by SNESSetFunction()
275: Output Parameter:
276: . F - function vector
277: */
278: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
279: {
280: AppCtx *user = (AppCtx*)ptr;
282: PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
283: PetscReal two = 2.0,lambda,hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc;
284: PetscScalar u_north,u_south,u_east,u_west,u_up,u_down,u,u_xx,u_yy,u_zz,***x,***f;
285: Vec localX;
288: DMGetLocalVector(user->da,&localX);
289: DMDAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,
290: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
292: lambda = user->param;
293: hx = 1.0/(PetscReal)(Mx-1);
294: hy = 1.0/(PetscReal)(My-1);
295: hz = 1.0/(PetscReal)(Mz-1);
296: sc = hx*hy*hz*lambda;
297: hxhzdhy = hx*hz/hy;
298: hyhzdhx = hy*hz/hx;
299: hxhydhz = hx*hy/hz;
301: /*
302: Scatter ghost points to local vector,using the 2-step process
303: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
304: By placing code between these two statements, computations can be
305: done while messages are in transition.
306: */
307: DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
308: DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
310: /*
311: Get pointers to vector data
312: */
313: DMDAVecGetArray(user->da,localX,&x);
314: DMDAVecGetArray(user->da,F,&f);
316: /*
317: Get local grid boundaries
318: */
319: DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);
321: /*
322: Compute function over the locally owned part of the grid
323: */
324: for (k=zs; k<zs+zm; k++) {
325: for (j=ys; j<ys+ym; j++) {
326: for (i=xs; i<xs+xm; i++) {
327: if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
328: f[k][j][i] = x[k][j][i];
329: } else {
330: u = x[k][j][i];
331: u_east = x[k][j][i+1];
332: u_west = x[k][j][i-1];
333: u_north = x[k][j+1][i];
334: u_south = x[k][j-1][i];
335: u_up = x[k+1][j][i];
336: u_down = x[k-1][j][i];
337: u_xx = (-u_east + two*u - u_west)*hyhzdhx;
338: u_yy = (-u_north + two*u - u_south)*hxhzdhy;
339: u_zz = (-u_up + two*u - u_down)*hxhydhz;
340: f[k][j][i] = u_xx + u_yy + u_zz - sc*PetscExpScalar(u);
341: }
342: }
343: }
344: }
346: /*
347: Restore vectors
348: */
349: DMDAVecRestoreArray(user->da,localX,&x);
350: DMDAVecRestoreArray(user->da,F,&f);
351: DMRestoreLocalVector(user->da,&localX);
352: PetscLogFlops(11.0*ym*xm);
353: return(0);
354: }
355: /* ------------------------------------------------------------------- */
358: /*
359: FormJacobian - Evaluates Jacobian matrix.
361: Input Parameters:
362: . snes - the SNES context
363: . x - input vector
364: . ptr - optional user-defined context, as set by SNESSetJacobian()
366: Output Parameters:
367: . A - Jacobian matrix
368: . B - optionally different preconditioning matrix
369: . flag - flag indicating matrix structure
371: */
372: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
373: {
374: AppCtx *user = (AppCtx*)ptr; /* user-defined application context */
375: Mat jac = *B; /* Jacobian matrix */
376: Vec localX;
378: PetscInt i,j,k,Mx,My,Mz;
379: MatStencil col[7],row;
380: PetscInt xs,ys,zs,xm,ym,zm;
381: PetscScalar lambda,v[7],hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc,***x;
385: DMGetLocalVector(user->da,&localX);
386: DMDAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,
387: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
389: lambda = user->param;
390: hx = 1.0/(PetscReal)(Mx-1);
391: hy = 1.0/(PetscReal)(My-1);
392: hz = 1.0/(PetscReal)(Mz-1);
393: sc = hx*hy*hz*lambda;
394: hxhzdhy = hx*hz/hy;
395: hyhzdhx = hy*hz/hx;
396: hxhydhz = hx*hy/hz;
398: /*
399: Scatter ghost points to local vector, using the 2-step process
400: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
401: By placing code between these two statements, computations can be
402: done while messages are in transition.
403: */
404: DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
405: DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
407: /*
408: Get pointer to vector data
409: */
410: DMDAVecGetArray(user->da,localX,&x);
412: /*
413: Get local grid boundaries
414: */
415: DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);
417: /*
418: Compute entries for the locally owned part of the Jacobian.
419: - Currently, all PETSc parallel matrix formats are partitioned by
420: contiguous chunks of rows across the processors.
421: - Each processor needs to insert only elements that it owns
422: locally (but any non-local elements will be sent to the
423: appropriate processor during matrix assembly).
424: - Here, we set all entries for a particular row at once.
425: - We can set matrix entries either using either
426: MatSetValuesLocal() or MatSetValues(), as discussed above.
427: */
428: for (k=zs; k<zs+zm; k++) {
429: for (j=ys; j<ys+ym; j++) {
430: for (i=xs; i<xs+xm; i++) {
431: row.k = k; row.j = j; row.i = i;
432: /* boundary points */
433: if (i == 0 || j == 0 || k == 0|| i == Mx-1 || j == My-1 || k == Mz-1) {
434: v[0] = 1.0;
435: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
436: } else {
437: /* interior grid points */
438: v[0] = -hxhydhz; col[0].k=k-1;col[0].j=j; col[0].i = i;
439: v[1] = -hxhzdhy; col[1].k=k; col[1].j=j-1;col[1].i = i;
440: v[2] = -hyhzdhx; col[2].k=k; col[2].j=j; col[2].i = i-1;
441: v[3] = 2.0*(hyhzdhx+hxhzdhy+hxhydhz)-sc*PetscExpScalar(x[k][j][i]);col[3].k=row.k;col[3].j=row.j;col[3].i = row.i;
442: v[4] = -hyhzdhx; col[4].k=k; col[4].j=j; col[4].i = i+1;
443: v[5] = -hxhzdhy; col[5].k=k; col[5].j=j+1;col[5].i = i;
444: v[6] = -hxhydhz; col[6].k=k+1;col[6].j=j; col[6].i = i;
445: MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);
446: }
447: }
448: }
449: }
450: DMDAVecRestoreArray(user->da,localX,&x);
451: DMRestoreLocalVector(user->da,&localX);
453: /*
454: Assemble matrix, using the 2-step process:
455: MatAssemblyBegin(), MatAssemblyEnd().
456: */
457: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
458: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
460: /*
461: Normally since the matrix has already been assembled above; this
462: would do nothing. But in the matrix free mode -snes_mf_operator
463: this tells the "matrix-free" matrix that a new linear system solve
464: is about to be done.
465: */
467: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
468: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
470: /*
471: Set flag to indicate that the Jacobian matrix retains an identical
472: nonzero structure throughout all nonlinear iterations (although the
473: values of the entries change). Thus, we can save some work in setting
474: up the preconditioner (e.g., no need to redo symbolic factorization for
475: ILU/ICC preconditioners).
476: - If the nonzero structure of the matrix is different during
477: successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
478: must be used instead. If you are unsure whether the matrix
479: structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
480: - Caution: If you specify SAME_NONZERO_PATTERN, PETSc
481: believes your assertion and does not check the structure
482: of the matrix. If you erroneously claim that the structure
483: is the same when it actually is not, the new preconditioner
484: will not function correctly. Thus, use this optimization
485: feature with caution!
486: */
487: *flag = SAME_NONZERO_PATTERN;
490: /*
491: Tell the matrix we will never add a new nonzero location to the
492: matrix. If we do, it will generate an error.
493: */
494: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
495: return(0);
496: }