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: * 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: */
33: #include <petsc.h>
35: /* z = psi(x,y) is the hemispherical obstacle, but made C^1 with "skirt" at r=r0 */
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: }
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. */
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: }
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*);
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;
82: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
84: 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);
89: DMSetFromOptions(da);
90: DMSetUp(da);
91: DMDASetUniformCoordinates(da,-2.0,2.0,-2.0,2.0,0.0,1.0);
93: DMCreateGlobalVector(da,&u);
94: VecSet(u,0.0);
96: SNESCreate(PETSC_COMM_WORLD,&snes);
97: SNESSetDM(snes,da);
98: SNESSetType(snes,SNESVINEWTONRSLS);
99: SNESVISetComputeVariableBounds(snes,&FormBounds);
100: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocal,NULL);
101: DMDASNESSetJacobianLocal(da,(DMDASNESJacobian)FormJacobianLocal,NULL);
102: SNESSetFromOptions(snes);
104: /* solve nonlinear system */
105: SNESSolve(snes,NULL,u);
106: VecDestroy(&u);
107: DMDestroy(&da);
108: /* DMDA after solve may be different, e.g. with -snes_grid_sequence */
109: SNESGetDM(snes,&da_after);
110: SNESGetSolution(snes,&u); /* do not destroy u */
111: DMDAGetLocalInfo(da_after,&info);
112: VecDuplicate(u,&u_exact);
113: FormExactSolution(&info,u_exact);
114: VecAXPY(u,-1.0,u_exact); /* u <-- u - u_exact */
115: VecNorm(u,NORM_1,&error1);
116: error1 /= (PetscReal)info.mx * (PetscReal)info.my; /* average error */
117: VecNorm(u,NORM_INFINITY,&errorinf);
118: 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);
119: VecDestroy(&u_exact);
120: SNESDestroy(&snes);
121: DMDestroy(&da);
122: PetscFinalize();
123: return ierr;
124: }
126: PetscErrorCode FormExactSolution(DMDALocalInfo *info, Vec u)
127: {
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: DMDAVecGetArray(info->da, u, &au);
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: DMDAVecRestoreArray(info->da, u, &au);
142: return 0;
143: }
145: PetscErrorCode FormBounds(SNES snes, Vec Xl, Vec Xu)
146: {
148: DM da;
149: DMDALocalInfo info;
150: PetscInt i, j;
151: PetscReal **aXl, dx, dy, x, y;
153: SNESGetDM(snes,&da);
154: DMDAGetLocalInfo(da,&info);
155: dx = 4.0 / (PetscReal)(info.mx-1);
156: dy = 4.0 / (PetscReal)(info.my-1);
157: DMDAVecGetArray(da, Xl, &aXl);
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: DMDAVecRestoreArray(da, Xl, &aXl);
166: VecSet(Xu,PETSC_INFINITY);
167: return 0;
168: }
170: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, PetscScalar **au, PetscScalar **af, void *user)
171: {
173: PetscInt i,j;
174: PetscReal dx,dy,x,y,ue,un,us,uw;
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: PetscLogFlops(12.0*info->ym*info->xm);
195: return(0);
196: }
198: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, PetscScalar **au, Mat A, Mat jac, void *user)
199: {
201: PetscInt i,j,n;
202: MatStencil col[5],row;
203: PetscReal v[5],dx,dy,oxx,oyy;
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: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
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: MatSetValuesStencil(jac,1,&row,n,col,v,INSERT_VALUES);
232: }
233: }
234: }
236: /* Assemble matrix, using the 2-step process: */
237: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
238: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
239: if (A != jac) {
240: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
241: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
242: }
243: PetscLogFlops(2.0*info->ym*info->xm);
244: return(0);
245: }
247: /*TEST
249: build:
250: requires: !complex
252: test:
253: suffix: 1
254: requires: !single
255: nsize: 1
256: args: -da_refine 1 -snes_monitor_short -snes_type vinewtonrsls
258: test:
259: suffix: 2
260: requires: !single
261: nsize: 2
262: args: -da_refine 1 -snes_monitor_short -snes_type vinewtonssls
264: test:
265: suffix: 3
266: requires: !single
267: nsize: 2
268: args: -snes_grid_sequence 2 -snes_vi_monitor -snes_type vinewtonrsls
270: test:
271: suffix: mg
272: requires: !single
273: nsize: 4
274: args: -snes_grid_sequence 3 -snes_converged_reason -pc_type mg
276: test:
277: suffix: 4
278: nsize: 1
279: args: -mat_is_symmetric
281: test:
282: suffix: 5
283: nsize: 1
284: args: -ksp_converged_reason -snes_fd_color
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
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
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
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
310: TEST*/