xref: /petsc/src/snes/tutorials/ex64.c (revision bfe80ac4a46d58cb7760074b25f5e81b2f541d8a) 
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";
4 
5 #include <petscdmplex.h>
6 #include <petscsnes.h>
7 #include <petscds.h>
8 
9 typedef struct {
10   PetscScalar mu;
11   PetscScalar lam;
12   PetscScalar alpha;
13   PetscScalar kappa;
14 } Parameter;
15 
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]);
19 
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 }
25 
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 }
31 
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];
36 
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 }
44 
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]);
49 
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 }
53 
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 }
58 
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]);
62 
63   J[0] = -1.0 / lam;
64 }
65 
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]);
69 
70   J[0] = 1.0 / mu;
71 }
72 
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]);
77 
78   J[0] = alpha / lam;
79 }
80 
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]);
85 
86   f[0] = alpha / lam * (-2 * alpha * u[uOff[2]] + u[uOff[1]]);
87 }
88 
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]);
93 
94   J[0] = alpha / lam;
95 }
96 
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]);
101 
102   J[0] = -2 * PetscSqr(alpha) / lam;
103 }
104 
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]);
109 
110   J[0] = 2 * PetscSqr(alpha) / lam;
111 }
112 
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 }
118 
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]);
122 
123   for (PetscInt d = 0; d < dim; ++d) J[d * dim + d] = -kappa;
124 }
125 
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]);
129 
130   for (PetscInt d = 0; d < dim; ++d) J[d * dim + d] = kappa;
131 }
132 
133 typedef struct {
134   PetscScalar E;
135   PetscScalar nu;
136   PetscScalar alpha;
137   PetscScalar kappa;
138 } AppCtx;
139 
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 }
155 
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 }
165 
166 static PetscErrorCode SetupEqn(DM dm, AppCtx *user)
167 {
168   PetscDS        ds;
169   DMLabel        label;
170   const PetscInt id = 1;
171   PetscScalar    constants[4];
172 
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));
190 
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));
194 
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));
200 
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 }
206 
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_"};
215 
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;
224 
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 }
240 
241 int main(int argc, char **argv)
242 {
243   SNES   snes;
244   DM     dm;
245   Vec    u;
246   AppCtx user;
247 
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));
256 
257   PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
258   PetscCall(SNESSetDM(snes, dm));
259   PetscCall(SNESSetType(snes, SNESKSPONLY));
260   PetscCall(SNESSetFromOptions(snes));
261 
262   /* Solve with random rhs */
263   PetscCall(DMCreateGlobalVector(dm, &u));
264   PetscCall(VecSetRandom(u, NULL));
265   PetscCall(SNESSolve(snes, NULL, u));
266 
267   PetscCall(VecDestroy(&u));
268   PetscCall(SNESDestroy(&snes));
269   PetscCall(DMDestroy(&dm));
270   PetscCall(PetscFinalize());
271   return 0;
272 }
273 
274 /*TEST
275 
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
282 
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
294 
295 TEST*/
296