1 static char help[] = "Poisson problem in 2d with finite elements and bound constraints.\n\ 2 This example is intended to test VI solvers.\n\n\n"; 3 4 /* 5 This is the ball obstacle problem, taken from the MS thesis ``Adaptive Mesh Refinement for Variational Inequalities'' 6 by Stefano Fochesatto, University of Alaska Fairbanks, 2025 7 This is the same VI problem as in src/snes/tutorials/ex9.c, which uses DMDA. The example 8 is also documented by Chapter 12 of E. Bueler, "PETSc for Partial Differential Equations", 9 SIAM Press 2021. 10 11 To visualize the solution, configure with petsc4py, pip install pyvista, and use 12 13 -potential_view pyvista -view_pyvista_warp 1. 14 15 To look at the error use 16 17 -snes_convergence_estimate -convest_num_refine 2 -convest_monitor -convest_error_view pyvista 18 19 and for the inactive residual and active set use 20 21 -snes_vi_monitor_residual pyvista -snes_vi_monitor_active pyvista 22 23 To see the convergence history use 24 25 -snes_vi_monitor -snes_converged_reason -convest_monitor 26 */ 27 28 #include <petscdmplex.h> 29 #include <petscsnes.h> 30 #include <petscds.h> 31 #include <petscbag.h> 32 #include <petscconvest.h> 33 34 typedef enum { 35 OBSTACLE_NONE, 36 OBSTACLE_BALL, 37 NUM_OBSTACLE_TYPES 38 } ObstacleType; 39 const char *obstacleTypes[NUM_OBSTACLE_TYPES + 1] = {"none", "ball", "unknown"}; 40 41 typedef struct { 42 PetscReal r_0; // Ball radius 43 PetscReal r_free; // Radius of the free boundary for the ball obstacle 44 PetscReal A; // Logarithmic coefficient in exact ball solution 45 PetscReal B; // Constant coefficient in exact ball solution 46 } Parameter; 47 48 typedef struct { 49 // Problem definition 50 ObstacleType obsType; // Type of obstacle 51 PetscBag bag; // Problem parameters 52 } AppCtx; 53 54 static PetscErrorCode obstacle_ball(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 55 { 56 Parameter *par = (Parameter *)ctx; 57 const PetscReal r_0 = par->r_0; 58 const PetscReal r = PetscSqrtReal(PetscSqr(x[0]) + PetscSqr(x[1])); 59 const PetscReal psi_0 = PetscSqrtReal(1. - PetscSqr(r_0)); 60 const PetscReal dpsi_0 = -r_0 / psi_0; 61 62 PetscFunctionBegin; 63 PetscCheck(dim == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Ball obstacle is only defined in 2D"); 64 if (r < r_0) u[0] = PetscSqrtReal(1.0 - PetscSqr(r)); 65 else u[0] = psi_0 + dpsi_0 * (r - r_0); 66 PetscFunctionReturn(PETSC_SUCCESS); 67 } 68 69 static PetscErrorCode exactSol_ball(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 70 { 71 Parameter *par = (Parameter *)ctx; 72 const PetscReal r_free = par->r_free; 73 const PetscReal A = par->A; 74 const PetscReal B = par->B; 75 const PetscReal r = PetscSqrtReal(PetscSqr(x[0]) + PetscSqr(x[1])); 76 77 PetscFunctionBegin; 78 PetscCheck(dim == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Ball obstacle is only defined in 2D"); 79 if (r < r_free) PetscCall(obstacle_ball(dim, time, x, Nc, u, ctx)); 80 else u[0] = -A * PetscLogReal(r) + B; 81 PetscFunctionReturn(PETSC_SUCCESS); 82 } 83 84 static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 85 { 86 for (PetscInt d = 0; d < dim; ++d) f1[d] = u_x[d]; 87 } 88 89 static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 90 { 91 for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = 1.0; 92 } 93 94 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 95 { 96 PetscInt obs = OBSTACLE_BALL; 97 98 options->obsType = OBSTACLE_BALL; 99 100 PetscFunctionBeginUser; 101 PetscOptionsBegin(comm, "", "Ball Obstacle Problem Options", "DMPLEX"); 102 PetscCall(PetscOptionsEList("-obs_type", "Type of obstacle", "ex34.c", obstacleTypes, NUM_OBSTACLE_TYPES, obstacleTypes[options->obsType], &obs, NULL)); 103 options->obsType = (ObstacleType)obs; 104 PetscOptionsEnd(); 105 PetscFunctionReturn(PETSC_SUCCESS); 106 } 107 108 static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx) 109 { 110 PetscBag bag; 111 Parameter *p; 112 113 PetscFunctionBeginUser; 114 /* setup PETSc parameter bag */ 115 PetscCall(PetscBagGetData(ctx->bag, (void **)&p)); 116 PetscCall(PetscBagSetName(ctx->bag, "par", "Obstacle Parameters")); 117 bag = ctx->bag; 118 PetscCall(PetscBagRegisterReal(bag, &p->r_0, 0.9, "r_0", "Ball radius, m")); 119 PetscCall(PetscBagRegisterReal(bag, &p->r_free, 0.697965148223374, "r_free", "Ball free boundary radius, m")); 120 PetscCall(PetscBagRegisterReal(bag, &p->A, 0.680259411891719, "A", "Logarithmic coefficient in exact ball solution")); 121 PetscCall(PetscBagRegisterReal(bag, &p->B, 0.471519893402112, "B", "Constant coefficient in exact ball solution")); 122 PetscCall(PetscBagSetFromOptions(bag)); 123 { 124 PetscViewer viewer; 125 PetscViewerFormat format; 126 PetscBool flg; 127 128 PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg)); 129 if (flg) { 130 PetscCall(PetscViewerPushFormat(viewer, format)); 131 PetscCall(PetscBagView(bag, viewer)); 132 PetscCall(PetscViewerFlush(viewer)); 133 PetscCall(PetscViewerPopFormat(viewer)); 134 PetscCall(PetscViewerDestroy(&viewer)); 135 } 136 } 137 PetscFunctionReturn(PETSC_SUCCESS); 138 } 139 140 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 141 { 142 PetscFunctionBeginUser; 143 PetscCall(DMCreate(comm, dm)); 144 PetscCall(DMSetType(*dm, DMPLEX)); 145 PetscCall(DMSetFromOptions(*dm)); 146 PetscCall(DMSetApplicationContext(*dm, user)); 147 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 148 PetscCall(DMGetCoordinatesLocalSetUp(*dm)); 149 PetscFunctionReturn(PETSC_SUCCESS); 150 } 151 152 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 153 { 154 PetscErrorCode (*exact)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 155 Parameter *param; 156 PetscDS ds; 157 PetscWeakForm wf; 158 DMLabel label; 159 PetscInt dim, id = 1; 160 void *ctx; 161 162 PetscFunctionBeginUser; 163 PetscCall(DMGetDS(dm, &ds)); 164 PetscCall(PetscDSGetWeakForm(ds, &wf)); 165 PetscCall(PetscDSGetSpatialDimension(ds, &dim)); 166 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 167 switch (user->obsType) { 168 case OBSTACLE_BALL: 169 PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u)); 170 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu)); 171 PetscCall(PetscDSSetLowerBound(ds, 0, obstacle_ball, param)); 172 exact = exactSol_ball; 173 break; 174 default: 175 SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid obstacle type: %s (%d)", obstacleTypes[PetscMin(user->obsType, NUM_OBSTACLE_TYPES)], user->obsType); 176 } 177 PetscCall(PetscBagGetData(user->bag, (void **)&ctx)); 178 PetscCall(PetscDSSetExactSolution(ds, 0, exact, ctx)); 179 PetscCall(DMGetLabel(dm, "marker", &label)); 180 if (label) PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))exact, NULL, ctx, NULL)); 181 /* Setup constants */ 182 { 183 PetscScalar constants[4]; 184 185 constants[0] = param->r_0; // Ball radius 186 constants[1] = param->r_free; // Radius of the free boundary for the ball obstacle 187 constants[2] = param->A; // Logarithmic coefficient in exact ball solution 188 constants[3] = param->B; // Constant coefficient in exact ball solution 189 PetscCall(PetscDSSetConstants(ds, 4, constants)); 190 } 191 PetscFunctionReturn(PETSC_SUCCESS); 192 } 193 194 PetscErrorCode SetupFE(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), void *ctx) 195 { 196 AppCtx *user = (AppCtx *)ctx; 197 DM cdm = dm; 198 PetscFE fe; 199 char prefix[PETSC_MAX_PATH_LEN]; 200 DMPolytopeType ct; 201 PetscInt dim, cStart; 202 203 PetscFunctionBegin; 204 /* Create finite element */ 205 PetscCall(DMGetDimension(dm, &dim)); 206 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL)); 207 PetscCall(DMPlexGetCellType(dm, cStart, &ct)); 208 PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name)); 209 PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, name ? prefix : NULL, -1, &fe)); 210 PetscCall(PetscObjectSetName((PetscObject)fe, name)); 211 // Set discretization and boundary conditions for each mesh 212 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 213 PetscCall(DMCreateDS(dm)); 214 PetscCall((*setup)(dm, user)); 215 while (cdm) { 216 PetscCall(DMCopyDisc(dm, cdm)); 217 PetscCall(DMGetCoarseDM(cdm, &cdm)); 218 } 219 PetscCall(PetscFEDestroy(&fe)); 220 PetscFunctionReturn(PETSC_SUCCESS); 221 } 222 223 int main(int argc, char **argv) 224 { 225 DM dm; /* Problem specification */ 226 SNES snes; /* Nonlinear solver */ 227 Vec u; /* Solutions */ 228 AppCtx user; /* User-defined work context */ 229 230 PetscFunctionBeginUser; 231 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 232 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 233 PetscCall(PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &user.bag)); 234 PetscCall(SetupParameters(PETSC_COMM_WORLD, &user)); 235 /* Primal system */ 236 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 237 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 238 PetscCall(SNESSetDM(snes, dm)); 239 PetscCall(SetupFE(dm, "potential", SetupPrimalProblem, &user)); 240 PetscCall(DMCreateGlobalVector(dm, &u)); 241 PetscCall(VecSet(u, 0.0)); 242 PetscCall(PetscObjectSetName((PetscObject)u, "potential")); 243 PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user)); 244 PetscCall(DMPlexSetSNESVariableBounds(dm, snes)); 245 246 PetscCall(SNESSetFromOptions(snes)); 247 PetscCall(DMSNESCheckFromOptions(snes, u)); 248 PetscCall(SNESSolve(snes, NULL, u)); 249 PetscCall(SNESGetSolution(snes, &u)); 250 PetscCall(VecViewFromOptions(u, NULL, "-potential_view")); 251 /* Cleanup */ 252 PetscCall(VecDestroy(&u)); 253 PetscCall(SNESDestroy(&snes)); 254 PetscCall(DMDestroy(&dm)); 255 PetscCall(PetscBagDestroy(&user.bag)); 256 PetscCall(PetscFinalize()); 257 return 0; 258 } 259 260 /*TEST 261 262 testset: 263 args: -dm_plex_box_lower -2.,-2. -dm_plex_box_upper 2.,2. -dm_plex_box_faces 20,20 \ 264 -potential_petscspace_degree 1 \ 265 -snes_type vinewtonrsls -snes_vi_zero_tolerance 1.0e-12 \ 266 -ksp_type preonly -pc_type lu 267 268 # Check the exact solution 269 test: 270 suffix: ball_0 271 requires: triangle 272 args: -dmsnes_check 273 274 # Check convergence 275 test: 276 suffix: ball_1 277 requires: triangle 278 args: -snes_convergence_estimate -convest_num_refine 2 279 280 # Check different size obstacle 281 test: 282 suffix: ball_2 283 requires: triangle 284 args: -r_0 0.4 285 286 # Check quadrilateral mesh 287 test: 288 suffix: ball_3 289 args: -dm_plex_simplex 0 -snes_convergence_estimate -convest_num_refine 2 290 291 TEST*/ 292