Actual source code: ex9.c
petsc-3.8.4 2018-03-24
1: static const char help[] = "Solves obstacle problem in 2D as a variational inequality.\n\
2: An elliptic problem with solution u constrained to be above a given function psi. \n\
3: Exact solution is known.\n";
5: /* Solve on a square R = {-2<x<2,-2<y<2}:
6: u_{xx} + u_{yy} = 0
7: on the set where membrane is above obstacle. Constraint is u(x,y) >= psi(x,y).
8: Here psi is the upper hemisphere of the unit ball. On the boundary of R
9: we have nonhomogenous Dirichlet boundary conditions coming from the exact solution.
11: Method is centered finite differences.
13: This example was contributed by Ed Bueler. The exact solution is known for the
14: given psi and boundary values in question. See
15: https://github.com/bueler/fem-code-challenge/blob/master/obstacleDOC.pdf?raw=true.
17: Example usage follows.
19: Get help:
20: ./ex9 -help
22: Monitor run:
23: ./ex9 -snes_converged_reason -snes_monitor -snes_vi_monitor
25: Use finite difference evaluation of Jacobian by coloring, instead of analytical:
26: ./ex9 -snes_fd_color
28: Graphical:
29: ./ex9 -snes_monitor_solution draw -draw_pause 1 -da_refine 2
31: Convergence evidence:
32: for M in 1 2 3 4 5; do mpiexec -n 4 ./ex9 -da_refine $M; done
33: */
35: #include <petscdm.h>
36: #include <petscdmda.h>
37: #include <petscsnes.h>
39: /* application context for obstacle problem solver */
40: typedef struct {
41: Vec psi, uexact;
42: } ObsCtx;
44: extern PetscErrorCode FormPsiAndExactSoln(DM);
45: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,ObsCtx*);
46: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,Mat,ObsCtx*);
48: int main(int argc,char **argv)
49: {
50: PetscErrorCode ierr;
51: ObsCtx user;
52: SNES snes;
53: DM da;
54: Vec u, /* solution */
55: Xu; /* upper bound */
56: DMDALocalInfo info;
57: PetscReal error1,errorinf;
59: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
61: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,11,11,/* default to 10x10 grid */
62: PETSC_DECIDE,PETSC_DECIDE, /* number of processors in each dimension */1,/* dof = 1 */1,/* s = 1; stencil extends out one cell */
63: NULL,NULL,/* do not specify processor decomposition */&da);
64: DMSetFromOptions(da);
65: DMSetUp(da);
66: DMCreateGlobalVector(da,&u);
67: VecDuplicate(u,&(user.uexact));
68: VecDuplicate(u,&(user.psi));
70: DMDASetUniformCoordinates(da,-2.0,2.0,-2.0,2.0,0.0,1.0);
71: DMSetApplicationContext(da,&user);
73: FormPsiAndExactSoln(da);
74: VecSet(u,0.0);
76: SNESCreate(PETSC_COMM_WORLD,&snes);
77: SNESSetDM(snes,da);
78: SNESSetApplicationContext(snes,&user);
79: SNESSetType(snes,SNESVINEWTONRSLS);
81: /* set upper and lower bound constraints for VI */
82: VecDuplicate(u,&Xu);
83: VecSet(Xu,PETSC_INFINITY);
84: SNESVISetVariableBounds(snes,user.psi,Xu);
85: VecDestroy(&Xu);
87: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);
88: DMDASNESSetJacobianLocal(da,(PetscErrorCode (*)(DMDALocalInfo*,void*,Mat,Mat,void*))FormJacobianLocal,&user);
90: SNESSetFromOptions(snes);
92: /* report on setup */
93: DMDAGetLocalInfo(da,&info);
94: PetscPrintf(PETSC_COMM_WORLD,"setup done: grid Mx,My = %D,%D with spacing dx,dy = %.4f,%.4f\n",
95: info.mx,info.my,(double)(4.0/(PetscReal)(info.mx-1)),(double)(4.0/(PetscReal)(info.my-1)));
97: /* solve nonlinear system */
98: SNESSolve(snes,NULL,u);
100: /* compare to exact */
101: VecAXPY(u,-1.0,user.uexact); /* u <- u - uexact */
102: VecNorm(u,NORM_1,&error1);
103: error1 /= (PetscReal)info.mx * (PetscReal)info.my;
104: VecNorm(u,NORM_INFINITY,&errorinf);
105: PetscPrintf(PETSC_COMM_WORLD,"errors: av |u-uexact| = %.3e |u-uexact|_inf = %.3e\n",(double)error1,(double)errorinf);
107: /* Free work space. */
108: VecDestroy(&u);
109: VecDestroy(&(user.psi));
110: VecDestroy(&(user.uexact));
112: SNESDestroy(&snes);
113: DMDestroy(&da);
114: PetscFinalize();
115: return 0;
116: }
119: PetscErrorCode FormPsiAndExactSoln(DM da) {
120: ObsCtx *user;
122: DMDALocalInfo info;
123: PetscInt i,j;
124: DM coordDA;
125: Vec coordinates;
126: DMDACoor2d **coords;
127: PetscReal **psi, **uexact, r;
128: const PetscReal afree = 0.69797, A = 0.68026, B = 0.47152;
131: DMGetApplicationContext(da,&user);
132: DMDAGetLocalInfo(da,&info);
134: DMGetCoordinateDM(da, &coordDA);
135: DMGetCoordinates(da, &coordinates);
137: DMDAVecGetArray(coordDA, coordinates, &coords);
138: DMDAVecGetArray(da, user->psi, &psi);
139: DMDAVecGetArray(da, user->uexact, &uexact);
140: for (j=info.ys; j<info.ys+info.ym; j++) {
141: for (i=info.xs; i<info.xs+info.xm; i++) {
142: r = PetscSqrtReal(PetscPowScalarInt(coords[j][i].x,2) + PetscPowScalarInt(coords[j][i].y,2));
143: if (r <= 1.0) psi[j][i] = PetscSqrtReal(1.0 - r * r);
144: else psi[j][i] = -1.0;
145: if (r <= afree) uexact[j][i] = psi[j][i]; /* on the obstacle */
146: else uexact[j][i] = - A * PetscLogReal(r) + B; /* solves the laplace eqn */
147: }
148: }
149: DMDAVecRestoreArray(da, user->psi, &psi);
150: DMDAVecRestoreArray(da, user->uexact, &uexact);
151: DMDAVecRestoreArray(coordDA, coordinates, &coords);
152: return(0);
153: }
156: /* FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch */
157: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,ObsCtx *user) {
159: PetscInt i,j;
160: PetscReal dx,dy,uxx,uyy;
161: PetscReal **uexact; /* used for boundary values only */
164: dx = 4.0 / (PetscReal)(info->mx-1);
165: dy = 4.0 / (PetscReal)(info->my-1);
167: DMDAVecGetArray(info->da, user->uexact, &uexact);
168: for (j=info->ys; j<info->ys+info->ym; j++) {
169: for (i=info->xs; i<info->xs+info->xm; i++) {
170: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
171: f[j][i] = 4.0*(x[j][i] - uexact[j][i]);
172: } else {
173: uxx = dy*(x[j][i-1] - 2.0 * x[j][i] + x[j][i+1]) / dx;
174: uyy = dx*(x[j-1][i] - 2.0 * x[j][i] + x[j+1][i]) / dy;
175: f[j][i] = -uxx - uyy;
176: }
177: }
178: }
179: DMDAVecRestoreArray(info->da, user->uexact, &uexact);
181: PetscLogFlops(10.0*info->ym*info->xm);
182: return(0);
183: }
186: /* FormJacobianLocal - Evaluates Jacobian matrix on local process patch */
187: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat A,Mat jac, ObsCtx *user)
188: {
190: PetscInt i,j;
191: MatStencil col[5],row;
192: PetscReal v[5],dx,dy,oxx,oyy;
195: dx = 4.0 / (PetscReal)(info->mx-1);
196: dy = 4.0 / (PetscReal)(info->my-1);
197: oxx = dy / dx;
198: oyy = dx / dy;
200: for (j=info->ys; j<info->ys+info->ym; j++) {
201: for (i=info->xs; i<info->xs+info->xm; i++) {
202: row.j = j; row.i = i;
203: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) { /* boundary */
204: v[0] = 4.0;
205: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
206: } else { /* interior grid points */
207: v[0] = -oyy; col[0].j = j - 1; col[0].i = i;
208: v[1] = -oxx; col[1].j = j; col[1].i = i - 1;
209: v[2] = 2.0 * (oxx + oyy); col[2].j = j; col[2].i = i;
210: v[3] = -oxx; col[3].j = j; col[3].i = i + 1;
211: v[4] = -oyy; col[4].j = j + 1; col[4].i = i;
212: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
213: }
214: }
215: }
217: /* Assemble matrix, using the 2-step process: */
218: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
219: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
220: if (A != jac) {
221: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
222: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
223: }
224: PetscLogFlops(2.0*info->ym*info->xm);
225: return(0);
226: }