1 static const char help[] = "Solves obstacle problem in 2D as a variational inequality\n\
2 or nonlinear complementarity problem. This is a form of the Laplace equation in\n\
3 which the solution u is constrained to be above a given function psi. In the\n\
4 problem here an exact solution is known.\n";
5
6 /* On a square S = {-2<x<2,-2<y<2}, the PDE
7 u_{xx} + u_{yy} = 0
8 is solved on the set where membrane is above obstacle (u(x,y) >= psi(x,y)).
9 Here psi is the upper hemisphere of the unit ball. On the boundary of S
10 we have Dirichlet boundary conditions from the exact solution. Uses centered
11 FD scheme. This example contributed by Ed Bueler.
12
13 Example usage:
14 * get help:
15 ./ex9 -help
16 * monitor run:
17 ./ex9 -da_refine 2 -snes_vi_monitor
18 * use other SNESVI type (default is SNESVINEWTONRSLS):
19 ./ex9 -da_refine 2 -snes_vi_monitor -snes_type vinewtonssls
20 * use FD evaluation of Jacobian by coloring, instead of analytical:
21 ./ex9 -da_refine 2 -snes_fd_color
22 * X windows visualizations:
23 ./ex9 -snes_monitor_solution draw -draw_pause 1 -da_refine 4
24 ./ex9 -snes_vi_monitor_residual -draw_pause 1 -da_refine 4
25 * serial convergence evidence:
26 for M in 3 4 5 6 7; do ./ex9 -snes_grid_sequence $M -pc_type mg; done
27 * parallel full-cycle multigrid from enlarged coarse mesh:
28 mpiexec -n 4 ./ex9 -da_grid_x 12 -da_grid_y 12 -snes_converged_reason -snes_grid_sequence 4 -pc_type mg
29 */
30
31 #include <petsc.h>
32
33 /* z = psi(x,y) is the hemispherical obstacle, but made C^1 with "skirt" at r=r0 */
psi(PetscReal x,PetscReal y)34 PetscReal psi(PetscReal x, PetscReal y)
35 {
36 const PetscReal r = x * x + y * y, r0 = 0.9, psi0 = PetscSqrtReal(1.0 - r0 * r0), dpsi0 = -r0 / psi0;
37 if (r <= r0) {
38 return PetscSqrtReal(1.0 - r);
39 } else {
40 return psi0 + dpsi0 * (r - r0);
41 }
42 }
43
44 /* This exact solution solves a 1D radial free-boundary problem for the
45 Laplace equation, on the interval 0 < r < 2, with above obstacle psi(x,y).
46 The Laplace equation applies where u(r) > psi(r),
47 u''(r) + r^-1 u'(r) = 0
48 with boundary conditions including free b.c.s at an unknown location r = a:
49 u(a) = psi(a), u'(a) = psi'(a), u(2) = 0
50 The solution is u(r) = - A log(r) + B on r > a. The boundary conditions
51 can then be reduced to a root-finding problem for a:
52 a^2 (log(2) - log(a)) = 1 - a^2
53 The solution is a = 0.697965148223374 (giving residual 1.5e-15). Then
54 A = a^2*(1-a^2)^(-0.5) and B = A*log(2) are as given below in the code. */
u_exact(PetscReal x,PetscReal y)55 PetscReal u_exact(PetscReal x, PetscReal y)
56 {
57 const PetscReal afree = 0.697965148223374, A = 0.680259411891719, B = 0.471519893402112;
58 PetscReal r;
59 r = PetscSqrtReal(x * x + y * y);
60 return (r <= afree) ? psi(x, y) /* active set; on the obstacle */
61 : -A * PetscLogReal(r) + B; /* solves laplace eqn */
62 }
63
64 extern PetscErrorCode FormExactSolution(DMDALocalInfo *, Vec);
65 extern PetscErrorCode FormBounds(SNES, Vec, Vec);
66 extern PetscErrorCode FormFunctionLocal(DMDALocalInfo *, PetscReal **, PetscReal **, void *);
67 extern PetscErrorCode FormJacobianLocal(DMDALocalInfo *, PetscReal **, Mat, Mat, void *);
68
main(int argc,char ** argv)69 int main(int argc, char **argv)
70 {
71 SNES snes;
72 DM da, da_after;
73 Vec u, u_exact;
74 DMDALocalInfo info;
75 PetscReal error1, errorinf;
76
77 PetscFunctionBeginUser;
78 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
79
80 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 5, 5, /* 5x5 coarse grid; override with -da_grid_x,_y */
81 PETSC_DECIDE, PETSC_DECIDE, 1, 1, /* dof=1 and s = 1 (stencil extends out one cell) */
82 NULL, NULL, &da));
83 PetscCall(DMSetFromOptions(da));
84 PetscCall(DMSetUp(da));
85 PetscCall(DMDASetUniformCoordinates(da, -2.0, 2.0, -2.0, 2.0, 0.0, 1.0));
86
87 PetscCall(DMCreateGlobalVector(da, &u));
88 PetscCall(VecSet(u, 0.0));
89
90 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
91 PetscCall(SNESSetDM(snes, da));
92 PetscCall(SNESSetType(snes, SNESVINEWTONRSLS));
93 PetscCall(SNESVISetComputeVariableBounds(snes, &FormBounds));
94 PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (DMDASNESFunctionFn *)FormFunctionLocal, NULL));
95 PetscCall(DMDASNESSetJacobianLocal(da, (DMDASNESJacobianFn *)FormJacobianLocal, NULL));
96 PetscCall(SNESSetFromOptions(snes));
97
98 /* solve nonlinear system */
99 PetscCall(SNESSolve(snes, NULL, u));
100 PetscCall(VecDestroy(&u));
101 PetscCall(DMDestroy(&da));
102 /* DMDA after solve may be different, e.g. with -snes_grid_sequence */
103 PetscCall(SNESGetDM(snes, &da_after));
104 PetscCall(SNESGetSolution(snes, &u)); /* do not destroy u */
105 PetscCall(DMDAGetLocalInfo(da_after, &info));
106 PetscCall(VecDuplicate(u, &u_exact));
107 PetscCall(FormExactSolution(&info, u_exact));
108 PetscCall(VecAXPY(u, -1.0, u_exact)); /* u <-- u - u_exact */
109 PetscCall(VecNorm(u, NORM_1, &error1));
110 error1 /= (PetscReal)info.mx * (PetscReal)info.my; /* average error */
111 PetscCall(VecNorm(u, NORM_INFINITY, &errorinf));
112 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "errors on %" PetscInt_FMT " x %" PetscInt_FMT " grid: av |u-uexact| = %.3e, |u-uexact|_inf = %.3e\n", info.mx, info.my, (double)error1, (double)errorinf));
113 PetscCall(VecDestroy(&u_exact));
114 PetscCall(SNESDestroy(&snes));
115 PetscCall(DMDestroy(&da));
116 PetscCall(PetscFinalize());
117 return 0;
118 }
119
FormExactSolution(DMDALocalInfo * info,Vec u)120 PetscErrorCode FormExactSolution(DMDALocalInfo *info, Vec u)
121 {
122 PetscInt i, j;
123 PetscReal **au, dx, dy, x, y;
124
125 PetscFunctionBeginUser;
126 dx = 4.0 / (PetscReal)(info->mx - 1);
127 dy = 4.0 / (PetscReal)(info->my - 1);
128 PetscCall(DMDAVecGetArray(info->da, u, &au));
129 for (j = info->ys; j < info->ys + info->ym; j++) {
130 y = -2.0 + j * dy;
131 for (i = info->xs; i < info->xs + info->xm; i++) {
132 x = -2.0 + i * dx;
133 au[j][i] = u_exact(x, y);
134 }
135 }
136 PetscCall(DMDAVecRestoreArray(info->da, u, &au));
137 PetscFunctionReturn(PETSC_SUCCESS);
138 }
139
FormBounds(SNES snes,Vec Xl,Vec Xu)140 PetscErrorCode FormBounds(SNES snes, Vec Xl, Vec Xu)
141 {
142 DM da;
143 DMDALocalInfo info;
144 PetscInt i, j;
145 PetscReal **aXl, dx, dy, x, y;
146
147 PetscFunctionBeginUser;
148 PetscCall(SNESGetDM(snes, &da));
149 PetscCall(DMDAGetLocalInfo(da, &info));
150 dx = 4.0 / (PetscReal)(info.mx - 1);
151 dy = 4.0 / (PetscReal)(info.my - 1);
152 PetscCall(DMDAVecGetArray(da, Xl, &aXl));
153 for (j = info.ys; j < info.ys + info.ym; j++) {
154 y = -2.0 + j * dy;
155 for (i = info.xs; i < info.xs + info.xm; i++) {
156 x = -2.0 + i * dx;
157 aXl[j][i] = psi(x, y);
158 }
159 }
160 PetscCall(DMDAVecRestoreArray(da, Xl, &aXl));
161 PetscCall(VecSet(Xu, PETSC_INFINITY));
162 PetscFunctionReturn(PETSC_SUCCESS);
163 }
164
FormFunctionLocal(DMDALocalInfo * info,PetscScalar ** au,PetscScalar ** af,void * user)165 PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, PetscScalar **au, PetscScalar **af, void *user)
166 {
167 PetscInt i, j;
168 PetscReal dx, dy, x, y, ue, un, us, uw;
169
170 PetscFunctionBeginUser;
171 dx = 4.0 / (PetscReal)(info->mx - 1);
172 dy = 4.0 / (PetscReal)(info->my - 1);
173 for (j = info->ys; j < info->ys + info->ym; j++) {
174 y = -2.0 + j * dy;
175 for (i = info->xs; i < info->xs + info->xm; i++) {
176 x = -2.0 + i * dx;
177 if (i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1) {
178 af[j][i] = 4.0 * (au[j][i] - u_exact(x, y));
179 } else {
180 uw = (i - 1 == 0) ? u_exact(x - dx, y) : au[j][i - 1];
181 ue = (i + 1 == info->mx - 1) ? u_exact(x + dx, y) : au[j][i + 1];
182 us = (j - 1 == 0) ? u_exact(x, y - dy) : au[j - 1][i];
183 un = (j + 1 == info->my - 1) ? u_exact(x, y + dy) : au[j + 1][i];
184 af[j][i] = -(dy / dx) * (uw - 2.0 * au[j][i] + ue) - (dx / dy) * (us - 2.0 * au[j][i] + un);
185 }
186 }
187 }
188 PetscCall(PetscLogFlops(12.0 * info->ym * info->xm));
189 PetscFunctionReturn(PETSC_SUCCESS);
190 }
191
FormJacobianLocal(DMDALocalInfo * info,PetscScalar ** au,Mat A,Mat jac,void * user)192 PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, PetscScalar **au, Mat A, Mat jac, void *user)
193 {
194 PetscInt i, j, n;
195 MatStencil col[5], row;
196 PetscReal v[5], dx, dy, oxx, oyy;
197
198 PetscFunctionBeginUser;
199 dx = 4.0 / (PetscReal)(info->mx - 1);
200 dy = 4.0 / (PetscReal)(info->my - 1);
201 oxx = dy / dx;
202 oyy = dx / dy;
203 for (j = info->ys; j < info->ys + info->ym; j++) {
204 for (i = info->xs; i < info->xs + info->xm; i++) {
205 row.j = j;
206 row.i = i;
207 if (i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1) { /* boundary */
208 v[0] = 4.0;
209 PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES));
210 } else { /* interior grid points */
211 v[0] = 2.0 * (oxx + oyy);
212 col[0].j = j;
213 col[0].i = i;
214 n = 1;
215 if (i - 1 > 0) {
216 v[n] = -oxx;
217 col[n].j = j;
218 col[n++].i = i - 1;
219 }
220 if (i + 1 < info->mx - 1) {
221 v[n] = -oxx;
222 col[n].j = j;
223 col[n++].i = i + 1;
224 }
225 if (j - 1 > 0) {
226 v[n] = -oyy;
227 col[n].j = j - 1;
228 col[n++].i = i;
229 }
230 if (j + 1 < info->my - 1) {
231 v[n] = -oyy;
232 col[n].j = j + 1;
233 col[n++].i = i;
234 }
235 PetscCall(MatSetValuesStencil(jac, 1, &row, n, col, v, INSERT_VALUES));
236 }
237 }
238 }
239
240 /* Assemble matrix, using the 2-step process: */
241 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
242 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
243 if (A != jac) {
244 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
245 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
246 }
247 PetscCall(PetscLogFlops(2.0 * info->ym * info->xm));
248 PetscFunctionReturn(PETSC_SUCCESS);
249 }
250
251 /*TEST
252
253 build:
254 requires: !complex
255
256 test:
257 suffix: 1
258 requires: !single
259 nsize: 1
260 args: -da_refine 1 -snes_monitor_short -snes_type vinewtonrsls
261
262 test:
263 suffix: 2
264 requires: !single
265 nsize: 2
266 args: -da_refine 1 -snes_monitor_short -snes_type vinewtonssls
267
268 test:
269 suffix: 3
270 requires: !single
271 nsize: 2
272 args: -snes_grid_sequence 2 -snes_vi_monitor -snes_type vinewtonrsls
273
274 test:
275 suffix: mg
276 requires: !single
277 nsize: 4
278 args: -snes_grid_sequence 3 -snes_converged_reason -pc_type mg
279
280 test:
281 suffix: 4
282 nsize: 1
283 args: -mat_is_symmetric
284
285 test:
286 suffix: 5
287 nsize: 1
288 args: -ksp_converged_reason -snes_fd_color
289
290 test:
291 suffix: 6
292 requires: !single
293 nsize: 2
294 args: -snes_grid_sequence 2 -pc_type mg -snes_monitor_short -ksp_converged_reason
295
296 test:
297 suffix: 7
298 nsize: 2
299 args: -da_refine 1 -snes_monitor_short -snes_type composite -snes_composite_type multiplicative -snes_composite_sneses vinewtonrsls,vinewtonssls -sub_0_snes_vi_monitor -sub_1_snes_vi_monitor
300 TODO: fix nasty memory leak in SNESCOMPOSITE
301
302 test:
303 suffix: 8
304 nsize: 2
305 args: -da_refine 1 -snes_monitor_short -snes_type composite -snes_composite_type additive -snes_composite_sneses vinewtonrsls -sub_0_snes_vi_monitor
306 TODO: fix nasty memory leak in SNESCOMPOSITE
307
308 test:
309 suffix: 9
310 nsize: 2
311 args: -da_refine 1 -snes_monitor_short -snes_type composite -snes_composite_type additiveoptimal -snes_composite_sneses vinewtonrsls -sub_0_snes_vi_monitor
312 TODO: fix nasty memory leak in SNESCOMPOSITE
313
314 TEST*/
315