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