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