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*/
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 newtonls -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 <petscdm.h>
58: #include <petscdmda.h>
59: #include <petscsnes.h>
61: /*
62: User-defined routines
63: */
64: extern PetscErrorCode FormMatrix(DM,Mat);
65: extern PetscErrorCode MyComputeFunction(SNES,Vec,Vec,void*);
66: extern PetscErrorCode MyComputeJacobian(SNES,Vec,Mat,Mat,void*);
67: extern PetscErrorCode NonlinearGS(SNES,Vec);
69: int main(int argc,char **argv)
70: {
71: SNES snes; /* nonlinear solver */
72: SNES psnes; /* nonlinear Gauss-Seidel approximate solver */
73: Vec x,b; /* solution vector */
74: PetscInt its; /* iterations for convergence */
75: DM da;
76: PetscBool use_ngs_as_npc = PETSC_FALSE; /* use the nonlinear Gauss-Seidel approximate solver */
78: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: Initialize program
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
82: PetscInitialize(&argc,&argv,(char*)0,help);
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Create nonlinear solver context
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87: SNESCreate(PETSC_COMM_WORLD,&snes);
89: PetscOptionsGetBool(NULL,NULL,"-use_ngs_as_npc",&use_ngs_as_npc,0);
91: if (use_ngs_as_npc) {
92: SNESGetNPC(snes,&psnes);
93: SNESSetType(psnes,SNESSHELL);
94: SNESShellSetSolve(psnes,NonlinearGS);
95: }
97: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98: Create distributed array (DMDA) to manage parallel grid and vectors
99: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
101: DMSetFromOptions(da);
102: DMSetUp(da);
103: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
104: SNESSetDM(snes,da);
105: if (use_ngs_as_npc) {
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: VecSet(b,1.0);
116: SNESSetFunction(snes,NULL,MyComputeFunction,NULL);
117: SNESSetJacobian(snes,NULL,NULL,MyComputeJacobian,NULL);
119: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: Customize nonlinear solver; set runtime options
121: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122: SNESSetFromOptions(snes);
124: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: Solve nonlinear system
126: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
127: SNESSolve(snes,b,x);
128: SNESGetIterationNumber(snes,&its);
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: Free work space. All PETSc objects should be destroyed when they
135: are no longer needed.
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137: VecDestroy(&x);
138: VecDestroy(&b);
139: SNESDestroy(&snes);
140: DMDestroy(&da);
141: PetscFinalize();
142: return 0;
143: }
145: /* ------------------------------------------------------------------- */
146: PetscErrorCode MyComputeFunction(SNES snes,Vec x,Vec F,void *ctx)
147: {
148: Mat J;
149: DM dm;
152: SNESGetDM(snes,&dm);
153: DMGetApplicationContext(dm,&J);
154: if (!J) {
155: DMSetMatType(dm,MATAIJ);
156: DMCreateMatrix(dm,&J);
157: MatSetDM(J, NULL);
158: FormMatrix(dm,J);
159: DMSetApplicationContext(dm,J);
160: DMSetApplicationContextDestroy(dm,(PetscErrorCode (*)(void**))MatDestroy);
161: }
162: MatMult(J,x,F);
163: return 0;
164: }
166: PetscErrorCode MyComputeJacobian(SNES snes,Vec x,Mat J,Mat Jp,void *ctx)
167: {
168: DM dm;
171: SNESGetDM(snes,&dm);
172: FormMatrix(dm,Jp);
173: return 0;
174: }
176: PetscErrorCode FormMatrix(DM da,Mat jac)
177: {
178: PetscInt i,j,nrows = 0;
179: MatStencil col[5],row,*rows;
180: PetscScalar v[5],hx,hy,hxdhy,hydhx;
181: DMDALocalInfo info;
184: DMDAGetLocalInfo(da,&info);
185: hx = 1.0/(PetscReal)(info.mx-1);
186: hy = 1.0/(PetscReal)(info.my-1);
187: hxdhy = hx/hy;
188: hydhx = hy/hx;
190: PetscMalloc1(info.ym*info.xm,&rows);
191: /*
192: Compute entries for the locally owned part of the Jacobian.
193: - Currently, all PETSc parallel matrix formats are partitioned by
194: contiguous chunks of rows across the processors.
195: - Each processor needs to insert only elements that it owns
196: locally (but any non-local elements will be sent to the
197: appropriate processor during matrix assembly).
198: - Here, we set all entries for a particular row at once.
199: - We can set matrix entries either using either
200: MatSetValuesLocal() or MatSetValues(), as discussed above.
201: */
202: for (j=info.ys; j<info.ys+info.ym; j++) {
203: for (i=info.xs; i<info.xs+info.xm; i++) {
204: row.j = j; row.i = i;
205: /* boundary points */
206: if (i == 0 || j == 0 || i == info.mx-1 || j == info.my-1) {
207: v[0] = 2.0*(hydhx + hxdhy);
208: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
209: rows[nrows].i = i;
210: rows[nrows++].j = j;
211: } else {
212: /* interior grid points */
213: v[0] = -hxdhy; col[0].j = j - 1; col[0].i = i;
214: v[1] = -hydhx; col[1].j = j; col[1].i = i-1;
215: v[2] = 2.0*(hydhx + hxdhy); col[2].j = row.j; col[2].i = row.i;
216: v[3] = -hydhx; col[3].j = j; col[3].i = i+1;
217: v[4] = -hxdhy; col[4].j = j + 1; col[4].i = i;
218: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
219: }
220: }
221: }
223: /*
224: Assemble matrix, using the 2-step process:
225: MatAssemblyBegin(), MatAssemblyEnd().
226: */
227: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
228: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
229: MatZeroRowsColumnsStencil(jac,nrows,rows,2.0*(hydhx + hxdhy),NULL,NULL);
230: PetscFree(rows);
231: /*
232: Tell the matrix we will never add a new nonzero location to the
233: matrix. If we do, it will generate an error.
234: */
235: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
236: return 0;
237: }
239: /* ------------------------------------------------------------------- */
240: /*
241: Applies some sweeps on nonlinear Gauss-Seidel on each process
243: */
244: PetscErrorCode NonlinearGS(SNES snes,Vec X)
245: {
246: PetscInt i,j,Mx,My,xs,ys,xm,ym,its,l;
247: PetscReal hx,hy,hxdhy,hydhx;
248: PetscScalar **x,F,J,u,uxx,uyy;
249: DM da;
250: Vec localX;
253: SNESGetTolerances(snes,NULL,NULL,NULL,&its,NULL);
254: SNESShellGetContext(snes,&da);
256: 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);
258: hx = 1.0/(PetscReal)(Mx-1);
259: hy = 1.0/(PetscReal)(My-1);
260: hxdhy = hx/hy;
261: hydhx = hy/hx;
263: DMGetLocalVector(da,&localX);
265: for (l=0; l<its; l++) {
267: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
268: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
269: /*
270: Get a pointer to vector data.
271: - For default PETSc vectors, VecGetArray() returns a pointer to
272: the data array. Otherwise, the routine is implementation dependent.
273: - You MUST call VecRestoreArray() when you no longer need access to
274: the array.
275: */
276: DMDAVecGetArray(da,localX,&x);
278: /*
279: Get local grid boundaries (for 2-dimensional DMDA):
280: xs, ys - starting grid indices (no ghost points)
281: xm, ym - widths of local grid (no ghost points)
283: */
284: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
286: for (j=ys; j<ys+ym; j++) {
287: for (i=xs; i<xs+xm; i++) {
288: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
289: /* boundary conditions are all zero Dirichlet */
290: x[j][i] = 0.0;
291: } else {
292: u = x[j][i];
293: uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
294: uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
295: F = uxx + uyy;
296: J = 2.0*(hydhx + hxdhy);
297: u = u - F/J;
299: x[j][i] = u;
300: }
301: }
302: }
304: /*
305: Restore vector
306: */
307: DMDAVecRestoreArray(da,localX,&x);
308: DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
309: DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
310: }
311: DMRestoreLocalVector(da,&localX);
312: return 0;
313: }
315: /*TEST
317: test:
318: args: -snes_monitor_short -snes_type nrichardson
319: requires: !single
321: test:
322: suffix: 2
323: args: -snes_monitor_short -ksp_monitor_short -ksp_type richardson -pc_type none -ksp_richardson_self_scale
324: requires: !single
326: test:
327: suffix: 3
328: args: -snes_monitor_short -snes_type ngmres
330: test:
331: suffix: 4
332: args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none
334: test:
335: suffix: 5
336: args: -snes_monitor_short -snes_type ncg
338: test:
339: suffix: 6
340: args: -snes_monitor_short -ksp_type cg -ksp_monitor_short -pc_type none
342: test:
343: suffix: 7
344: 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
345: requires: !single
347: test:
348: suffix: 8
349: 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
351: test:
352: suffix: 9
353: args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none -snes_type newtontrdc
355: test:
356: suffix: 10
357: args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none -snes_type newtontrdc -snes_trdc_use_cauchy false
359: TEST*/