1 static char help[] = "Time dependent Biot Poroelasticity problem with finite elements.\n\
2 We solve three field, quasi-static poroelasticity in a rectangular\n\
3 domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\
4 Contributed by: Robert Walker <rwalker6@buffalo.edu>\n\n\n";
5
6 #include <petscdmplex.h>
7 #include <petscsnes.h>
8 #include <petscts.h>
9 #include <petscds.h>
10 #include <petscbag.h>
11
12 #include <petsc/private/tsimpl.h>
13
14 /* This presentation of poroelasticity is taken from
15
16 @book{Cheng2016,
17 title={Poroelasticity},
18 author={Cheng, Alexander H-D},
19 volume={27},
20 year={2016},
21 publisher={Springer}
22 }
23
24 For visualization, use
25
26 -dm_view hdf5:${PETSC_DIR}/sol.h5 -monitor_solution hdf5:${PETSC_DIR}/sol.h5::append
27
28 The weak form would then be, using test function $(v, q, \tau)$,
29
30 (q, \frac{1}{M} \frac{dp}{dt}) + (q, \alpha \frac{d\varepsilon}{dt}) + (\nabla q, \kappa \nabla p) = (q, g)
31 -(\nabla v, 2 G \epsilon) - (\nabla\cdot v, \frac{2 G \nu}{1 - 2\nu} \varepsilon) + \alpha (\nabla\cdot v, p) = (v, f)
32 (\tau, \nabla \cdot u - \varepsilon) = 0
33 */
34
35 typedef enum {SOL_QUADRATIC_LINEAR, SOL_QUADRATIC_TRIG, SOL_TRIG_LINEAR, SOL_TERZAGHI, SOL_MANDEL, SOL_CRYER, NUM_SOLUTION_TYPES} SolutionType;
36 const char *solutionTypes[NUM_SOLUTION_TYPES+1] = {"quadratic_linear", "quadratic_trig", "trig_linear", "terzaghi", "mandel", "cryer", "unknown"};
37
38 typedef struct {
39 PetscScalar mu; /* shear modulus */
40 PetscScalar K_u; /* undrained bulk modulus */
41 PetscScalar alpha; /* Biot effective stress coefficient */
42 PetscScalar M; /* Biot modulus */
43 PetscScalar k; /* (isotropic) permeability */
44 PetscScalar mu_f; /* fluid dynamic viscosity */
45 PetscReal zmax; /* depth maximum extent */
46 PetscReal zmin; /* depth minimum extent */
47 PetscReal ymax; /* vertical maximum extent */
48 PetscReal ymin; /* vertical minimum extent */
49 PetscReal xmax; /* horizontal maximum extent */
50 PetscReal xmin; /* horizontal minimum extent */
51 PetscScalar P_0; /* magnitude of vertical stress */
52 } Parameter;
53
54 typedef struct {
55 /* Domain and mesh definition */
56 char dmType[256]; /* DM type for the solve */
57 PetscInt dim; /* The topological mesh dimension */
58 PetscBool simplex; /* Simplicial mesh */
59 PetscReal refLimit; /* Refine mesh with generator */
60 /* Problem definition */
61 SolutionType solType; /* Type of exact solution */
62 PetscBag bag; /* Problem parameters */
63 PetscReal t_r; /* Relaxation time: 4 L^2 / c */
64 PetscReal dtInitial; /* Override the choice for first timestep */
65 /* Exact solution terms */
66 PetscInt niter; /* Number of series term iterations in exact solutions */
67 PetscReal eps; /* Precision value for root finding */
68 PetscReal *zeroArray; /* Array of root locations */
69 } AppCtx;
70
zero(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)71 static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
72 {
73 PetscInt c;
74 for (c = 0; c < Nc; ++c) u[c] = 0.0;
75 return 0;
76 }
77
78 /* Quadratic space and linear time solution
79
80 2D:
81 u = x^2
82 v = y^2 - 2xy
83 p = (x + y) t
84 e = 2y
85 f = <2 G, 4 G + 2 \lambda > - <alpha t, alpha t>
86 g = 0
87 \epsilon = / 2x -y \
88 \ -y 2y - 2x /
89 Tr(\epsilon) = e = div u = 2y
90 div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
91 = 2 G < 2-1, 2 > + \lambda <0, 2> - alpha <t, t>
92 = <2 G, 4 G + 2 \lambda> - <alpha t, alpha t>
93 \frac{1}{M} \frac{dp}{dt} + \alpha \frac{d\varepsilon}{dt} - \nabla \cdot \kappa \nabla p
94 = \frac{1}{M} \frac{dp}{dt} + \kappa \Delta p
95 = (x + y)/M
96
97 3D:
98 u = x^2
99 v = y^2 - 2xy
100 w = z^2 - 2yz
101 p = (x + y + z) t
102 e = 2z
103 f = <2 G, 4 G + 2 \lambda > - <alpha t, alpha t, alpha t>
104 g = 0
105 \varepsilon = / 2x -y 0 \
106 | -y 2y - 2x -z |
107 \ 0 -z 2z - 2y/
108 Tr(\varepsilon) = div u = 2z
109 div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
110 = 2 G < 2-1, 2-1, 2 > + \lambda <0, 0, 2> - alpha <t, t, t>
111 = <2 G, 2G, 4 G + 2 \lambda> - <alpha t, alpha t, alpha t>
112 */
quadratic_u(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)113 static PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
114 {
115 PetscInt d;
116
117 for (d = 0; d < dim; ++d) {
118 u[d] = PetscSqr(x[d]) - (d > 0 ? 2.0 * x[d-1] * x[d] : 0.0);
119 }
120 return 0;
121 }
122
linear_eps(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)123 static PetscErrorCode linear_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
124 {
125 u[0] = 2.0*x[dim-1];
126 return 0;
127 }
128
linear_linear_p(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)129 static PetscErrorCode linear_linear_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
130 {
131 PetscReal sum = 0.0;
132 PetscInt d;
133
134 for (d = 0; d < dim; ++d) sum += x[d];
135 u[0] = sum*time;
136 return 0;
137 }
138
linear_linear_p_t(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)139 static PetscErrorCode linear_linear_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
140 {
141 PetscReal sum = 0.0;
142 PetscInt d;
143
144 for (d = 0; d < dim; ++d) sum += x[d];
145 u[0] = sum;
146 return 0;
147 }
148
f0_quadratic_linear_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 f0[])149 static void f0_quadratic_linear_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
150 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
151 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
152 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
153 {
154 const PetscReal G = PetscRealPart(constants[0]);
155 const PetscReal K_u = PetscRealPart(constants[1]);
156 const PetscReal alpha = PetscRealPart(constants[2]);
157 const PetscReal M = PetscRealPart(constants[3]);
158 const PetscReal K_d = K_u - alpha*alpha*M;
159 const PetscReal lambda = K_d - (2.0 * G) / 3.0;
160 PetscInt d;
161
162 for (d = 0; d < dim-1; ++d) {
163 f0[d] -= 2.0*G - alpha*t;
164 }
165 f0[dim-1] -= 2.0*lambda + 4.0*G - alpha*t;
166 }
167
f0_quadratic_linear_p(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 f0[])168 static void f0_quadratic_linear_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
169 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
170 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
171 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
172 {
173 const PetscReal alpha = PetscRealPart(constants[2]);
174 const PetscReal M = PetscRealPart(constants[3]);
175 PetscReal sum = 0.0;
176 PetscInt d;
177
178 for (d = 0; d < dim; ++d) sum += x[d];
179 f0[0] += u_t ? alpha*u_t[uOff[1]] : 0.0;
180 f0[0] += u_t ? u_t[uOff[2]]/M : 0.0;
181 f0[0] -= sum/M;
182 }
183
184 /* Quadratic space and trigonometric time solution
185
186 2D:
187 u = x^2
188 v = y^2 - 2xy
189 p = (x + y) cos(t)
190 e = 2y
191 f = <2 G, 4 G + 2 \lambda > - <alpha cos(t), alpha cos(t)>
192 g = 0
193 \epsilon = / 2x -y \
194 \ -y 2y - 2x /
195 Tr(\epsilon) = e = div u = 2y
196 div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
197 = 2 G < 2-1, 2 > + \lambda <0, 2> - alpha <cos(t), cos(t)>
198 = <2 G, 4 G + 2 \lambda> - <alpha cos(t), alpha cos(t)>
199 \frac{1}{M} \frac{dp}{dt} + \alpha \frac{d\varepsilon}{dt} - \nabla \cdot \kappa \nabla p
200 = \frac{1}{M} \frac{dp}{dt} + \kappa \Delta p
201 = -(x + y)/M sin(t)
202
203 3D:
204 u = x^2
205 v = y^2 - 2xy
206 w = z^2 - 2yz
207 p = (x + y + z) cos(t)
208 e = 2z
209 f = <2 G, 4 G + 2 \lambda > - <alpha cos(t), alpha cos(t), alpha cos(t)>
210 g = 0
211 \varepsilon = / 2x -y 0 \
212 | -y 2y - 2x -z |
213 \ 0 -z 2z - 2y/
214 Tr(\varepsilon) = div u = 2z
215 div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
216 = 2 G < 2-1, 2-1, 2 > + \lambda <0, 0, 2> - alpha <cos(t), cos(t), cos(t)>
217 = <2 G, 2G, 4 G + 2 \lambda> - <alpha cos(t), alpha cos(t), alpha cos(t)>
218 */
linear_trig_p(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)219 static PetscErrorCode linear_trig_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
220 {
221 PetscReal sum = 0.0;
222 PetscInt d;
223
224 for (d = 0; d < dim; ++d) sum += x[d];
225 u[0] = sum*PetscCosReal(time);
226 return 0;
227 }
228
linear_trig_p_t(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)229 static PetscErrorCode linear_trig_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
230 {
231 PetscReal sum = 0.0;
232 PetscInt d;
233
234 for (d = 0; d < dim; ++d) sum += x[d];
235 u[0] = -sum*PetscSinReal(time);
236 return 0;
237 }
238
f0_quadratic_trig_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 f0[])239 static void f0_quadratic_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
240 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
241 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
242 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
243 {
244 const PetscReal G = PetscRealPart(constants[0]);
245 const PetscReal K_u = PetscRealPart(constants[1]);
246 const PetscReal alpha = PetscRealPart(constants[2]);
247 const PetscReal M = PetscRealPart(constants[3]);
248 const PetscReal K_d = K_u - alpha*alpha*M;
249 const PetscReal lambda = K_d - (2.0 * G) / 3.0;
250 PetscInt d;
251
252 for (d = 0; d < dim-1; ++d) {
253 f0[d] -= 2.0*G - alpha*PetscCosReal(t);
254 }
255 f0[dim-1] -= 2.0*lambda + 4.0*G - alpha*PetscCosReal(t);
256 }
257
f0_quadratic_trig_p(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 f0[])258 static void f0_quadratic_trig_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
259 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
260 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
261 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
262 {
263 const PetscReal alpha = PetscRealPart(constants[2]);
264 const PetscReal M = PetscRealPart(constants[3]);
265 PetscReal sum = 0.0;
266 PetscInt d;
267
268 for (d = 0; d < dim; ++d) sum += x[d];
269
270 f0[0] += u_t ? alpha*u_t[uOff[1]] : 0.0;
271 f0[0] += u_t ? u_t[uOff[2]]/M : 0.0;
272 f0[0] += PetscSinReal(t)*sum/M;
273 }
274
275 /* Trigonometric space and linear time solution
276
277 u = sin(2 pi x)
278 v = sin(2 pi y) - 2xy
279 \varepsilon = / 2 pi cos(2 pi x) -y \
280 \ -y 2 pi cos(2 pi y) - 2x /
281 Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
282 div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
283 = \lambda \partial_j 2 pi (cos(2 pi x) + cos(2 pi y)) + 2\mu < -4 pi^2 sin(2 pi x) - 1, -4 pi^2 sin(2 pi y) >
284 = \lambda < -4 pi^2 sin(2 pi x) - 2, -4 pi^2 sin(2 pi y) > + \mu < -8 pi^2 sin(2 pi x) - 2, -8 pi^2 sin(2 pi y) >
285
286 2D:
287 u = sin(2 pi x)
288 v = sin(2 pi y) - 2xy
289 p = (cos(2 pi x) + cos(2 pi y)) t
290 e = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
291 f = < -4 pi^2 sin(2 pi x) (2 G + lambda) - (2 G - 2 lambda), -4 pi^2 sin(2 pi y) (2G + lambda) > + 2 pi alpha t <sin(2 pi x), sin(2 pi y)>
292 g = 0
293 \varepsilon = / 2 pi cos(2 pi x) -y \
294 \ -y 2 pi cos(2 pi y) - 2x /
295 Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
296 div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
297 = 2 G < -4 pi^2 sin(2 pi x) - 1, -4 pi^2 sin(2 pi y) > + \lambda <-4 pi^2 sin(2 pi x) - 2, -4 pi^2 sin(2 pi y)> - alpha <-2 pi sin(2 pi x) t, -2 pi sin(2 pi y) t>
298 = < -4 pi^2 sin(2 pi x) (2 G + lambda) - (2 G + 2 lambda), -4 pi^2 sin(2 pi y) (2G + lambda) > + 2 pi alpha t <sin(2 pi x), sin(2 pi y)>
299 \frac{1}{M} \frac{dp}{dt} + \alpha \frac{d\varepsilon}{dt} - \nabla \cdot \kappa \nabla p
300 = \frac{1}{M} \frac{dp}{dt} + \kappa \Delta p
301 = (cos(2 pi x) + cos(2 pi y))/M - 4 pi^2 \kappa (cos(2 pi x) + cos(2 pi y)) t
302
303 3D:
304 u = sin(2 pi x)
305 v = sin(2 pi y) - 2xy
306 v = sin(2 pi y) - 2yz
307 p = (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) t
308 e = 2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2y
309 f = < -4 pi^2 sin(2 pi x) (2 G + lambda) - (2 G + 2 lambda), -4 pi^2 sin(2 pi y) (2 G + lambda) - (2 G + 2 lambda), -4 pi^2 sin(2 pi z) (2G + lambda) > + 2 pi alpha t <sin(2 pi x), sin(2 pi y), , sin(2 pi z)>
310 g = 0
311 \varepsilon = / 2 pi cos(2 pi x) -y 0 \
312 | -y 2 pi cos(2 pi y) - 2x -z |
313 \ 0 -z 2 pi cos(2 pi z) - 2y /
314 Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2 y
315 div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
316 = 2 G < -4 pi^2 sin(2 pi x) - 1, -4 pi^2 sin(2 pi y) - 1, -4 pi^2 sin(2 pi z) > + \lambda <-4 pi^2 sin(2 pi x) - 2, 4 pi^2 sin(2 pi y) - 2, -4 pi^2 sin(2 pi y)> - alpha <-2 pi sin(2 pi x) t, -2 pi sin(2 pi y) t, -2 pi sin(2 pi z) t>
317 = < -4 pi^2 sin(2 pi x) (2 G + lambda) - (2 G + 2 lambda), -4 pi^2 sin(2 pi y) (2 G + lambda) - (2 G + 2 lambda), -4 pi^2 sin(2 pi z) (2G + lambda) > + 2 pi alpha t <sin(2 pi x), sin(2 pi y), , sin(2 pi z)>
318 \frac{1}{M} \frac{dp}{dt} + \alpha \frac{d\varepsilon}{dt} - \nabla \cdot \kappa \nabla p
319 = \frac{1}{M} \frac{dp}{dt} + \kappa \Delta p
320 = (cos(2 pi x) + cos(2 pi y) + cos(2 pi z))/M - 4 pi^2 \kappa (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) t
321 */
trig_u(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)322 static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
323 {
324 PetscInt d;
325
326 for (d = 0; d < dim; ++d) {
327 u[d] = PetscSinReal(2.*PETSC_PI*x[d]) - (d > 0 ? 2.0 * x[d-1] * x[d] : 0.0);
328 }
329 return 0;
330 }
331
trig_eps(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)332 static PetscErrorCode trig_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
333 {
334 PetscReal sum = 0.0;
335 PetscInt d;
336
337 for (d = 0; d < dim; ++d) sum += 2.*PETSC_PI*PetscCosReal(2.*PETSC_PI*x[d]) - (d < dim-1 ? 2.*x[d] : 0.0);
338 u[0] = sum;
339 return 0;
340 }
341
trig_linear_p(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)342 static PetscErrorCode trig_linear_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
343 {
344 PetscReal sum = 0.0;
345 PetscInt d;
346
347 for (d = 0; d < dim; ++d) sum += PetscCosReal(2.*PETSC_PI*x[d]);
348 u[0] = sum*time;
349 return 0;
350 }
351
trig_linear_p_t(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)352 static PetscErrorCode trig_linear_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
353 {
354 PetscReal sum = 0.0;
355 PetscInt d;
356
357 for (d = 0; d < dim; ++d) sum += PetscCosReal(2.*PETSC_PI*x[d]);
358 u[0] = sum;
359 return 0;
360 }
361
f0_trig_linear_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 f0[])362 static void f0_trig_linear_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
363 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
364 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
365 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
366 {
367 const PetscReal G = PetscRealPart(constants[0]);
368 const PetscReal K_u = PetscRealPart(constants[1]);
369 const PetscReal alpha = PetscRealPart(constants[2]);
370 const PetscReal M = PetscRealPart(constants[3]);
371 const PetscReal K_d = K_u - alpha*alpha*M;
372 const PetscReal lambda = K_d - (2.0 * G) / 3.0;
373 PetscInt d;
374
375 for (d = 0; d < dim-1; ++d) {
376 f0[d] += PetscSqr(2.*PETSC_PI)*PetscSinReal(2.*PETSC_PI*x[d])*(2.*G + lambda) + 2.0*(G + lambda) - 2.*PETSC_PI*alpha*PetscSinReal(2.*PETSC_PI*x[d])*t;
377 }
378 f0[dim-1] += PetscSqr(2.*PETSC_PI)*PetscSinReal(2.*PETSC_PI*x[dim-1])*(2.*G + lambda) - 2.*PETSC_PI*alpha*PetscSinReal(2.*PETSC_PI*x[dim-1])*t;
379 }
380
f0_trig_linear_p(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 f0[])381 static void f0_trig_linear_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
382 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
383 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
384 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
385 {
386 const PetscReal alpha = PetscRealPart(constants[2]);
387 const PetscReal M = PetscRealPart(constants[3]);
388 const PetscReal kappa = PetscRealPart(constants[4]);
389 PetscReal sum = 0.0;
390 PetscInt d;
391
392 for (d = 0; d < dim; ++d) sum += PetscCosReal(2.*PETSC_PI*x[d]);
393 f0[0] += u_t ? alpha*u_t[uOff[1]] : 0.0;
394 f0[0] += u_t ? u_t[uOff[2]]/M : 0.0;
395 f0[0] -= sum/M - 4*PetscSqr(PETSC_PI)*kappa*sum*t;
396 }
397
398 /* Terzaghi Solutions */
399 /* The analytical solutions given here are drawn from chapter 7, section 3, */
400 /* "One-Dimensional Consolidation Problem," from Poroelasticity, by Cheng. */
terzaghi_drainage_pressure(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)401 static PetscErrorCode terzaghi_drainage_pressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
402 {
403 AppCtx *user = (AppCtx *) ctx;
404 Parameter *param;
405 PetscErrorCode ierr;
406
407 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
408 if (time <= 0.0) {
409 PetscScalar alpha = param->alpha; /* - */
410 PetscScalar K_u = param->K_u; /* Pa */
411 PetscScalar M = param->M; /* Pa */
412 PetscScalar G = param->mu; /* Pa */
413 PetscScalar P_0 = param->P_0; /* Pa */
414 PetscScalar K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */
415 PetscScalar eta = (3.0*alpha*G) / (3.0*K_d + 4.0*G); /* -, Cheng (B.11) */
416 PetscScalar S = (3.0*K_u + 4.0*G) / (M*(3.0*K_d + 4.0*G)); /* Pa^{-1}, Cheng (B.14) */
417
418 u[0] = ((P_0*eta) / (G*S));
419 } else {
420 u[0] = 0.0;
421 }
422 return 0;
423 }
424
terzaghi_initial_u(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)425 static PetscErrorCode terzaghi_initial_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
426 {
427 AppCtx *user = (AppCtx *) ctx;
428 Parameter *param;
429 PetscErrorCode ierr;
430
431 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
432 {
433 PetscScalar K_u = param->K_u; /* Pa */
434 PetscScalar G = param->mu; /* Pa */
435 PetscScalar P_0 = param->P_0; /* Pa */
436 PetscReal L = param->ymax - param->ymin; /* m */
437 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); /* -, Cheng (B.9) */
438 PetscReal zstar = x[1] / L; /* - */
439
440 u[0] = 0.0;
441 u[1] = ((P_0*L*(1.0 - 2.0*nu_u)) / (2.0*G*(1.0 - nu_u))) * (1.0 - zstar);
442 }
443 return 0;
444 }
445
terzaghi_initial_eps(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)446 static PetscErrorCode terzaghi_initial_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
447 {
448 AppCtx *user = (AppCtx *) ctx;
449 Parameter *param;
450 PetscErrorCode ierr;
451
452 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
453 {
454 PetscScalar K_u = param->K_u; /* Pa */
455 PetscScalar G = param->mu; /* Pa */
456 PetscScalar P_0 = param->P_0; /* Pa */
457 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); /* -, Cheng (B.9) */
458
459 u[0] = -(P_0*(1.0 - 2.0*nu_u)) / (2.0*G*(1.0 - nu_u));
460 }
461 return 0;
462 }
463
terzaghi_2d_u(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)464 static PetscErrorCode terzaghi_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
465 {
466 AppCtx *user = (AppCtx *) ctx;
467 Parameter *param;
468 PetscErrorCode ierr;
469
470 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
471 if (time < 0.0) {
472 ierr = terzaghi_initial_u(dim, time, x, Nc, u, ctx);CHKERRQ(ierr);
473 } else {
474 PetscScalar alpha = param->alpha; /* - */
475 PetscScalar K_u = param->K_u; /* Pa */
476 PetscScalar M = param->M; /* Pa */
477 PetscScalar G = param->mu; /* Pa */
478 PetscScalar P_0 = param->P_0; /* Pa */
479 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
480 PetscReal L = param->ymax - param->ymin; /* m */
481 PetscInt N = user->niter, m;
482
483 PetscScalar K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */
484 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G)); /* -, Cheng (B.8) */
485 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); /* -, Cheng (B.9) */
486 PetscScalar S = (3.0*K_u + 4.0*G) / (M*(3.0*K_d + 4.0*G)); /* Pa^{-1}, Cheng (B.14) */
487 PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
488
489 PetscReal zstar = x[1] / L; /* - */
490 PetscReal tstar = PetscRealPart(c*time) / PetscSqr(2.0*L); /* - */
491 PetscScalar F2 = 0.0;
492
493 for (m = 1; m < 2*N+1; ++m) {
494 if (m%2 == 1) {
495 F2 += (8.0 / PetscSqr(m*PETSC_PI)) * PetscCosReal(0.5*m*PETSC_PI*zstar) * (1.0 - PetscExpReal(-PetscSqr(m*PETSC_PI)*tstar));
496 }
497 }
498 u[0] = 0.0;
499 u[1] = ((P_0*L*(1.0 - 2.0*nu_u)) / (2.0*G*(1.0 - nu_u))) * (1.0 - zstar) + ((P_0*L*(nu_u - nu)) / (2.0*G*(1.0 - nu_u)*(1.0 - nu)))*F2; /* m */
500 }
501 return 0;
502 }
503
terzaghi_2d_eps(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)504 static PetscErrorCode terzaghi_2d_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
505 {
506 AppCtx *user = (AppCtx *) ctx;
507 Parameter *param;
508 PetscErrorCode ierr;
509
510 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
511 if (time < 0.0) {
512 ierr = terzaghi_initial_eps(dim, time, x, Nc, u, ctx);CHKERRQ(ierr);
513 } else {
514 PetscScalar alpha = param->alpha; /* - */
515 PetscScalar K_u = param->K_u; /* Pa */
516 PetscScalar M = param->M; /* Pa */
517 PetscScalar G = param->mu; /* Pa */
518 PetscScalar P_0 = param->P_0; /* Pa */
519 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
520 PetscReal L = param->ymax - param->ymin; /* m */
521 PetscInt N = user->niter, m;
522
523 PetscScalar K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */
524 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G)); /* -, Cheng (B.8) */
525 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); /* -, Cheng (B.9) */
526 PetscScalar S = (3.0*K_u + 4.0*G) / (M*(3.0*K_d + 4.0*G)); /* Pa^{-1}, Cheng (B.14) */
527 PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
528
529 PetscReal zstar = x[1] / L; /* - */
530 PetscReal tstar = PetscRealPart(c*time) / PetscSqr(2.0*L); /* - */
531 PetscScalar F2_z = 0.0;
532
533 for (m = 1; m < 2*N+1; ++m) {
534 if (m%2 == 1) {
535 F2_z += (-4.0 / (m*PETSC_PI*L)) * PetscSinReal(0.5*m*PETSC_PI*zstar) * (1.0 - PetscExpReal(-PetscSqr(m*PETSC_PI)*tstar));
536 }
537 }
538 u[0] = -((P_0*L*(1.0 - 2.0*nu_u)) / (2.0*G*(1.0 - nu_u)*L)) + ((P_0*L*(nu_u - nu)) / (2.0*G*(1.0 - nu_u)*(1.0 - nu)))*F2_z; /* - */
539 }
540 return 0;
541 }
542
543 // Pressure
terzaghi_2d_p(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)544 static PetscErrorCode terzaghi_2d_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
545 {
546 AppCtx *user = (AppCtx *) ctx;
547 Parameter *param;
548 PetscErrorCode ierr;
549
550 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
551 if (time <= 0.0) {
552 ierr = terzaghi_drainage_pressure(dim, time, x, Nc, u, ctx);CHKERRQ(ierr);
553 } else {
554 PetscScalar alpha = param->alpha; /* - */
555 PetscScalar K_u = param->K_u; /* Pa */
556 PetscScalar M = param->M; /* Pa */
557 PetscScalar G = param->mu; /* Pa */
558 PetscScalar P_0 = param->P_0; /* Pa */
559 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
560 PetscReal L = param->ymax - param->ymin; /* m */
561 PetscInt N = user->niter, m;
562
563 PetscScalar K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */
564 PetscScalar eta = (3.0*alpha*G) / (3.0*K_d + 4.0*G); /* -, Cheng (B.11) */
565 PetscScalar S = (3.0*K_u + 4.0*G) / (M*(3.0*K_d + 4.0*G)); /* Pa^{-1}, Cheng (B.14) */
566 PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
567
568 PetscReal zstar = x[1] / L; /* - */
569 PetscReal tstar = PetscRealPart(c*time) / PetscSqr(2.0*L); /* - */
570 PetscScalar F1 = 0.0;
571
572 if (PetscAbsScalar((1/M + (alpha*eta)/G) - S) > 1.0e-10) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "S %g != check %g", S, (1/M + (alpha*eta)/G));
573
574 for (m = 1; m < 2*N+1; ++m) {
575 if (m%2 == 1) {
576 F1 += (4.0 / (m*PETSC_PI)) * PetscSinReal(0.5*m*PETSC_PI*zstar) * PetscExpReal(-PetscSqr(m*PETSC_PI)*tstar);
577 }
578 }
579 u[0] = ((P_0*eta) / (G*S)) * F1; /* Pa */
580 }
581 return 0;
582 }
583
terzaghi_2d_u_t(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)584 static PetscErrorCode terzaghi_2d_u_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
585 {
586 AppCtx *user = (AppCtx *) ctx;
587 Parameter *param;
588 PetscErrorCode ierr;
589
590 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
591 if (time <= 0.0) {
592 u[0] = 0.0;
593 u[1] = 0.0;
594 } else {
595 PetscScalar alpha = param->alpha; /* - */
596 PetscScalar K_u = param->K_u; /* Pa */
597 PetscScalar M = param->M; /* Pa */
598 PetscScalar G = param->mu; /* Pa */
599 PetscScalar P_0 = param->P_0; /* Pa */
600 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
601 PetscReal L = param->ymax - param->ymin; /* m */
602 PetscInt N = user->niter, m;
603
604 PetscScalar K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */
605 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G)); /* -, Cheng (B.8) */
606 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); /* -, Cheng (B.9) */
607 PetscScalar S = (3.0*K_u + 4.0*G) / (M*(3.0*K_d + 4.0*G)); /* Pa^{-1}, Cheng (B.14) */
608 PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
609
610 PetscReal zstar = x[1] / L; /* - */
611 PetscReal tstar = PetscRealPart(c*time) / PetscSqr(2.0*L); /* - */
612 PetscScalar F2_t = 0.0;
613
614 for (m = 1; m < 2*N+1; ++m) {
615 if (m%2 == 1) {
616 F2_t += (2.0*c / PetscSqr(L)) * PetscCosReal(0.5*m*PETSC_PI*zstar) * PetscExpReal(-PetscSqr(m*PETSC_PI)*tstar);
617 }
618 }
619 u[0] = 0.0;
620 u[1] = ((P_0*L*(nu_u - nu)) / (2.0*G*(1.0 - nu_u)*(1.0 - nu)))*F2_t; /* m / s */
621 }
622 return 0;
623 }
624
terzaghi_2d_eps_t(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)625 static PetscErrorCode terzaghi_2d_eps_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
626 {
627 AppCtx *user = (AppCtx *) ctx;
628 Parameter *param;
629 PetscErrorCode ierr;
630
631 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
632 if (time <= 0.0) {
633 u[0] = 0.0;
634 } else {
635 PetscScalar alpha = param->alpha; /* - */
636 PetscScalar K_u = param->K_u; /* Pa */
637 PetscScalar M = param->M; /* Pa */
638 PetscScalar G = param->mu; /* Pa */
639 PetscScalar P_0 = param->P_0; /* Pa */
640 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
641 PetscReal L = param->ymax - param->ymin; /* m */
642 PetscInt N = user->niter, m;
643
644 PetscScalar K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */
645 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G)); /* -, Cheng (B.8) */
646 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); /* -, Cheng (B.9) */
647 PetscScalar S = (3.0*K_u + 4.0*G) / (M*(3.0*K_d + 4.0*G)); /* Pa^{-1}, Cheng (B.14) */
648 PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
649
650 PetscReal zstar = x[1] / L; /* - */
651 PetscReal tstar = PetscRealPart(c*time) / PetscSqr(2.0*L); /* - */
652 PetscScalar F2_zt = 0.0;
653
654 for (m = 1; m < 2*N+1; ++m) {
655 if (m%2 == 1) {
656 F2_zt += ((-m*PETSC_PI*c) / (L*L*L)) * PetscSinReal(0.5*m*PETSC_PI*zstar) * PetscExpReal(-PetscSqr(m*PETSC_PI)*tstar);
657 }
658 }
659 u[0] = ((P_0*L*(nu_u - nu)) / (2.0*G*(1.0 - nu_u)*(1.0 - nu)))*F2_zt; /* 1 / s */
660 }
661 return 0;
662 }
663
terzaghi_2d_p_t(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)664 static PetscErrorCode terzaghi_2d_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
665 {
666
667 AppCtx *user = (AppCtx *) ctx;
668 Parameter *param;
669 PetscErrorCode ierr;
670
671 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
672 if (time <= 0.0) {
673 PetscScalar alpha = param->alpha; /* - */
674 PetscScalar K_u = param->K_u; /* Pa */
675 PetscScalar M = param->M; /* Pa */
676 PetscScalar G = param->mu; /* Pa */
677 PetscScalar P_0 = param->P_0; /* Pa */
678 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
679 PetscReal L = param->ymax - param->ymin; /* m */
680
681 PetscScalar K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */
682 PetscScalar eta = (3.0*alpha*G) / (3.0*K_d + 4.0*G); /* -, Cheng (B.11) */
683 PetscScalar S = (3.0*K_u + 4.0*G) / (M*(3.0*K_d + 4.0*G)); /* Pa^{-1}, Cheng (B.14) */
684 PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
685
686 u[0] = -((P_0*eta) / (G*S)) * PetscSqr(0*PETSC_PI)*c / PetscSqr(2.0*L); /* Pa / s */
687 } else {
688 PetscScalar alpha = param->alpha; /* - */
689 PetscScalar K_u = param->K_u; /* Pa */
690 PetscScalar M = param->M; /* Pa */
691 PetscScalar G = param->mu; /* Pa */
692 PetscScalar P_0 = param->P_0; /* Pa */
693 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
694 PetscReal L = param->ymax - param->ymin; /* m */
695 PetscInt N = user->niter, m;
696
697 PetscScalar K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */
698 PetscScalar eta = (3.0*alpha*G) / (3.0*K_d + 4.0*G); /* -, Cheng (B.11) */
699 PetscScalar S = (3.0*K_u + 4.0*G) / (M*(3.0*K_d + 4.0*G)); /* Pa^{-1}, Cheng (B.14) */
700 PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
701
702 PetscReal zstar = x[1] / L; /* - */
703 PetscReal tstar = PetscRealPart(c*time) / PetscSqr(2.0*L); /* - */
704 PetscScalar F1_t = 0.0;
705 PetscScalar F1_zz = 0.0;
706
707 if (PetscAbsScalar((1/M + (alpha*eta)/G) - S) > 1.0e-10) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "S %g != check %g", S, (1/M + (alpha*eta)/G));
708
709 for (m = 1; m < 2*N+1; ++m) {
710 if (m%2 == 1) {
711 F1_t += ((-m*PETSC_PI*c) / PetscSqr(L)) * PetscSinReal(0.5*m*PETSC_PI*zstar) * PetscExpReal(-PetscSqr(m*PETSC_PI)*tstar);
712 F1_zz += (-m*PETSC_PI / PetscSqr(L)) * PetscSinReal(0.5*m*PETSC_PI*zstar) * PetscExpReal(-PetscSqr(m*PETSC_PI)*tstar);
713 }
714 }
715 u[0] = ((P_0*eta) / (G*S)) * F1_t; /* Pa / s */
716 }
717 return 0;
718 }
719
720 /* Mandel Solutions */
mandel_drainage_pressure(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)721 static PetscErrorCode mandel_drainage_pressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
722 {
723 AppCtx *user = (AppCtx *) ctx;
724 Parameter *param;
725 PetscErrorCode ierr;
726
727 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
728 if (time <= 0.0) {
729 PetscScalar alpha = param->alpha; /* - */
730 PetscScalar K_u = param->K_u; /* Pa */
731 PetscScalar M = param->M; /* Pa */
732 PetscScalar G = param->mu; /* Pa */
733 PetscScalar P_0 = param->P_0; /* Pa */
734 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
735 PetscReal a = 0.5*(param->xmax - param->xmin); /* m */
736 PetscInt N = user->niter, n;
737
738 PetscScalar K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */
739 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); /* -, Cheng (B.9) */
740 PetscScalar B = alpha*M / K_u; /* -, Cheng (B.12) */
741 PetscScalar S = (3.0*K_u + 4.0*G) / (M*(3.0*K_d + 4.0*G)); /* Pa^{-1}, Cheng (B.14) */
742 PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
743
744 PetscScalar A1 = 3.0 / (B * (1.0 + nu_u));
745 PetscReal aa = 0.0;
746 PetscReal p = 0.0;
747 PetscReal time = 0.0;
748
749 for (n = 1; n < N+1; ++n) {
750 aa = user->zeroArray[n-1];
751 p += (PetscSinReal(aa) / (aa - PetscSinReal(aa)*PetscCosReal(aa))) * (PetscCosReal( (aa*x[0]) / a) - PetscCosReal(aa)) * PetscExpReal(-1.0*(aa*aa * PetscRealPart(c) * time)/(a*a));
752 }
753 u[0] = ((2.0 * P_0) / (a*A1)) * p;
754 } else {
755 u[0] = 0.0;
756 }
757 return 0;
758 }
759
mandel_initial_u(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)760 static PetscErrorCode mandel_initial_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
761 {
762 AppCtx *user = (AppCtx *) ctx;
763 Parameter *param;
764 PetscErrorCode ierr;
765
766 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
767 {
768 PetscScalar alpha = param->alpha; /* - */
769 PetscScalar K_u = param->K_u; /* Pa */
770 PetscScalar M = param->M; /* Pa */
771 PetscScalar G = param->mu; /* Pa */
772 PetscScalar P_0 = param->P_0; /* Pa */
773 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
774 PetscScalar a = 0.5*(param->xmax - param->xmin); /* m */
775 PetscInt N = user->niter, n;
776
777 PetscScalar K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */
778 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G)); /* -, Cheng (B.8) */
779 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); /* -, Cheng (B.9) */
780 PetscScalar S = (3.0*K_u + 4.0*G) / (M*(3.0*K_d + 4.0*G)); /* Pa^{-1}, Cheng (B.14) */
781 PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
782
783 PetscScalar A_s = 0.0;
784 PetscScalar B_s = 0.0;
785 PetscScalar time = 0.0;
786 PetscScalar alpha_n = 0.0;
787
788 for (n = 1; n < N+1; ++n) {
789 alpha_n = user->zeroArray[n-1];
790 A_s += ((PetscSinReal(alpha_n) * PetscCosReal(alpha_n)) / (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscExpReal(-1*(alpha_n*alpha_n*c*time)/(a*a));
791 B_s += (PetscCosReal(alpha_n) / (alpha_n - PetscSinReal(alpha_n)*PetscCosReal(alpha_n))) * PetscSinReal( (alpha_n * x[0])/a) * PetscExpReal(-1*(alpha_n*alpha_n*c*time)/(a*a));
792 }
793 u[0] = ((P_0*nu)/(2.0*G*a) - (P_0*nu_u)/(G*a) * A_s)* x[0] + P_0/G * B_s;
794 u[1] = (-1*(P_0*(1.0-nu))/(2*G*a) + (P_0*(1-nu_u))/(G*a) * A_s)*x[1];
795 }
796 return 0;
797 }
798
mandel_initial_eps(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)799 static PetscErrorCode mandel_initial_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
800 {
801 AppCtx *user = (AppCtx *) ctx;
802 Parameter *param;
803 PetscErrorCode ierr;
804
805 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
806 {
807 PetscScalar alpha = param->alpha; /* - */
808 PetscScalar K_u = param->K_u; /* Pa */
809 PetscScalar M = param->M; /* Pa */
810 PetscScalar G = param->mu; /* Pa */
811 PetscScalar P_0 = param->P_0; /* Pa */
812 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
813 PetscReal a = 0.5*(param->xmax - param->xmin); /* m */
814 PetscInt N = user->niter, n;
815
816 PetscScalar K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */
817 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G)); /* -, Cheng (B.8) */
818 PetscScalar S = (3.0*K_u + 4.0*G) / (M*(3.0*K_d + 4.0*G)); /* Pa^{-1}, Cheng (B.14) */
819 PetscReal c = PetscRealPart(kappa / S); /* m^2 / s, Cheng (B.16) */
820
821 PetscReal aa = 0.0;
822 PetscReal eps_A = 0.0;
823 PetscReal eps_B = 0.0;
824 PetscReal eps_C = 0.0;
825 PetscReal time = 0.0;
826
827 for (n = 1; n < N+1; ++n) {
828 aa = user->zeroArray[n-1];
829 eps_A += (aa * PetscExpReal( (-1.0*aa*aa*c*time)/(a*a))*PetscCosReal(aa)*PetscCosReal( (aa*x[0])/a)) / (a * (aa - PetscSinReal(aa)*PetscCosReal(aa)));
830 eps_B += ( PetscExpReal( (-1.0*aa*aa*c*time)/(a*a))*PetscSinReal(aa)*PetscCosReal(aa)) / (aa - PetscSinReal(aa)*PetscCosReal(aa));
831 eps_C += ( PetscExpReal( (-1.0*aa*aa*c*time)/(aa*aa))*PetscSinReal(aa)*PetscCosReal(aa)) / (aa - PetscSinReal(aa)*PetscCosReal(aa));
832 }
833 u[0] = (P_0/G)*eps_A + ( (P_0*nu)/(2.0*G*a)) - eps_B/(G*a) - (P_0*(1-nu))/(2*G*a) + eps_C/(G*a);
834 }
835 return 0;
836 }
837
838 // Displacement
mandel_2d_u(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)839 static PetscErrorCode mandel_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
840 {
841
842 Parameter *param;
843 PetscErrorCode ierr;
844
845 AppCtx *user = (AppCtx *) ctx;
846
847 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
848 if (time <= 0.0) {
849 ierr = mandel_initial_u(dim, time, x, Nc, u, ctx);CHKERRQ(ierr);
850 } else {
851 PetscInt NITER = user->niter;
852 PetscScalar alpha = param->alpha;
853 PetscScalar K_u = param->K_u;
854 PetscScalar M = param->M;
855 PetscScalar G = param->mu;
856 PetscScalar k = param->k;
857 PetscScalar mu_f = param->mu_f;
858 PetscScalar F = param->P_0;
859
860 PetscScalar K_d = K_u - alpha*alpha*M;
861 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G));
862 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G));
863 PetscScalar kappa = k / mu_f;
864 PetscReal a = (param->xmax - param->xmin) / 2.0;
865 PetscReal c = PetscRealPart(((2.0*kappa*G) * (1.0 - nu) * (nu_u - nu)) / ( alpha*alpha * (1.0 - 2.0*nu) * (1.0 - nu_u)));
866
867 // Series term
868 PetscScalar A_x = 0.0;
869 PetscScalar B_x = 0.0;
870
871 for (PetscInt n=1; n < NITER+1; n++) {
872 PetscReal alpha_n = user->zeroArray[n-1];
873
874 A_x += ( (PetscSinReal(alpha_n) * PetscCosReal(alpha_n)) / (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscExpReal( -1*(alpha_n*alpha_n*c*time)/(a*a));
875 B_x += ( PetscCosReal(alpha_n) / (alpha_n - PetscSinReal(alpha_n)*PetscCosReal(alpha_n))) * PetscSinReal( (alpha_n * x[0])/a) * PetscExpReal( -1*(alpha_n*alpha_n*c*time)/(a*a));
876 }
877 u[0] = ((F*nu)/(2.0*G*a) - (F*nu_u)/(G*a) * A_x)* x[0] + F/G * B_x;
878 u[1] = (-1*(F*(1.0-nu))/(2*G*a) + (F*(1-nu_u))/(G*a) * A_x)*x[1];
879 }
880 return 0;
881 }
882
883 // Trace strain
mandel_2d_eps(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)884 static PetscErrorCode mandel_2d_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
885 {
886
887 Parameter *param;
888 PetscErrorCode ierr;
889
890 AppCtx *user = (AppCtx *) ctx;
891
892 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
893 if (time <= 0.0) {
894 ierr = mandel_initial_eps(dim, time, x, Nc, u, ctx);CHKERRQ(ierr);
895 } else {
896 PetscInt NITER = user->niter;
897 PetscScalar alpha = param->alpha;
898 PetscScalar K_u = param->K_u;
899 PetscScalar M = param->M;
900 PetscScalar G = param->mu;
901 PetscScalar k = param->k;
902 PetscScalar mu_f = param->mu_f;
903 PetscScalar F = param->P_0;
904
905 PetscScalar K_d = K_u - alpha*alpha*M;
906 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G));
907 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G));
908 PetscScalar kappa = k / mu_f;
909 //const PetscScalar B = (alpha*M)/(K_d + alpha*alpha * M);
910
911 //const PetscScalar b = (YMAX - YMIN) / 2.0;
912 PetscScalar a = (param->xmax - param->xmin) / 2.0;
913 PetscReal c = PetscRealPart(((2.0*kappa*G) * (1.0 - nu) * (nu_u - nu)) / (alpha*alpha * (1.0 - 2.0*nu) * (1.0 - nu_u)));
914
915 // Series term
916 PetscScalar eps_A = 0.0;
917 PetscScalar eps_B = 0.0;
918 PetscScalar eps_C = 0.0;
919
920 for (PetscInt n=1; n < NITER+1; n++)
921 {
922 PetscReal aa = user->zeroArray[n-1];
923
924 eps_A += (aa * PetscExpReal( (-1.0*aa*aa*c*time)/(a*a))*PetscCosReal(aa)*PetscCosReal( (aa*x[0])/a)) / (a * (aa - PetscSinReal(aa)*PetscCosReal(aa)));
925
926 eps_B += ( PetscExpReal( (-1.0*aa*aa*c*time)/(a*a))*PetscSinReal(aa)*PetscCosReal(aa)) / (aa - PetscSinReal(aa)*PetscCosReal(aa));
927
928 eps_C += ( PetscExpReal( (-1.0*aa*aa*c*time)/(aa*aa))*PetscSinReal(aa)*PetscCosReal(aa)) / (aa - PetscSinReal(aa)*PetscCosReal(aa));
929 }
930
931 u[0] = (F/G)*eps_A + ( (F*nu)/(2.0*G*a)) - eps_B/(G*a) - (F*(1-nu))/(2*G*a) + eps_C/(G*a);
932 }
933 return 0;
934
935 }
936
937 // Pressure
mandel_2d_p(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)938 static PetscErrorCode mandel_2d_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
939 {
940
941 Parameter *param;
942 PetscErrorCode ierr;
943
944 AppCtx *user = (AppCtx *) ctx;
945
946 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
947 if (time <= 0.0) {
948 ierr = mandel_drainage_pressure(dim, time, x, Nc, u, ctx);CHKERRQ(ierr);
949 } else {
950 PetscInt NITER = user->niter;
951
952 PetscScalar alpha = param->alpha;
953 PetscScalar K_u = param->K_u;
954 PetscScalar M = param->M;
955 PetscScalar G = param->mu;
956 PetscScalar k = param->k;
957 PetscScalar mu_f = param->mu_f;
958 PetscScalar F = param->P_0;
959
960 PetscScalar K_d = K_u - alpha*alpha*M;
961 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G));
962 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G));
963 PetscScalar kappa = k / mu_f;
964 PetscScalar B = (alpha*M)/(K_d + alpha*alpha * M);
965
966 PetscReal a = (param->xmax - param->xmin) / 2.0;
967 PetscReal c = PetscRealPart(((2.0*kappa*G) * (1.0 - nu) * (nu_u - nu)) / (alpha*alpha * (1.0 - 2.0*nu) * (1.0 - nu_u)));
968 PetscScalar A1 = 3.0 / (B * (1.0 + nu_u));
969 //PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu);
970
971 // Series term
972 PetscScalar aa = 0.0;
973 PetscScalar p = 0.0;
974
975 for (PetscInt n=1; n < NITER+1; n++)
976 {
977 aa = user->zeroArray[n-1];
978 p += (PetscSinReal(aa)/ (aa - PetscSinReal(aa)*PetscCosReal(aa))) * (PetscCosReal( (aa*x[0]) / a) - PetscCosReal(aa)) * PetscExpReal(-1.0*(aa*aa * c * time)/(a*a));
979 }
980 u[0] = ((2.0 * F) / (a*A1)) * p;
981 }
982 return 0;
983 }
984
985 // Time derivative of displacement
mandel_2d_u_t(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)986 static PetscErrorCode mandel_2d_u_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
987 {
988
989 Parameter *param;
990 PetscErrorCode ierr;
991
992 AppCtx *user = (AppCtx *) ctx;
993
994 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
995
996 PetscInt NITER = user->niter;
997 PetscScalar alpha = param->alpha;
998 PetscScalar K_u = param->K_u;
999 PetscScalar M = param->M;
1000 PetscScalar G = param->mu;
1001 PetscScalar F = param->P_0;
1002
1003 PetscScalar K_d = K_u - alpha*alpha*M;
1004 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G));
1005 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G));
1006 PetscScalar kappa = param->k / param->mu_f;
1007 PetscReal a = (param->xmax - param->xmin) / 2.0;
1008 PetscReal c = PetscRealPart(((2.0*kappa*G) * (1.0 - nu) * (nu_u - nu)) / (alpha*alpha * (1.0 - 2.0*nu) * (1.0 - nu_u)));
1009
1010 // Series term
1011 PetscScalar A_s_t = 0.0;
1012 PetscScalar B_s_t = 0.0;
1013
1014 for (PetscInt n=1; n < NITER+1; n++)
1015 {
1016 PetscReal alpha_n = user->zeroArray[n-1];
1017
1018 A_s_t += (-1.0*alpha_n*alpha_n*c*PetscExpReal( (-1.0*alpha_n*alpha_n*time)/(a*a))*PetscSinReal( (alpha_n*x[0])/a) * PetscCosReal(alpha_n)) / ( a*a*(alpha_n - PetscSinReal(alpha_n)*PetscCosReal(alpha_n)));
1019 B_s_t += (-1.0*alpha_n*alpha_n*c*PetscExpReal( (-1.0*alpha_n*alpha_n*time)/(a*a))*PetscSinReal( alpha_n) * PetscCosReal(alpha_n)) / ( a*a*(alpha_n - PetscSinReal(alpha_n)*PetscCosReal(alpha_n)));
1020 }
1021
1022 u[0] = (F/G)*A_s_t - ( (F*nu_u*x[0])/(G*a))*B_s_t;
1023 u[1] = ( (F*x[1]*(1 - nu_u)) / (G*a))*B_s_t;
1024
1025 return 0;
1026 }
1027
1028 // Time derivative of trace strain
mandel_2d_eps_t(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)1029 static PetscErrorCode mandel_2d_eps_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1030 {
1031
1032 Parameter *param;
1033 PetscErrorCode ierr;
1034
1035 AppCtx *user = (AppCtx *) ctx;
1036
1037 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
1038
1039 PetscInt NITER = user->niter;
1040 PetscScalar alpha = param->alpha;
1041 PetscScalar K_u = param->K_u;
1042 PetscScalar M = param->M;
1043 PetscScalar G = param->mu;
1044 PetscScalar k = param->k;
1045 PetscScalar mu_f = param->mu_f;
1046 PetscScalar F = param->P_0;
1047
1048 PetscScalar K_d = K_u - alpha*alpha*M;
1049 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G));
1050 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G));
1051 PetscScalar kappa = k / mu_f;
1052 //const PetscScalar B = (alpha*M)/(K_d + alpha*alpha * M);
1053
1054 //const PetscScalar b = (YMAX - YMIN) / 2.0;
1055 PetscReal a = (param->xmax - param->xmin) / 2.0;
1056 PetscReal c = PetscRealPart(((2.0*kappa*G) * (1.0 - nu) * (nu_u - nu)) / (alpha*alpha * (1.0 - 2.0*nu) * (1.0 - nu_u)));
1057
1058 // Series term
1059 PetscScalar eps_As = 0.0;
1060 PetscScalar eps_Bs = 0.0;
1061 PetscScalar eps_Cs = 0.0;
1062
1063 for (PetscInt n=1; n < NITER+1; n++)
1064 {
1065 PetscReal alpha_n = user->zeroArray[n-1];
1066
1067 eps_As += (-1.0*alpha_n*alpha_n*alpha_n*c*PetscExpReal( (-1.0*alpha_n*alpha_n*c*time)/(a*a))*PetscCosReal(alpha_n)*PetscCosReal( (alpha_n*x[0])/a)) / ( alpha_n*alpha_n*alpha_n*(alpha_n - PetscSinReal(alpha_n)*PetscCosReal(alpha_n)));
1068 eps_Bs += (-1.0*alpha_n*alpha_n*c*PetscExpReal( (-1.0*alpha_n*alpha_n*c*time)/(a*a))*PetscSinReal(alpha_n)*PetscCosReal(alpha_n)) / (alpha_n*alpha_n * (alpha_n - PetscSinReal(alpha_n)*PetscCosReal(alpha_n)));
1069 eps_Cs += (-1.0*alpha_n*alpha_n*c*PetscExpReal( (-1.0*alpha_n*alpha_n*c*time)/(a*a))*PetscSinReal(alpha_n)*PetscCosReal(alpha_n)) / (alpha_n*alpha_n * (alpha_n - PetscSinReal(alpha_n)*PetscCosReal(alpha_n)));
1070 }
1071
1072 u[0] = (F/G)*eps_As - ( (F*nu_u)/(G*a))*eps_Bs + ( (F*(1-nu_u))/(G*a))*eps_Cs;
1073 return 0;
1074
1075 }
1076
1077 // Time derivative of pressure
mandel_2d_p_t(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)1078 static PetscErrorCode mandel_2d_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1079 {
1080
1081 Parameter *param;
1082 PetscErrorCode ierr;
1083
1084 AppCtx *user = (AppCtx *) ctx;
1085
1086 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
1087
1088 PetscInt NITER = user->niter;
1089
1090 PetscScalar alpha = param->alpha;
1091 PetscScalar K_u = param->K_u;
1092 PetscScalar M = param->M;
1093 PetscScalar G = param->mu;
1094 PetscScalar k = param->k;
1095 PetscScalar mu_f = param->mu_f;
1096 PetscScalar F = param->P_0;
1097
1098 PetscScalar K_d = K_u - alpha*alpha*M;
1099 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G));
1100 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G));
1101 PetscScalar kappa = k / mu_f;
1102
1103 PetscReal a = (param->xmax - param->xmin) / 2.0;
1104 PetscReal c = PetscRealPart(((2.0*kappa*G) * (1.0 - nu) * (nu_u - nu)) / (alpha*alpha * (1.0 - 2.0*nu) * (1.0 - nu_u)));
1105 //PetscScalar A1 = 3.0 / (B * (1.0 + nu_u));
1106 //PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu);
1107
1108 // Series term
1109 PetscScalar P_s = 0.0;
1110
1111 for (PetscInt n=1; n < NITER+1; n++)
1112 {
1113 PetscReal alpha_n = user->zeroArray[n-1];
1114
1115 P_s += (-1.0*alpha_n*alpha_n*c*( -1.0*PetscCosReal(alpha_n) + PetscCosReal( (alpha_n*x[0])/a))*PetscExpReal( (-1.0*alpha_n*alpha_n*c*time)/(a*a))*PetscSinReal(alpha_n)) / ( a*a*(alpha_n - PetscSinReal(alpha_n)*PetscCosReal(alpha_n)));
1116 }
1117 u[0] = ( (2.0*F*(-2.0*nu + 3.0*nu_u))/(3.0*a*alpha*(1.0 - 2.0*nu)));
1118
1119 return 0;
1120 }
1121
1122 /* Cryer Solutions */
cryer_drainage_pressure(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)1123 static PetscErrorCode cryer_drainage_pressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1124 {
1125 AppCtx *user = (AppCtx *) ctx;
1126 Parameter *param;
1127 PetscErrorCode ierr;
1128
1129 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
1130 if (time <= 0.0) {
1131 PetscScalar alpha = param->alpha; /* - */
1132 PetscScalar K_u = param->K_u; /* Pa */
1133 PetscScalar M = param->M; /* Pa */
1134 PetscScalar P_0 = param->P_0; /* Pa */
1135 PetscScalar B = alpha*M / K_u; /* -, Cheng (B.12) */
1136
1137 u[0] = P_0*B;
1138 } else {
1139 u[0] = 0.0;
1140 }
1141 return 0;
1142 }
1143
cryer_initial_u(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)1144 static PetscErrorCode cryer_initial_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1145 {
1146 AppCtx *user = (AppCtx *) ctx;
1147 Parameter *param;
1148 PetscErrorCode ierr;
1149
1150 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
1151 {
1152 PetscScalar K_u = param->K_u; /* Pa */
1153 PetscScalar G = param->mu; /* Pa */
1154 PetscScalar P_0 = param->P_0; /* Pa */
1155 PetscReal R_0 = param->ymax; /* m */
1156 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); /* -, Cheng (B.9) */
1157
1158 PetscScalar u_0 = -P_0*R_0*(1. - 2.*nu_u) / (2.*G*(1. + nu_u)); /* Cheng (7.407) */
1159 PetscReal u_sc = PetscRealPart(u_0)/R_0;
1160
1161 u[0] = u_sc * x[0];
1162 u[1] = u_sc * x[1];
1163 u[2] = u_sc * x[2];
1164 }
1165 return 0;
1166 }
1167
cryer_initial_eps(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)1168 static PetscErrorCode cryer_initial_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1169 {
1170 AppCtx *user = (AppCtx *) ctx;
1171 Parameter *param;
1172 PetscErrorCode ierr;
1173
1174 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
1175 {
1176 PetscScalar K_u = param->K_u; /* Pa */
1177 PetscScalar G = param->mu; /* Pa */
1178 PetscScalar P_0 = param->P_0; /* Pa */
1179 PetscReal R_0 = param->ymax; /* m */
1180 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); /* -, Cheng (B.9) */
1181
1182 PetscScalar u_0 = -P_0*R_0*(1. - 2.*nu_u) / (2.*G*(1. + nu_u)); /* Cheng (7.407) */
1183 PetscReal u_sc = PetscRealPart(u_0)/R_0;
1184
1185 /* div R = 1/R^2 d/dR R^2 R = 3 */
1186 u[0] = 3.*u_sc;
1187 u[1] = 3.*u_sc;
1188 u[2] = 3.*u_sc;
1189 }
1190 return 0;
1191 }
1192
1193 // Displacement
cryer_3d_u(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)1194 static PetscErrorCode cryer_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1195 {
1196 AppCtx *user = (AppCtx *) ctx;
1197 Parameter *param;
1198 PetscErrorCode ierr;
1199
1200 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
1201 if (time <= 0.0) {
1202 ierr = cryer_initial_u(dim, time, x, Nc, u, ctx);CHKERRQ(ierr);
1203 } else {
1204 PetscScalar alpha = param->alpha; /* - */
1205 PetscScalar K_u = param->K_u; /* Pa */
1206 PetscScalar M = param->M; /* Pa */
1207 PetscScalar G = param->mu; /* Pa */
1208 PetscScalar P_0 = param->P_0; /* Pa */
1209 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
1210 PetscReal R_0 = param->ymax; /* m */
1211 PetscInt N = user->niter, n;
1212
1213 PetscScalar K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */
1214 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G)); /* -, Cheng (B.8) */
1215 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); /* -, Cheng (B.9) */
1216 PetscScalar S = (3.0*K_u + 4.0*G) / (M*(3.0*K_d + 4.0*G)); /* Pa^{-1}, Cheng (B.14) */
1217 PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
1218 PetscScalar u_inf = -P_0*R_0*(1. - 2.*nu) / (2.*G*(1. + nu)); /* m, Cheng (7.388) */
1219
1220 PetscReal R = PetscSqrtReal(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
1221 PetscReal R_star = R/R_0;
1222 PetscReal tstar = PetscRealPart(c*time) / PetscSqr(R_0); /* - */
1223 PetscReal A_n = 0.0;
1224 PetscScalar u_sc;
1225
1226 for (n = 1; n < N+1; ++n) {
1227 const PetscReal x_n = user->zeroArray[n-1];
1228 const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu)*PetscSqr(1 + nu_u)*x_n - 18.0*(1 + nu)*(nu_u - nu)*(1 - nu_u));
1229
1230 /* m , Cheng (7.404) */
1231 A_n += PetscRealPart(
1232 (12.0*(1.0 + nu)*(nu_u - nu))/((1.0 - 2.0*nu)*E_n*PetscSqr(R_star)*x_n*PetscSinReal(PetscSqrtReal(x_n))) *
1233 (3.0*(nu_u - nu) * (PetscSinReal(R_star * PetscSqrtReal(x_n)) - R_star*PetscSqrtReal(x_n)*PetscCosReal(R_star * PetscSqrtReal(x_n)))
1234 + (1.0 - nu)*(1.0 - 2.0*nu)*PetscPowRealInt(R_star, 3)*x_n*PetscSinReal(PetscSqrtReal(x_n))) * PetscExpReal(-x_n * tstar));
1235 }
1236 u_sc = PetscRealPart(u_inf) * (R_star - A_n);
1237 u[0] = u_sc * x[0] / R;
1238 u[1] = u_sc * x[1] / R;
1239 u[2] = u_sc * x[2] / R;
1240 }
1241 return 0;
1242 }
1243
1244 // Volumetric Strain
cryer_3d_eps(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)1245 static PetscErrorCode cryer_3d_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1246 {
1247 AppCtx *user = (AppCtx *) ctx;
1248 Parameter *param;
1249 PetscErrorCode ierr;
1250
1251 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
1252 if (time <= 0.0) {
1253 ierr = cryer_initial_eps(dim, time, x, Nc, u, ctx);CHKERRQ(ierr);
1254 } else {
1255 PetscScalar alpha = param->alpha; /* - */
1256 PetscScalar K_u = param->K_u; /* Pa */
1257 PetscScalar M = param->M; /* Pa */
1258 PetscScalar G = param->mu; /* Pa */
1259 PetscScalar P_0 = param->P_0; /* Pa */
1260 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
1261 PetscReal R_0 = param->ymax; /* m */
1262 PetscInt N = user->niter, n;
1263
1264 PetscScalar K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */
1265 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G)); /* -, Cheng (B.8) */
1266 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); /* -, Cheng (B.9) */
1267 PetscScalar S = (3.0*K_u + 4.0*G) / (M*(3.0*K_d + 4.0*G)); /* Pa^{-1}, Cheng (B.14) */
1268 PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
1269 PetscScalar u_inf = -P_0*R_0*(1. - 2.*nu) / (2.*G*(1. + nu)); /* m, Cheng (7.388) */
1270
1271 PetscReal R = PetscSqrtReal(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
1272 PetscReal R_star = R/R_0;
1273 PetscReal tstar = PetscRealPart(c*time) / PetscSqr(R_0); /* - */
1274 PetscReal divA_n = 0.0;
1275
1276 if (R_star < PETSC_SMALL) {
1277 for (n = 1; n < N+1; ++n) {
1278 const PetscReal x_n = user->zeroArray[n-1];
1279 const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu)*PetscSqr(1 + nu_u)*x_n - 18.0*(1 + nu)*(nu_u - nu)*(1 - nu_u));
1280
1281 divA_n += PetscRealPart(
1282 (12.0*(1.0 + nu)*(nu_u - nu))/((1.0 - 2.0*nu)*E_n*PetscSqr(R_star)*x_n*PetscSinReal(PetscSqrtReal(x_n))) *
1283 (3.0*(nu_u - nu) * PetscSqrtReal(x_n) * ((2.0 + PetscSqr(R_star*PetscSqrtReal(x_n))) - 2.0*PetscCosReal(R_star * PetscSqrtReal(x_n)))
1284 + 5.0 * (1.0 - nu)*(1.0 - 2.0*nu)*PetscPowRealInt(R_star, 2)*x_n*PetscSinReal(PetscSqrtReal(x_n))) * PetscExpReal(-x_n * tstar));
1285 }
1286 } else {
1287 for (n = 1; n < N+1; ++n) {
1288 const PetscReal x_n = user->zeroArray[n-1];
1289 const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu)*PetscSqr(1 + nu_u)*x_n - 18.0*(1 + nu)*(nu_u - nu)*(1 - nu_u));
1290
1291 divA_n += PetscRealPart(
1292 (12.0*(1.0 + nu)*(nu_u - nu))/((1.0 - 2.0*nu)*E_n*PetscSqr(R_star)*x_n*PetscSinReal(PetscSqrtReal(x_n))) *
1293 (3.0*(nu_u - nu) * PetscSqrtReal(x_n) * ((2.0/(R_star*PetscSqrtReal(x_n)) + R_star*PetscSqrtReal(x_n))*PetscSinReal(R_star * PetscSqrtReal(x_n)) - 2.0*PetscCosReal(R_star * PetscSqrtReal(x_n)))
1294 + 5.0 * (1.0 - nu)*(1.0 - 2.0*nu)*PetscPowRealInt(R_star, 2)*x_n*PetscSinReal(PetscSqrtReal(x_n))) * PetscExpReal(-x_n * tstar));
1295 }
1296 }
1297 if (PetscAbsReal(divA_n) > 1e3) PetscPrintf(PETSC_COMM_SELF, "(%g, %g, %g) divA_n: %g\n", x[0], x[1], x[2], divA_n);
1298 u[0] = PetscRealPart(u_inf)/R_0 * (3.0 - divA_n);
1299 }
1300 return 0;
1301 }
1302
1303 // Pressure
cryer_3d_p(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * ctx)1304 static PetscErrorCode cryer_3d_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1305 {
1306 AppCtx *user = (AppCtx *) ctx;
1307 Parameter *param;
1308 PetscErrorCode ierr;
1309
1310 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
1311 if (time <= 0.0) {
1312 ierr = cryer_drainage_pressure(dim, time, x, Nc, u, ctx);CHKERRQ(ierr);
1313 } else {
1314 PetscScalar alpha = param->alpha; /* - */
1315 PetscScalar K_u = param->K_u; /* Pa */
1316 PetscScalar M = param->M; /* Pa */
1317 PetscScalar G = param->mu; /* Pa */
1318 PetscScalar P_0 = param->P_0; /* Pa */
1319 PetscReal R_0 = param->ymax; /* m */
1320 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
1321 PetscInt N = user->niter, n;
1322
1323 PetscScalar K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */
1324 PetscScalar eta = (3.0*alpha*G) / (3.0*K_d + 4.0*G); /* -, Cheng (B.11) */
1325 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G)); /* -, Cheng (B.8) */
1326 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); /* -, Cheng (B.9) */
1327 PetscScalar S = (3.0*K_u + 4.0*G) / (M*(3.0*K_d + 4.0*G)); /* Pa^{-1}, Cheng (B.14) */
1328 PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
1329 PetscScalar R = PetscSqrtReal(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
1330
1331 PetscScalar R_star = R / R_0;
1332 PetscScalar t_star = PetscRealPart(c * time) / PetscSqr(R_0);
1333 PetscReal A_x = 0.0;
1334
1335 for (n = 1; n < N+1; ++n) {
1336 const PetscReal x_n = user->zeroArray[n-1];
1337 const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu)*PetscSqr(1 + nu_u)*x_n - 18.0*(1 + nu)*(nu_u - nu)*(1 - nu_u));
1338
1339 A_x += PetscRealPart(((18.0*PetscSqr(nu_u - nu)) / (eta * E_n)) * (PetscSinReal(R_star * PetscSqrtReal(x_n)) / (R_star * PetscSinReal(PetscSqrtReal(x_n))) - 1.0) * PetscExpReal(-x_n * t_star)); /* Cheng (7.395) */
1340 }
1341 u[0] = P_0 * A_x;
1342 }
1343 return 0;
1344 }
1345
1346 /* Boundary Kernels */
f0_terzaghi_bd_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[],const PetscReal n[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[])1347 static void f0_terzaghi_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1348 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1349 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1350 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1351 {
1352 const PetscReal P = PetscRealPart(constants[5]);
1353
1354 f0[0] = 0.0;
1355 f0[1] = P;
1356 }
1357
f0_mandel_bd_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[],const PetscReal n[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[])1358 static void f0_mandel_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1359 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1360 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1361 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1362 {
1363 // Uniform stress distribution
1364 /* PetscScalar xmax = 0.5;
1365 PetscScalar xmin = -0.5;
1366 PetscScalar ymax = 0.5;
1367 PetscScalar ymin = -0.5;
1368 PetscScalar P = constants[5];
1369 PetscScalar aL = (xmax - xmin) / 2.0;
1370 PetscScalar sigma_zz = -1.0*P / aL; */
1371
1372 // Analytical (parabolic) stress distribution
1373 PetscReal a1, a2, am;
1374 PetscReal y1, y2, ym;
1375
1376 PetscInt NITER = 500;
1377 PetscReal EPS = 0.000001;
1378 PetscReal zeroArray[500]; /* NITER */
1379 PetscReal xmax = 1.0;
1380 PetscReal xmin = 0.0;
1381 PetscReal ymax = 0.1;
1382 PetscReal ymin = 0.0;
1383 PetscReal lower[2], upper[2];
1384
1385 lower[0] = xmin - (xmax - xmin) / 2.0;
1386 lower[1] = ymin - (ymax - ymin) / 2.0;
1387 upper[0] = xmax - (xmax - xmin) / 2.0;
1388 upper[1] = ymax - (ymax - ymin) / 2.0;
1389
1390 xmin = lower[0];
1391 ymin = lower[1];
1392 xmax = upper[0];
1393 ymax = upper[1];
1394
1395 PetscScalar G = constants[0];
1396 PetscScalar K_u = constants[1];
1397 PetscScalar alpha = constants[2];
1398 PetscScalar M = constants[3];
1399 PetscScalar kappa = constants[4];
1400 PetscScalar F = constants[5];
1401
1402 PetscScalar K_d = K_u - alpha*alpha*M;
1403 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G));
1404 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G));
1405 PetscReal aL = (xmax - xmin) / 2.0;
1406 PetscReal c = PetscRealPart(((2.0*kappa*G) * (1.0 - nu) * (nu_u - nu)) / (alpha*alpha * (1.0 - 2.0*nu) * (1.0 - nu_u)));
1407 PetscScalar B = (3.0 * (nu_u - nu)) / ( alpha * (1.0 - 2.0*nu) * (1.0 + nu_u));
1408 PetscScalar A1 = 3.0 / (B * (1.0 + nu_u));
1409 PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu);
1410
1411 // Generate zero values
1412 for (PetscInt i=1; i < NITER+1; i++)
1413 {
1414 a1 = ((PetscReal) i - 1.0) * PETSC_PI * PETSC_PI / 4.0 + EPS;
1415 a2 = a1 + PETSC_PI/2;
1416 for (PetscInt j=0; j<NITER; j++)
1417 {
1418 y1 = PetscTanReal(a1) - PetscRealPart(A1/A2)*a1;
1419 y2 = PetscTanReal(a2) - PetscRealPart(A1/A2)*a2;
1420 am = (a1 + a2)/2.0;
1421 ym = PetscTanReal(am) - PetscRealPart(A1/A2)*am;
1422 if ((ym*y1) > 0)
1423 {
1424 a1 = am;
1425 } else {
1426 a2 = am;
1427 }
1428 if (PetscAbsReal(y2) < EPS)
1429 {
1430 am = a2;
1431 }
1432 }
1433 zeroArray[i-1] = am;
1434 }
1435
1436 // Solution for sigma_zz
1437 PetscScalar A_x = 0.0;
1438 PetscScalar B_x = 0.0;
1439
1440 for (PetscInt n=1; n < NITER+1; n++)
1441 {
1442 PetscReal alpha_n = zeroArray[n-1];
1443
1444 A_x += ( PetscSinReal(alpha_n) / (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscCosReal( (alpha_n * x[0]) / aL) * PetscExpReal( -1.0*( (alpha_n*alpha_n*c*t)/(aL*aL)));
1445 B_x += ( (PetscSinReal(alpha_n) * PetscCosReal(alpha_n))/(alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscExpReal( -1.0*( (alpha_n*alpha_n*c*t)/(aL*aL)));
1446 }
1447
1448 PetscScalar sigma_zz = -1.0*(F/aL) - ((2.0*F)/aL) * (A2/A1) * A_x + ((2.0*F)/aL) * B_x;
1449
1450
1451 if (x[1] == ymax) {
1452 f0[1] += sigma_zz;
1453 } else if (x[1] == ymin) {
1454 f0[1] -= sigma_zz;
1455 }
1456 }
1457
f0_cryer_bd_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[],const PetscReal n[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[])1458 static void f0_cryer_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1459 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1460 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1461 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1462 {
1463 const PetscReal P_0 = PetscRealPart(constants[5]);
1464 const PetscReal R = PetscSqrtReal(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
1465 PetscInt d;
1466
1467 for (d = 0; d < dim; ++d) f0[d] = -P_0*n[d];
1468 //PetscPrintf(PETSC_COMM_SELF, "R: %g P_0: %g n: (%g, %g, %g) hat n (%g, %g, %g)\n", R, P_0, n[0], n[1], n[2], x[0]/R, x[1]/R, x[2]/R);
1469 for (d = 0; d < dim; ++d) if (PetscAbsReal(n[d] - x[d]/R) > 1.0) PetscPrintf(PETSC_COMM_SELF, "WTF? R: %g P_0: %g n: (%g, %g, %g) hat n (%g, %g, %g)\n", R, P_0, n[0], n[1], n[2], x[0]/R, x[1]/R, x[2]/R);
1470 //for (d = 0; d < dim; ++d) f0[d] = -P_0*x[d]/R;
1471 }
1472
1473 /* Standard Kernels - Residual */
1474 /* f0_e */
f0_epsilon(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 f0[])1475 static void f0_epsilon(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1476 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1477 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1478 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1479 {
1480 PetscInt d;
1481
1482 for (d = 0; d < dim; ++d) {
1483 f0[0] += u_x[d*dim+d];
1484 }
1485 f0[0] -= u[uOff[1]];
1486 }
1487
1488 /* f0_p */
f0_p(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 f0[])1489 static void f0_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1490 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1491 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1492 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1493 {
1494 const PetscReal alpha = PetscRealPart(constants[2]);
1495 const PetscReal M = PetscRealPart(constants[3]);
1496
1497 f0[0] += alpha*u_t[uOff[1]];
1498 f0[0] += u_t[uOff[2]]/M;
1499 }
1500
1501 /* f1_u */
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[])1502 static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1503 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1504 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1505 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
1506 {
1507 const PetscInt Nc = dim;
1508 const PetscReal G = PetscRealPart(constants[0]);
1509 const PetscReal K_u = PetscRealPart(constants[1]);
1510 const PetscReal alpha = PetscRealPart(constants[2]);
1511 const PetscReal M = PetscRealPart(constants[3]);
1512 const PetscReal K_d = K_u - alpha*alpha*M;
1513 const PetscReal lambda = K_d - (2.0 * G) / 3.0;
1514 PetscInt c, d;
1515
1516 for (c = 0; c < Nc; ++c)
1517 {
1518 for (d = 0; d < dim; ++d)
1519 {
1520 f1[c*dim+d] -= G*(u_x[c*dim+d] + u_x[d*dim+c]);
1521 }
1522 f1[c*dim+c] -= lambda*u[uOff[1]];
1523 f1[c*dim+c] += alpha*u[uOff[2]];
1524 }
1525 }
1526
1527 /* f1_p */
f1_p(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[])1528 static void f1_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1529 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1530 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1531 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
1532 {
1533 const PetscReal kappa = PetscRealPart(constants[4]);
1534 PetscInt d;
1535
1536 for (d = 0; d < dim; ++d) {
1537 f1[d] += kappa*u_x[uOff_x[2]+d];
1538 }
1539 }
1540
1541 /*
1542 \partial_df \phi_fc g_{fc,gc,df,dg} \partial_dg \phi_gc
1543
1544 \partial_df \phi_fc \lambda \delta_{fc,df} \sum_gc \partial_dg \phi_gc \delta_{gc,dg}
1545 = \partial_fc \phi_fc \sum_gc \partial_gc \phi_gc
1546 */
1547
1548
1549 /* Standard Kernels - Jacobian */
1550 /* g0_ee */
g0_ee(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 g0[])1551 static void g0_ee(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1552 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1553 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1554 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1555 {
1556 g0[0] = -1.0;
1557 }
1558
1559 /* g0_pe */
g0_pe(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 g0[])1560 static void g0_pe(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1561 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1562 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1563 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1564 {
1565 const PetscReal alpha = PetscRealPart(constants[2]);
1566
1567 g0[0] = u_tShift*alpha;
1568 }
1569
1570 /* g0_pp */
g0_pp(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 g0[])1571 static void g0_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1572 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1573 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1574 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1575 {
1576 const PetscReal M = PetscRealPart(constants[3]);
1577
1578 g0[0] = u_tShift/M;
1579 }
1580
1581 /* g1_eu */
g1_eu(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 g1[])1582 static void g1_eu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1583 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1584 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1585 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
1586 {
1587 PetscInt d;
1588 for (d = 0; d < dim; ++d) g1[d*dim+d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */
1589 }
1590
1591 /* g2_ue */
g2_ue(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 g2[])1592 static void g2_ue(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1593 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1594 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1595 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
1596 {
1597 const PetscReal G = PetscRealPart(constants[0]);
1598 const PetscReal K_u = PetscRealPart(constants[1]);
1599 const PetscReal alpha = PetscRealPart(constants[2]);
1600 const PetscReal M = PetscRealPart(constants[3]);
1601 const PetscReal K_d = K_u - alpha*alpha*M;
1602 const PetscReal lambda = K_d - (2.0 * G) / 3.0;
1603 PetscInt d;
1604
1605 for (d = 0; d < dim; ++d) {
1606 g2[d*dim + d] -= lambda;
1607 }
1608 }
1609 /* g2_up */
g2_up(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 g2[])1610 static void g2_up(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1611 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1612 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1613 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
1614 {
1615 const PetscReal alpha = PetscRealPart(constants[2]);
1616 PetscInt d;
1617
1618 for (d = 0; d < dim; ++d) {
1619 g2[d*dim + d] += alpha;
1620 }
1621 }
1622
1623 /* g3_uu */
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[])1624 static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1625 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1626 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1627 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
1628 {
1629 const PetscInt Nc = dim;
1630 const PetscReal G = PetscRealPart(constants[0]);
1631 PetscInt c, d;
1632
1633 for (c = 0; c < Nc; ++c) {
1634 for (d = 0; d < dim; ++d) {
1635 g3[((c*Nc + c)*dim + d)*dim + d] -= G;
1636 g3[((c*Nc + d)*dim + d)*dim + c] -= G;
1637 }
1638 }
1639 }
1640
1641 /* g3_pp */
g3_pp(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[])1642 static void g3_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1643 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1644 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1645 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
1646 {
1647 const PetscReal kappa = PetscRealPart(constants[4]);
1648 PetscInt d;
1649
1650 for (d = 0; d < dim; ++d) g3[d*dim+d] += kappa;
1651 }
1652
ProcessOptions(MPI_Comm comm,AppCtx * options)1653 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
1654 {
1655 PetscInt sol;
1656 PetscErrorCode ierr;
1657
1658 PetscFunctionBeginUser;
1659 options->dim = 2;
1660 options->simplex = PETSC_TRUE;
1661 options->refLimit = -1.0;
1662 options->solType = SOL_QUADRATIC_TRIG;
1663 options->niter = 500;
1664 options->eps = PETSC_SMALL;
1665 options->dtInitial = -1.0;
1666 ierr = PetscStrncpy(options->dmType, DMPLEX, 256);CHKERRQ(ierr);
1667
1668 ierr = PetscOptionsBegin(comm, "", "Biot Poroelasticity Options", "DMPLEX");CHKERRQ(ierr);
1669 ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex53.c", options->dim, &options->dim, NULL);CHKERRQ(ierr);
1670 ierr = PetscOptionsInt("-niter", "Number of series term iterations in exact solutions", "ex53.c", options->niter, &options->niter, NULL);CHKERRQ(ierr);
1671 ierr = PetscOptionsBool("-simplex", "Simplicial (true) or tensor (false) mesh", "ex53.c", options->simplex, &options->simplex, NULL);CHKERRQ(ierr);
1672 ierr = PetscOptionsReal("-ref_limit", "Maximum cell volume for refined mesh", "ex53.c", options->refLimit, &options->refLimit, NULL);CHKERRQ(ierr);
1673 sol = options->solType;
1674 ierr = PetscOptionsEList("-sol_type", "Type of exact solution", "ex53.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL);CHKERRQ(ierr);
1675 options->solType = (SolutionType) sol;
1676 ierr = PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex53.c", DMList, options->dmType, options->dmType, 256, NULL);CHKERRQ(ierr);
1677 ierr = PetscOptionsReal("-eps", "Precision value for root finding", "ex53.c", options->eps, &options->eps, NULL);CHKERRQ(ierr);
1678 ierr = PetscOptionsReal("-dt_initial", "Override the initial timestep", "ex53.c", options->dtInitial, &options->dtInitial, NULL);CHKERRQ(ierr);
1679
1680 // Wrap up loose ends
1681 if (options->solType == SOL_CRYER) {
1682 options->dim = 3;
1683 }
1684
1685 ierr = PetscOptionsEnd();CHKERRQ(ierr);
1686 PetscFunctionReturn(0);
1687 }
1688
mandelZeros(MPI_Comm comm,AppCtx * ctx,Parameter * param)1689 static PetscErrorCode mandelZeros(MPI_Comm comm, AppCtx *ctx, Parameter *param)
1690 {
1691 //PetscBag bag;
1692 PetscReal a1, a2, am;
1693 PetscReal y1, y2, ym;
1694
1695 PetscFunctionBeginUser;
1696 //ierr = PetscBagGetData(ctx->bag, (void **) ¶m);CHKERRQ(ierr);
1697 PetscInt NITER = ctx->niter;
1698 PetscReal EPS = ctx->eps;
1699 //const PetscScalar YMAX = param->ymax;
1700 //const PetscScalar YMIN = param->ymin;
1701 PetscScalar alpha = param->alpha;
1702 PetscScalar K_u = param->K_u;
1703 PetscScalar M = param->M;
1704 PetscScalar G = param->mu;
1705 //const PetscScalar k = param->k;
1706 //const PetscScalar mu_f = param->mu_f;
1707 //const PetscScalar P_0 = param->P_0;
1708
1709 PetscScalar K_d = K_u - alpha*alpha*M;
1710 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G));
1711 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G));
1712 //const PetscScalar kappa = k / mu_f;
1713
1714 // Generate zero values
1715 for (PetscInt i=1; i < NITER+1; i++)
1716 {
1717 a1 = ((PetscReal) i - 1.0) * PETSC_PI * PETSC_PI / 4.0 + EPS;
1718 a2 = a1 + PETSC_PI/2;
1719 am = a1;
1720 for (PetscInt j=0; j<NITER; j++)
1721 {
1722 y1 = PetscTanReal(a1) - PetscRealPart((1.0 - nu)/(nu_u - nu))*a1;
1723 y2 = PetscTanReal(a2) - PetscRealPart((1.0 - nu)/(nu_u - nu))*a2;
1724 am = (a1 + a2)/2.0;
1725 ym = PetscTanReal(am) - PetscRealPart((1.0 - nu)/(nu_u - nu))*am;
1726 if ((ym*y1) > 0)
1727 {
1728 a1 = am;
1729 } else {
1730 a2 = am;
1731 }
1732 if (PetscAbsReal(y2) < EPS)
1733 {
1734 am = a2;
1735 }
1736 }
1737 ctx->zeroArray[i-1] = am;
1738 }
1739 PetscFunctionReturn(0);
1740 }
1741
CryerFunction(PetscReal nu_u,PetscReal nu,PetscReal x)1742 static PetscReal CryerFunction(PetscReal nu_u, PetscReal nu, PetscReal x)
1743 {
1744 return PetscTanReal(PetscSqrtReal(x))*(6.0*(nu_u - nu) - (1.0 - nu)*(1.0 + nu_u)*x) - (6.0*(nu_u - nu)*PetscSqrtReal(x));
1745 }
1746
cryerZeros(MPI_Comm comm,AppCtx * ctx,Parameter * param)1747 static PetscErrorCode cryerZeros(MPI_Comm comm, AppCtx *ctx, Parameter *param)
1748 {
1749 PetscReal alpha = PetscRealPart(param->alpha); /* - */
1750 PetscReal K_u = PetscRealPart(param->K_u); /* Pa */
1751 PetscReal M = PetscRealPart(param->M); /* Pa */
1752 PetscReal G = PetscRealPart(param->mu); /* Pa */
1753 PetscInt N = ctx->niter, n;
1754
1755 PetscReal K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */
1756 PetscReal nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G)); /* -, Cheng (B.8) */
1757 PetscReal nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); /* -, Cheng (B.9) */
1758
1759 PetscFunctionBeginUser;
1760 for (n = 1; n < N+1; ++n) {
1761 PetscReal tol = PetscPowReal(n, 1.5)*ctx->eps;
1762 PetscReal a1 = 0., a2 = 0., am = 0.;
1763 PetscReal y1, y2, ym;
1764 PetscInt j, k = n-1;
1765
1766 y1 = y2 = 1.;
1767 while (y1*y2 > 0) {
1768 ++k;
1769 a1 = PetscSqr(n*PETSC_PI) - k*PETSC_PI;
1770 a2 = PetscSqr(n*PETSC_PI) + k*PETSC_PI;
1771 y1 = CryerFunction(nu_u, nu, a1);
1772 y2 = CryerFunction(nu_u, nu, a2);
1773 }
1774 for (j = 0; j < 50000; ++j) {
1775 y1 = CryerFunction(nu_u, nu, a1);
1776 y2 = CryerFunction(nu_u, nu, a2);
1777 if (y1*y2 > 0) SETERRQ5(comm, PETSC_ERR_PLIB, "Invalid root finding initialization for root %D, (%g, %g)--(%g, %g)", n, a1, y1, a2, y2);
1778 am = (a1 + a2) / 2.0;
1779 ym = CryerFunction(nu_u, nu, am);
1780 if ((ym * y1) < 0) a2 = am;
1781 else a1 = am;
1782 if (PetscAbsScalar(ym) < tol) break;
1783 }
1784 if (PetscAbsScalar(ym) >= tol) SETERRQ2(comm, PETSC_ERR_PLIB, "Root finding did not converge for root %D (%g)", n, PetscAbsScalar(ym));
1785 ctx->zeroArray[n-1] = am;
1786 }
1787 PetscFunctionReturn(0);
1788 }
1789
SetupParameters(MPI_Comm comm,AppCtx * ctx)1790 static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
1791 {
1792 PetscBag bag;
1793 Parameter *p;
1794 PetscErrorCode ierr;
1795
1796 PetscFunctionBeginUser;
1797 /* setup PETSc parameter bag */
1798 ierr = PetscBagGetData(ctx->bag,(void**)&p);CHKERRQ(ierr);
1799 ierr = PetscBagSetName(ctx->bag,"par","Poroelastic Parameters");CHKERRQ(ierr);
1800 bag = ctx->bag;
1801 if (ctx->solType == SOL_TERZAGHI) {
1802 // Realistic values - Terzaghi
1803 ierr = PetscBagRegisterScalar(bag, &p->mu, 3.0, "mu", "Shear Modulus, Pa");CHKERRQ(ierr);
1804 ierr = PetscBagRegisterScalar(bag, &p->K_u, 9.76, "K_u", "Undrained Bulk Modulus, Pa");CHKERRQ(ierr);
1805 ierr = PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -");CHKERRQ(ierr);
1806 ierr = PetscBagRegisterScalar(bag, &p->M, 16.0, "M", "Biot Modulus, Pa");CHKERRQ(ierr);
1807 ierr = PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2");CHKERRQ(ierr);
1808 ierr = PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s");CHKERRQ(ierr);
1809 ierr = PetscBagRegisterScalar(bag, &p->zmax, 1.0, "zmax", "Depth Maximum Extent, m");CHKERRQ(ierr);
1810 ierr = PetscBagRegisterScalar(bag, &p->zmin, 0.0, "zmin", "Depth Minimum Extent, m");CHKERRQ(ierr);
1811 ierr = PetscBagRegisterScalar(bag, &p->ymax, 10.0, "ymax", "Vertical Maximum Extent, m");CHKERRQ(ierr);
1812 ierr = PetscBagRegisterScalar(bag, &p->ymin, 0.0, "ymin", "Vertical Minimum Extent, m");CHKERRQ(ierr);
1813 ierr = PetscBagRegisterScalar(bag, &p->xmax, 10.0, "xmax", "Horizontal Maximum Extent, m");CHKERRQ(ierr);
1814 ierr = PetscBagRegisterScalar(bag, &p->xmin, 0.0, "xmin", "Horizontal Minimum Extent, m");CHKERRQ(ierr);
1815 ierr = PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa");CHKERRQ(ierr);
1816 } else if (ctx->solType == SOL_MANDEL) {
1817 // Realistic values - Mandel
1818 ierr = PetscBagRegisterScalar(bag, &p->mu, 0.75, "mu", "Shear Modulus, Pa");CHKERRQ(ierr);
1819 ierr = PetscBagRegisterScalar(bag, &p->K_u, 2.6941176470588233, "K_u", "Undrained Bulk Modulus, Pa");CHKERRQ(ierr);
1820 ierr = PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -");CHKERRQ(ierr);
1821 ierr = PetscBagRegisterScalar(bag, &p->M, 4.705882352941176, "M", "Biot Modulus, Pa");CHKERRQ(ierr);
1822 ierr = PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2");CHKERRQ(ierr);
1823 ierr = PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s");CHKERRQ(ierr);
1824 ierr = PetscBagRegisterScalar(bag, &p->zmax, 1.0, "zmax", "Depth Maximum Extent, m");CHKERRQ(ierr);
1825 ierr = PetscBagRegisterScalar(bag, &p->zmin, 0.0, "zmin", "Depth Minimum Extent, m");CHKERRQ(ierr);
1826 ierr = PetscBagRegisterScalar(bag, &p->ymax, 0.25, "ymax", "Vertical Maximum Extent, m");CHKERRQ(ierr);
1827 ierr = PetscBagRegisterScalar(bag, &p->ymin, 0.0, "ymin", "Vertical Minimum Extent, m");CHKERRQ(ierr);
1828 ierr = PetscBagRegisterScalar(bag, &p->xmax, 1.0, "xmax", "Horizontal Maximum Extent, m");CHKERRQ(ierr);
1829 ierr = PetscBagRegisterScalar(bag, &p->xmin, 0.0, "xmin", "Horizontal Minimum Extent, m");CHKERRQ(ierr);
1830 ierr = PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa");CHKERRQ(ierr);
1831 } else if (ctx->solType == SOL_CRYER) {
1832 // Realistic values - Mandel
1833 ierr = PetscBagRegisterScalar(bag, &p->mu, 0.75, "mu", "Shear Modulus, Pa");CHKERRQ(ierr);
1834 ierr = PetscBagRegisterScalar(bag, &p->K_u, 2.6941176470588233, "K_u", "Undrained Bulk Modulus, Pa");CHKERRQ(ierr);
1835 ierr = PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -");CHKERRQ(ierr);
1836 ierr = PetscBagRegisterScalar(bag, &p->M, 4.705882352941176, "M", "Biot Modulus, Pa");CHKERRQ(ierr);
1837 ierr = PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2");CHKERRQ(ierr);
1838 ierr = PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s");CHKERRQ(ierr);
1839 ierr = PetscBagRegisterScalar(bag, &p->zmax, 1.0, "zmax", "Depth Maximum Extent, m");CHKERRQ(ierr);
1840 ierr = PetscBagRegisterScalar(bag, &p->zmin, 0.0, "zmin", "Depth Minimum Extent, m");CHKERRQ(ierr);
1841 ierr = PetscBagRegisterScalar(bag, &p->ymax, 1.0, "ymax", "Vertical Maximum Extent, m");CHKERRQ(ierr);
1842 ierr = PetscBagRegisterScalar(bag, &p->ymin, 0.0, "ymin", "Vertical Minimum Extent, m");CHKERRQ(ierr);
1843 ierr = PetscBagRegisterScalar(bag, &p->xmax, 1.0, "xmax", "Horizontal Maximum Extent, m");CHKERRQ(ierr);
1844 ierr = PetscBagRegisterScalar(bag, &p->xmin, 0.0, "xmin", "Horizontal Minimum Extent, m");CHKERRQ(ierr);
1845 ierr = PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa");CHKERRQ(ierr);
1846 } else {
1847 // Nonsense values
1848 ierr = PetscBagRegisterScalar(bag, &p->mu, 1.0, "mu", "Shear Modulus, Pa");CHKERRQ(ierr);
1849 ierr = PetscBagRegisterScalar(bag, &p->K_u, 1.0, "K_u", "Undrained Bulk Modulus, Pa");CHKERRQ(ierr);
1850 ierr = PetscBagRegisterScalar(bag, &p->alpha, 1.0, "alpha", "Biot Effective Stress Coefficient, -");CHKERRQ(ierr);
1851 ierr = PetscBagRegisterScalar(bag, &p->M, 1.0, "M", "Biot Modulus, Pa");CHKERRQ(ierr);
1852 ierr = PetscBagRegisterScalar(bag, &p->k, 1.0, "k", "Isotropic Permeability, m**2");CHKERRQ(ierr);
1853 ierr = PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s");CHKERRQ(ierr);
1854 ierr = PetscBagRegisterScalar(bag, &p->zmax, 1.0, "zmax", "Depth Maximum Extent, m");CHKERRQ(ierr);
1855 ierr = PetscBagRegisterScalar(bag, &p->zmin, 0.0, "zmin", "Depth Minimum Extent, m");CHKERRQ(ierr);
1856 ierr = PetscBagRegisterScalar(bag, &p->ymax, 1.0, "ymax", "Vertical Maximum Extent, m");CHKERRQ(ierr);
1857 ierr = PetscBagRegisterScalar(bag, &p->ymin, 0.0, "ymin", "Vertical Minimum Extent, m");CHKERRQ(ierr);
1858 ierr = PetscBagRegisterScalar(bag, &p->xmax, 1.0, "xmax", "Horizontal Maximum Extent, m");CHKERRQ(ierr);
1859 ierr = PetscBagRegisterScalar(bag, &p->xmin, 0.0, "xmin", "Horizontal Minimum Extent, m");CHKERRQ(ierr);
1860 ierr = PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa");CHKERRQ(ierr);
1861 }
1862 ierr = PetscBagSetFromOptions(bag);CHKERRQ(ierr);
1863 {
1864 PetscScalar K_d = p->K_u - p->alpha*p->alpha*p->M;
1865 PetscScalar nu_u = (3.0*p->K_u - 2.0*p->mu) / (2.0*(3.0*p->K_u + p->mu));
1866 PetscScalar nu = (3.0*K_d - 2.0*p->mu) / (2.0*(3.0*K_d + p->mu));
1867 PetscScalar S = (3.0*p->K_u + 4.0*p->mu) / (p->M*(3.0*K_d + 4.0*p->mu));
1868 PetscReal c = PetscRealPart((p->k/p->mu_f) / S);
1869
1870 PetscViewer viewer;
1871 PetscViewerFormat format;
1872 PetscBool flg;
1873
1874 switch (ctx->solType) {
1875 case SOL_QUADRATIC_LINEAR:
1876 case SOL_QUADRATIC_TRIG:
1877 case SOL_TRIG_LINEAR: ctx->t_r = PetscSqr(p->xmax - p->xmin)/c; break;
1878 case SOL_TERZAGHI: ctx->t_r = PetscSqr(2.0*(p->ymax - p->ymin))/c; break;
1879 case SOL_MANDEL: ctx->t_r = PetscSqr(2.0*(p->ymax - p->ymin))/c; break;
1880 case SOL_CRYER: ctx->t_r = PetscSqr(p->ymax)/c; break;
1881 default: SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%D)", solutionTypes[PetscMin(ctx->solType, NUM_SOLUTION_TYPES)], ctx->solType);
1882 }
1883 ierr = PetscOptionsGetViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg);CHKERRQ(ierr);
1884 if (flg) {
1885 ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr);
1886 ierr = PetscBagView(bag, viewer);CHKERRQ(ierr);
1887 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
1888 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
1889 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
1890 ierr = PetscPrintf(comm, " Max displacement: %g %g\n", p->P_0*(p->ymax - p->ymin)*(1. - 2.*nu_u)/(2.*p->mu*(1. - nu_u)), p->P_0*(p->ymax - p->ymin)*(1. - 2.*nu)/(2.*p->mu*(1. - nu)));
1891 ierr = PetscPrintf(comm, " Relaxation time: %g\n", ctx->t_r);
1892 }
1893 }
1894 PetscFunctionReturn(0);
1895 }
1896
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)1897 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
1898 {
1899 Parameter *param;
1900 PetscBool flg;
1901 PetscErrorCode ierr;
1902
1903 PetscFunctionBeginUser;
1904 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
1905 if (user->solType == SOL_CRYER) {
1906 DM rdm;
1907
1908 if (!user->simplex) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot create ball with cubic cells");
1909 if (param->xmin != 0.0 || param->ymin != 0.0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Cannot shift center of ball to (%g, %g)", param->xmin, param->ymin);
1910 if (param->xmax != param->ymax) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Cannot radius of ball must be equal in x and y: %g != %g", param->xmax, param->ymax);
1911 ierr = DMPlexCreateBallMesh(comm, user->dim, param->xmax, dm);CHKERRQ(ierr);
1912
1913 ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
1914 ierr = DMPlexSetRefinementLimit(*dm, user->refLimit);CHKERRQ(ierr);
1915 ierr = DMRefine(*dm, comm, &rdm);CHKERRQ(ierr);
1916 if (rdm) {
1917 ierr = DMDestroy(dm);CHKERRQ(ierr);
1918 *dm = rdm;
1919 }
1920 ierr = DMPlexSetRefinementUniform(*dm, PETSC_TRUE);CHKERRQ(ierr);
1921 } else if (user->solType == SOL_MANDEL) {
1922 PetscReal lower[2], upper[2];
1923
1924 lower[0] = param->xmin - (param->xmax - param->xmin) / 2.0;
1925 lower[1] = param->ymin - (param->ymax - param->ymin) / 2.0;
1926 upper[0] = param->xmax - (param->xmax - param->xmin) / 2.0;
1927 upper[1] = param->ymax - (param->ymax - param->ymin) / 2.0;
1928 //reset min / max values for mandel
1929 param->xmin = lower[0];
1930 param->ymin = lower[1];
1931 param->xmax = upper[0];
1932 param->ymax = upper[1];
1933 ierr = DMPlexCreateBoxMesh(comm, user->dim, user->simplex, NULL, lower, upper, NULL, PETSC_TRUE, dm);CHKERRQ(ierr);
1934 } else {
1935 Parameter *param;
1936 PetscReal lower[3], upper[3];
1937
1938 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
1939 lower[0] = param->xmin;
1940 lower[1] = param->ymin;
1941 lower[2] = param->zmin;
1942 upper[0] = param->xmax;
1943 upper[1] = param->ymax;
1944 upper[2] = param->zmax;
1945 ierr = DMPlexCreateBoxMesh(comm, user->dim, user->simplex, NULL, lower, upper, NULL, PETSC_TRUE, dm);CHKERRQ(ierr);
1946 }
1947 {
1948 DM pdm = NULL;
1949 PetscPartitioner part;
1950
1951 ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr);
1952 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr);
1953 ierr = DMPlexDistribute(*dm, 0, NULL, &pdm);CHKERRQ(ierr);
1954 if (pdm) {
1955 ierr = DMDestroy(dm);CHKERRQ(ierr);
1956 *dm = pdm;
1957 }
1958 }
1959 ierr = PetscStrcmp(user->dmType, DMPLEX, &flg);CHKERRQ(ierr);
1960 if (flg) {
1961 DM ndm;
1962
1963 ierr = DMConvert(*dm, user->dmType, &ndm);CHKERRQ(ierr);
1964 if (ndm) {
1965 ierr = DMDestroy(dm);CHKERRQ(ierr);
1966 *dm = ndm;
1967 }
1968 }
1969 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr);
1970
1971 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr);
1972 ierr = DMSetApplicationContext(*dm, user);CHKERRQ(ierr);
1973 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
1974 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
1975 PetscFunctionReturn(0);
1976 }
1977
SetupPrimalProblem(DM dm,AppCtx * user)1978 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
1979 {
1980 PetscErrorCode (*exact[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
1981 PetscErrorCode (*exact_t[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
1982 PetscDS prob;
1983 Parameter *param;
1984 PetscInt id_mandel[2];
1985 PetscInt comp[1];
1986 PetscInt comp_mandel[2];
1987 PetscInt dim, id, f;
1988 PetscErrorCode ierr;
1989
1990 PetscFunctionBeginUser;
1991 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
1992 ierr = PetscDSGetSpatialDimension(prob, &dim);CHKERRQ(ierr);
1993 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr);
1994 exact_t[0] = exact_t[1] = exact_t[2] = zero;
1995
1996 /* Setup Problem Formulation and Boundary Conditions */
1997 switch (user->solType) {
1998 case SOL_QUADRATIC_LINEAR:
1999 ierr = PetscDSSetResidual(prob, 0, f0_quadratic_linear_u, f1_u);CHKERRQ(ierr);
2000 ierr = PetscDSSetResidual(prob, 1, f0_epsilon, NULL);CHKERRQ(ierr);
2001 ierr = PetscDSSetResidual(prob, 2, f0_quadratic_linear_p, f1_p);CHKERRQ(ierr);
2002 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr);
2003 ierr = PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_ue, NULL);CHKERRQ(ierr);
2004 ierr = PetscDSSetJacobian(prob, 0, 2, NULL, NULL, g2_up, NULL);CHKERRQ(ierr);
2005 ierr = PetscDSSetJacobian(prob, 1, 0, NULL, g1_eu, NULL, NULL);CHKERRQ(ierr);
2006 ierr = PetscDSSetJacobian(prob, 1, 1, g0_ee, NULL, NULL, NULL);CHKERRQ(ierr);
2007 ierr = PetscDSSetJacobian(prob, 2, 1, g0_pe, NULL, NULL, NULL);CHKERRQ(ierr);
2008 ierr = PetscDSSetJacobian(prob, 2, 2, g0_pp, NULL, NULL, g3_pp);CHKERRQ(ierr);
2009 exact[0] = quadratic_u;
2010 exact[1] = linear_eps;
2011 exact[2] = linear_linear_p;
2012 exact_t[2] = linear_linear_p_t;
2013
2014 id = 1;
2015 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", "marker", 0, 0, NULL, (void (*)(void)) exact[0], NULL, 1, &id, user);CHKERRQ(ierr);
2016 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", "marker", 2, 0, NULL, (void (*)(void)) exact[2], (void (*)(void)) exact_t[2], 1, &id, user);CHKERRQ(ierr);
2017 break;
2018 case SOL_TRIG_LINEAR:
2019 ierr = PetscDSSetResidual(prob, 0, f0_trig_linear_u, f1_u);CHKERRQ(ierr);
2020 ierr = PetscDSSetResidual(prob, 1, f0_epsilon, NULL);CHKERRQ(ierr);
2021 ierr = PetscDSSetResidual(prob, 2, f0_trig_linear_p, f1_p);CHKERRQ(ierr);
2022 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr);
2023 ierr = PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_ue, NULL);CHKERRQ(ierr);
2024 ierr = PetscDSSetJacobian(prob, 0, 2, NULL, NULL, g2_up, NULL);CHKERRQ(ierr);
2025 ierr = PetscDSSetJacobian(prob, 1, 0, NULL, g1_eu, NULL, NULL);CHKERRQ(ierr);
2026 ierr = PetscDSSetJacobian(prob, 1, 1, g0_ee, NULL, NULL, NULL);CHKERRQ(ierr);
2027 ierr = PetscDSSetJacobian(prob, 2, 1, g0_pe, NULL, NULL, NULL);CHKERRQ(ierr);
2028 ierr = PetscDSSetJacobian(prob, 2, 2, g0_pp, NULL, NULL, g3_pp);CHKERRQ(ierr);
2029 exact[0] = trig_u;
2030 exact[1] = trig_eps;
2031 exact[2] = trig_linear_p;
2032 exact_t[2] = trig_linear_p_t;
2033
2034 id = 1;
2035 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", "marker", 0, 0, NULL, (void (*)(void)) exact[0], NULL, 1, &id, user);CHKERRQ(ierr);
2036 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", "marker", 2, 0, NULL, (void (*)(void)) exact[2], (void (*)(void)) exact_t[2], 1, &id, user);CHKERRQ(ierr);
2037 break;
2038 case SOL_QUADRATIC_TRIG:
2039 ierr = PetscDSSetResidual(prob, 0, f0_quadratic_trig_u, f1_u);CHKERRQ(ierr);
2040 ierr = PetscDSSetResidual(prob, 1, f0_epsilon, NULL);CHKERRQ(ierr);
2041 ierr = PetscDSSetResidual(prob, 2, f0_quadratic_trig_p, f1_p);CHKERRQ(ierr);
2042 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr);
2043 ierr = PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_ue, NULL);CHKERRQ(ierr);
2044 ierr = PetscDSSetJacobian(prob, 0, 2, NULL, NULL, g2_up, NULL);CHKERRQ(ierr);
2045 ierr = PetscDSSetJacobian(prob, 1, 0, NULL, g1_eu, NULL, NULL);CHKERRQ(ierr);
2046 ierr = PetscDSSetJacobian(prob, 1, 1, g0_ee, NULL, NULL, NULL);CHKERRQ(ierr);
2047 ierr = PetscDSSetJacobian(prob, 2, 1, g0_pe, NULL, NULL, NULL);CHKERRQ(ierr);
2048 ierr = PetscDSSetJacobian(prob, 2, 2, g0_pp, NULL, NULL, g3_pp);CHKERRQ(ierr);
2049 exact[0] = quadratic_u;
2050 exact[1] = linear_eps;
2051 exact[2] = linear_trig_p;
2052 exact_t[2] = linear_trig_p_t;
2053
2054 id = 1;
2055 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", "marker", 0, 0, NULL, (void (*)(void)) exact[0], NULL, 1, &id, user);CHKERRQ(ierr);
2056 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", "marker", 2, 0, NULL, (void (*)(void)) exact[2], (void (*)(void)) exact_t[2], 1, &id, user);CHKERRQ(ierr);
2057 break;
2058 case SOL_TERZAGHI:
2059 ierr = PetscDSSetResidual(prob, 0, NULL, f1_u);CHKERRQ(ierr);
2060 ierr = PetscDSSetResidual(prob, 1, f0_epsilon, NULL);CHKERRQ(ierr);
2061 ierr = PetscDSSetResidual(prob, 2, f0_p, f1_p);CHKERRQ(ierr);
2062 ierr = PetscDSSetBdResidual(prob, 0, f0_terzaghi_bd_u, NULL);CHKERRQ(ierr);
2063 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr);
2064 ierr = PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_ue, NULL);CHKERRQ(ierr);
2065 ierr = PetscDSSetJacobian(prob, 0, 2, NULL, NULL, g2_up, NULL);CHKERRQ(ierr);
2066 ierr = PetscDSSetJacobian(prob, 1, 0, NULL, g1_eu, NULL, NULL);CHKERRQ(ierr);
2067 ierr = PetscDSSetJacobian(prob, 1, 1, g0_ee, NULL, NULL, NULL);CHKERRQ(ierr);
2068 ierr = PetscDSSetJacobian(prob, 2, 1, g0_pe, NULL, NULL, NULL);CHKERRQ(ierr);
2069 ierr = PetscDSSetJacobian(prob, 2, 2, g0_pp, NULL, NULL, g3_pp);CHKERRQ(ierr);
2070
2071 exact[0] = terzaghi_2d_u;
2072 exact[1] = terzaghi_2d_eps;
2073 exact[2] = terzaghi_2d_p;
2074 exact_t[0] = terzaghi_2d_u_t;
2075 exact_t[1] = terzaghi_2d_eps_t;
2076 exact_t[2] = terzaghi_2d_p_t;
2077
2078 id = 1;
2079 ierr = DMAddBoundary(dm, DM_BC_NATURAL, "vertical stress", "marker", 0, 0, NULL, NULL, NULL, 1, &id, user);CHKERRQ(ierr);
2080 id = 3;
2081 comp[0] = 1;
2082 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed base", "marker", 0, 1, comp, (void (*)(void)) zero, NULL, 1, &id, user);CHKERRQ(ierr);
2083 id = 2;
2084 comp[0] = 0;
2085 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed side", "marker", 0, 1, comp, (void (*)(void)) zero, NULL, 1, &id, user);CHKERRQ(ierr);
2086 id = 4;
2087 comp[0] = 0;
2088 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed side", "marker", 0, 1, comp, (void (*)(void)) zero, NULL, 1, &id, user);CHKERRQ(ierr);
2089 id = 1;
2090 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", "marker", 2, 0, NULL, (void (*)(void)) terzaghi_drainage_pressure, NULL, 1, &id, user);CHKERRQ(ierr);
2091 break;
2092 case SOL_MANDEL:
2093 ierr = PetscDSSetResidual(prob, 0, NULL, f1_u);CHKERRQ(ierr);
2094 ierr = PetscDSSetResidual(prob, 1, f0_epsilon, NULL);CHKERRQ(ierr);
2095 ierr = PetscDSSetResidual(prob, 2, f0_p, f1_p);CHKERRQ(ierr);
2096 ierr = PetscDSSetBdResidual(prob, 0, f0_mandel_bd_u, NULL);CHKERRQ(ierr);
2097 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr);
2098 ierr = PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_ue, NULL);CHKERRQ(ierr);
2099 ierr = PetscDSSetJacobian(prob, 0, 2, NULL, NULL, g2_up, NULL);CHKERRQ(ierr);
2100 ierr = PetscDSSetJacobian(prob, 1, 0, NULL, g1_eu, NULL, NULL);CHKERRQ(ierr);
2101 ierr = PetscDSSetJacobian(prob, 1, 1, g0_ee, NULL, NULL, NULL);CHKERRQ(ierr);
2102 ierr = PetscDSSetJacobian(prob, 2, 1, g0_pe, NULL, NULL, NULL);CHKERRQ(ierr);
2103 ierr = PetscDSSetJacobian(prob, 2, 2, g0_pp, NULL, NULL, g3_pp);CHKERRQ(ierr);
2104
2105 ierr = mandelZeros(PETSC_COMM_WORLD, user, param);CHKERRQ(ierr);
2106
2107 exact[0] = mandel_2d_u;
2108 exact[1] = mandel_2d_eps;
2109 exact[2] = mandel_2d_p;
2110 exact_t[0] = mandel_2d_u_t;
2111 exact_t[1] = mandel_2d_eps_t;
2112 exact_t[2] = mandel_2d_p_t;
2113
2114 id_mandel[0] = 3;
2115 id_mandel[1] = 1;
2116 //comp[0] = 1;
2117 comp_mandel[0] = 0;
2118 comp_mandel[1] = 1;
2119 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "vertical stress", "marker", 0, 2, comp_mandel, (void (*)(void)) mandel_2d_u, NULL, 2, id_mandel, user);CHKERRQ(ierr);
2120 //ierr = DMAddBoundary(dm, DM_BC_NATURAL, "vertical stress", "marker", 0, 1, comp, NULL, 2, id_mandel, user);CHKERRQ(ierr);
2121 //ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed base", "marker", 0, 1, comp, (void (*)(void)) zero, 2, id_mandel, user);CHKERRQ(ierr);
2122
2123 id_mandel[0] = 2;
2124 id_mandel[1] = 4;
2125 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", "marker", 2, 0, NULL, (void (*)(void)) zero, NULL, 2, id_mandel, user);CHKERRQ(ierr);
2126 break;
2127 case SOL_CRYER:
2128 ierr = PetscDSSetResidual(prob, 0, NULL, f1_u);CHKERRQ(ierr);
2129 ierr = PetscDSSetResidual(prob, 1, f0_epsilon, NULL);CHKERRQ(ierr);
2130 ierr = PetscDSSetResidual(prob, 2, f0_p, f1_p);CHKERRQ(ierr);
2131 ierr = PetscDSSetBdResidual(prob, 0, f0_cryer_bd_u, NULL);CHKERRQ(ierr);
2132 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr);
2133 ierr = PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_ue, NULL);CHKERRQ(ierr);
2134 ierr = PetscDSSetJacobian(prob, 0, 2, NULL, NULL, g2_up, NULL);CHKERRQ(ierr);
2135 ierr = PetscDSSetJacobian(prob, 1, 0, NULL, g1_eu, NULL, NULL);CHKERRQ(ierr);
2136 ierr = PetscDSSetJacobian(prob, 1, 1, g0_ee, NULL, NULL, NULL);CHKERRQ(ierr);
2137 ierr = PetscDSSetJacobian(prob, 2, 1, g0_pe, NULL, NULL, NULL);CHKERRQ(ierr);
2138 ierr = PetscDSSetJacobian(prob, 2, 2, g0_pp, NULL, NULL, g3_pp);CHKERRQ(ierr);
2139
2140 ierr = cryerZeros(PETSC_COMM_WORLD, user, param);CHKERRQ(ierr);
2141
2142 exact[0] = cryer_3d_u;
2143 exact[1] = cryer_3d_eps;
2144 exact[2] = cryer_3d_p;
2145
2146 id = 1;
2147 ierr = DMAddBoundary(dm, DM_BC_NATURAL, "normal stress", "marker", 0, 0, NULL, NULL, NULL, 1, &id, user);CHKERRQ(ierr);
2148 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", "marker", 2, 0, NULL, (void (*)(void)) cryer_drainage_pressure, NULL, 1, &id, user);CHKERRQ(ierr);
2149 break;
2150 default: SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%D)", solutionTypes[PetscMin(user->solType, NUM_SOLUTION_TYPES)], user->solType);
2151 }
2152 for (f = 0; f < 3; ++f) {
2153 ierr = PetscDSSetExactSolution(prob, f, exact[f], user);CHKERRQ(ierr);
2154 ierr = PetscDSSetExactSolutionTimeDerivative(prob, f, exact_t[f], user);CHKERRQ(ierr);
2155 }
2156
2157 /* Setup constants */
2158 {
2159 PetscScalar constants[6];
2160 constants[0] = param->mu; /* shear modulus, Pa */
2161 constants[1] = param->K_u; /* undrained bulk modulus, Pa */
2162 constants[2] = param->alpha; /* Biot effective stress coefficient, - */
2163 constants[3] = param->M; /* Biot modulus, Pa */
2164 constants[4] = param->k/param->mu_f; /* Darcy coefficient, m**2 / Pa*s */
2165 constants[5] = param->P_0; /* Magnitude of Vertical Stress, Pa */
2166 ierr = PetscDSSetConstants(prob, 6, constants);CHKERRQ(ierr);
2167 }
2168 PetscFunctionReturn(0);
2169 }
2170
CreateElasticityNullSpace(DM dm,PetscInt origField,PetscInt field,MatNullSpace * nullspace)2171 static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace)
2172 {
2173 PetscErrorCode ierr;
2174
2175 PetscFunctionBegin;
2176 ierr = DMPlexCreateRigidBody(dm, origField, nullspace);CHKERRQ(ierr);
2177 PetscFunctionReturn(0);
2178 }
2179
SetupFE(DM dm,PetscBool simplex,PetscInt Nf,PetscInt Nc[],const char * name[],PetscErrorCode (* setup)(DM,AppCtx *),void * ctx)2180 static PetscErrorCode SetupFE(DM dm, PetscBool simplex, PetscInt Nf, PetscInt Nc[], const char *name[], PetscErrorCode (*setup)(DM, AppCtx *), void *ctx)
2181 {
2182 AppCtx *user = (AppCtx *) ctx;
2183 DM cdm = dm;
2184 PetscFE fe;
2185 PetscQuadrature q = NULL;
2186 char prefix[PETSC_MAX_PATH_LEN];
2187 PetscInt dim, f;
2188 PetscErrorCode ierr;
2189
2190 PetscFunctionBegin;
2191 /* Create finite element */
2192 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
2193 for (f = 0; f < Nf; ++f) {
2194 ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name[f]);CHKERRQ(ierr);
2195 ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc[f], simplex, name[f] ? prefix : NULL, -1, &fe);CHKERRQ(ierr);
2196 ierr = PetscObjectSetName((PetscObject) fe, name[f]);CHKERRQ(ierr);
2197 if (!q) {ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);}
2198 ierr = PetscFESetQuadrature(fe, q);CHKERRQ(ierr);
2199 ierr = DMSetField(dm, f, NULL, (PetscObject) fe);CHKERRQ(ierr);
2200 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
2201 }
2202 ierr = DMCreateDS(dm);CHKERRQ(ierr);
2203 ierr = (*setup)(dm, user);CHKERRQ(ierr);
2204 while (cdm) {
2205 ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr);
2206 if (0) {ierr = DMSetNearNullSpaceConstructor(cdm, 0, CreateElasticityNullSpace);CHKERRQ(ierr);}
2207 /* TODO: Check whether the boundary of coarse meshes is marked */
2208 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr);
2209 }
2210 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
2211 PetscFunctionReturn(0);
2212 }
2213
SetInitialConditions(TS ts,Vec u)2214 static PetscErrorCode SetInitialConditions(TS ts, Vec u)
2215 {
2216 DM dm;
2217 PetscReal t;
2218 PetscErrorCode ierr;
2219
2220 PetscFunctionBegin;
2221 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
2222 ierr = TSGetTime(ts, &t);CHKERRQ(ierr);
2223 if (t <= 0.0) {
2224 PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *);
2225 void *ctxs[3];
2226 AppCtx *ctx;
2227
2228 ierr = DMGetApplicationContext(dm, (void **) &ctx);CHKERRQ(ierr);
2229 switch (ctx->solType) {
2230 case SOL_TERZAGHI:
2231 funcs[0] = terzaghi_initial_u; ctxs[0] = ctx;
2232 funcs[1] = terzaghi_initial_eps; ctxs[1] = ctx;
2233 funcs[2] = terzaghi_drainage_pressure; ctxs[2] = ctx;
2234 ierr = DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u);CHKERRQ(ierr);
2235 break;
2236 case SOL_MANDEL:
2237 funcs[0] = mandel_initial_u; ctxs[0] = ctx;
2238 funcs[1] = mandel_initial_eps; ctxs[1] = ctx;
2239 funcs[2] = mandel_drainage_pressure; ctxs[2] = ctx;
2240 ierr = DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u);CHKERRQ(ierr);
2241 break;
2242 case SOL_CRYER:
2243 funcs[0] = cryer_initial_u; ctxs[0] = ctx;
2244 funcs[1] = cryer_initial_eps; ctxs[1] = ctx;
2245 funcs[2] = cryer_drainage_pressure; ctxs[2] = ctx;
2246 ierr = DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u);CHKERRQ(ierr);
2247 break;
2248 default:
2249 ierr = DMComputeExactSolution(dm, t, u, NULL);CHKERRQ(ierr);
2250 }
2251 } else {
2252 ierr = DMComputeExactSolution(dm, t, u, NULL);CHKERRQ(ierr);
2253 }
2254 PetscFunctionReturn(0);
2255 }
2256
2257 /* Need to create Viewer each time because HDF5 can get corrupted */
SolutionMonitor(TS ts,PetscInt steps,PetscReal time,Vec u,void * mctx)2258 static PetscErrorCode SolutionMonitor(TS ts, PetscInt steps, PetscReal time, Vec u, void *mctx)
2259 {
2260 DM dm;
2261 Vec exact;
2262 PetscViewer viewer;
2263 PetscViewerFormat format;
2264 PetscOptions options;
2265 const char *prefix;
2266 PetscErrorCode ierr;
2267
2268 PetscFunctionBegin;
2269 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
2270 ierr = PetscObjectGetOptions((PetscObject) ts, &options);CHKERRQ(ierr);
2271 ierr = PetscObjectGetOptionsPrefix((PetscObject) ts, &prefix);CHKERRQ(ierr);
2272 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject) ts), options, prefix, "-monitor_solution", &viewer, &format, NULL);CHKERRQ(ierr);
2273 ierr = DMGetGlobalVector(dm, &exact);CHKERRQ(ierr);
2274 ierr = DMComputeExactSolution(dm, time, exact, NULL);CHKERRQ(ierr);
2275 ierr = DMSetOutputSequenceNumber(dm, steps, time);CHKERRQ(ierr);
2276 ierr = VecView(exact, viewer);CHKERRQ(ierr);
2277 ierr = VecView(u, viewer);CHKERRQ(ierr);
2278 ierr = DMRestoreGlobalVector(dm, &exact);CHKERRQ(ierr);
2279 {
2280 PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
2281 void **ectxs;
2282 PetscReal *err;
2283 PetscInt Nf, f;
2284
2285 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr);
2286 ierr = PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err);CHKERRQ(ierr);
2287 {
2288 PetscInt Nds, s;
2289
2290 ierr = DMGetNumDS(dm, &Nds);CHKERRQ(ierr);
2291 for (s = 0; s < Nds; ++s) {
2292 PetscDS ds;
2293 DMLabel label;
2294 IS fieldIS;
2295 const PetscInt *fields;
2296 PetscInt dsNf, f;
2297
2298 ierr = DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds);CHKERRQ(ierr);
2299 ierr = PetscDSGetNumFields(ds, &dsNf);CHKERRQ(ierr);
2300 ierr = ISGetIndices(fieldIS, &fields);CHKERRQ(ierr);
2301 for (f = 0; f < dsNf; ++f) {
2302 const PetscInt field = fields[f];
2303 ierr = PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]);CHKERRQ(ierr);
2304 }
2305 ierr = ISRestoreIndices(fieldIS, &fields);CHKERRQ(ierr);
2306 }
2307 }
2308 ierr = DMComputeL2FieldDiff(dm, time, exacts, ectxs, u, err);CHKERRQ(ierr);
2309 ierr = PetscPrintf(PetscObjectComm((PetscObject) ts), "Time: %g L_2 Error: [", time);CHKERRQ(ierr);
2310 for (f = 0; f < Nf; ++f) {
2311 if (f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) ts), ", ");CHKERRQ(ierr);}
2312 ierr = PetscPrintf(PetscObjectComm((PetscObject) ts), "%g", (double) err[f]);CHKERRQ(ierr);
2313 }
2314 ierr = PetscPrintf(PetscObjectComm((PetscObject) ts), "]\n");CHKERRQ(ierr);
2315 ierr = PetscFree3(exacts, ectxs, err);CHKERRQ(ierr);
2316 }
2317 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
2318 PetscFunctionReturn(0);
2319 }
2320
SetupMonitor(TS ts,AppCtx * ctx)2321 static PetscErrorCode SetupMonitor(TS ts, AppCtx *ctx)
2322 {
2323 PetscViewer viewer;
2324 PetscViewerFormat format;
2325 PetscOptions options;
2326 const char *prefix;
2327 PetscBool flg;
2328 PetscErrorCode ierr;
2329
2330 PetscFunctionBegin;
2331 ierr = PetscObjectGetOptions((PetscObject) ts, &options);CHKERRQ(ierr);
2332 ierr = PetscObjectGetOptionsPrefix((PetscObject) ts, &prefix);CHKERRQ(ierr);
2333 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject) ts), options, prefix, "-monitor_solution", &viewer, &format, &flg);CHKERRQ(ierr);
2334 if (flg) {ierr = TSMonitorSet(ts, SolutionMonitor, ctx, NULL);CHKERRQ(ierr);}
2335 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
2336 PetscFunctionReturn(0);
2337 }
2338
TSAdaptChoose_Terzaghi(TSAdapt adapt,TS ts,PetscReal h,PetscInt * next_sc,PetscReal * next_h,PetscBool * accept,PetscReal * wlte,PetscReal * wltea,PetscReal * wlter)2339 static PetscErrorCode TSAdaptChoose_Terzaghi(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept, PetscReal *wlte, PetscReal *wltea, PetscReal *wlter)
2340 {
2341 static PetscReal dtTarget = -1.0;
2342 PetscReal dtInitial;
2343 DM dm;
2344 AppCtx *ctx;
2345 PetscInt step;
2346 PetscErrorCode ierr;
2347
2348 PetscFunctionBegin;
2349 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
2350 ierr = DMGetApplicationContext(dm, (void **) &ctx);CHKERRQ(ierr);
2351 ierr = TSGetStepNumber(ts, &step);CHKERRQ(ierr);
2352 dtInitial = ctx->dtInitial < 0.0 ? 1.0e-4*ctx->t_r : ctx->dtInitial;
2353 if (!step) {
2354 if (PetscAbsReal(dtInitial - h) > PETSC_SMALL) {
2355 *accept = PETSC_FALSE;
2356 *next_h = dtInitial;
2357 dtTarget = h;
2358 } else {
2359 *accept = PETSC_TRUE;
2360 *next_h = dtTarget < 0.0 ? dtInitial : dtTarget;
2361 dtTarget = -1.0;
2362 }
2363 } else {
2364 *accept = PETSC_TRUE;
2365 *next_h = h;
2366 }
2367 *next_sc = 0; /* Reuse the same order scheme */
2368 *wlte = -1; /* Weighted local truncation error was not evaluated */
2369 *wltea = -1; /* Weighted absolute local truncation error was not evaluated */
2370 *wlter = -1; /* Weighted relative local truncation error was not evaluated */
2371 PetscFunctionReturn(0);
2372 }
2373
main(int argc,char ** argv)2374 int main(int argc, char **argv)
2375 {
2376 AppCtx ctx; /* User-defined work context */
2377 DM dm; /* Problem specification */
2378 TS ts; /* Time Series / Nonlinear solver */
2379 Vec u; /* Solutions */
2380 const char *name[3] = {"displacement", "tracestrain", "pressure"};
2381 PetscReal t;
2382 PetscInt Nc[3];
2383 PetscErrorCode ierr;
2384
2385 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
2386 ierr = ProcessOptions(PETSC_COMM_WORLD, &ctx);CHKERRQ(ierr);
2387 ierr = PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &ctx.bag);CHKERRQ(ierr);
2388 ierr = PetscMalloc1(ctx.niter, &ctx.zeroArray);CHKERRQ(ierr);
2389 ierr = SetupParameters(PETSC_COMM_WORLD, &ctx);CHKERRQ(ierr);
2390 /* Primal System */
2391 ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr);
2392 ierr = CreateMesh(PETSC_COMM_WORLD, &ctx, &dm);CHKERRQ(ierr);
2393 ierr = DMSetApplicationContext(dm, &ctx);CHKERRQ(ierr);
2394 ierr = TSSetDM(ts, dm);CHKERRQ(ierr);
2395
2396 Nc[0] = ctx.dim;
2397 Nc[1] = 1;
2398 Nc[2] = 1;
2399
2400 ierr = SetupFE(dm, ctx.simplex, 3, Nc, name, SetupPrimalProblem, &ctx);CHKERRQ(ierr);
2401 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr);
2402 ierr = DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);CHKERRQ(ierr);
2403 ierr = DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);CHKERRQ(ierr);
2404 ierr = DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);CHKERRQ(ierr);
2405 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);
2406 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
2407 ierr = TSSetComputeInitialCondition(ts, SetInitialConditions);CHKERRQ(ierr);
2408 ierr = SetupMonitor(ts, &ctx);CHKERRQ(ierr);
2409
2410 if (ctx.solType != SOL_QUADRATIC_TRIG) {
2411 TSAdapt adapt;
2412
2413 ierr = TSGetAdapt(ts, &adapt);CHKERRQ(ierr);
2414 adapt->ops->choose = TSAdaptChoose_Terzaghi;
2415 }
2416 if (ctx.solType == SOL_CRYER) {
2417 Mat J;
2418 MatNullSpace sp;
2419
2420 ierr = TSSetUp(ts);CHKERRQ(ierr);
2421 ierr = TSGetIJacobian(ts, &J, NULL, NULL, NULL);CHKERRQ(ierr);
2422 ierr = DMPlexCreateRigidBody(dm, 0, &sp);CHKERRQ(ierr);
2423 ierr = MatSetNullSpace(J, sp);CHKERRQ(ierr);
2424 ierr = MatNullSpaceDestroy(&sp);CHKERRQ(ierr);
2425 }
2426 ierr = TSGetTime(ts, &t);CHKERRQ(ierr);
2427 ierr = DMSetOutputSequenceNumber(dm, 0, t);CHKERRQ(ierr);
2428 ierr = DMTSCheckFromOptions(ts, u);CHKERRQ(ierr);
2429 ierr = SetInitialConditions(ts, u);CHKERRQ(ierr);
2430 ierr = PetscObjectSetName((PetscObject) u, "solution");CHKERRQ(ierr);
2431 ierr = TSSolve(ts, u);CHKERRQ(ierr);
2432 ierr = DMTSCheckFromOptions(ts, u);CHKERRQ(ierr);
2433 ierr = TSGetSolution(ts, &u);CHKERRQ(ierr);
2434 ierr = VecViewFromOptions(u, NULL, "-sol_vec_view");CHKERRQ(ierr);
2435
2436 /* Cleanup */
2437 ierr = VecDestroy(&u);CHKERRQ(ierr);
2438 ierr = TSDestroy(&ts);CHKERRQ(ierr);
2439 ierr = DMDestroy(&dm);CHKERRQ(ierr);
2440 ierr = PetscBagDestroy(&ctx.bag);CHKERRQ(ierr);
2441 ierr = PetscFree(ctx.zeroArray);CHKERRQ(ierr);
2442 ierr = PetscFinalize();
2443 return ierr;
2444 }
2445
2446 /*TEST
2447
2448 test:
2449 suffix: 2d_quad_linear
2450 requires: triangle
2451 args: -sol_type quadratic_linear -dm_refine 2 \
2452 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2453 -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme
2454
2455 test:
2456 suffix: 3d_quad_linear
2457 requires: ctetgen
2458 args: -dim 3 -sol_type quadratic_linear -dm_refine 1 \
2459 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2460 -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme
2461
2462 test:
2463 suffix: 2d_trig_linear
2464 requires: triangle
2465 args: -sol_type trig_linear -dm_refine 1 \
2466 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2467 -dmts_check .0001 -ts_max_steps 5 -ts_dt 0.00001 -ts_monitor_extreme
2468
2469 test:
2470 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9, 2.1, 1.8]
2471 suffix: 2d_trig_linear_sconv
2472 requires: triangle
2473 args: -sol_type trig_linear -dm_refine 1 \
2474 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2475 -convest_num_refine 1 -ts_convergence_estimate -ts_convergence_temporal 0 -ts_max_steps 1 -ts_dt 0.00001 -pc_type lu
2476
2477 test:
2478 suffix: 3d_trig_linear
2479 requires: ctetgen
2480 args: -dim 3 -sol_type trig_linear -dm_refine 1 \
2481 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2482 -dmts_check .0001 -ts_max_steps 2 -ts_monitor_extreme
2483
2484 test:
2485 # -dm_refine 1 -convest_num_refine 2 gets L_2 convergence rate: [2.0, 2.1, 1.9]
2486 suffix: 3d_trig_linear_sconv
2487 requires: ctetgen
2488 args: -dim 3 -sol_type trig_linear -dm_refine 1 \
2489 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2490 -convest_num_refine 1 -ts_convergence_estimate -ts_convergence_temporal 0 -ts_max_steps 1 -pc_type lu
2491
2492 test:
2493 suffix: 2d_quad_trig
2494 requires: triangle
2495 args: -sol_type quadratic_trig -dm_refine 2 \
2496 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2497 -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme
2498
2499 test:
2500 # Using -dm_refine 4 gets the convergence rates to [0.95, 0.97, 0.90]
2501 suffix: 2d_quad_trig_tconv
2502 requires: triangle
2503 args: -sol_type quadratic_trig -dm_refine 1 \
2504 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2505 -convest_num_refine 3 -ts_convergence_estimate -ts_max_steps 5 -pc_type lu
2506
2507 test:
2508 suffix: 3d_quad_trig
2509 requires: ctetgen
2510 args: -dim 3 -sol_type quadratic_trig -dm_refine 1 \
2511 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2512 -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme
2513
2514 test:
2515 # Using -dm_refine 2 -convest_num_refine 3 gets the convergence rates to [1.0, 1.0, 1.0]
2516 suffix: 3d_quad_trig_tconv
2517 requires: ctetgen
2518 args: -dim 3 -sol_type quadratic_trig -dm_refine 1 \
2519 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2520 -convest_num_refine 1 -ts_convergence_estimate -ts_max_steps 5 -pc_type lu
2521
2522 test:
2523 suffix: 2d_terzaghi
2524 requires: triangle
2525 args: -sol_type terzaghi -dm_plex_separate_marker -dm_plex_box_faces 1,8 -simplex 0 \
2526 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 -niter 16000 \
2527 -ts_dt 0.0028666667 -ts_max_steps 2 -ts_monitor -dmts_check .0001 -pc_type lu
2528
2529 test:
2530 # -dm_plex_box_faces 1,64 -ts_max_steps 4 -convest_num_refine 3 gives L_2 convergence rate: [1.1, 1.1, 1.1]
2531 suffix: 2d_terzaghi_tconv
2532 requires: triangle
2533 args: -sol_type terzaghi -dm_plex_separate_marker -dm_plex_box_faces 1,8 -simplex 0 \
2534 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 -niter 16000 \
2535 -ts_dt 0.023 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1 -pc_type lu
2536
2537 test:
2538 # -dm_plex_box_faces 1,16 -convest_num_refine 4 gives L_2 convergence rate: [1.7, 1.2, 1.1]
2539 # if we add -displacement_petscspace_degree 3 -tracestrain_petscspace_degree 2 -pressure_petscspace_degree 2, we get [2.1, 1.6, 1.5], so I think we lose an order
2540 suffix: 2d_terzaghi_sconv
2541 requires: triangle
2542 args: -sol_type terzaghi -dm_plex_separate_marker -dm_plex_box_faces 1,8 -simplex 0 \
2543 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 -niter 16000 \
2544 -ts_dt 1e-5 -dt_initial 1e-5 -ts_max_steps 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 -pc_type lu
2545
2546 test:
2547 suffix: 2d_mandel
2548 requires: triangle
2549 args: -sol_type mandel -dm_plex_separate_marker -dm_refine 1 \
2550 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2551 -ts_dt 0.0028666667 -ts_max_steps 2 -ts_monitor -dmts_check .0001 -pc_type lu
2552
2553 test:
2554 # -dm_refine 5 -ts_max_steps 4 -convest_num_refine 3 gives L_2 convergence rate: [0.26, -0.0058, 0.26]
2555 suffix: 2d_mandel_tconv
2556 requires: triangle
2557 args: -sol_type mandel -dm_plex_separate_marker -dm_refine 1 \
2558 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2559 -ts_dt 0.023 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1 -pc_type lu
2560
2561 test:
2562 suffix: 3d_cryer
2563 requires: ctetgen !complex
2564 args: -sol_type cryer \
2565 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2566 -ts_dt 0.0028666667 -ts_max_time 0.014333 -ts_max_steps 2 -dmts_check .0001 -pc_type lu -pc_factor_shift_type nonzero
2567
2568 test:
2569 # Displacement and Pressure converge. The analytic expression for trace strain is inaccurate at the origin
2570 # -bd_dm_refine 3 -ref_limit 0.00666667 -ts_max_steps 5 -convest_num_refine 2 gives L_2 convergence rate: [0.47, -0.43, 1.5]
2571 suffix: 3d_cryer_tconv
2572 requires: ctetgen !complex
2573 args: -sol_type cryer -bd_dm_refine 1 -ref_limit 0.00666667 \
2574 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2575 -ts_dt 0.023 -ts_max_time 0.092 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1 -pc_type lu -pc_factor_shift_type nonzero
2576
2577 TEST*/
2578