Actual source code: ex53.c
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";
6: #include <petscdmplex.h>
7: #include <petscsnes.h>
8: #include <petscts.h>
9: #include <petscds.h>
10: #include <petscbag.h>
12: #include <petsc/private/tsimpl.h>
14: /* This presentation of poroelasticity is taken from
16: @book{Cheng2016,
17: title={Poroelasticity},
18: author={Cheng, Alexander H-D},
19: volume={27},
20: year={2016},
21: publisher={Springer}
22: }
24: For visualization, use
26: -dm_view hdf5:${PETSC_DIR}/sol.h5 -monitor_solution hdf5:${PETSC_DIR}/sol.h5::append
28: The weak form would then be, using test function $(v, q, \tau)$,
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: */
35: typedef enum {
36: SOL_QUADRATIC_LINEAR,
37: SOL_QUADRATIC_TRIG,
38: SOL_TRIG_LINEAR,
39: SOL_TERZAGHI,
40: SOL_MANDEL,
41: SOL_CRYER,
42: NUM_SOLUTION_TYPES
43: } SolutionType;
44: const char *solutionTypes[NUM_SOLUTION_TYPES + 1] = {"quadratic_linear", "quadratic_trig", "trig_linear", "terzaghi", "mandel", "cryer", "unknown"};
46: typedef struct {
47: PetscScalar mu; /* shear modulus */
48: PetscScalar K_u; /* undrained bulk modulus */
49: PetscScalar alpha; /* Biot effective stress coefficient */
50: PetscScalar M; /* Biot modulus */
51: PetscScalar k; /* (isotropic) permeability */
52: PetscScalar mu_f; /* fluid dynamic viscosity */
53: PetscScalar P_0; /* magnitude of vertical stress */
54: } Parameter;
56: typedef struct {
57: /* Domain and mesh definition */
58: PetscReal xmin[3]; /* Lower left bottom corner of bounding box */
59: PetscReal xmax[3]; /* Upper right top corner of bounding box */
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;
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 PETSC_SUCCESS;
76: }
78: /* Quadratic space and linear time solution
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
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: */
113: static PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
114: {
115: PetscInt d;
117: for (d = 0; d < dim; ++d) u[d] = PetscSqr(x[d]) - (d > 0 ? 2.0 * x[d - 1] * x[d] : 0.0);
118: return PETSC_SUCCESS;
119: }
121: static PetscErrorCode linear_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
122: {
123: u[0] = 2.0 * x[dim - 1];
124: return PETSC_SUCCESS;
125: }
127: static PetscErrorCode linear_linear_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
128: {
129: PetscReal sum = 0.0;
130: PetscInt d;
132: for (d = 0; d < dim; ++d) sum += x[d];
133: u[0] = sum * time;
134: return PETSC_SUCCESS;
135: }
137: static PetscErrorCode linear_linear_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
138: {
139: PetscReal sum = 0.0;
140: PetscInt d;
142: for (d = 0; d < dim; ++d) sum += x[d];
143: u[0] = sum;
144: return PETSC_SUCCESS;
145: }
147: static void 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[])
148: {
149: const PetscReal G = PetscRealPart(constants[0]);
150: const PetscReal K_u = PetscRealPart(constants[1]);
151: const PetscReal alpha = PetscRealPart(constants[2]);
152: const PetscReal M = PetscRealPart(constants[3]);
153: const PetscReal K_d = K_u - alpha * alpha * M;
154: const PetscReal lambda = K_d - (2.0 * G) / 3.0;
155: PetscInt d;
157: for (d = 0; d < dim - 1; ++d) f0[d] -= 2.0 * G - alpha * t;
158: f0[dim - 1] -= 2.0 * lambda + 4.0 * G - alpha * t;
159: }
161: static void 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[])
162: {
163: const PetscReal alpha = PetscRealPart(constants[2]);
164: const PetscReal M = PetscRealPart(constants[3]);
165: PetscReal sum = 0.0;
166: PetscInt d;
168: for (d = 0; d < dim; ++d) sum += x[d];
169: f0[0] += u_t ? alpha * u_t[uOff[1]] : 0.0;
170: f0[0] += u_t ? u_t[uOff[2]] / M : 0.0;
171: f0[0] -= sum / M;
172: }
174: /* Quadratic space and trigonometric time solution
176: 2D:
177: u = x^2
178: v = y^2 - 2xy
179: p = (x + y) cos(t)
180: e = 2y
181: f = <2 G, 4 G + 2 \lambda > - <alpha cos(t), alpha cos(t)>
182: g = 0
183: \epsilon = / 2x -y \
184: \ -y 2y - 2x /
185: Tr(\epsilon) = e = div u = 2y
186: div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
187: = 2 G < 2-1, 2 > + \lambda <0, 2> - alpha <cos(t), cos(t)>
188: = <2 G, 4 G + 2 \lambda> - <alpha cos(t), alpha cos(t)>
189: \frac{1}{M} \frac{dp}{dt} + \alpha \frac{d\varepsilon}{dt} - \nabla \cdot \kappa \nabla p
190: = \frac{1}{M} \frac{dp}{dt} + \kappa \Delta p
191: = -(x + y)/M sin(t)
193: 3D:
194: u = x^2
195: v = y^2 - 2xy
196: w = z^2 - 2yz
197: p = (x + y + z) cos(t)
198: e = 2z
199: f = <2 G, 4 G + 2 \lambda > - <alpha cos(t), alpha cos(t), alpha cos(t)>
200: g = 0
201: \varepsilon = / 2x -y 0 \
202: | -y 2y - 2x -z |
203: \ 0 -z 2z - 2y/
204: Tr(\varepsilon) = div u = 2z
205: div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
206: = 2 G < 2-1, 2-1, 2 > + \lambda <0, 0, 2> - alpha <cos(t), cos(t), cos(t)>
207: = <2 G, 2G, 4 G + 2 \lambda> - <alpha cos(t), alpha cos(t), alpha cos(t)>
208: */
209: static PetscErrorCode linear_trig_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
210: {
211: PetscReal sum = 0.0;
212: PetscInt d;
214: for (d = 0; d < dim; ++d) sum += x[d];
215: u[0] = sum * PetscCosReal(time);
216: return PETSC_SUCCESS;
217: }
219: static PetscErrorCode linear_trig_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
220: {
221: PetscReal sum = 0.0;
222: PetscInt d;
224: for (d = 0; d < dim; ++d) sum += x[d];
225: u[0] = -sum * PetscSinReal(time);
226: return PETSC_SUCCESS;
227: }
229: static void 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[])
230: {
231: const PetscReal G = PetscRealPart(constants[0]);
232: const PetscReal K_u = PetscRealPart(constants[1]);
233: const PetscReal alpha = PetscRealPart(constants[2]);
234: const PetscReal M = PetscRealPart(constants[3]);
235: const PetscReal K_d = K_u - alpha * alpha * M;
236: const PetscReal lambda = K_d - (2.0 * G) / 3.0;
237: PetscInt d;
239: for (d = 0; d < dim - 1; ++d) f0[d] -= 2.0 * G - alpha * PetscCosReal(t);
240: f0[dim - 1] -= 2.0 * lambda + 4.0 * G - alpha * PetscCosReal(t);
241: }
243: static void 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[])
244: {
245: const PetscReal alpha = PetscRealPart(constants[2]);
246: const PetscReal M = PetscRealPart(constants[3]);
247: PetscReal sum = 0.0;
248: PetscInt d;
250: for (d = 0; d < dim; ++d) sum += x[d];
252: f0[0] += u_t ? alpha * u_t[uOff[1]] : 0.0;
253: f0[0] += u_t ? u_t[uOff[2]] / M : 0.0;
254: f0[0] += PetscSinReal(t) * sum / M;
255: }
257: /* Trigonometric space and linear time solution
259: u = sin(2 pi x)
260: v = sin(2 pi y) - 2xy
261: \varepsilon = / 2 pi cos(2 pi x) -y \
262: \ -y 2 pi cos(2 pi y) - 2x /
263: Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
264: div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
265: = \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) >
266: = \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) >
268: 2D:
269: u = sin(2 pi x)
270: v = sin(2 pi y) - 2xy
271: p = (cos(2 pi x) + cos(2 pi y)) t
272: e = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
273: 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)>
274: g = 0
275: \varepsilon = / 2 pi cos(2 pi x) -y \
276: \ -y 2 pi cos(2 pi y) - 2x /
277: Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
278: div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
279: = 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>
280: = < -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)>
281: \frac{1}{M} \frac{dp}{dt} + \alpha \frac{d\varepsilon}{dt} - \nabla \cdot \kappa \nabla p
282: = \frac{1}{M} \frac{dp}{dt} + \kappa \Delta p
283: = (cos(2 pi x) + cos(2 pi y))/M - 4 pi^2 \kappa (cos(2 pi x) + cos(2 pi y)) t
285: 3D:
286: u = sin(2 pi x)
287: v = sin(2 pi y) - 2xy
288: v = sin(2 pi y) - 2yz
289: p = (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) t
290: e = 2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2y
291: 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)>
292: g = 0
293: \varepsilon = / 2 pi cos(2 pi x) -y 0 \
294: | -y 2 pi cos(2 pi y) - 2x -z |
295: \ 0 -z 2 pi cos(2 pi z) - 2y /
296: Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2 y
297: div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
298: = 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>
299: = < -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)>
300: \frac{1}{M} \frac{dp}{dt} + \alpha \frac{d\varepsilon}{dt} - \nabla \cdot \kappa \nabla p
301: = \frac{1}{M} \frac{dp}{dt} + \kappa \Delta p
302: = (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
303: */
304: static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
305: {
306: PetscInt d;
308: for (d = 0; d < dim; ++d) u[d] = PetscSinReal(2. * PETSC_PI * x[d]) - (d > 0 ? 2.0 * x[d - 1] * x[d] : 0.0);
309: return PETSC_SUCCESS;
310: }
312: static PetscErrorCode trig_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
313: {
314: PetscReal sum = 0.0;
315: PetscInt d;
317: for (d = 0; d < dim; ++d) sum += 2. * PETSC_PI * PetscCosReal(2. * PETSC_PI * x[d]) - (d < dim - 1 ? 2. * x[d] : 0.0);
318: u[0] = sum;
319: return PETSC_SUCCESS;
320: }
322: static PetscErrorCode trig_linear_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
323: {
324: PetscReal sum = 0.0;
325: PetscInt d;
327: for (d = 0; d < dim; ++d) sum += PetscCosReal(2. * PETSC_PI * x[d]);
328: u[0] = sum * time;
329: return PETSC_SUCCESS;
330: }
332: static PetscErrorCode trig_linear_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
333: {
334: PetscReal sum = 0.0;
335: PetscInt d;
337: for (d = 0; d < dim; ++d) sum += PetscCosReal(2. * PETSC_PI * x[d]);
338: u[0] = sum;
339: return PETSC_SUCCESS;
340: }
342: static void 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[])
343: {
344: const PetscReal G = PetscRealPart(constants[0]);
345: const PetscReal K_u = PetscRealPart(constants[1]);
346: const PetscReal alpha = PetscRealPart(constants[2]);
347: const PetscReal M = PetscRealPart(constants[3]);
348: const PetscReal K_d = K_u - alpha * alpha * M;
349: const PetscReal lambda = K_d - (2.0 * G) / 3.0;
350: PetscInt d;
352: for (d = 0; d < dim - 1; ++d) 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;
353: 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;
354: }
356: static void 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[])
357: {
358: const PetscReal alpha = PetscRealPart(constants[2]);
359: const PetscReal M = PetscRealPart(constants[3]);
360: const PetscReal kappa = PetscRealPart(constants[4]);
361: PetscReal sum = 0.0;
362: PetscInt d;
364: for (d = 0; d < dim; ++d) sum += PetscCosReal(2. * PETSC_PI * x[d]);
365: f0[0] += u_t ? alpha * u_t[uOff[1]] : 0.0;
366: f0[0] += u_t ? u_t[uOff[2]] / M : 0.0;
367: f0[0] -= sum / M - 4 * PetscSqr(PETSC_PI) * kappa * sum * t;
368: }
370: /* Terzaghi Solutions */
371: /* The analytical solutions given here are drawn from chapter 7, section 3, */
372: /* "One-Dimensional Consolidation Problem," from Poroelasticity, by Cheng. */
373: static PetscErrorCode terzaghi_drainage_pressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
374: {
375: AppCtx *user = (AppCtx *)ctx;
376: Parameter *param;
378: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
379: if (time <= 0.0) {
380: PetscScalar alpha = param->alpha; /* - */
381: PetscScalar K_u = param->K_u; /* Pa */
382: PetscScalar M = param->M; /* Pa */
383: PetscScalar G = param->mu; /* Pa */
384: PetscScalar P_0 = param->P_0; /* Pa */
385: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
386: PetscScalar eta = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G); /* -, Cheng (B.11) */
387: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
389: u[0] = ((P_0 * eta) / (G * S));
390: } else {
391: u[0] = 0.0;
392: }
393: return PETSC_SUCCESS;
394: }
396: static PetscErrorCode terzaghi_initial_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
397: {
398: AppCtx *user = (AppCtx *)ctx;
399: Parameter *param;
401: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
402: {
403: PetscScalar K_u = param->K_u; /* Pa */
404: PetscScalar G = param->mu; /* Pa */
405: PetscScalar P_0 = param->P_0; /* Pa */
406: PetscReal L = user->xmax[1] - user->xmin[1]; /* m */
407: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
408: PetscReal zstar = x[1] / L; /* - */
410: u[0] = 0.0;
411: u[1] = ((P_0 * L * (1.0 - 2.0 * nu_u)) / (2.0 * G * (1.0 - nu_u))) * (1.0 - zstar);
412: }
413: return PETSC_SUCCESS;
414: }
416: static PetscErrorCode terzaghi_initial_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
417: {
418: AppCtx *user = (AppCtx *)ctx;
419: Parameter *param;
421: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
422: {
423: PetscScalar K_u = param->K_u; /* Pa */
424: PetscScalar G = param->mu; /* Pa */
425: PetscScalar P_0 = param->P_0; /* Pa */
426: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
428: u[0] = -(P_0 * (1.0 - 2.0 * nu_u)) / (2.0 * G * (1.0 - nu_u));
429: }
430: return PETSC_SUCCESS;
431: }
433: static PetscErrorCode terzaghi_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
434: {
435: AppCtx *user = (AppCtx *)ctx;
436: Parameter *param;
438: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
439: if (time < 0.0) {
440: PetscCall(terzaghi_initial_u(dim, time, x, Nc, u, ctx));
441: } else {
442: PetscScalar alpha = param->alpha; /* - */
443: PetscScalar K_u = param->K_u; /* Pa */
444: PetscScalar M = param->M; /* Pa */
445: PetscScalar G = param->mu; /* Pa */
446: PetscScalar P_0 = param->P_0; /* Pa */
447: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
448: PetscReal L = user->xmax[1] - user->xmin[1]; /* m */
449: PetscInt N = user->niter, m;
451: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
452: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
453: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
454: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
455: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
457: PetscReal zstar = x[1] / L; /* - */
458: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
459: PetscScalar F2 = 0.0;
461: for (m = 1; m < 2 * N + 1; ++m) {
462: if (m % 2 == 1) F2 += (8.0 / PetscSqr(m * PETSC_PI)) * PetscCosReal(0.5 * m * PETSC_PI * zstar) * (1.0 - PetscExpReal(-PetscSqr(m * PETSC_PI) * tstar));
463: }
464: u[0] = 0.0;
465: 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 */
466: }
467: return PETSC_SUCCESS;
468: }
470: static PetscErrorCode terzaghi_2d_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
471: {
472: AppCtx *user = (AppCtx *)ctx;
473: Parameter *param;
475: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
476: if (time < 0.0) {
477: PetscCall(terzaghi_initial_eps(dim, time, x, Nc, u, ctx));
478: } else {
479: PetscScalar alpha = param->alpha; /* - */
480: PetscScalar K_u = param->K_u; /* Pa */
481: PetscScalar M = param->M; /* Pa */
482: PetscScalar G = param->mu; /* Pa */
483: PetscScalar P_0 = param->P_0; /* Pa */
484: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
485: PetscReal L = user->xmax[1] - user->xmin[1]; /* m */
486: PetscInt N = user->niter, m;
488: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
489: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
490: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
491: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
492: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
494: PetscReal zstar = x[1] / L; /* - */
495: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
496: PetscScalar F2_z = 0.0;
498: for (m = 1; m < 2 * N + 1; ++m) {
499: if (m % 2 == 1) F2_z += (-4.0 / (m * PETSC_PI * L)) * PetscSinReal(0.5 * m * PETSC_PI * zstar) * (1.0 - PetscExpReal(-PetscSqr(m * PETSC_PI) * tstar));
500: }
501: 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; /* - */
502: }
503: return PETSC_SUCCESS;
504: }
506: // Pressure
507: static PetscErrorCode terzaghi_2d_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
508: {
509: AppCtx *user = (AppCtx *)ctx;
510: Parameter *param;
512: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
513: if (time <= 0.0) {
514: PetscCall(terzaghi_drainage_pressure(dim, time, x, Nc, u, ctx));
515: } else {
516: PetscScalar alpha = param->alpha; /* - */
517: PetscScalar K_u = param->K_u; /* Pa */
518: PetscScalar M = param->M; /* Pa */
519: PetscScalar G = param->mu; /* Pa */
520: PetscScalar P_0 = param->P_0; /* Pa */
521: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
522: PetscReal L = user->xmax[1] - user->xmin[1]; /* m */
523: PetscInt N = user->niter, m;
525: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
526: PetscScalar eta = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G); /* -, Cheng (B.11) */
527: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
528: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
530: PetscReal zstar = x[1] / L; /* - */
531: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
532: PetscScalar F1 = 0.0;
534: PetscCheck(PetscAbsScalar((1 / M + (alpha * eta) / G) - S) <= 1.0e-10, PETSC_COMM_SELF, PETSC_ERR_PLIB, "S %g != check %g", (double)PetscAbsScalar(S), (double)PetscAbsScalar(1 / M + (alpha * eta) / G));
536: for (m = 1; m < 2 * N + 1; ++m) {
537: if (m % 2 == 1) F1 += (4.0 / (m * PETSC_PI)) * PetscSinReal(0.5 * m * PETSC_PI * zstar) * PetscExpReal(-PetscSqr(m * PETSC_PI) * tstar);
538: }
539: u[0] = ((P_0 * eta) / (G * S)) * F1; /* Pa */
540: }
541: return PETSC_SUCCESS;
542: }
544: static PetscErrorCode terzaghi_2d_u_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
545: {
546: AppCtx *user = (AppCtx *)ctx;
547: Parameter *param;
549: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
550: if (time <= 0.0) {
551: u[0] = 0.0;
552: u[1] = 0.0;
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 = user->xmax[1] - user->xmin[1]; /* m */
561: PetscInt N = user->niter, m;
563: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
564: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
565: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
566: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
567: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
569: PetscReal zstar = x[1] / L; /* - */
570: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
571: PetscScalar F2_t = 0.0;
573: for (m = 1; m < 2 * N + 1; ++m) {
574: if (m % 2 == 1) F2_t += (2.0 * c / PetscSqr(L)) * PetscCosReal(0.5 * m * PETSC_PI * zstar) * PetscExpReal(-PetscSqr(m * PETSC_PI) * tstar);
575: }
576: u[0] = 0.0;
577: u[1] = ((P_0 * L * (nu_u - nu)) / (2.0 * G * (1.0 - nu_u) * (1.0 - nu))) * F2_t; /* m / s */
578: }
579: return PETSC_SUCCESS;
580: }
582: static PetscErrorCode terzaghi_2d_eps_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
583: {
584: AppCtx *user = (AppCtx *)ctx;
585: Parameter *param;
587: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
588: if (time <= 0.0) {
589: u[0] = 0.0;
590: } else {
591: PetscScalar alpha = param->alpha; /* - */
592: PetscScalar K_u = param->K_u; /* Pa */
593: PetscScalar M = param->M; /* Pa */
594: PetscScalar G = param->mu; /* Pa */
595: PetscScalar P_0 = param->P_0; /* Pa */
596: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
597: PetscReal L = user->xmax[1] - user->xmin[1]; /* m */
598: PetscInt N = user->niter, m;
600: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
601: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
602: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
603: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
604: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
606: PetscReal zstar = x[1] / L; /* - */
607: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
608: PetscScalar F2_zt = 0.0;
610: for (m = 1; m < 2 * N + 1; ++m) {
611: if (m % 2 == 1) F2_zt += ((-m * PETSC_PI * c) / (L * L * L)) * PetscSinReal(0.5 * m * PETSC_PI * zstar) * PetscExpReal(-PetscSqr(m * PETSC_PI) * tstar);
612: }
613: u[0] = ((P_0 * L * (nu_u - nu)) / (2.0 * G * (1.0 - nu_u) * (1.0 - nu))) * F2_zt; /* 1 / s */
614: }
615: return PETSC_SUCCESS;
616: }
618: static PetscErrorCode terzaghi_2d_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
619: {
620: AppCtx *user = (AppCtx *)ctx;
621: Parameter *param;
623: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
624: if (time <= 0.0) {
625: PetscScalar alpha = param->alpha; /* - */
626: PetscScalar K_u = param->K_u; /* Pa */
627: PetscScalar M = param->M; /* Pa */
628: PetscScalar G = param->mu; /* Pa */
629: PetscScalar P_0 = param->P_0; /* Pa */
630: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
631: PetscReal L = user->xmax[1] - user->xmin[1]; /* m */
633: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
634: PetscScalar eta = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G); /* -, Cheng (B.11) */
635: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
636: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
638: u[0] = -((P_0 * eta) / (G * S)) * PetscSqr(0 * PETSC_PI) * c / PetscSqr(2.0 * L); /* Pa / s */
639: } else {
640: PetscScalar alpha = param->alpha; /* - */
641: PetscScalar K_u = param->K_u; /* Pa */
642: PetscScalar M = param->M; /* Pa */
643: PetscScalar G = param->mu; /* Pa */
644: PetscScalar P_0 = param->P_0; /* Pa */
645: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
646: PetscReal L = user->xmax[1] - user->xmin[1]; /* m */
647: PetscInt N = user->niter, m;
649: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
650: PetscScalar eta = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G); /* -, Cheng (B.11) */
651: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
652: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
654: PetscReal zstar = x[1] / L; /* - */
655: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
656: PetscScalar F1_t = 0.0;
658: PetscCheck(PetscAbsScalar((1 / M + (alpha * eta) / G) - S) <= 1.0e-10, PETSC_COMM_SELF, PETSC_ERR_PLIB, "S %g != check %g", (double)PetscAbsScalar(S), (double)PetscAbsScalar(1 / M + (alpha * eta) / G));
660: for (m = 1; m < 2 * N + 1; ++m) {
661: if (m % 2 == 1) F1_t += ((-m * PETSC_PI * c) / PetscSqr(L)) * PetscSinReal(0.5 * m * PETSC_PI * zstar) * PetscExpReal(-PetscSqr(m * PETSC_PI) * tstar);
662: }
663: u[0] = ((P_0 * eta) / (G * S)) * F1_t; /* Pa / s */
664: }
665: return PETSC_SUCCESS;
666: }
668: /* Mandel Solutions */
669: static PetscErrorCode mandel_drainage_pressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
670: {
671: AppCtx *user = (AppCtx *)ctx;
672: Parameter *param;
674: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
675: if (time <= 0.0) {
676: PetscScalar alpha = param->alpha; /* - */
677: PetscScalar K_u = param->K_u; /* Pa */
678: PetscScalar M = param->M; /* Pa */
679: PetscScalar G = param->mu; /* Pa */
680: PetscScalar P_0 = param->P_0; /* Pa */
681: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
682: PetscReal a = 0.5 * (user->xmax[0] - user->xmin[0]); /* m */
683: PetscInt N = user->niter, n;
685: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
686: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
687: PetscScalar B = alpha * M / K_u; /* -, Cheng (B.12) */
688: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
689: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
691: PetscScalar A1 = 3.0 / (B * (1.0 + nu_u));
692: PetscReal aa = 0.0;
693: PetscReal p = 0.0;
694: PetscReal time = 0.0;
696: for (n = 1; n < N + 1; ++n) {
697: aa = user->zeroArray[n - 1];
698: 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));
699: }
700: u[0] = ((2.0 * P_0) / (a * A1)) * p;
701: } else {
702: u[0] = 0.0;
703: }
704: return PETSC_SUCCESS;
705: }
707: static PetscErrorCode mandel_initial_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
708: {
709: AppCtx *user = (AppCtx *)ctx;
710: Parameter *param;
712: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
713: {
714: PetscScalar alpha = param->alpha; /* - */
715: PetscScalar K_u = param->K_u; /* Pa */
716: PetscScalar M = param->M; /* Pa */
717: PetscScalar G = param->mu; /* Pa */
718: PetscScalar P_0 = param->P_0; /* Pa */
719: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
720: PetscReal a = 0.5 * (user->xmax[0] - user->xmin[0]); /* m */
721: PetscInt N = user->niter, n;
723: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
724: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
725: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
726: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
727: PetscReal c = PetscRealPart(kappa / S); /* m^2 / s, Cheng (B.16) */
729: PetscReal A_s = 0.0;
730: PetscReal B_s = 0.0;
731: for (n = 1; n < N + 1; ++n) {
732: PetscReal alpha_n = user->zeroArray[n - 1];
733: 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));
734: 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));
735: }
736: u[0] = ((P_0 * nu) / (2.0 * G * a) - (P_0 * nu_u) / (G * a) * A_s) * x[0] + P_0 / G * B_s;
737: u[1] = (-1 * (P_0 * (1.0 - nu)) / (2 * G * a) + (P_0 * (1 - nu_u)) / (G * a) * A_s) * x[1];
738: }
739: return PETSC_SUCCESS;
740: }
742: static PetscErrorCode mandel_initial_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
743: {
744: AppCtx *user = (AppCtx *)ctx;
745: Parameter *param;
747: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
748: {
749: PetscScalar alpha = param->alpha; /* - */
750: PetscScalar K_u = param->K_u; /* Pa */
751: PetscScalar M = param->M; /* Pa */
752: PetscScalar G = param->mu; /* Pa */
753: PetscScalar P_0 = param->P_0; /* Pa */
754: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
755: PetscReal a = 0.5 * (user->xmax[0] - user->xmin[0]); /* m */
756: PetscInt N = user->niter, n;
758: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
759: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
760: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
761: PetscReal c = PetscRealPart(kappa / S); /* m^2 / s, Cheng (B.16) */
763: PetscReal aa = 0.0;
764: PetscReal eps_A = 0.0;
765: PetscReal eps_B = 0.0;
766: PetscReal eps_C = 0.0;
767: PetscReal time = 0.0;
769: for (n = 1; n < N + 1; ++n) {
770: aa = user->zeroArray[n - 1];
771: 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)));
772: eps_B += (PetscExpReal((-1.0 * aa * aa * c * time) / (a * a)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa));
773: eps_C += (PetscExpReal((-1.0 * aa * aa * c * time) / (aa * aa)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa));
774: }
775: 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);
776: }
777: return PETSC_SUCCESS;
778: }
780: // Displacement
781: static PetscErrorCode mandel_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
782: {
783: Parameter *param;
785: AppCtx *user = (AppCtx *)ctx;
787: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
788: if (time <= 0.0) {
789: PetscCall(mandel_initial_u(dim, time, x, Nc, u, ctx));
790: } else {
791: PetscInt NITER = user->niter;
792: PetscScalar alpha = param->alpha;
793: PetscScalar K_u = param->K_u;
794: PetscScalar M = param->M;
795: PetscScalar G = param->mu;
796: PetscScalar k = param->k;
797: PetscScalar mu_f = param->mu_f;
798: PetscScalar F = param->P_0;
800: PetscScalar K_d = K_u - alpha * alpha * M;
801: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
802: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
803: PetscScalar kappa = k / mu_f;
804: PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0;
805: PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u)));
807: // Series term
808: PetscScalar A_x = 0.0;
809: PetscScalar B_x = 0.0;
811: for (PetscInt n = 1; n < NITER + 1; n++) {
812: PetscReal alpha_n = user->zeroArray[n - 1];
814: 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));
815: 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));
816: }
817: u[0] = ((F * nu) / (2.0 * G * a) - (F * nu_u) / (G * a) * A_x) * x[0] + F / G * B_x;
818: u[1] = (-1 * (F * (1.0 - nu)) / (2 * G * a) + (F * (1 - nu_u)) / (G * a) * A_x) * x[1];
819: }
820: return PETSC_SUCCESS;
821: }
823: // Trace strain
824: static PetscErrorCode mandel_2d_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
825: {
826: Parameter *param;
828: AppCtx *user = (AppCtx *)ctx;
830: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
831: if (time <= 0.0) {
832: PetscCall(mandel_initial_eps(dim, time, x, Nc, u, ctx));
833: } else {
834: PetscInt NITER = user->niter;
835: PetscScalar alpha = param->alpha;
836: PetscScalar K_u = param->K_u;
837: PetscScalar M = param->M;
838: PetscScalar G = param->mu;
839: PetscScalar k = param->k;
840: PetscScalar mu_f = param->mu_f;
841: PetscScalar F = param->P_0;
843: PetscScalar K_d = K_u - alpha * alpha * M;
844: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
845: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
846: PetscScalar kappa = k / mu_f;
847: //const PetscScalar B = (alpha*M)/(K_d + alpha*alpha * M);
849: //const PetscScalar b = (YMAX - YMIN) / 2.0;
850: PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0;
851: PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u)));
853: // Series term
854: PetscReal eps_A = 0.0;
855: PetscReal eps_B = 0.0;
856: PetscReal eps_C = 0.0;
858: for (PetscInt n = 1; n < NITER + 1; n++) {
859: PetscReal aa = user->zeroArray[n - 1];
861: 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)));
863: eps_B += (PetscExpReal((-1.0 * aa * aa * c * time) / (a * a)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa));
865: eps_C += (PetscExpReal((-1.0 * aa * aa * c * time) / (aa * aa)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa));
866: }
868: 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);
869: }
870: return PETSC_SUCCESS;
871: }
873: // Pressure
874: static PetscErrorCode mandel_2d_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
875: {
876: Parameter *param;
878: AppCtx *user = (AppCtx *)ctx;
880: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
881: if (time <= 0.0) {
882: PetscCall(mandel_drainage_pressure(dim, time, x, Nc, u, ctx));
883: } else {
884: PetscInt NITER = user->niter;
886: PetscScalar alpha = param->alpha;
887: PetscScalar K_u = param->K_u;
888: PetscScalar M = param->M;
889: PetscScalar G = param->mu;
890: PetscScalar k = param->k;
891: PetscScalar mu_f = param->mu_f;
892: PetscScalar F = param->P_0;
894: PetscScalar K_d = K_u - alpha * alpha * M;
895: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
896: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
897: PetscScalar kappa = k / mu_f;
898: PetscScalar B = (alpha * M) / (K_d + alpha * alpha * M);
900: PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0;
901: PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u)));
902: PetscScalar A1 = 3.0 / (B * (1.0 + nu_u));
903: //PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu);
905: // Series term
906: PetscReal p = 0.0;
908: for (PetscInt n = 1; n < NITER + 1; n++) {
909: PetscReal aa = user->zeroArray[n - 1];
910: p += (PetscSinReal(aa) / (aa - PetscSinReal(aa) * PetscCosReal(aa))) * (PetscCosReal((aa * x[0]) / a) - PetscCosReal(aa)) * PetscExpReal(-1.0 * (aa * aa * c * time) / (a * a));
911: }
912: u[0] = ((2.0 * F) / (a * A1)) * p;
913: }
914: return PETSC_SUCCESS;
915: }
917: // Time derivative of displacement
918: static PetscErrorCode mandel_2d_u_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
919: {
920: Parameter *param;
922: AppCtx *user = (AppCtx *)ctx;
924: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
926: PetscInt NITER = user->niter;
927: PetscScalar alpha = param->alpha;
928: PetscScalar K_u = param->K_u;
929: PetscScalar M = param->M;
930: PetscScalar G = param->mu;
931: PetscScalar F = param->P_0;
933: PetscScalar K_d = K_u - alpha * alpha * M;
934: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
935: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
936: PetscScalar kappa = param->k / param->mu_f;
937: PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0;
938: PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u)));
940: // Series term
941: PetscScalar A_s_t = 0.0;
942: PetscScalar B_s_t = 0.0;
944: for (PetscInt n = 1; n < NITER + 1; n++) {
945: PetscReal alpha_n = user->zeroArray[n - 1];
947: 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)));
948: 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)));
949: }
951: u[0] = (F / G) * A_s_t - ((F * nu_u * x[0]) / (G * a)) * B_s_t;
952: u[1] = ((F * x[1] * (1 - nu_u)) / (G * a)) * B_s_t;
954: return PETSC_SUCCESS;
955: }
957: // Time derivative of trace strain
958: static PetscErrorCode mandel_2d_eps_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
959: {
960: Parameter *param;
962: AppCtx *user = (AppCtx *)ctx;
964: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
966: PetscInt NITER = user->niter;
967: PetscScalar alpha = param->alpha;
968: PetscScalar K_u = param->K_u;
969: PetscScalar M = param->M;
970: PetscScalar G = param->mu;
971: PetscScalar k = param->k;
972: PetscScalar mu_f = param->mu_f;
973: PetscScalar F = param->P_0;
975: PetscScalar K_d = K_u - alpha * alpha * M;
976: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
977: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
978: PetscScalar kappa = k / mu_f;
979: //const PetscScalar B = (alpha*M)/(K_d + alpha*alpha * M);
981: //const PetscScalar b = (YMAX - YMIN) / 2.0;
982: PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0;
983: PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u)));
985: // Series term
986: PetscScalar eps_As = 0.0;
987: PetscScalar eps_Bs = 0.0;
988: PetscScalar eps_Cs = 0.0;
990: for (PetscInt n = 1; n < NITER + 1; n++) {
991: PetscReal alpha_n = user->zeroArray[n - 1];
993: 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)));
994: 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)));
995: 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)));
996: }
998: u[0] = (F / G) * eps_As - ((F * nu_u) / (G * a)) * eps_Bs + ((F * (1 - nu_u)) / (G * a)) * eps_Cs;
999: return PETSC_SUCCESS;
1000: }
1002: // Time derivative of pressure
1003: static PetscErrorCode mandel_2d_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1004: {
1005: Parameter *param;
1007: AppCtx *user = (AppCtx *)ctx;
1009: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1011: PetscScalar alpha = param->alpha;
1012: PetscScalar K_u = param->K_u;
1013: PetscScalar M = param->M;
1014: PetscScalar G = param->mu;
1015: PetscScalar F = param->P_0;
1017: PetscScalar K_d = K_u - alpha * alpha * M;
1018: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
1019: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
1021: PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0;
1022: //PetscScalar A1 = 3.0 / (B * (1.0 + nu_u));
1023: //PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu);
1025: u[0] = ((2.0 * F * (-2.0 * nu + 3.0 * nu_u)) / (3.0 * a * alpha * (1.0 - 2.0 * nu)));
1027: return PETSC_SUCCESS;
1028: }
1030: /* Cryer Solutions */
1031: static PetscErrorCode cryer_drainage_pressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1032: {
1033: AppCtx *user = (AppCtx *)ctx;
1034: Parameter *param;
1036: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1037: if (time <= 0.0) {
1038: PetscScalar alpha = param->alpha; /* - */
1039: PetscScalar K_u = param->K_u; /* Pa */
1040: PetscScalar M = param->M; /* Pa */
1041: PetscScalar P_0 = param->P_0; /* Pa */
1042: PetscScalar B = alpha * M / K_u; /* -, Cheng (B.12) */
1044: u[0] = P_0 * B;
1045: } else {
1046: u[0] = 0.0;
1047: }
1048: return PETSC_SUCCESS;
1049: }
1051: static PetscErrorCode cryer_initial_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1052: {
1053: AppCtx *user = (AppCtx *)ctx;
1054: Parameter *param;
1056: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1057: {
1058: PetscScalar K_u = param->K_u; /* Pa */
1059: PetscScalar G = param->mu; /* Pa */
1060: PetscScalar P_0 = param->P_0; /* Pa */
1061: PetscReal R_0 = user->xmax[1]; /* m */
1062: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
1064: PetscScalar u_0 = -P_0 * R_0 * (1. - 2. * nu_u) / (2. * G * (1. + nu_u)); /* Cheng (7.407) */
1065: PetscReal u_sc = PetscRealPart(u_0) / R_0;
1067: u[0] = u_sc * x[0];
1068: u[1] = u_sc * x[1];
1069: u[2] = u_sc * x[2];
1070: }
1071: return PETSC_SUCCESS;
1072: }
1074: static PetscErrorCode cryer_initial_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1075: {
1076: AppCtx *user = (AppCtx *)ctx;
1077: Parameter *param;
1079: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1080: {
1081: PetscScalar K_u = param->K_u; /* Pa */
1082: PetscScalar G = param->mu; /* Pa */
1083: PetscScalar P_0 = param->P_0; /* Pa */
1084: PetscReal R_0 = user->xmax[1]; /* m */
1085: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
1087: PetscScalar u_0 = -P_0 * R_0 * (1. - 2. * nu_u) / (2. * G * (1. + nu_u)); /* Cheng (7.407) */
1088: PetscReal u_sc = PetscRealPart(u_0) / R_0;
1090: /* div R = 1/R^2 d/dR R^2 R = 3 */
1091: u[0] = 3. * u_sc;
1092: u[1] = 3. * u_sc;
1093: u[2] = 3. * u_sc;
1094: }
1095: return PETSC_SUCCESS;
1096: }
1098: // Displacement
1099: static PetscErrorCode cryer_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1100: {
1101: AppCtx *user = (AppCtx *)ctx;
1102: Parameter *param;
1104: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1105: if (time <= 0.0) {
1106: PetscCall(cryer_initial_u(dim, time, x, Nc, u, ctx));
1107: } else {
1108: PetscScalar alpha = param->alpha; /* - */
1109: PetscScalar K_u = param->K_u; /* Pa */
1110: PetscScalar M = param->M; /* Pa */
1111: PetscScalar G = param->mu; /* Pa */
1112: PetscScalar P_0 = param->P_0; /* Pa */
1113: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
1114: PetscReal R_0 = user->xmax[1]; /* m */
1115: PetscInt N = user->niter, n;
1117: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
1118: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
1119: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
1120: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
1121: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
1122: PetscScalar u_inf = -P_0 * R_0 * (1. - 2. * nu) / (2. * G * (1. + nu)); /* m, Cheng (7.388) */
1124: PetscReal R = PetscSqrtReal(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
1125: PetscReal R_star = R / R_0;
1126: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(R_0); /* - */
1127: PetscReal A_n = 0.0;
1128: PetscScalar u_sc;
1130: for (n = 1; n < N + 1; ++n) {
1131: const PetscReal x_n = user->zeroArray[n - 1];
1132: const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u));
1134: /* m , Cheng (7.404) */
1135: if (R_star != 0) {
1136: A_n += PetscRealPart((12.0 * (1.0 + nu) * (nu_u - nu)) / ((1.0 - 2.0 * nu) * E_n * PetscSqr(R_star) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * (3.0 * (nu_u - nu) * (PetscSinReal(R_star * PetscSqrtReal(x_n)) - R_star * PetscSqrtReal(x_n) * PetscCosReal(R_star * PetscSqrtReal(x_n))) + (1.0 - nu) * (1.0 - 2.0 * nu) * PetscPowRealInt(R_star, 3) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * PetscExpReal(-x_n * tstar));
1137: }
1138: }
1139: if (R_star != 0) u_sc = PetscRealPart(u_inf) * (R_star - A_n) / R;
1140: else u_sc = PetscRealPart(u_inf) / R_0;
1141: u[0] = u_sc * x[0];
1142: u[1] = u_sc * x[1];
1143: u[2] = u_sc * x[2];
1144: }
1145: return PETSC_SUCCESS;
1146: }
1148: // Volumetric Strain
1149: static PetscErrorCode cryer_3d_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1150: {
1151: AppCtx *user = (AppCtx *)ctx;
1152: Parameter *param;
1154: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1155: if (time <= 0.0) {
1156: PetscCall(cryer_initial_eps(dim, time, x, Nc, u, ctx));
1157: } else {
1158: PetscScalar alpha = param->alpha; /* - */
1159: PetscScalar K_u = param->K_u; /* Pa */
1160: PetscScalar M = param->M; /* Pa */
1161: PetscScalar G = param->mu; /* Pa */
1162: PetscScalar P_0 = param->P_0; /* Pa */
1163: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
1164: PetscReal R_0 = user->xmax[1]; /* m */
1165: PetscInt N = user->niter, n;
1167: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
1168: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
1169: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
1170: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
1171: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
1172: PetscScalar u_inf = -P_0 * R_0 * (1. - 2. * nu) / (2. * G * (1. + nu)); /* m, Cheng (7.388) */
1174: PetscReal R = PetscSqrtReal(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
1175: PetscReal R_star = R / R_0;
1176: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(R_0); /* - */
1177: PetscReal divA_n = 0.0;
1179: if (R_star < PETSC_SMALL) {
1180: for (n = 1; n < N + 1; ++n) {
1181: const PetscReal x_n = user->zeroArray[n - 1];
1182: const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u));
1184: divA_n += PetscRealPart((12.0 * (1.0 + nu) * (nu_u - nu)) / ((1.0 - 2.0 * nu) * E_n * PetscSqr(R_star) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * (3.0 * (nu_u - nu) * PetscSqrtReal(x_n) * ((2.0 + PetscSqr(R_star * PetscSqrtReal(x_n))) - 2.0 * PetscCosReal(R_star * PetscSqrtReal(x_n))) + 5.0 * (1.0 - nu) * (1.0 - 2.0 * nu) * PetscPowRealInt(R_star, 2) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * PetscExpReal(-x_n * tstar));
1185: }
1186: } else {
1187: for (n = 1; n < N + 1; ++n) {
1188: const PetscReal x_n = user->zeroArray[n - 1];
1189: const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u));
1191: divA_n += PetscRealPart((12.0 * (1.0 + nu) * (nu_u - nu)) / ((1.0 - 2.0 * nu) * E_n * PetscSqr(R_star) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * (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))) + 5.0 * (1.0 - nu) * (1.0 - 2.0 * nu) * PetscPowRealInt(R_star, 2) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * PetscExpReal(-x_n * tstar));
1192: }
1193: }
1194: if (PetscAbsReal(divA_n) > 1e3) PetscCall(PetscPrintf(PETSC_COMM_SELF, "(%g, %g, %g) divA_n: %g\n", (double)x[0], (double)x[1], (double)x[2], (double)divA_n));
1195: u[0] = PetscRealPart(u_inf) / R_0 * (3.0 - divA_n);
1196: }
1197: return PETSC_SUCCESS;
1198: }
1200: // Pressure
1201: static PetscErrorCode cryer_3d_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1202: {
1203: AppCtx *user = (AppCtx *)ctx;
1204: Parameter *param;
1206: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1207: if (time <= 0.0) {
1208: PetscCall(cryer_drainage_pressure(dim, time, x, Nc, u, ctx));
1209: } else {
1210: PetscScalar alpha = param->alpha; /* - */
1211: PetscScalar K_u = param->K_u; /* Pa */
1212: PetscScalar M = param->M; /* Pa */
1213: PetscScalar G = param->mu; /* Pa */
1214: PetscScalar P_0 = param->P_0; /* Pa */
1215: PetscReal R_0 = user->xmax[1]; /* m */
1216: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
1217: PetscInt N = user->niter, n;
1219: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
1220: PetscScalar eta = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G); /* -, Cheng (B.11) */
1221: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
1222: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
1223: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
1224: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
1225: PetscReal R = PetscSqrtReal(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
1227: PetscReal R_star = R / R_0;
1228: PetscReal t_star = PetscRealPart(c * time) / PetscSqr(R_0);
1229: PetscReal A_x = 0.0;
1231: for (n = 1; n < N + 1; ++n) {
1232: const PetscReal x_n = user->zeroArray[n - 1];
1233: const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u));
1235: 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) */
1236: }
1237: u[0] = P_0 * A_x;
1238: }
1239: return PETSC_SUCCESS;
1240: }
1242: /* Boundary Kernels */
1243: static void 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[])
1244: {
1245: const PetscReal P = PetscRealPart(constants[5]);
1247: f0[0] = 0.0;
1248: f0[1] = P;
1249: }
1251: #if 0
1252: static void f0_mandel_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1253: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1254: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1255: PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1256: {
1257: // Uniform stress distribution
1258: /* PetscScalar xmax = 0.5;
1259: PetscScalar xmin = -0.5;
1260: PetscScalar ymax = 0.5;
1261: PetscScalar ymin = -0.5;
1262: PetscScalar P = constants[5];
1263: PetscScalar aL = (xmax - xmin) / 2.0;
1264: PetscScalar sigma_zz = -1.0*P / aL; */
1266: // Analytical (parabolic) stress distribution
1267: PetscReal a1, a2, am;
1268: PetscReal y1, y2, ym;
1270: PetscInt NITER = 500;
1271: PetscReal EPS = 0.000001;
1272: PetscReal zeroArray[500]; /* NITER */
1273: PetscReal xmax = 1.0;
1274: PetscReal xmin = 0.0;
1275: PetscReal ymax = 0.1;
1276: PetscReal ymin = 0.0;
1277: PetscReal lower[2], upper[2];
1279: lower[0] = xmin - (xmax - xmin) / 2.0;
1280: lower[1] = ymin - (ymax - ymin) / 2.0;
1281: upper[0] = xmax - (xmax - xmin) / 2.0;
1282: upper[1] = ymax - (ymax - ymin) / 2.0;
1284: xmin = lower[0];
1285: ymin = lower[1];
1286: xmax = upper[0];
1287: ymax = upper[1];
1289: PetscScalar G = constants[0];
1290: PetscScalar K_u = constants[1];
1291: PetscScalar alpha = constants[2];
1292: PetscScalar M = constants[3];
1293: PetscScalar kappa = constants[4];
1294: PetscScalar F = constants[5];
1296: PetscScalar K_d = K_u - alpha*alpha*M;
1297: PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G));
1298: PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G));
1299: PetscReal aL = (xmax - xmin) / 2.0;
1300: PetscReal c = PetscRealPart(((2.0*kappa*G) * (1.0 - nu) * (nu_u - nu)) / (alpha*alpha * (1.0 - 2.0*nu) * (1.0 - nu_u)));
1301: PetscScalar B = (3.0 * (nu_u - nu)) / ( alpha * (1.0 - 2.0*nu) * (1.0 + nu_u));
1302: PetscScalar A1 = 3.0 / (B * (1.0 + nu_u));
1303: PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu);
1305: // Generate zero values
1306: for (PetscInt i=1; i < NITER+1; i++)
1307: {
1308: a1 = ((PetscReal) i - 1.0) * PETSC_PI * PETSC_PI / 4.0 + EPS;
1309: a2 = a1 + PETSC_PI/2;
1310: for (PetscInt j=0; j<NITER; j++)
1311: {
1312: y1 = PetscTanReal(a1) - PetscRealPart(A1/A2)*a1;
1313: y2 = PetscTanReal(a2) - PetscRealPart(A1/A2)*a2;
1314: am = (a1 + a2)/2.0;
1315: ym = PetscTanReal(am) - PetscRealPart(A1/A2)*am;
1316: if ((ym*y1) > 0)
1317: {
1318: a1 = am;
1319: } else {
1320: a2 = am;
1321: }
1322: if (PetscAbsReal(y2) < EPS)
1323: {
1324: am = a2;
1325: }
1326: }
1327: zeroArray[i-1] = am;
1328: }
1330: // Solution for sigma_zz
1331: PetscScalar A_x = 0.0;
1332: PetscScalar B_x = 0.0;
1334: for (PetscInt n=1; n < NITER+1; n++)
1335: {
1336: PetscReal alpha_n = zeroArray[n-1];
1338: 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)));
1339: 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)));
1340: }
1342: PetscScalar sigma_zz = -1.0*(F/aL) - ((2.0*F)/aL) * (A2/A1) * A_x + ((2.0*F)/aL) * B_x;
1344: if (x[1] == ymax) {
1345: f0[1] += sigma_zz;
1346: } else if (x[1] == ymin) {
1347: f0[1] -= sigma_zz;
1348: }
1349: }
1350: #endif
1352: static void 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[])
1353: {
1354: const PetscReal P_0 = PetscRealPart(constants[5]);
1355: PetscInt d;
1357: for (d = 0; d < dim; ++d) f0[d] = -P_0 * n[d];
1358: }
1360: /* Standard Kernels - Residual */
1361: /* f0_e */
1362: static void 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[])
1363: {
1364: PetscInt d;
1366: for (d = 0; d < dim; ++d) f0[0] += u_x[d * dim + d];
1367: f0[0] -= u[uOff[1]];
1368: }
1370: /* f0_p */
1371: static void 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[])
1372: {
1373: const PetscReal alpha = PetscRealPart(constants[2]);
1374: const PetscReal M = PetscRealPart(constants[3]);
1376: f0[0] += alpha * u_t[uOff[1]];
1377: f0[0] += u_t[uOff[2]] / M;
1378: if (f0[0] != f0[0]) abort();
1379: }
1381: /* f1_u */
1382: static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
1383: {
1384: const PetscInt Nc = dim;
1385: const PetscReal G = PetscRealPart(constants[0]);
1386: const PetscReal K_u = PetscRealPart(constants[1]);
1387: const PetscReal alpha = PetscRealPart(constants[2]);
1388: const PetscReal M = PetscRealPart(constants[3]);
1389: const PetscReal K_d = K_u - alpha * alpha * M;
1390: const PetscReal lambda = K_d - (2.0 * G) / 3.0;
1391: PetscInt c, d;
1393: for (c = 0; c < Nc; ++c) {
1394: for (d = 0; d < dim; ++d) f1[c * dim + d] -= G * (u_x[c * dim + d] + u_x[d * dim + c]);
1395: f1[c * dim + c] -= lambda * u[uOff[1]];
1396: f1[c * dim + c] += alpha * u[uOff[2]];
1397: }
1398: }
1400: /* f1_p */
1401: static void 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[])
1402: {
1403: const PetscReal kappa = PetscRealPart(constants[4]);
1404: PetscInt d;
1406: for (d = 0; d < dim; ++d) f1[d] += kappa * u_x[uOff_x[2] + d];
1407: }
1409: /*
1410: \partial_df \phi_fc g_{fc,gc,df,dg} \partial_dg \phi_gc
1412: \partial_df \phi_fc \lambda \delta_{fc,df} \sum_gc \partial_dg \phi_gc \delta_{gc,dg}
1413: = \partial_fc \phi_fc \sum_gc \partial_gc \phi_gc
1414: */
1416: /* Standard Kernels - Jacobian */
1417: /* g0_ee */
1418: static void 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[])
1419: {
1420: g0[0] = -1.0;
1421: }
1423: /* g0_pe */
1424: static void 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[])
1425: {
1426: const PetscReal alpha = PetscRealPart(constants[2]);
1428: g0[0] = u_tShift * alpha;
1429: }
1431: /* g0_pp */
1432: static void 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[])
1433: {
1434: const PetscReal M = PetscRealPart(constants[3]);
1436: g0[0] = u_tShift / M;
1437: }
1439: /* g1_eu */
1440: static void 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[])
1441: {
1442: PetscInt d;
1443: for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */
1444: }
1446: /* g2_ue */
1447: static void 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[])
1448: {
1449: const PetscReal G = PetscRealPart(constants[0]);
1450: const PetscReal K_u = PetscRealPart(constants[1]);
1451: const PetscReal alpha = PetscRealPart(constants[2]);
1452: const PetscReal M = PetscRealPart(constants[3]);
1453: const PetscReal K_d = K_u - alpha * alpha * M;
1454: const PetscReal lambda = K_d - (2.0 * G) / 3.0;
1455: PetscInt d;
1457: for (d = 0; d < dim; ++d) g2[d * dim + d] -= lambda;
1458: }
1459: /* g2_up */
1460: static void 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[])
1461: {
1462: const PetscReal alpha = PetscRealPart(constants[2]);
1463: PetscInt d;
1465: for (d = 0; d < dim; ++d) g2[d * dim + d] += alpha;
1466: }
1468: /* g3_uu */
1469: static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
1470: {
1471: const PetscInt Nc = dim;
1472: const PetscReal G = PetscRealPart(constants[0]);
1473: PetscInt c, d;
1475: for (c = 0; c < Nc; ++c) {
1476: for (d = 0; d < dim; ++d) {
1477: g3[((c * Nc + c) * dim + d) * dim + d] -= G;
1478: g3[((c * Nc + d) * dim + d) * dim + c] -= G;
1479: }
1480: }
1481: }
1483: /* g3_pp */
1484: static void 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[])
1485: {
1486: const PetscReal kappa = PetscRealPart(constants[4]);
1487: PetscInt d;
1489: for (d = 0; d < dim; ++d) g3[d * dim + d] += kappa;
1490: }
1492: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
1493: {
1494: PetscInt sol;
1496: PetscFunctionBeginUser;
1497: options->solType = SOL_QUADRATIC_TRIG;
1498: options->niter = 500;
1499: options->eps = PETSC_SMALL;
1500: options->dtInitial = -1.0;
1501: PetscOptionsBegin(comm, "", "Biot Poroelasticity Options", "DMPLEX");
1502: PetscCall(PetscOptionsInt("-niter", "Number of series term iterations in exact solutions", "ex53.c", options->niter, &options->niter, NULL));
1503: sol = options->solType;
1504: PetscCall(PetscOptionsEList("-sol_type", "Type of exact solution", "ex53.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL));
1505: options->solType = (SolutionType)sol;
1506: PetscCall(PetscOptionsReal("-eps", "Precision value for root finding", "ex53.c", options->eps, &options->eps, NULL));
1507: PetscCall(PetscOptionsReal("-dt_initial", "Override the initial timestep", "ex53.c", options->dtInitial, &options->dtInitial, NULL));
1508: PetscOptionsEnd();
1509: PetscFunctionReturn(PETSC_SUCCESS);
1510: }
1512: static PetscErrorCode mandelZeros(MPI_Comm comm, AppCtx *ctx, Parameter *param)
1513: {
1514: //PetscBag bag;
1515: PetscReal a1, a2, am;
1516: PetscReal y1, y2, ym;
1518: PetscFunctionBeginUser;
1519: //PetscCall(PetscBagGetData(ctx->bag, (void **) ¶m));
1520: PetscInt NITER = ctx->niter;
1521: PetscReal EPS = ctx->eps;
1522: //const PetscScalar YMAX = param->ymax;
1523: //const PetscScalar YMIN = param->ymin;
1524: PetscScalar alpha = param->alpha;
1525: PetscScalar K_u = param->K_u;
1526: PetscScalar M = param->M;
1527: PetscScalar G = param->mu;
1528: //const PetscScalar k = param->k;
1529: //const PetscScalar mu_f = param->mu_f;
1530: //const PetscScalar P_0 = param->P_0;
1532: PetscScalar K_d = K_u - alpha * alpha * M;
1533: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
1534: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
1535: //const PetscScalar kappa = k / mu_f;
1537: // Generate zero values
1538: for (PetscInt i = 1; i < NITER + 1; i++) {
1539: a1 = ((PetscReal)i - 1.0) * PETSC_PI * PETSC_PI / 4.0 + EPS;
1540: a2 = a1 + PETSC_PI / 2;
1541: am = a1;
1542: for (PetscInt j = 0; j < NITER; j++) {
1543: y1 = PetscTanReal(a1) - PetscRealPart((1.0 - nu) / (nu_u - nu)) * a1;
1544: y2 = PetscTanReal(a2) - PetscRealPart((1.0 - nu) / (nu_u - nu)) * a2;
1545: am = (a1 + a2) / 2.0;
1546: ym = PetscTanReal(am) - PetscRealPart((1.0 - nu) / (nu_u - nu)) * am;
1547: if ((ym * y1) > 0) {
1548: a1 = am;
1549: } else {
1550: a2 = am;
1551: }
1552: if (PetscAbsReal(y2) < EPS) am = a2;
1553: }
1554: ctx->zeroArray[i - 1] = am;
1555: }
1556: PetscFunctionReturn(PETSC_SUCCESS);
1557: }
1559: static PetscReal CryerFunction(PetscReal nu_u, PetscReal nu, PetscReal x)
1560: {
1561: return PetscTanReal(PetscSqrtReal(x)) * (6.0 * (nu_u - nu) - (1.0 - nu) * (1.0 + nu_u) * x) - (6.0 * (nu_u - nu) * PetscSqrtReal(x));
1562: }
1564: static PetscErrorCode cryerZeros(MPI_Comm comm, AppCtx *ctx, Parameter *param)
1565: {
1566: PetscReal alpha = PetscRealPart(param->alpha); /* - */
1567: PetscReal K_u = PetscRealPart(param->K_u); /* Pa */
1568: PetscReal M = PetscRealPart(param->M); /* Pa */
1569: PetscReal G = PetscRealPart(param->mu); /* Pa */
1570: PetscInt N = ctx->niter, n;
1572: PetscReal K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
1573: PetscReal nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
1574: PetscReal nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
1576: PetscFunctionBeginUser;
1577: for (n = 1; n < N + 1; ++n) {
1578: PetscReal tol = PetscPowReal(n, 1.5) * ctx->eps;
1579: PetscReal a1 = 0., a2 = 0., am = 0.;
1580: PetscReal y1, y2, ym;
1581: PetscInt j, k = n - 1;
1583: y1 = y2 = 1.;
1584: while (y1 * y2 > 0) {
1585: ++k;
1586: a1 = PetscSqr(n * PETSC_PI) - k * PETSC_PI;
1587: a2 = PetscSqr(n * PETSC_PI) + k * PETSC_PI;
1588: y1 = CryerFunction(nu_u, nu, a1);
1589: y2 = CryerFunction(nu_u, nu, a2);
1590: }
1591: for (j = 0; j < 50000; ++j) {
1592: y1 = CryerFunction(nu_u, nu, a1);
1593: y2 = CryerFunction(nu_u, nu, a2);
1594: PetscCheck(y1 * y2 <= 0, comm, PETSC_ERR_PLIB, "Invalid root finding initialization for root %" PetscInt_FMT ", (%g, %g)--(%g, %g)", n, (double)a1, (double)y1, (double)a2, (double)y2);
1595: am = (a1 + a2) / 2.0;
1596: ym = CryerFunction(nu_u, nu, am);
1597: if ((ym * y1) < 0) a2 = am;
1598: else a1 = am;
1599: if (PetscAbsReal(ym) < tol) break;
1600: }
1601: PetscCheck(PetscAbsReal(ym) < tol, comm, PETSC_ERR_PLIB, "Root finding did not converge for root %" PetscInt_FMT " (%g)", n, (double)PetscAbsReal(ym));
1602: ctx->zeroArray[n - 1] = am;
1603: }
1604: PetscFunctionReturn(PETSC_SUCCESS);
1605: }
1607: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
1608: {
1609: PetscBag bag;
1610: Parameter *p;
1612: PetscFunctionBeginUser;
1613: /* setup PETSc parameter bag */
1614: PetscCall(PetscBagGetData(ctx->bag, (void **)&p));
1615: PetscCall(PetscBagSetName(ctx->bag, "par", "Poroelastic Parameters"));
1616: bag = ctx->bag;
1617: if (ctx->solType == SOL_TERZAGHI) {
1618: // Realistic values - Terzaghi
1619: PetscCall(PetscBagRegisterScalar(bag, &p->mu, 3.0, "mu", "Shear Modulus, Pa"));
1620: PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 9.76, "K_u", "Undrained Bulk Modulus, Pa"));
1621: PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -"));
1622: PetscCall(PetscBagRegisterScalar(bag, &p->M, 16.0, "M", "Biot Modulus, Pa"));
1623: PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2"));
1624: PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s"));
1625: PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa"));
1626: } else if (ctx->solType == SOL_MANDEL) {
1627: // Realistic values - Mandel
1628: PetscCall(PetscBagRegisterScalar(bag, &p->mu, 0.75, "mu", "Shear Modulus, Pa"));
1629: PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 2.6941176470588233, "K_u", "Undrained Bulk Modulus, Pa"));
1630: PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -"));
1631: PetscCall(PetscBagRegisterScalar(bag, &p->M, 4.705882352941176, "M", "Biot Modulus, Pa"));
1632: PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2"));
1633: PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s"));
1634: PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa"));
1635: } else if (ctx->solType == SOL_CRYER) {
1636: // Realistic values - Mandel
1637: PetscCall(PetscBagRegisterScalar(bag, &p->mu, 0.75, "mu", "Shear Modulus, Pa"));
1638: PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 2.6941176470588233, "K_u", "Undrained Bulk Modulus, Pa"));
1639: PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -"));
1640: PetscCall(PetscBagRegisterScalar(bag, &p->M, 4.705882352941176, "M", "Biot Modulus, Pa"));
1641: PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2"));
1642: PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s"));
1643: PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa"));
1644: } else {
1645: // Nonsense values
1646: PetscCall(PetscBagRegisterScalar(bag, &p->mu, 1.0, "mu", "Shear Modulus, Pa"));
1647: PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 1.0, "K_u", "Undrained Bulk Modulus, Pa"));
1648: PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 1.0, "alpha", "Biot Effective Stress Coefficient, -"));
1649: PetscCall(PetscBagRegisterScalar(bag, &p->M, 1.0, "M", "Biot Modulus, Pa"));
1650: PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.0, "k", "Isotropic Permeability, m**2"));
1651: PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s"));
1652: PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa"));
1653: }
1654: PetscCall(PetscBagSetFromOptions(bag));
1655: {
1656: PetscScalar K_d = p->K_u - p->alpha * p->alpha * p->M;
1657: PetscScalar nu_u = (3.0 * p->K_u - 2.0 * p->mu) / (2.0 * (3.0 * p->K_u + p->mu));
1658: PetscScalar nu = (3.0 * K_d - 2.0 * p->mu) / (2.0 * (3.0 * K_d + p->mu));
1659: PetscScalar S = (3.0 * p->K_u + 4.0 * p->mu) / (p->M * (3.0 * K_d + 4.0 * p->mu));
1660: PetscReal c = PetscRealPart((p->k / p->mu_f) / S);
1662: PetscViewer viewer;
1663: PetscViewerFormat format;
1664: PetscBool flg;
1666: switch (ctx->solType) {
1667: case SOL_QUADRATIC_LINEAR:
1668: case SOL_QUADRATIC_TRIG:
1669: case SOL_TRIG_LINEAR:
1670: ctx->t_r = PetscSqr(ctx->xmax[0] - ctx->xmin[0]) / c;
1671: break;
1672: case SOL_TERZAGHI:
1673: ctx->t_r = PetscSqr(2.0 * (ctx->xmax[1] - ctx->xmin[1])) / c;
1674: break;
1675: case SOL_MANDEL:
1676: ctx->t_r = PetscSqr(2.0 * (ctx->xmax[1] - ctx->xmin[1])) / c;
1677: break;
1678: case SOL_CRYER:
1679: ctx->t_r = PetscSqr(ctx->xmax[1]) / c;
1680: break;
1681: default:
1682: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%d)", solutionTypes[PetscMin(ctx->solType, NUM_SOLUTION_TYPES)], ctx->solType);
1683: }
1684: PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg));
1685: if (flg) {
1686: PetscCall(PetscViewerPushFormat(viewer, format));
1687: PetscCall(PetscBagView(bag, viewer));
1688: PetscCall(PetscViewerFlush(viewer));
1689: PetscCall(PetscViewerPopFormat(viewer));
1690: PetscCall(PetscViewerDestroy(&viewer));
1691: PetscCall(PetscPrintf(comm, " Max displacement: %g %g\n", (double)PetscRealPart(p->P_0 * (ctx->xmax[1] - ctx->xmin[1]) * (1. - 2. * nu_u) / (2. * p->mu * (1. - nu_u))), (double)PetscRealPart(p->P_0 * (ctx->xmax[1] - ctx->xmin[1]) * (1. - 2. * nu) / (2. * p->mu * (1. - nu)))));
1692: PetscCall(PetscPrintf(comm, " Relaxation time: %g\n", (double)ctx->t_r));
1693: }
1694: }
1695: PetscFunctionReturn(PETSC_SUCCESS);
1696: }
1698: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
1699: {
1700: PetscFunctionBeginUser;
1701: PetscCall(DMCreate(comm, dm));
1702: PetscCall(DMSetType(*dm, DMPLEX));
1703: PetscCall(DMSetFromOptions(*dm));
1704: PetscCall(DMSetApplicationContext(*dm, user));
1705: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
1706: PetscCall(DMGetBoundingBox(*dm, user->xmin, user->xmax));
1707: PetscFunctionReturn(PETSC_SUCCESS);
1708: }
1710: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
1711: {
1712: PetscErrorCode (*exact[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
1713: PetscErrorCode (*exact_t[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
1714: PetscDS ds;
1715: DMLabel label;
1716: PetscWeakForm wf;
1717: Parameter *param;
1718: PetscInt id_mandel[2];
1719: PetscInt comp[1];
1720: PetscInt comp_mandel[2];
1721: PetscInt dim, id, bd, f;
1723: PetscFunctionBeginUser;
1724: PetscCall(DMGetLabel(dm, "marker", &label));
1725: PetscCall(DMGetDS(dm, &ds));
1726: PetscCall(PetscDSGetSpatialDimension(ds, &dim));
1727: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1728: exact_t[0] = exact_t[1] = exact_t[2] = zero;
1730: /* Setup Problem Formulation and Boundary Conditions */
1731: switch (user->solType) {
1732: case SOL_QUADRATIC_LINEAR:
1733: PetscCall(PetscDSSetResidual(ds, 0, f0_quadratic_linear_u, f1_u));
1734: PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1735: PetscCall(PetscDSSetResidual(ds, 2, f0_quadratic_linear_p, f1_p));
1736: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1737: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1738: PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1739: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1740: PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1741: PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1742: PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));
1743: exact[0] = quadratic_u;
1744: exact[1] = linear_eps;
1745: exact[2] = linear_linear_p;
1746: exact_t[2] = linear_linear_p_t;
1748: id = 1;
1749: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", label, 1, &id, 0, 0, NULL, (void (*)(void))exact[0], NULL, user, NULL));
1750: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", label, 1, &id, 2, 0, NULL, (void (*)(void))exact[2], (void (*)(void))exact_t[2], user, NULL));
1751: break;
1752: case SOL_TRIG_LINEAR:
1753: PetscCall(PetscDSSetResidual(ds, 0, f0_trig_linear_u, f1_u));
1754: PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1755: PetscCall(PetscDSSetResidual(ds, 2, f0_trig_linear_p, f1_p));
1756: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1757: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1758: PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1759: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1760: PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1761: PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1762: PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));
1763: exact[0] = trig_u;
1764: exact[1] = trig_eps;
1765: exact[2] = trig_linear_p;
1766: exact_t[2] = trig_linear_p_t;
1768: id = 1;
1769: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", label, 1, &id, 0, 0, NULL, (void (*)(void))exact[0], NULL, user, NULL));
1770: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", label, 1, &id, 2, 0, NULL, (void (*)(void))exact[2], (void (*)(void))exact_t[2], user, NULL));
1771: break;
1772: case SOL_QUADRATIC_TRIG:
1773: PetscCall(PetscDSSetResidual(ds, 0, f0_quadratic_trig_u, f1_u));
1774: PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1775: PetscCall(PetscDSSetResidual(ds, 2, f0_quadratic_trig_p, f1_p));
1776: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1777: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1778: PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1779: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1780: PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1781: PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1782: PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));
1783: exact[0] = quadratic_u;
1784: exact[1] = linear_eps;
1785: exact[2] = linear_trig_p;
1786: exact_t[2] = linear_trig_p_t;
1788: id = 1;
1789: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", label, 1, &id, 0, 0, NULL, (void (*)(void))exact[0], NULL, user, NULL));
1790: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", label, 1, &id, 2, 0, NULL, (void (*)(void))exact[2], (void (*)(void))exact_t[2], user, NULL));
1791: break;
1792: case SOL_TERZAGHI:
1793: PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u));
1794: PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1795: PetscCall(PetscDSSetResidual(ds, 2, f0_p, f1_p));
1796: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1797: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1798: PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1799: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1800: PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1801: PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1802: PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));
1804: exact[0] = terzaghi_2d_u;
1805: exact[1] = terzaghi_2d_eps;
1806: exact[2] = terzaghi_2d_p;
1807: exact_t[0] = terzaghi_2d_u_t;
1808: exact_t[1] = terzaghi_2d_eps_t;
1809: exact_t[2] = terzaghi_2d_p_t;
1811: id = 1;
1812: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "vertical stress", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
1813: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1814: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_terzaghi_bd_u, 0, NULL));
1816: id = 3;
1817: comp[0] = 1;
1818: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed base", label, 1, &id, 0, 1, comp, (void (*)(void))zero, NULL, user, NULL));
1819: id = 2;
1820: comp[0] = 0;
1821: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed side", label, 1, &id, 0, 1, comp, (void (*)(void))zero, NULL, user, NULL));
1822: id = 4;
1823: comp[0] = 0;
1824: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed side", label, 1, &id, 0, 1, comp, (void (*)(void))zero, NULL, user, NULL));
1825: id = 1;
1826: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", label, 1, &id, 2, 0, NULL, (void (*)(void))terzaghi_drainage_pressure, NULL, user, NULL));
1827: break;
1828: case SOL_MANDEL:
1829: PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u));
1830: PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1831: PetscCall(PetscDSSetResidual(ds, 2, f0_p, f1_p));
1832: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1833: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1834: PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1835: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1836: PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1837: PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1838: PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));
1840: PetscCall(mandelZeros(PETSC_COMM_WORLD, user, param));
1842: exact[0] = mandel_2d_u;
1843: exact[1] = mandel_2d_eps;
1844: exact[2] = mandel_2d_p;
1845: exact_t[0] = mandel_2d_u_t;
1846: exact_t[1] = mandel_2d_eps_t;
1847: exact_t[2] = mandel_2d_p_t;
1849: id_mandel[0] = 3;
1850: id_mandel[1] = 1;
1851: //comp[0] = 1;
1852: comp_mandel[0] = 0;
1853: comp_mandel[1] = 1;
1854: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "vertical stress", label, 2, id_mandel, 0, 2, comp_mandel, (void (*)(void))mandel_2d_u, NULL, user, NULL));
1855: //PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "vertical stress", "marker", 0, 1, comp, NULL, 2, id_mandel, user));
1856: //PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed base", "marker", 0, 1, comp, (void (*)(void)) zero, 2, id_mandel, user));
1857: //PetscCall(PetscDSSetBdResidual(ds, 0, f0_mandel_bd_u, NULL));
1859: id_mandel[0] = 2;
1860: id_mandel[1] = 4;
1861: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", label, 2, id_mandel, 2, 0, NULL, (void (*)(void))zero, NULL, user, NULL));
1862: break;
1863: case SOL_CRYER:
1864: PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u));
1865: PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1866: PetscCall(PetscDSSetResidual(ds, 2, f0_p, f1_p));
1867: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1868: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1869: PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1870: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1871: PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1872: PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1873: PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));
1875: PetscCall(cryerZeros(PETSC_COMM_WORLD, user, param));
1877: exact[0] = cryer_3d_u;
1878: exact[1] = cryer_3d_eps;
1879: exact[2] = cryer_3d_p;
1881: id = 1;
1882: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "normal stress", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
1883: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1884: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_cryer_bd_u, 0, NULL));
1886: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", label, 1, &id, 2, 0, NULL, (void (*)(void))cryer_drainage_pressure, NULL, user, NULL));
1887: break;
1888: default:
1889: SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%d)", solutionTypes[PetscMin(user->solType, NUM_SOLUTION_TYPES)], user->solType);
1890: }
1891: for (f = 0; f < 3; ++f) {
1892: PetscCall(PetscDSSetExactSolution(ds, f, exact[f], user));
1893: PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, f, exact_t[f], user));
1894: }
1896: /* Setup constants */
1897: {
1898: PetscScalar constants[6];
1899: constants[0] = param->mu; /* shear modulus, Pa */
1900: constants[1] = param->K_u; /* undrained bulk modulus, Pa */
1901: constants[2] = param->alpha; /* Biot effective stress coefficient, - */
1902: constants[3] = param->M; /* Biot modulus, Pa */
1903: constants[4] = param->k / param->mu_f; /* Darcy coefficient, m**2 / Pa*s */
1904: constants[5] = param->P_0; /* Magnitude of Vertical Stress, Pa */
1905: PetscCall(PetscDSSetConstants(ds, 6, constants));
1906: }
1907: PetscFunctionReturn(PETSC_SUCCESS);
1908: }
1910: static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace)
1911: {
1912: PetscFunctionBeginUser;
1913: PetscCall(DMPlexCreateRigidBody(dm, origField, nullspace));
1914: PetscFunctionReturn(PETSC_SUCCESS);
1915: }
1917: static PetscErrorCode SetupFE(DM dm, PetscInt Nf, PetscInt Nc[], const char *name[], PetscErrorCode (*setup)(DM, AppCtx *), void *ctx)
1918: {
1919: AppCtx *user = (AppCtx *)ctx;
1920: DM cdm = dm;
1921: PetscFE fe;
1922: PetscQuadrature q = NULL, fq = NULL;
1923: char prefix[PETSC_MAX_PATH_LEN];
1924: PetscInt dim, f;
1925: PetscBool simplex;
1927: PetscFunctionBeginUser;
1928: /* Create finite element */
1929: PetscCall(DMGetDimension(dm, &dim));
1930: PetscCall(DMPlexIsSimplex(dm, &simplex));
1931: for (f = 0; f < Nf; ++f) {
1932: PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name[f]));
1933: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc[f], simplex, name[f] ? prefix : NULL, -1, &fe));
1934: PetscCall(PetscObjectSetName((PetscObject)fe, name[f]));
1935: if (!q) PetscCall(PetscFEGetQuadrature(fe, &q));
1936: if (!fq) PetscCall(PetscFEGetFaceQuadrature(fe, &fq));
1937: PetscCall(PetscFESetQuadrature(fe, q));
1938: PetscCall(PetscFESetFaceQuadrature(fe, fq));
1939: PetscCall(DMSetField(dm, f, NULL, (PetscObject)fe));
1940: PetscCall(PetscFEDestroy(&fe));
1941: }
1942: PetscCall(DMCreateDS(dm));
1943: PetscCall((*setup)(dm, user));
1944: while (cdm) {
1945: PetscCall(DMCopyDisc(dm, cdm));
1946: if (0) PetscCall(DMSetNearNullSpaceConstructor(cdm, 0, CreateElasticityNullSpace));
1947: /* TODO: Check whether the boundary of coarse meshes is marked */
1948: PetscCall(DMGetCoarseDM(cdm, &cdm));
1949: }
1950: PetscCall(PetscFEDestroy(&fe));
1951: PetscFunctionReturn(PETSC_SUCCESS);
1952: }
1954: static PetscErrorCode SetInitialConditions(TS ts, Vec u)
1955: {
1956: DM dm;
1957: PetscReal t;
1959: PetscFunctionBeginUser;
1960: PetscCall(TSGetDM(ts, &dm));
1961: PetscCall(TSGetTime(ts, &t));
1962: if (t <= 0.0) {
1963: PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
1964: void *ctxs[3];
1965: AppCtx *ctx;
1967: PetscCall(DMGetApplicationContext(dm, &ctx));
1968: switch (ctx->solType) {
1969: case SOL_TERZAGHI:
1970: funcs[0] = terzaghi_initial_u;
1971: ctxs[0] = ctx;
1972: funcs[1] = terzaghi_initial_eps;
1973: ctxs[1] = ctx;
1974: funcs[2] = terzaghi_drainage_pressure;
1975: ctxs[2] = ctx;
1976: PetscCall(DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u));
1977: break;
1978: case SOL_MANDEL:
1979: funcs[0] = mandel_initial_u;
1980: ctxs[0] = ctx;
1981: funcs[1] = mandel_initial_eps;
1982: ctxs[1] = ctx;
1983: funcs[2] = mandel_drainage_pressure;
1984: ctxs[2] = ctx;
1985: PetscCall(DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u));
1986: break;
1987: case SOL_CRYER:
1988: funcs[0] = cryer_initial_u;
1989: ctxs[0] = ctx;
1990: funcs[1] = cryer_initial_eps;
1991: ctxs[1] = ctx;
1992: funcs[2] = cryer_drainage_pressure;
1993: ctxs[2] = ctx;
1994: PetscCall(DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u));
1995: break;
1996: default:
1997: PetscCall(DMComputeExactSolution(dm, t, u, NULL));
1998: }
1999: } else {
2000: PetscCall(DMComputeExactSolution(dm, t, u, NULL));
2001: }
2002: PetscFunctionReturn(PETSC_SUCCESS);
2003: }
2005: /* Need to create Viewer each time because HDF5 can get corrupted */
2006: static PetscErrorCode SolutionMonitor(TS ts, PetscInt steps, PetscReal time, Vec u, void *mctx)
2007: {
2008: DM dm;
2009: Vec exact;
2010: PetscViewer viewer;
2011: PetscViewerFormat format;
2012: PetscOptions options;
2013: const char *prefix;
2015: PetscFunctionBeginUser;
2016: PetscCall(TSGetDM(ts, &dm));
2017: PetscCall(PetscObjectGetOptions((PetscObject)ts, &options));
2018: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, &prefix));
2019: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), options, prefix, "-monitor_solution", &viewer, &format, NULL));
2020: PetscCall(DMGetGlobalVector(dm, &exact));
2021: PetscCall(DMComputeExactSolution(dm, time, exact, NULL));
2022: PetscCall(DMSetOutputSequenceNumber(dm, steps, time));
2023: PetscCall(VecView(exact, viewer));
2024: PetscCall(VecView(u, viewer));
2025: PetscCall(DMRestoreGlobalVector(dm, &exact));
2026: {
2027: PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
2028: void **ectxs;
2029: PetscReal *err;
2030: PetscInt Nf, f;
2032: PetscCall(DMGetNumFields(dm, &Nf));
2033: PetscCall(PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err));
2034: {
2035: PetscInt Nds, s;
2037: PetscCall(DMGetNumDS(dm, &Nds));
2038: for (s = 0; s < Nds; ++s) {
2039: PetscDS ds;
2040: DMLabel label;
2041: IS fieldIS;
2042: const PetscInt *fields;
2043: PetscInt dsNf, f;
2045: PetscCall(DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds, NULL));
2046: PetscCall(PetscDSGetNumFields(ds, &dsNf));
2047: PetscCall(ISGetIndices(fieldIS, &fields));
2048: for (f = 0; f < dsNf; ++f) {
2049: const PetscInt field = fields[f];
2050: PetscCall(PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]));
2051: }
2052: PetscCall(ISRestoreIndices(fieldIS, &fields));
2053: }
2054: }
2055: PetscCall(DMComputeL2FieldDiff(dm, time, exacts, ectxs, u, err));
2056: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "Time: %g L_2 Error: [", (double)time));
2057: for (f = 0; f < Nf; ++f) {
2058: if (f) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), ", "));
2059: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%g", (double)err[f]));
2060: }
2061: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "]\n"));
2062: PetscCall(PetscFree3(exacts, ectxs, err));
2063: }
2064: PetscCall(PetscViewerDestroy(&viewer));
2065: PetscFunctionReturn(PETSC_SUCCESS);
2066: }
2068: static PetscErrorCode SetupMonitor(TS ts, AppCtx *ctx)
2069: {
2070: PetscViewer viewer;
2071: PetscViewerFormat format;
2072: PetscOptions options;
2073: const char *prefix;
2074: PetscBool flg;
2076: PetscFunctionBeginUser;
2077: PetscCall(PetscObjectGetOptions((PetscObject)ts, &options));
2078: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, &prefix));
2079: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), options, prefix, "-monitor_solution", &viewer, &format, &flg));
2080: if (flg) PetscCall(TSMonitorSet(ts, SolutionMonitor, ctx, NULL));
2081: PetscCall(PetscViewerDestroy(&viewer));
2082: PetscFunctionReturn(PETSC_SUCCESS);
2083: }
2085: static PetscErrorCode TSAdaptChoose_Terzaghi(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept, PetscReal *wlte, PetscReal *wltea, PetscReal *wlter)
2086: {
2087: static PetscReal dtTarget = -1.0;
2088: PetscReal dtInitial;
2089: DM dm;
2090: AppCtx *ctx;
2091: PetscInt step;
2093: PetscFunctionBeginUser;
2094: PetscCall(TSGetDM(ts, &dm));
2095: PetscCall(DMGetApplicationContext(dm, &ctx));
2096: PetscCall(TSGetStepNumber(ts, &step));
2097: dtInitial = ctx->dtInitial < 0.0 ? 1.0e-4 * ctx->t_r : ctx->dtInitial;
2098: if (!step) {
2099: if (PetscAbsReal(dtInitial - h) > PETSC_SMALL) {
2100: *accept = PETSC_FALSE;
2101: *next_h = dtInitial;
2102: dtTarget = h;
2103: } else {
2104: *accept = PETSC_TRUE;
2105: *next_h = dtTarget < 0.0 ? dtInitial : dtTarget;
2106: dtTarget = -1.0;
2107: }
2108: } else {
2109: *accept = PETSC_TRUE;
2110: *next_h = h;
2111: }
2112: *next_sc = 0; /* Reuse the same order scheme */
2113: *wlte = -1; /* Weighted local truncation error was not evaluated */
2114: *wltea = -1; /* Weighted absolute local truncation error was not evaluated */
2115: *wlter = -1; /* Weighted relative local truncation error was not evaluated */
2116: PetscFunctionReturn(PETSC_SUCCESS);
2117: }
2119: int main(int argc, char **argv)
2120: {
2121: AppCtx ctx; /* User-defined work context */
2122: DM dm; /* Problem specification */
2123: TS ts; /* Time Series / Nonlinear solver */
2124: Vec u; /* Solutions */
2125: const char *name[3] = {"displacement", "tracestrain", "pressure"};
2126: PetscReal t;
2127: PetscInt dim, Nc[3];
2129: PetscFunctionBeginUser;
2130: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
2131: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
2132: PetscCall(PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &ctx.bag));
2133: PetscCall(PetscMalloc1(ctx.niter, &ctx.zeroArray));
2134: PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm));
2135: PetscCall(SetupParameters(PETSC_COMM_WORLD, &ctx));
2136: /* Primal System */
2137: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
2138: PetscCall(DMSetApplicationContext(dm, &ctx));
2139: PetscCall(TSSetDM(ts, dm));
2141: PetscCall(DMGetDimension(dm, &dim));
2142: Nc[0] = dim;
2143: Nc[1] = 1;
2144: Nc[2] = 1;
2146: PetscCall(SetupFE(dm, 3, Nc, name, SetupPrimalProblem, &ctx));
2147: PetscCall(DMCreateGlobalVector(dm, &u));
2148: PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx));
2149: PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx));
2150: PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx));
2151: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
2152: PetscCall(TSSetFromOptions(ts));
2153: PetscCall(TSSetComputeInitialCondition(ts, SetInitialConditions));
2154: PetscCall(SetupMonitor(ts, &ctx));
2156: if (ctx.solType != SOL_QUADRATIC_TRIG) {
2157: TSAdapt adapt;
2159: PetscCall(TSGetAdapt(ts, &adapt));
2160: adapt->ops->choose = TSAdaptChoose_Terzaghi;
2161: }
2162: if (ctx.solType == SOL_CRYER) {
2163: Mat J;
2164: MatNullSpace sp;
2166: PetscCall(TSSetUp(ts));
2167: PetscCall(TSGetIJacobian(ts, &J, NULL, NULL, NULL));
2168: PetscCall(DMPlexCreateRigidBody(dm, 0, &sp));
2169: PetscCall(MatSetNullSpace(J, sp));
2170: PetscCall(MatNullSpaceDestroy(&sp));
2171: }
2172: PetscCall(TSGetTime(ts, &t));
2173: PetscCall(DMSetOutputSequenceNumber(dm, 0, t));
2174: PetscCall(DMTSCheckFromOptions(ts, u));
2175: PetscCall(SetInitialConditions(ts, u));
2176: PetscCall(PetscObjectSetName((PetscObject)u, "solution"));
2177: PetscCall(TSSolve(ts, u));
2178: PetscCall(DMTSCheckFromOptions(ts, u));
2179: PetscCall(TSGetSolution(ts, &u));
2180: PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));
2182: /* Cleanup */
2183: PetscCall(VecDestroy(&u));
2184: PetscCall(TSDestroy(&ts));
2185: PetscCall(DMDestroy(&dm));
2186: PetscCall(PetscBagDestroy(&ctx.bag));
2187: PetscCall(PetscFree(ctx.zeroArray));
2188: PetscCall(PetscFinalize());
2189: return 0;
2190: }
2192: /*TEST
2194: test:
2195: suffix: 2d_quad_linear
2196: requires: triangle
2197: args: -sol_type quadratic_linear -dm_refine 2 \
2198: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2199: -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme
2201: test:
2202: suffix: 3d_quad_linear
2203: requires: ctetgen
2204: args: -dm_plex_dim 3 -sol_type quadratic_linear -dm_refine 1 \
2205: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2206: -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme
2208: test:
2209: suffix: 2d_trig_linear
2210: requires: triangle
2211: args: -sol_type trig_linear -dm_refine 1 \
2212: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2213: -dmts_check .0001 -ts_max_steps 5 -ts_dt 0.00001 -ts_monitor_extreme
2215: test:
2216: # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9, 2.1, 1.8]
2217: suffix: 2d_trig_linear_sconv
2218: requires: triangle
2219: args: -sol_type trig_linear -dm_refine 1 \
2220: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2221: -convest_num_refine 1 -ts_convergence_estimate -ts_convergence_temporal 0 -ts_max_steps 1 -ts_dt 0.00001 -pc_type lu
2223: test:
2224: suffix: 3d_trig_linear
2225: requires: ctetgen
2226: args: -dm_plex_dim 3 -sol_type trig_linear -dm_refine 1 \
2227: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2228: -dmts_check .0001 -ts_max_steps 2 -ts_monitor_extreme
2230: test:
2231: # -dm_refine 1 -convest_num_refine 2 gets L_2 convergence rate: [2.0, 2.1, 1.9]
2232: suffix: 3d_trig_linear_sconv
2233: requires: ctetgen
2234: args: -dm_plex_dim 3 -sol_type trig_linear -dm_refine 1 \
2235: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2236: -convest_num_refine 1 -ts_convergence_estimate -ts_convergence_temporal 0 -ts_max_steps 1 -pc_type lu
2238: test:
2239: suffix: 2d_quad_trig
2240: requires: triangle
2241: args: -sol_type quadratic_trig -dm_refine 2 \
2242: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2243: -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme
2245: test:
2246: # Using -dm_refine 4 gets the convergence rates to [0.95, 0.97, 0.90]
2247: suffix: 2d_quad_trig_tconv
2248: requires: triangle
2249: args: -sol_type quadratic_trig -dm_refine 1 \
2250: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2251: -convest_num_refine 3 -ts_convergence_estimate -ts_max_steps 5 -pc_type lu
2253: test:
2254: suffix: 3d_quad_trig
2255: requires: ctetgen
2256: args: -dm_plex_dim 3 -sol_type quadratic_trig -dm_refine 1 \
2257: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2258: -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme
2260: test:
2261: # Using -dm_refine 2 -convest_num_refine 3 gets the convergence rates to [1.0, 1.0, 1.0]
2262: suffix: 3d_quad_trig_tconv
2263: requires: ctetgen
2264: args: -dm_plex_dim 3 -sol_type quadratic_trig -dm_refine 1 \
2265: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2266: -convest_num_refine 1 -ts_convergence_estimate -ts_max_steps 5 -pc_type lu
2268: testset:
2269: args: -sol_type terzaghi -dm_plex_simplex 0 -dm_plex_box_faces 1,8 -dm_plex_box_lower 0,0 -dm_plex_box_upper 10,10 -dm_plex_separate_marker \
2270: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 -niter 16000 \
2271: -pc_type lu
2273: test:
2274: suffix: 2d_terzaghi
2275: requires: double
2276: args: -ts_dt 0.0028666667 -ts_max_steps 2 -ts_monitor -dmts_check .0001
2278: test:
2279: # -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]
2280: suffix: 2d_terzaghi_tconv
2281: args: -ts_dt 0.023 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1
2283: test:
2284: # -dm_plex_box_faces 1,16 -convest_num_refine 4 gives L_2 convergence rate: [1.7, 1.2, 1.1]
2285: # if we add -displacement_petscspace_degree 3 -tracestrain_petscspace_degree 2 -pressure_petscspace_degree 2, we get [2.1, 1.6, 1.5]
2286: suffix: 2d_terzaghi_sconv
2287: args: -ts_dt 1e-5 -dt_initial 1e-5 -ts_max_steps 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1
2289: testset:
2290: args: -sol_type mandel -dm_plex_simplex 0 -dm_plex_box_lower -0.5,-0.125 -dm_plex_box_upper 0.5,0.125 -dm_plex_separate_marker -dm_refine 1 \
2291: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2292: -pc_type lu
2294: test:
2295: suffix: 2d_mandel
2296: requires: double
2297: args: -ts_dt 0.0028666667 -ts_max_steps 2 -ts_monitor -dmts_check .0001
2299: test:
2300: # -dm_refine 3 -ts_max_steps 4 -convest_num_refine 3 gives L_2 convergence rate: [1.6, 0.93, 1.2]
2301: suffix: 2d_mandel_sconv
2302: args: -ts_dt 1e-5 -dt_initial 1e-5 -ts_max_steps 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1
2304: test:
2305: # -dm_refine 5 -ts_max_steps 4 -convest_num_refine 3 gives L_2 convergence rate: [0.26, -0.0058, 0.26]
2306: suffix: 2d_mandel_tconv
2307: args: -ts_dt 0.023 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1
2309: testset:
2310: requires: ctetgen !complex
2311: args: -sol_type cryer -dm_plex_dim 3 -dm_plex_shape ball \
2312: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1
2314: test:
2315: suffix: 3d_cryer
2316: args: -ts_dt 0.0028666667 -ts_max_time 0.014333 -ts_max_steps 2 -dmts_check .0001 \
2317: -pc_type svd
2319: test:
2320: # -bd_dm_refine 3 -dm_refine_volume_limit_pre 0.004 -convest_num_refine 2 gives L_2 convergence rate: []
2321: suffix: 3d_cryer_sconv
2322: args: -bd_dm_refine 1 -dm_refine_volume_limit_pre 0.00666667 \
2323: -ts_dt 1e-5 -dt_initial 1e-5 -ts_max_steps 2 \
2324: -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
2325: -pc_type lu -pc_factor_shift_type nonzero
2327: test:
2328: # Displacement and Pressure converge. The analytic expression for trace strain is inaccurate at the origin
2329: # -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]
2330: suffix: 3d_cryer_tconv
2331: args: -bd_dm_refine 1 -dm_refine_volume_limit_pre 0.00666667 \
2332: -ts_dt 0.023 -ts_max_time 0.092 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1 \
2333: -pc_type lu -pc_factor_shift_type nonzero
2335: TEST*/