Actual source code: ex35.c
petsc-3.3-p7 2013-05-11
1: static const char help[] = "-Laplacian u = b as a nonlinear problem.\n\n";
3: /*T
4: Concepts: SNES^parallel Bratu example
5: Concepts: DMDA^using distributed arrays;
6: Concepts: IS coloirng types;
7: Processors: n
8: T*/
10: /*
12: The linear and nonlinear versions of these should give almost identical results on this problem
14: Richardson
15: Nonlinear:
16: -snes_rtol 1.e-12 -snes_monitor -snes_type nrichardson -snes_linesearch_monitor
18: Linear:
19: -snes_rtol 1.e-12 -snes_monitor -ksp_rtol 1.e-12 -ksp_monitor -ksp_type richardson -pc_type none -ksp_richardson_self_scale -info
21: GMRES
22: Nonlinear:
23: -snes_rtol 1.e-12 -snes_monitor -snes_type ngmres
25: Linear:
26: -snes_rtol 1.e-12 -snes_monitor -ksp_type gmres -ksp_monitor -ksp_rtol 1.e-12 -pc_type none
28: CG
29: Nonlinear:
30: -snes_rtol 1.e-12 -snes_monitor -snes_type ncg -snes_linesearch_monitor
32: Linear:
33: -snes_rtol 1.e-12 -snes_monitor -ksp_type cg -ksp_monitor -ksp_rtol 1.e-12 -pc_type none
35: Multigrid
36: Linear:
37: 1 level:
38: -snes_rtol 1.e-12 -snes_monitor -pc_type mg -mg_levels_ksp_type richardson -mg_levels_pc_type none -mg_levels_ksp_monitor
39: -mg_levels_ksp_richardson_self_scale -ksp_type richardson -ksp_monitor -ksp_rtol 1.e-12 -ksp_monitor_true_residual
41: n levels:
42: -da_refine n
44: Nonlinear:
45: 1 level:
46: -snes_rtol 1.e-12 -snes_monitor -snes_type fas -fas_levels_snes_monitor
48: n levels:
49: -da_refine n -fas_coarse_snes_type ls -fas_coarse_pc_type lu -fas_coarse_ksp_type preonly
51: */
53: /*
54: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
55: Include "petscsnes.h" so that we can use SNES solvers. Note that this
56: */
57: #include <petscdmda.h>
58: #include <petscsnes.h>
60: /*
61: User-defined routines
62: */
63: extern PetscErrorCode FormMatrix(DM,Mat);
64: extern PetscErrorCode MyDMComputeFunction(DM,Vec,Vec);
65: extern PetscErrorCode MyDMComputeJacobian(DM,Vec,Mat,Mat,MatStructure*);
66: extern PetscErrorCode NonlinearGS(SNES,Vec);
70: int main(int argc,char **argv)
71: {
72: SNES snes; /* nonlinear solver */
73: SNES psnes; /* nonlinear Gauss-Seidel approximate solver */
74: Vec x,b; /* solution vector */
75: PetscInt its; /* iterations for convergence */
76: PetscErrorCode ierr;
77: DM da;
78: PetscBool use_ngs = PETSC_FALSE; /* use the nonlinear Gauss-Seidel approximate solver */
80: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: Initialize program
82: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84: PetscInitialize(&argc,&argv,(char *)0,help);
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Create nonlinear solver context
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: SNESCreate(PETSC_COMM_WORLD,&snes);
91: PetscOptionsGetBool(PETSC_NULL,"-use_ngs",&use_ngs,0);
93: if (use_ngs) {
94: SNESGetPC(snes,&psnes);
95: SNESSetType(psnes,SNESSHELL);
96: SNESShellSetSolve(psnes,NonlinearGS);
97: }
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Create distributed array (DMDA) to manage parallel grid and vectors
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: 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);
103: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
104: SNESSetDM(snes,da);
105: if (use_ngs) {
106: SNESShellSetContext(psnes,da);
107: }
108: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: Extract global vectors from DMDA; then duplicate for remaining
110: vectors that are the same types
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: DMCreateGlobalVector(da,&x);
113: DMCreateGlobalVector(da,&b);
114: VecSetRandom(b,PETSC_NULL);
116: DMSetFunction(da,MyDMComputeFunction);
117: DMSetJacobian(da,MyDMComputeJacobian);
120: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: Customize nonlinear solver; set runtime options
122: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123: SNESSetFromOptions(snes);
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Solve nonlinear system
127: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128: SNESSolve(snes,b,x);
129: SNESGetIterationNumber(snes,&its);
131: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Free work space. All PETSc objects should be destroyed when they
136: are no longer needed.
137: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: VecDestroy(&x);
139: VecDestroy(&b);
140: SNESDestroy(&snes);
141: DMDestroy(&da);
142: PetscFinalize();
144: return(0);
145: }
147: /* ------------------------------------------------------------------- */
150: PetscErrorCode MyDMComputeFunction(DM dm,Vec x,Vec F)
151: {
153: Mat J;
156: DMGetApplicationContext(dm,&J);
157: if (!J) {
158: DMCreateMatrix(dm,MATAIJ,&J);
159: PetscObjectCompose((PetscObject)J,"DM",(PetscObject)PETSC_NULL);
160: FormMatrix(dm,J);
161: DMSetApplicationContext(dm,J);
162: DMSetApplicationContextDestroy(dm,(PetscErrorCode (*)(void**))MatDestroy);
163: }
164: MatMult(J,x,F);
165: return(0);
166: }
170: PetscErrorCode MyDMComputeJacobian(DM dm,Vec x,Mat J,Mat Jp,MatStructure *str)
171: {
175: FormMatrix(dm,Jp);
176: return(0);
177: }
181: PetscErrorCode FormMatrix(DM da,Mat jac)
182: {
184: PetscInt i,j,nrows = 0;
185: MatStencil col[5],row,*rows;
186: PetscScalar v[5],hx,hy,hxdhy,hydhx;
187: DMDALocalInfo info;
190: DMDAGetLocalInfo(da,&info);
191: hx = 1.0/(PetscReal)(info.mx-1);
192: hy = 1.0/(PetscReal)(info.my-1);
193: hxdhy = hx/hy;
194: hydhx = hy/hx;
196: PetscMalloc(info.ym*info.xm*sizeof(MatStencil),&rows);
197: /*
198: Compute entries for the locally owned part of the Jacobian.
199: - Currently, all PETSc parallel matrix formats are partitioned by
200: contiguous chunks of rows across the processors.
201: - Each processor needs to insert only elements that it owns
202: locally (but any non-local elements will be sent to the
203: appropriate processor during matrix assembly).
204: - Here, we set all entries for a particular row at once.
205: - We can set matrix entries either using either
206: MatSetValuesLocal() or MatSetValues(), as discussed above.
207: */
208: for (j=info.ys; j<info.ys+info.ym; j++) {
209: for (i=info.xs; i<info.xs+info.xm; i++) {
210: row.j = j; row.i = i;
211: /* boundary points */
212: if (i == 0 || j == 0 || i == info.mx-1 || j == info.my-1) {
213: v[0] = 2.0*(hydhx + hxdhy);
214: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
215: rows[nrows].i = i;
216: rows[nrows++].j = j;
217: } else {
218: /* interior grid points */
219: v[0] = -hxdhy; col[0].j = j - 1; col[0].i = i;
220: v[1] = -hydhx; col[1].j = j; col[1].i = i-1;
221: v[2] = 2.0*(hydhx + hxdhy); col[2].j = row.j; col[2].i = row.i;
222: v[3] = -hydhx; col[3].j = j; col[3].i = i+1;
223: v[4] = -hxdhy; col[4].j = j + 1; col[4].i = i;
224: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
225: }
226: }
227: }
229: /*
230: Assemble matrix, using the 2-step process:
231: MatAssemblyBegin(), MatAssemblyEnd().
232: */
233: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
234: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
235: MatZeroRowsColumnsStencil(jac,nrows,rows,2.0*(hydhx + hxdhy),PETSC_NULL,PETSC_NULL);
236: PetscFree(rows);
237: /*
238: Tell the matrix we will never add a new nonzero location to the
239: matrix. If we do, it will generate an error.
240: */
241: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
242: return(0);
243: }
247: /* ------------------------------------------------------------------- */
250: /*
251: Applies some sweeps on nonlinear Gauss-Seidel on each process
253: */
254: PetscErrorCode NonlinearGS(SNES snes,Vec X)
255: {
256: PetscInt i,j,Mx,My,xs,ys,xm,ym,its,l;
258: PetscReal hx,hy,hxdhy,hydhx;
259: PetscScalar **x,F,J,u,uxx,uyy;
260: DM da;
261: Vec localX;
264: SNESGetTolerances(snes,PETSC_NULL,PETSC_NULL,PETSC_NULL,&its,PETSC_NULL);
265: SNESShellGetContext(snes,(void**)&da);
267: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
268: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
270: hx = 1.0/(PetscReal)(Mx-1);
271: hy = 1.0/(PetscReal)(My-1);
272: hxdhy = hx/hy;
273: hydhx = hy/hx;
276: DMGetLocalVector(da,&localX);
278: for (l=0; l<its; l++) {
280: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
281: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
282: /*
283: Get a pointer to vector data.
284: - For default PETSc vectors, VecGetArray() returns a pointer to
285: the data array. Otherwise, the routine is implementation dependent.
286: - You MUST call VecRestoreArray() when you no longer need access to
287: the array.
288: */
289: DMDAVecGetArray(da,localX,&x);
291: /*
292: Get local grid boundaries (for 2-dimensional DMDA):
293: xs, ys - starting grid indices (no ghost points)
294: xm, ym - widths of local grid (no ghost points)
295:
296: */
297: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
298:
299: for (j=ys; j<ys+ym; j++) {
300: for (i=xs; i<xs+xm; i++) {
301: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
302: /* boundary conditions are all zero Dirichlet */
303: x[j][i] = 0.0;
304: } else {
305: u = x[j][i];
307: uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
308: uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
309: F = uxx + uyy;
310: J = 2.0*(hydhx + hxdhy);
311: u = u - F/J;
312:
313: x[j][i] = u;
314: }
315: }
316: }
317:
318: /*
319: Restore vector
320: */
321: DMDAVecRestoreArray(da,localX,&x);
322: DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
323: DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
324: }
325: DMRestoreLocalVector(da,&localX);
326: return(0);
327: }