Actual source code: ex35.c
1: static const char help[] = "-Laplacian u = b as a nonlinear problem.\n\n";
3: /*
5: The linear and nonlinear versions of these should give almost identical results on this problem
7: Richardson
8: Nonlinear:
9: -snes_rtol 1.e-12 -snes_monitor -snes_type nrichardson -snes_linesearch_monitor
11: Linear:
12: -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
14: GMRES
15: Nonlinear:
16: -snes_rtol 1.e-12 -snes_monitor -snes_type ngmres
18: Linear:
19: -snes_rtol 1.e-12 -snes_monitor -ksp_type gmres -ksp_monitor -ksp_rtol 1.e-12 -pc_type none
21: CG
22: Nonlinear:
23: -snes_rtol 1.e-12 -snes_monitor -snes_type ncg -snes_linesearch_monitor
25: Linear:
26: -snes_rtol 1.e-12 -snes_monitor -ksp_type cg -ksp_monitor -ksp_rtol 1.e-12 -pc_type none
28: Multigrid
29: Linear:
30: 1 level:
31: -snes_rtol 1.e-12 -snes_monitor -pc_type mg -mg_levels_ksp_type richardson -mg_levels_pc_type none -mg_levels_ksp_monitor
32: -mg_levels_ksp_richardson_self_scale -ksp_type richardson -ksp_monitor -ksp_rtol 1.e-12 -ksp_monitor_true_residual
34: n levels:
35: -da_refine n
37: Nonlinear:
38: 1 level:
39: -snes_rtol 1.e-12 -snes_monitor -snes_type fas -fas_levels_snes_monitor
41: n levels:
42: -da_refine n -fas_coarse_snes_type newtonls -fas_coarse_pc_type lu -fas_coarse_ksp_type preonly
44: */
46: /*
47: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
48: Include "petscsnes.h" so that we can use SNES solvers. Note that this
49: */
50: #include <petscdm.h>
51: #include <petscdmda.h>
52: #include <petscsnes.h>
54: /*
55: User-defined routines
56: */
57: extern PetscErrorCode FormMatrix(DM, Mat);
58: extern PetscErrorCode MyComputeFunction(SNES, Vec, Vec, void *);
59: extern PetscErrorCode MyComputeJacobian(SNES, Vec, Mat, Mat, void *);
60: extern PetscErrorCode NonlinearGS(SNES, Vec);
62: int main(int argc, char **argv)
63: {
64: SNES snes; /* nonlinear solver */
65: SNES psnes; /* nonlinear Gauss-Seidel approximate solver */
66: Vec x, b; /* solution vector */
67: PetscInt its; /* iterations for convergence */
68: DM da;
69: PetscBool use_ngs_as_npc = PETSC_FALSE; /* use the nonlinear Gauss-Seidel approximate solver */
71: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: Initialize program
73: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: PetscInitialize(&argc, &argv, (char *)0, help);
78: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: Create nonlinear solver context
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81: SNESCreate(PETSC_COMM_WORLD, &snes);
83: PetscOptionsGetBool(NULL, NULL, "-use_ngs_as_npc", &use_ngs_as_npc, 0);
85: if (use_ngs_as_npc) {
86: SNESGetNPC(snes, &psnes);
87: SNESSetType(psnes, SNESSHELL);
88: SNESShellSetSolve(psnes, NonlinearGS);
89: }
91: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: Create distributed array (DMDA) to manage parallel grid and vectors
93: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da);
95: DMSetFromOptions(da);
96: DMSetUp(da);
97: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
98: SNESSetDM(snes, da);
99: if (use_ngs_as_npc) SNESShellSetContext(psnes, da);
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Extract global vectors from DMDA; then duplicate for remaining
102: vectors that are the same types
103: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: DMCreateGlobalVector(da, &x);
105: DMCreateGlobalVector(da, &b);
106: VecSet(b, 1.0);
108: SNESSetFunction(snes, NULL, MyComputeFunction, NULL);
109: SNESSetJacobian(snes, NULL, NULL, MyComputeJacobian, NULL);
111: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: Customize nonlinear solver; set runtime options
113: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114: SNESSetFromOptions(snes);
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Solve nonlinear system
118: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: SNESSolve(snes, b, x);
120: SNESGetIterationNumber(snes, &its);
122: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Free work space. All PETSc objects should be destroyed when they
127: are no longer needed.
128: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129: VecDestroy(&x);
130: VecDestroy(&b);
131: SNESDestroy(&snes);
132: DMDestroy(&da);
133: PetscFinalize();
134: return 0;
135: }
137: /* ------------------------------------------------------------------- */
138: PetscErrorCode MyComputeFunction(SNES snes, Vec x, Vec F, void *ctx)
139: {
140: Mat J;
141: DM dm;
144: SNESGetDM(snes, &dm);
145: DMGetApplicationContext(dm, &J);
146: if (!J) {
147: DMSetMatType(dm, MATAIJ);
148: DMCreateMatrix(dm, &J);
149: MatSetDM(J, NULL);
150: FormMatrix(dm, J);
151: DMSetApplicationContext(dm, J);
152: DMSetApplicationContextDestroy(dm, (PetscErrorCode(*)(void **))MatDestroy);
153: }
154: MatMult(J, x, F);
155: return 0;
156: }
158: PetscErrorCode MyComputeJacobian(SNES snes, Vec x, Mat J, Mat Jp, void *ctx)
159: {
160: DM dm;
163: SNESGetDM(snes, &dm);
164: FormMatrix(dm, Jp);
165: return 0;
166: }
168: PetscErrorCode FormMatrix(DM da, Mat jac)
169: {
170: PetscInt i, j, nrows = 0;
171: MatStencil col[5], row, *rows;
172: PetscScalar v[5], hx, hy, hxdhy, hydhx;
173: DMDALocalInfo info;
176: DMDAGetLocalInfo(da, &info);
177: hx = 1.0 / (PetscReal)(info.mx - 1);
178: hy = 1.0 / (PetscReal)(info.my - 1);
179: hxdhy = hx / hy;
180: hydhx = hy / hx;
182: PetscMalloc1(info.ym * info.xm, &rows);
183: /*
184: Compute entries for the locally owned part of the Jacobian.
185: - Currently, all PETSc parallel matrix formats are partitioned by
186: contiguous chunks of rows across the processors.
187: - Each processor needs to insert only elements that it owns
188: locally (but any non-local elements will be sent to the
189: appropriate processor during matrix assembly).
190: - Here, we set all entries for a particular row at once.
191: - We can set matrix entries either using either
192: MatSetValuesLocal() or MatSetValues(), as discussed above.
193: */
194: for (j = info.ys; j < info.ys + info.ym; j++) {
195: for (i = info.xs; i < info.xs + info.xm; i++) {
196: row.j = j;
197: row.i = i;
198: /* boundary points */
199: if (i == 0 || j == 0 || i == info.mx - 1 || j == info.my - 1) {
200: v[0] = 2.0 * (hydhx + hxdhy);
201: MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES);
202: rows[nrows].i = i;
203: rows[nrows++].j = j;
204: } else {
205: /* interior grid points */
206: v[0] = -hxdhy;
207: col[0].j = j - 1;
208: col[0].i = i;
209: v[1] = -hydhx;
210: col[1].j = j;
211: col[1].i = i - 1;
212: v[2] = 2.0 * (hydhx + hxdhy);
213: col[2].j = row.j;
214: col[2].i = row.i;
215: v[3] = -hydhx;
216: col[3].j = j;
217: col[3].i = i + 1;
218: v[4] = -hxdhy;
219: col[4].j = j + 1;
220: col[4].i = i;
221: MatSetValuesStencil(jac, 1, &row, 5, col, v, INSERT_VALUES);
222: }
223: }
224: }
226: /*
227: Assemble matrix, using the 2-step process:
228: MatAssemblyBegin(), MatAssemblyEnd().
229: */
230: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
231: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
232: MatZeroRowsColumnsStencil(jac, nrows, rows, 2.0 * (hydhx + hxdhy), NULL, NULL);
233: PetscFree(rows);
234: /*
235: Tell the matrix we will never add a new nonzero location to the
236: matrix. If we do, it will generate an error.
237: */
238: MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
239: return 0;
240: }
242: /* ------------------------------------------------------------------- */
243: /*
244: Applies some sweeps on nonlinear Gauss-Seidel on each process
246: */
247: PetscErrorCode NonlinearGS(SNES snes, Vec X)
248: {
249: PetscInt i, j, Mx, My, xs, ys, xm, ym, its, l;
250: PetscReal hx, hy, hxdhy, hydhx;
251: PetscScalar **x, F, J, u, uxx, uyy;
252: DM da;
253: Vec localX;
256: SNESGetTolerances(snes, NULL, NULL, NULL, &its, NULL);
257: SNESShellGetContext(snes, &da);
259: 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);
261: hx = 1.0 / (PetscReal)(Mx - 1);
262: hy = 1.0 / (PetscReal)(My - 1);
263: hxdhy = hx / hy;
264: hydhx = hy / hx;
266: DMGetLocalVector(da, &localX);
268: for (l = 0; l < its; l++) {
269: DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
270: DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
271: /*
272: Get a pointer to vector data.
273: - For default PETSc vectors, VecGetArray() returns a pointer to
274: the data array. Otherwise, the routine is implementation dependent.
275: - You MUST call VecRestoreArray() when you no longer need access to
276: the array.
277: */
278: DMDAVecGetArray(da, localX, &x);
280: /*
281: Get local grid boundaries (for 2-dimensional DMDA):
282: xs, ys - starting grid indices (no ghost points)
283: xm, ym - widths of local grid (no ghost points)
285: */
286: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
288: for (j = ys; j < ys + ym; j++) {
289: for (i = xs; i < xs + xm; i++) {
290: if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) {
291: /* boundary conditions are all zero Dirichlet */
292: x[j][i] = 0.0;
293: } else {
294: u = x[j][i];
295: uxx = (2.0 * u - x[j][i - 1] - x[j][i + 1]) * hydhx;
296: uyy = (2.0 * u - x[j - 1][i] - x[j + 1][i]) * hxdhy;
297: F = uxx + uyy;
298: J = 2.0 * (hydhx + hxdhy);
299: u = u - F / J;
301: x[j][i] = u;
302: }
303: }
304: }
306: /*
307: Restore vector
308: */
309: DMDAVecRestoreArray(da, localX, &x);
310: DMLocalToGlobalBegin(da, localX, INSERT_VALUES, X);
311: DMLocalToGlobalEnd(da, localX, INSERT_VALUES, X);
312: }
313: DMRestoreLocalVector(da, &localX);
314: return 0;
315: }
317: /*TEST
319: test:
320: args: -snes_monitor_short -snes_type nrichardson
321: requires: !single
323: test:
324: suffix: 2
325: args: -snes_monitor_short -ksp_monitor_short -ksp_type richardson -pc_type none -ksp_richardson_self_scale
326: requires: !single
328: test:
329: suffix: 3
330: args: -snes_monitor_short -snes_type ngmres
332: test:
333: suffix: 4
334: args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none
336: test:
337: suffix: 5
338: args: -snes_monitor_short -snes_type ncg
340: test:
341: suffix: 6
342: args: -snes_monitor_short -ksp_type cg -ksp_monitor_short -pc_type none
344: test:
345: suffix: 7
346: 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
347: requires: !single
349: test:
350: suffix: 8
351: 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
353: test:
354: suffix: 9
355: args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none -snes_type newtontrdc
357: test:
358: suffix: 10
359: args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none -snes_type newtontrdc -snes_trdc_use_cauchy false
361: TEST*/