Actual source code: ex64.c
1: static char help[] = "Biot consolidation model discretized with finite elements,\n\
2: using a parallel unstructured mesh (DMPLEX) to represent the domain.\n\
3: We follow https://arxiv.org/pdf/1507.03199.\n\n\n";
5: #include <petscdmplex.h>
6: #include <petscsnes.h>
7: #include <petscds.h>
9: typedef struct {
10: PetscScalar mu;
11: PetscScalar lam;
12: PetscScalar alpha;
13: PetscScalar kappa;
14: } Parameter;
16: static void u_1(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 f[])
17: {
18: const PetscReal mu = PetscRealPart(constants[0]);
20: for (PetscInt c = 0; c < dim; ++c) {
21: for (PetscInt d = 0; d < dim; ++d) f[c * dim + d] = mu * (u_x[c * dim + d] + u_x[d * dim + c]);
22: f[c * dim + c] -= u[uOff[1]];
23: }
24: }
26: /* Jfunction_testtrial */
27: static void Ju_1_u1p0(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 J[])
28: {
29: for (PetscInt d = 0; d < dim; ++d) J[d * dim + d] = -1.0;
30: }
32: static void Ju_1_u1u1(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 J[])
33: {
34: const PetscReal mu = PetscRealPart(constants[0]);
35: const PetscInt Nc = uOff[1] - uOff[0];
37: for (PetscInt c = 0; c < Nc; ++c) {
38: for (PetscInt d = 0; d < dim; ++d) {
39: J[((c * Nc + c) * dim + d) * dim + d] += mu;
40: J[((c * Nc + d) * dim + d) * dim + c] += mu;
41: }
42: }
43: }
45: static void p_0(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 f[])
46: {
47: const PetscReal lam = PetscRealPart(constants[1]);
48: const PetscReal alpha = PetscRealPart(constants[2]);
50: f[0] = (-u[uOff[1]] + alpha * u[uOff[2]]) / lam;
51: for (PetscInt d = 0; d < dim; ++d) f[0] -= u_x[d * dim + d];
52: }
54: static void Jp_0_p0u1(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 J[])
55: {
56: for (PetscInt d = 0; d < dim; ++d) J[d * dim + d] = -1.0;
57: }
59: static void Jp_0_p0p0(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 J[])
60: {
61: const PetscReal lam = PetscRealPart(constants[1]);
63: J[0] = -1.0 / lam;
64: }
66: static void pJp_0_p0p0(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 J[])
67: {
68: const PetscReal mu = PetscRealPart(constants[0]);
70: J[0] = 1.0 / mu;
71: }
73: static void Jp_0_p0w0(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 J[])
74: {
75: const PetscReal lam = PetscRealPart(constants[1]);
76: const PetscReal alpha = PetscRealPart(constants[2]);
78: J[0] = alpha / lam;
79: }
81: static void w_0(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 f[])
82: {
83: const PetscReal lam = PetscRealPart(constants[1]);
84: const PetscReal alpha = PetscRealPart(constants[2]);
86: f[0] = alpha / lam * (-2 * alpha * u[uOff[2]] + u[uOff[1]]);
87: }
89: static void Jw_0_w0p0(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 J[])
90: {
91: const PetscReal lam = PetscRealPart(constants[1]);
92: const PetscReal alpha = PetscRealPart(constants[2]);
94: J[0] = alpha / lam;
95: }
97: static void Jw_0_w0w0(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 J[])
98: {
99: const PetscReal lam = PetscRealPart(constants[1]);
100: const PetscReal alpha = PetscRealPart(constants[2]);
102: J[0] = -2 * PetscSqr(alpha) / lam;
103: }
105: static void pJw_0_w0w0(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 J[])
106: {
107: const PetscReal lam = PetscRealPart(constants[1]);
108: const PetscReal alpha = PetscRealPart(constants[2]);
110: J[0] = 2 * PetscSqr(alpha) / lam;
111: }
113: static void w_1(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 f[])
114: {
115: const PetscReal kappa = PetscRealPart(constants[3]);
116: for (PetscInt d = 0; d < dim; ++d) f[d] = -kappa * u_x[uOff_x[2] + d];
117: }
119: static void Jw_1_w1w1(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 J[])
120: {
121: const PetscReal kappa = PetscRealPart(constants[3]);
123: for (PetscInt d = 0; d < dim; ++d) J[d * dim + d] = -kappa;
124: }
126: static void pJw_1_w1w1(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 J[])
127: {
128: const PetscReal kappa = PetscRealPart(constants[3]);
130: for (PetscInt d = 0; d < dim; ++d) J[d * dim + d] = kappa;
131: }
133: typedef struct {
134: PetscScalar E;
135: PetscScalar nu;
136: PetscScalar alpha;
137: PetscScalar kappa;
138: } AppCtx;
140: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *user)
141: {
142: PetscFunctionBeginUser;
143: user->E = 1.0;
144: user->nu = 0.3;
145: user->alpha = 0.5;
146: user->kappa = 1.0;
147: PetscOptionsBegin(comm, "", "Biot Problem Options", "DMPLEX");
148: PetscCall(PetscOptionsScalar("-E", "Young modulus", NULL, user->E, &user->E, NULL));
149: PetscCall(PetscOptionsScalar("-nu", "Poisson ratio", NULL, user->nu, &user->nu, NULL));
150: PetscCall(PetscOptionsScalar("-alpha", "Alpha", NULL, user->alpha, &user->alpha, NULL));
151: PetscCall(PetscOptionsScalar("-kappa", "kappa", NULL, user->kappa, &user->kappa, NULL));
152: PetscOptionsEnd();
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
157: {
158: PetscFunctionBeginUser;
159: PetscCall(DMCreate(comm, dm));
160: PetscCall(DMSetType(*dm, DMPLEX));
161: PetscCall(DMSetFromOptions(*dm));
162: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: static PetscErrorCode SetupEqn(DM dm, AppCtx *user)
167: {
168: PetscDS ds;
169: DMLabel label;
170: const PetscInt id = 1;
171: PetscScalar constants[4];
173: PetscFunctionBeginUser;
174: PetscCall(DMGetDS(dm, &ds));
175: PetscCall(PetscDSSetResidual(ds, 0, NULL, u_1));
176: PetscCall(PetscDSSetResidual(ds, 1, p_0, NULL));
177: PetscCall(PetscDSSetResidual(ds, 2, w_0, w_1));
178: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, Ju_1_u1u1));
179: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, Ju_1_u1p0, NULL));
180: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, Jp_0_p0u1, NULL, NULL));
181: PetscCall(PetscDSSetJacobian(ds, 1, 1, Jp_0_p0p0, NULL, NULL, NULL));
182: PetscCall(PetscDSSetJacobian(ds, 1, 2, Jp_0_p0w0, NULL, NULL, NULL));
183: PetscCall(PetscDSSetJacobian(ds, 2, 1, Jw_0_w0p0, NULL, NULL, NULL));
184: PetscCall(PetscDSSetJacobian(ds, 2, 2, Jw_0_w0w0, NULL, NULL, Jw_1_w1w1));
185: PetscCall(PetscDSSetJacobianPreconditioner(ds, 0, 0, NULL, NULL, NULL, Ju_1_u1u1));
186: PetscCall(PetscDSSetJacobianPreconditioner(ds, 0, 1, NULL, NULL, Ju_1_u1p0, NULL));
187: PetscCall(PetscDSSetJacobianPreconditioner(ds, 1, 0, NULL, Jp_0_p0u1, NULL, NULL));
188: PetscCall(PetscDSSetJacobianPreconditioner(ds, 1, 1, pJp_0_p0p0, NULL, NULL, NULL));
189: PetscCall(PetscDSSetJacobianPreconditioner(ds, 2, 2, pJw_0_w0w0, NULL, NULL, pJw_1_w1w1));
191: PetscCall(DMGetLabel(dm, "marker", &label));
192: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall-d", label, 1, &id, 0, 0, NULL, NULL, NULL, NULL, NULL));
193: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall-p", label, 1, &id, 2, 0, NULL, NULL, NULL, NULL, NULL));
195: constants[0] = user->E / (2.0 * (1.0 + user->nu));
196: constants[1] = user->nu * user->E / ((1.0 + user->nu) * (1.0 - 2.0 * user->nu));
197: constants[2] = user->alpha;
198: constants[3] = user->kappa;
199: PetscCall(PetscDSSetConstants(ds, 4, constants));
201: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "E = %g, nu = %g\n", (double)PetscRealPart(user->E), (double)PetscRealPart(user->nu)));
202: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "mu = %g, lambda = %g\n", (double)PetscRealPart(constants[0]), (double)PetscRealPart(constants[1])));
203: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "alpha = %g, kappa = %g\n", (double)PetscRealPart(constants[2]), (double)PetscRealPart(constants[3])));
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: static PetscErrorCode SetupProblem(DM dm, PetscErrorCode (*setupEqn)(DM, AppCtx *), AppCtx *user)
208: {
209: DM cdm = dm;
210: PetscQuadrature q = NULL;
211: PetscBool simplex;
212: PetscInt dim, Nf = 3, f, Nc[3];
213: const char *name[3] = {"displacement", "totalpressure", "pressure"};
214: const char *prefix[3] = {"displacement_", "totalpressure_", "pressure_"};
216: PetscFunctionBegin;
217: PetscCall(DMGetDimension(dm, &dim));
218: PetscCall(DMPlexIsSimplex(dm, &simplex));
219: Nc[0] = dim;
220: Nc[1] = 1;
221: Nc[2] = 1;
222: for (f = 0; f < Nf; ++f) {
223: PetscFE fe;
225: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc[f], simplex, prefix[f], -1, &fe));
226: PetscCall(PetscObjectSetName((PetscObject)fe, name[f]));
227: if (!q) PetscCall(PetscFEGetQuadrature(fe, &q));
228: PetscCall(PetscFESetQuadrature(fe, q));
229: PetscCall(DMSetField(dm, f, NULL, (PetscObject)fe));
230: PetscCall(PetscFEDestroy(&fe));
231: }
232: PetscCall(DMCreateDS(dm));
233: PetscCall((*setupEqn)(dm, user));
234: while (cdm) {
235: PetscCall(DMCopyDisc(dm, cdm));
236: PetscCall(DMGetCoarseDM(cdm, &cdm));
237: }
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: int main(int argc, char **argv)
242: {
243: SNES snes;
244: DM dm;
245: Vec u;
246: AppCtx user;
248: PetscFunctionBeginUser;
249: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
250: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
251: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
252: PetscCall(SetupProblem(dm, SetupEqn, &user));
253: PetscCall(DMPlexCreateClosureIndex(dm, NULL));
254: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
255: PetscCall(DMSetApplicationContext(dm, &user));
257: PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
258: PetscCall(SNESSetDM(snes, dm));
259: PetscCall(SNESSetType(snes, SNESKSPONLY));
260: PetscCall(SNESSetFromOptions(snes));
262: /* Solve with random rhs */
263: PetscCall(DMCreateGlobalVector(dm, &u));
264: PetscCall(VecSetRandom(u, NULL));
265: PetscCall(SNESSolve(snes, NULL, u));
267: PetscCall(VecDestroy(&u));
268: PetscCall(SNESDestroy(&snes));
269: PetscCall(DMDestroy(&dm));
270: PetscCall(PetscFinalize());
271: return 0;
272: }
274: /*TEST
276: test:
277: suffix: 2d_p2_p1_p2
278: args: -displacement_petscspace_degree 2 -totalpressure_petscspace_degree 1 -pressure_petscspace_degree 2 \
279: -dm_plex_box_faces 5,5 -dm_plex_separate_marker -dm_plex_simplex 0 \
280: -snes_error_if_not_converged -ksp_error_if_not_converged \
281: -ksp_type minres -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_totalpressure_pc_type jacobi -fieldsplit_displacement_pc_type gamg -fieldsplit_pressure_pc_type gamg
283: test:
284: nsize: 4
285: suffix: 2d_p2_p1_p2_fetidp
286: args: -displacement_petscspace_degree 2 -totalpressure_petscspace_degree 1 -pressure_petscspace_degree 2 \
287: -petscpartitioner_type simple -dm_mat_type is -dm_plex_box_faces 3,3 -dm_refine 2 -dm_plex_separate_marker -dm_plex_simplex 0 \
288: -snes_error_if_not_converged -ksp_error_if_not_converged \
289: -ksp_type fetidp -ksp_fetidp_saddlepoint -ksp_fetidp_pressure_field 1,2 -ksp_fetidp_pressure_schur -fetidp_ksp_type cg -fetidp_ksp_norm_type natural \
290: -fetidp_fieldsplit_lag_ksp_type preonly \
291: -fetidp_bddc_pc_bddc_detect_disconnected -fetidp_bddc_pc_bddc_corner_selection -fetidp_bddc_pc_bddc_symmetric -fetidp_bddc_pc_bddc_coarse_redundant_pc_type cholesky \
292: -fetidp_pc_discrete_harmonic 1 -fetidp_harmonic_pc_type cholesky \
293: -fetidp_fieldsplit_p_pc_type bddc -fetidp_fieldsplit_p_ksp_type preonly
295: TEST*/