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