Actual source code: ex9.c
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";
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.
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: */
31: #include <petsc.h>
33: /* z = psi(x,y) is the hemispherical obstacle, but made C^1 with "skirt" at r=r0 */
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: }
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. */
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: }
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 *);
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;
77: PetscFunctionBeginUser;
78: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
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));
87: PetscCall(DMCreateGlobalVector(da, &u));
89: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
90: PetscCall(SNESSetDM(snes, da));
91: PetscCall(SNESSetType(snes, SNESVINEWTONRSLS));
92: PetscCall(SNESVISetComputeVariableBounds(snes, &FormBounds));
93: PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (DMDASNESFunctionFn *)FormFunctionLocal, NULL));
94: PetscCall(DMDASNESSetJacobianLocal(da, (DMDASNESJacobianFn *)FormJacobianLocal, NULL));
95: PetscCall(SNESSetFromOptions(snes));
97: /* solve nonlinear system */
98: PetscCall(SNESSolve(snes, NULL, u));
99: PetscCall(VecDestroy(&u));
100: PetscCall(DMDestroy(&da));
101: /* DMDA after solve may be different, e.g. with -snes_grid_sequence */
102: PetscCall(SNESGetDM(snes, &da_after));
103: PetscCall(SNESGetSolution(snes, &u)); /* do not destroy u */
104: PetscCall(DMDAGetLocalInfo(da_after, &info));
105: PetscCall(VecDuplicate(u, &u_exact));
106: PetscCall(FormExactSolution(&info, u_exact));
107: PetscCall(VecAXPY(u, -1.0, u_exact)); /* u <-- u - u_exact */
108: PetscCall(VecNorm(u, NORM_1, &error1));
109: error1 /= (PetscReal)info.mx * (PetscReal)info.my; /* average error */
110: PetscCall(VecNorm(u, NORM_INFINITY, &errorinf));
111: 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));
112: PetscCall(VecDestroy(&u_exact));
113: PetscCall(SNESDestroy(&snes));
114: PetscCall(DMDestroy(&da));
115: PetscCall(PetscFinalize());
116: return 0;
117: }
119: PetscErrorCode FormExactSolution(DMDALocalInfo *info, Vec u)
120: {
121: PetscInt i, j;
122: PetscReal **au, dx, dy, x, y;
124: PetscFunctionBeginUser;
125: dx = 4.0 / (PetscReal)(info->mx - 1);
126: dy = 4.0 / (PetscReal)(info->my - 1);
127: PetscCall(DMDAVecGetArray(info->da, u, &au));
128: for (j = info->ys; j < info->ys + info->ym; j++) {
129: y = -2.0 + j * dy;
130: for (i = info->xs; i < info->xs + info->xm; i++) {
131: x = -2.0 + i * dx;
132: au[j][i] = u_exact(x, y);
133: }
134: }
135: PetscCall(DMDAVecRestoreArray(info->da, u, &au));
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: PetscErrorCode FormBounds(SNES snes, Vec Xl, Vec Xu)
140: {
141: DM da;
142: DMDALocalInfo info;
143: PetscInt i, j;
144: PetscReal **aXl, dx, dy, x, y;
146: PetscFunctionBeginUser;
147: PetscCall(SNESGetDM(snes, &da));
148: PetscCall(DMDAGetLocalInfo(da, &info));
149: dx = 4.0 / (PetscReal)(info.mx - 1);
150: dy = 4.0 / (PetscReal)(info.my - 1);
151: PetscCall(DMDAVecGetArray(da, Xl, &aXl));
152: for (j = info.ys; j < info.ys + info.ym; j++) {
153: y = -2.0 + j * dy;
154: for (i = info.xs; i < info.xs + info.xm; i++) {
155: x = -2.0 + i * dx;
156: aXl[j][i] = psi(x, y);
157: }
158: }
159: PetscCall(DMDAVecRestoreArray(da, Xl, &aXl));
160: PetscCall(VecSet(Xu, PETSC_INFINITY));
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, PetscScalar **au, PetscScalar **af, void *user)
165: {
166: PetscInt i, j;
167: PetscReal dx, dy, x, y, ue, un, us, uw;
169: PetscFunctionBeginUser;
170: dx = 4.0 / (PetscReal)(info->mx - 1);
171: dy = 4.0 / (PetscReal)(info->my - 1);
172: for (j = info->ys; j < info->ys + info->ym; j++) {
173: y = -2.0 + j * dy;
174: for (i = info->xs; i < info->xs + info->xm; i++) {
175: x = -2.0 + i * dx;
176: if (i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1) {
177: af[j][i] = 4.0 * (au[j][i] - u_exact(x, y));
178: } else {
179: uw = (i - 1 == 0) ? u_exact(x - dx, y) : au[j][i - 1];
180: ue = (i + 1 == info->mx - 1) ? u_exact(x + dx, y) : au[j][i + 1];
181: us = (j - 1 == 0) ? u_exact(x, y - dy) : au[j - 1][i];
182: un = (j + 1 == info->my - 1) ? u_exact(x, y + dy) : au[j + 1][i];
183: af[j][i] = -(dy / dx) * (uw - 2.0 * au[j][i] + ue) - (dx / dy) * (us - 2.0 * au[j][i] + un);
184: }
185: }
186: }
187: PetscCall(PetscLogFlops(12.0 * info->ym * info->xm));
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, PetscScalar **au, Mat A, Mat jac, void *user)
192: {
193: PetscInt i, j, n;
194: MatStencil col[5], row;
195: PetscReal v[5], dx, dy, oxx, oyy;
197: PetscFunctionBeginUser;
198: dx = 4.0 / (PetscReal)(info->mx - 1);
199: dy = 4.0 / (PetscReal)(info->my - 1);
200: oxx = dy / dx;
201: oyy = dx / dy;
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;
205: row.i = i;
206: if (i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1) { /* boundary */
207: v[0] = 4.0;
208: PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES));
209: } else { /* interior grid points */
210: v[0] = 2.0 * (oxx + oyy);
211: col[0].j = j;
212: col[0].i = i;
213: n = 1;
214: if (i - 1 > 0) {
215: v[n] = -oxx;
216: col[n].j = j;
217: col[n++].i = i - 1;
218: }
219: if (i + 1 < info->mx - 1) {
220: v[n] = -oxx;
221: col[n].j = j;
222: col[n++].i = i + 1;
223: }
224: if (j - 1 > 0) {
225: v[n] = -oyy;
226: col[n].j = j - 1;
227: col[n++].i = i;
228: }
229: if (j + 1 < info->my - 1) {
230: v[n] = -oyy;
231: col[n].j = j + 1;
232: col[n++].i = i;
233: }
234: PetscCall(MatSetValuesStencil(jac, 1, &row, n, col, v, INSERT_VALUES));
235: }
236: }
237: }
239: /* Assemble matrix, using the 2-step process: */
240: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
241: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
242: if (A != jac) {
243: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
244: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
245: }
246: PetscCall(PetscLogFlops(2.0 * info->ym * info->xm));
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: /*TEST
252: build:
253: requires: !complex
255: test:
256: suffix: 1
257: requires: !single
258: nsize: 1
259: args: -da_refine 1 -snes_monitor_short -snes_type vinewtonrsls
261: test:
262: suffix: 2
263: requires: !single
264: nsize: 2
265: args: -da_refine 1 -snes_monitor_short -snes_type vinewtonssls
267: test:
268: suffix: 3
269: requires: !single
270: nsize: 2
271: args: -snes_grid_sequence 2 -snes_vi_monitor -snes_type vinewtonrsls
273: test:
274: suffix: mg
275: requires: !single
276: nsize: 4
277: args: -snes_grid_sequence 3 -snes_converged_reason -pc_type mg
279: test:
280: suffix: 4
281: nsize: 1
282: args: -mat_is_symmetric
284: test:
285: suffix: 5
286: nsize: 1
287: args: -ksp_converged_reason -snes_fd_color
289: test:
290: suffix: 6
291: requires: !single
292: nsize: 2
293: args: -snes_grid_sequence 2 -pc_type mg -snes_monitor_short -ksp_converged_reason
295: test:
296: suffix: 7
297: nsize: 2
298: 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
299: TODO: fix nasty memory leak in SNESCOMPOSITE
301: test:
302: suffix: 8
303: nsize: 2
304: args: -da_refine 1 -snes_monitor_short -snes_type composite -snes_composite_type additive -snes_composite_sneses vinewtonrsls -sub_0_snes_vi_monitor
305: TODO: fix nasty memory leak in SNESCOMPOSITE
307: test:
308: suffix: 9
309: nsize: 2
310: args: -da_refine 1 -snes_monitor_short -snes_type composite -snes_composite_type additiveoptimal -snes_composite_sneses vinewtonrsls -sub_0_snes_vi_monitor
311: TODO: fix nasty memory leak in SNESCOMPOSITE
313: TEST*/