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
obstacle_ball(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)54 static PetscErrorCode obstacle_ball(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx 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
exactSol_ball(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)69 static PetscErrorCode exactSol_ball(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx 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
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[])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
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[])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
ProcessOptions(MPI_Comm comm,AppCtx * options)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
SetupParameters(MPI_Comm comm,AppCtx * ctx)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, &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
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)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
SetupPrimalProblem(DM dm,AppCtx * user)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, ¶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, &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, (PetscVoidFn *)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
SetupFE(DM dm,const char name[],PetscErrorCode (* setup)(DM,AppCtx *),PetscCtx ctx)194 PetscErrorCode SetupFE(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), PetscCtx 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
main(int argc,char ** argv)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 output_file: output/empty.out
286
287 # Check quadrilateral mesh
288 test:
289 suffix: ball_3
290 args: -dm_plex_simplex 0 -snes_convergence_estimate -convest_num_refine 2
291
292 TEST*/
293