Actual source code: ex34.c
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";
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.
11: To visualize the solution, configure with petsc4py, pip install pyvista, and use
13: -potential_view pyvista -view_pyvista_warp 1.
15: To look at the error use
17: -snes_convergence_estimate -convest_num_refine 2 -convest_monitor -convest_error_view pyvista
19: and for the inactive residual and active set use
21: -snes_vi_monitor_residual pyvista -snes_vi_monitor_active pyvista
23: To see the convergence history use
25: -snes_vi_monitor -snes_converged_reason -convest_monitor
26: */
28: #include <petscdmplex.h>
29: #include <petscsnes.h>
30: #include <petscds.h>
31: #include <petscbag.h>
32: #include <petscconvest.h>
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"};
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;
48: typedef struct {
49: // Problem definition
50: ObstacleType obsType; // Type of obstacle
51: PetscBag bag; // Problem parameters
52: } AppCtx;
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;
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: }
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]));
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: }
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: }
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: }
94: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
95: {
96: PetscInt obs = OBSTACLE_BALL;
98: options->obsType = OBSTACLE_BALL;
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: }
108: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
109: {
110: PetscBag bag;
111: Parameter *p;
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;
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: }
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: }
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;
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];
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: }
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;
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: }
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 */
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));
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: }
260: /*TEST
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
268: # Check the exact solution
269: test:
270: suffix: ball_0
271: requires: triangle
272: args: -dmsnes_check
274: # Check convergence
275: test:
276: suffix: ball_1
277: requires: triangle
278: args: -snes_convergence_estimate -convest_num_refine 2
280: # Check different size obstacle
281: test:
282: suffix: ball_2
283: requires: triangle
284: args: -r_0 0.4
286: # Check quadrilateral mesh
287: test:
288: suffix: ball_3
289: args: -dm_plex_simplex 0 -snes_convergence_estimate -convest_num_refine 2
291: TEST*/