Actual source code: ex9.c

petsc-3.9.4 2018-09-11
Report Typos and Errors
  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: }



230: /*TEST

232:    build:
233:       requires: !complex !single

235:    test:
236:       nsize: 2
237:       args: -da_refine 1 -snes_monitor_short -snes_type vinewtonrsls

239:    test:
240:       suffix: 2
241:       nsize: 2
242:       args: -da_refine 1 -snes_monitor_short -snes_type vinewtonssls

244:    test:
245:       suffix: 3
246:       nsize: 2
247:       args: -da_refine 1 -snes_monitor_short -snes_type vinewtonrsls

249:    test:
250:       suffix: 4
251:       nsize: 2
252:       args: -da_refine 1 -snes_monitor_short -snes_type vinewtonssls

254:    test:
255:       suffix: 5
256:       args: -snes_fd_color

258:    test:
259:       suffix: 5_2
260:       nsize: 4
261:       args: -snes_fd_color
262:       output_file: output/ex9_5.out

264:    test:
265:       suffix: 6
266:       args: -da_refine 3 -pc_type mg -snes_monitor_short -ksp_monitor_short

268:    test:
269:       suffix: 7
270:       nsize: 2
271:       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
272:       TODO: fix nasty memory leak in SNESCOMPOSITE

274:    test:
275:       suffix: 8
276:       nsize: 2
277:       args: -da_refine 1 -snes_monitor_short -snes_type composite -snes_composite_type additive -snes_composite_sneses vinewtonrsls -sub_0_snes_vi_monitor
278:       TODO: fix nasty memory leak in SNESCOMPOSITE

280:    test:
281:       suffix: 9
282:       nsize: 2
283:       args: -da_refine 1 -snes_monitor_short -snes_type composite -snes_composite_type additiveoptimal -snes_composite_sneses vinewtonrsls -sub_0_snes_vi_monitor
284:       TODO: fix nasty memory leak in SNESCOMPOSITE

286: TEST*/