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