Actual source code: ex35.c
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*/
12: /*
14: The linear and nonlinear versions of these should give almost identical results on this problem
16: Richardson
17: Nonlinear:
18: -snes_rtol 1.e-12 -snes_monitor -snes_type nrichardson -snes_linesearch_monitor
20: Linear:
21: -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
23: GMRES
24: Nonlinear:
25: -snes_rtol 1.e-12 -snes_monitor -snes_type ngmres
27: Linear:
28: -snes_rtol 1.e-12 -snes_monitor -ksp_type gmres -ksp_monitor -ksp_rtol 1.e-12 -pc_type none
30: CG
31: Nonlinear:
32: -snes_rtol 1.e-12 -snes_monitor -snes_type ncg -snes_linesearch_monitor
34: Linear:
35: -snes_rtol 1.e-12 -snes_monitor -ksp_type cg -ksp_monitor -ksp_rtol 1.e-12 -pc_type none
37: Multigrid
38: Linear:
39: 1 level:
40: -snes_rtol 1.e-12 -snes_monitor -pc_type mg -mg_levels_ksp_type richardson -mg_levels_pc_type none -mg_levels_ksp_monitor
41: -mg_levels_ksp_richardson_self_scale -ksp_type richardson -ksp_monitor -ksp_rtol 1.e-12 -ksp_monitor_true_residual
43: n levels:
44: -da_refine n
46: Nonlinear:
47: 1 level:
48: -snes_rtol 1.e-12 -snes_monitor -snes_type fas -fas_levels_snes_monitor
50: n levels:
51: -da_refine n -fas_coarse_snes_type newtonls -fas_coarse_pc_type lu -fas_coarse_ksp_type preonly
53: */
55: /*
56: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
57: Include "petscsnes.h" so that we can use SNES solvers. Note that this
58: */
59: #include <petscdm.h>
60: #include <petscdmda.h>
61: #include <petscsnes.h>
63: /*
64: User-defined routines
65: */
66: extern PetscErrorCode FormMatrix(DM,Mat);
67: extern PetscErrorCode MyComputeFunction(SNES,Vec,Vec,void*);
68: extern PetscErrorCode MyComputeJacobian(SNES,Vec,Mat,Mat,void*);
69: extern PetscErrorCode NonlinearGS(SNES,Vec);
71: int main(int argc,char **argv)
72: {
73: SNES snes; /* nonlinear solver */
74: SNES psnes; /* nonlinear Gauss-Seidel approximate solver */
75: Vec x,b; /* solution vector */
76: PetscInt its; /* iterations for convergence */
78: DM da;
79: PetscBool use_ngs_as_npc = PETSC_FALSE; /* use the nonlinear Gauss-Seidel approximate solver */
81: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82: Initialize program
83: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
87: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: Create nonlinear solver context
89: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90: SNESCreate(PETSC_COMM_WORLD,&snes);
92: PetscOptionsGetBool(NULL,NULL,"-use_ngs_as_npc",&use_ngs_as_npc,0);
94: if (use_ngs_as_npc) {
95: SNESGetNPC(snes,&psnes);
96: SNESSetType(psnes,SNESSHELL);
97: SNESShellSetSolve(psnes,NonlinearGS);
98: }
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Create distributed array (DMDA) to manage parallel grid and vectors
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
104: DMSetFromOptions(da);
105: DMSetUp(da);
106: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
107: SNESSetDM(snes,da);
108: if (use_ngs_as_npc) {
109: SNESShellSetContext(psnes,da);
110: }
111: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: Extract global vectors from DMDA; then duplicate for remaining
113: vectors that are the same types
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: DMCreateGlobalVector(da,&x);
116: DMCreateGlobalVector(da,&b);
117: VecSet(b,1.0);
119: SNESSetFunction(snes,NULL,MyComputeFunction,NULL);
120: SNESSetJacobian(snes,NULL,NULL,MyComputeJacobian,NULL);
122: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: Customize nonlinear solver; set runtime options
124: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: SNESSetFromOptions(snes);
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Solve nonlinear system
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: SNESSolve(snes,b,x);
131: SNESGetIterationNumber(snes,&its);
133: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: Free work space. All PETSc objects should be destroyed when they
138: are no longer needed.
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: VecDestroy(&x);
141: VecDestroy(&b);
142: SNESDestroy(&snes);
143: DMDestroy(&da);
144: PetscFinalize();
145: return ierr;
146: }
148: /* ------------------------------------------------------------------- */
149: PetscErrorCode MyComputeFunction(SNES snes,Vec x,Vec F,void *ctx)
150: {
152: Mat J;
153: DM dm;
156: SNESGetDM(snes,&dm);
157: DMGetApplicationContext(dm,&J);
158: if (!J) {
159: DMSetMatType(dm,MATAIJ);
160: DMCreateMatrix(dm,&J);
161: MatSetDM(J, NULL);
162: FormMatrix(dm,J);
163: DMSetApplicationContext(dm,J);
164: DMSetApplicationContextDestroy(dm,(PetscErrorCode (*)(void**))MatDestroy);
165: }
166: MatMult(J,x,F);
167: return(0);
168: }
170: PetscErrorCode MyComputeJacobian(SNES snes,Vec x,Mat J,Mat Jp,void *ctx)
171: {
173: DM dm;
176: SNESGetDM(snes,&dm);
177: FormMatrix(dm,Jp);
178: return(0);
179: }
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: PetscMalloc1(info.ym*info.xm,&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),NULL,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: /* ------------------------------------------------------------------- */
248: /*
249: Applies some sweeps on nonlinear Gauss-Seidel on each process
251: */
252: PetscErrorCode NonlinearGS(SNES snes,Vec X)
253: {
254: PetscInt i,j,Mx,My,xs,ys,xm,ym,its,l;
256: PetscReal hx,hy,hxdhy,hydhx;
257: PetscScalar **x,F,J,u,uxx,uyy;
258: DM da;
259: Vec localX;
262: SNESGetTolerances(snes,NULL,NULL,NULL,&its,NULL);
263: SNESShellGetContext(snes,(void**)&da);
265: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
267: hx = 1.0/(PetscReal)(Mx-1);
268: hy = 1.0/(PetscReal)(My-1);
269: hxdhy = hx/hy;
270: hydhx = hy/hx;
273: DMGetLocalVector(da,&localX);
275: for (l=0; l<its; l++) {
277: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
278: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
279: /*
280: Get a pointer to vector data.
281: - For default PETSc vectors, VecGetArray() returns a pointer to
282: the data array. Otherwise, the routine is implementation dependent.
283: - You MUST call VecRestoreArray() when you no longer need access to
284: the array.
285: */
286: DMDAVecGetArray(da,localX,&x);
288: /*
289: Get local grid boundaries (for 2-dimensional DMDA):
290: xs, ys - starting grid indices (no ghost points)
291: xm, ym - widths of local grid (no ghost points)
293: */
294: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
296: for (j=ys; j<ys+ym; j++) {
297: for (i=xs; i<xs+xm; i++) {
298: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
299: /* boundary conditions are all zero Dirichlet */
300: x[j][i] = 0.0;
301: } else {
302: u = x[j][i];
303: uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
304: uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
305: F = uxx + uyy;
306: J = 2.0*(hydhx + hxdhy);
307: u = u - F/J;
309: x[j][i] = u;
310: }
311: }
312: }
314: /*
315: Restore vector
316: */
317: DMDAVecRestoreArray(da,localX,&x);
318: DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
319: DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
320: }
321: DMRestoreLocalVector(da,&localX);
322: return(0);
323: }
326: /*TEST
328: test:
329: args: -snes_monitor_short -snes_type nrichardson
330: requires: !single
332: test:
333: suffix: 2
334: args: -snes_monitor_short -ksp_monitor_short -ksp_type richardson -pc_type none -ksp_richardson_self_scale
335: requires: !single
337: test:
338: suffix: 3
339: args: -snes_monitor_short -snes_type ngmres
341: test:
342: suffix: 4
343: args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none
345: test:
346: suffix: 5
347: args: -snes_monitor_short -snes_type ncg
349: test:
350: suffix: 6
351: args: -snes_monitor_short -ksp_type cg -ksp_monitor_short -pc_type none
353: test:
354: suffix: 7
355: args: -da_refine 2 -snes_monitor_short -pc_type mg -mg_levels_ksp_type richardson -mg_levels_pc_type none -mg_levels_ksp_monitor_short -mg_levels_ksp_richardson_self_scale -ksp_type richardson -ksp_monitor_short
356: requires: !single
358: test:
359: suffix: 8
360: args: -da_refine 2 -snes_monitor_short -snes_type fas -fas_levels_snes_monitor_short -fas_coarse_snes_type newtonls -fas_coarse_pc_type lu -fas_coarse_ksp_type preonly -snes_type fas -snes_rtol 1.e-5
362: TEST*/