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 */ 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. */ 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 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 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 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 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 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