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, PetscCtx ctx)
72: {
73: for (PetscInt c = 0; c < Nc; ++c) u[c] = 0.0;
74: return PETSC_SUCCESS;
75: }
77: /* Quadratic space and linear time solution
79: 2D:
80: u = x^2
81: v = y^2 - 2xy
82: p = (x + y) t
83: e = 2y
84: f = <2 G, 4 G + 2 \lambda > - <alpha t, alpha t>
85: g = 0
86: \epsilon = / 2x -y \
87: \ -y 2y - 2x /
88: Tr(\epsilon) = e = div u = 2y
89: div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
90: = 2 G < 2-1, 2 > + \lambda <0, 2> - alpha <t, t>
91: = <2 G, 4 G + 2 \lambda> - <alpha t, alpha t>
92: \frac{1}{M} \frac{dp}{dt} + \alpha \frac{d\varepsilon}{dt} - \nabla \cdot \kappa \nabla p
93: = \frac{1}{M} \frac{dp}{dt} + \kappa \Delta p
94: = (x + y)/M
96: 3D:
97: u = x^2
98: v = y^2 - 2xy
99: w = z^2 - 2yz
100: p = (x + y + z) t
101: e = 2z
102: f = <2 G, 4 G + 2 \lambda > - <alpha t, alpha t, alpha t>
103: g = 0
104: \varepsilon = / 2x -y 0 \
105: | -y 2y - 2x -z |
106: \ 0 -z 2z - 2y/
107: Tr(\varepsilon) = div u = 2z
108: div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
109: = 2 G < 2-1, 2-1, 2 > + \lambda <0, 0, 2> - alpha <t, t, t>
110: = <2 G, 2G, 4 G + 2 \lambda> - <alpha t, alpha t, alpha t>
111: */
112: static PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
113: {
114: PetscInt d;
116: for (d = 0; d < dim; ++d) u[d] = PetscSqr(x[d]) - (d > 0 ? 2.0 * x[d - 1] * x[d] : 0.0);
117: return PETSC_SUCCESS;
118: }
120: static PetscErrorCode linear_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
121: {
122: u[0] = 2.0 * x[dim - 1];
123: return PETSC_SUCCESS;
124: }
126: static PetscErrorCode linear_linear_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
127: {
128: PetscReal sum = 0.0;
129: PetscInt d;
131: for (d = 0; d < dim; ++d) sum += x[d];
132: u[0] = sum * time;
133: return PETSC_SUCCESS;
134: }
136: static PetscErrorCode linear_linear_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
137: {
138: PetscReal sum = 0.0;
139: PetscInt d;
141: for (d = 0; d < dim; ++d) sum += x[d];
142: u[0] = sum;
143: return PETSC_SUCCESS;
144: }
146: 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[])
147: {
148: const PetscReal G = PetscRealPart(constants[0]);
149: const PetscReal K_u = PetscRealPart(constants[1]);
150: const PetscReal alpha = PetscRealPart(constants[2]);
151: const PetscReal M = PetscRealPart(constants[3]);
152: const PetscReal K_d = K_u - alpha * alpha * M;
153: const PetscReal lambda = K_d - (2.0 * G) / 3.0;
154: PetscInt d;
156: for (d = 0; d < dim - 1; ++d) f0[d] -= 2.0 * G - alpha * t;
157: f0[dim - 1] -= 2.0 * lambda + 4.0 * G - alpha * t;
158: }
160: 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[])
161: {
162: const PetscReal alpha = PetscRealPart(constants[2]);
163: const PetscReal M = PetscRealPart(constants[3]);
164: PetscReal sum = 0.0;
165: PetscInt d;
167: for (d = 0; d < dim; ++d) sum += x[d];
168: f0[0] += u_t ? alpha * u_t[uOff[1]] : 0.0;
169: f0[0] += u_t ? u_t[uOff[2]] / M : 0.0;
170: f0[0] -= sum / M;
171: }
173: /* Quadratic space and trigonometric time solution
175: 2D:
176: u = x^2
177: v = y^2 - 2xy
178: p = (x + y) cos(t)
179: e = 2y
180: f = <2 G, 4 G + 2 \lambda > - <alpha cos(t), alpha cos(t)>
181: g = 0
182: \epsilon = / 2x -y \
183: \ -y 2y - 2x /
184: Tr(\epsilon) = e = div u = 2y
185: div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
186: = 2 G < 2-1, 2 > + \lambda <0, 2> - alpha <cos(t), cos(t)>
187: = <2 G, 4 G + 2 \lambda> - <alpha cos(t), alpha cos(t)>
188: \frac{1}{M} \frac{dp}{dt} + \alpha \frac{d\varepsilon}{dt} - \nabla \cdot \kappa \nabla p
189: = \frac{1}{M} \frac{dp}{dt} + \kappa \Delta p
190: = -(x + y)/M sin(t)
192: 3D:
193: u = x^2
194: v = y^2 - 2xy
195: w = z^2 - 2yz
196: p = (x + y + z) cos(t)
197: e = 2z
198: f = <2 G, 4 G + 2 \lambda > - <alpha cos(t), alpha cos(t), alpha cos(t)>
199: g = 0
200: \varepsilon = / 2x -y 0 \
201: | -y 2y - 2x -z |
202: \ 0 -z 2z - 2y/
203: Tr(\varepsilon) = div u = 2z
204: div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
205: = 2 G < 2-1, 2-1, 2 > + \lambda <0, 0, 2> - alpha <cos(t), cos(t), cos(t)>
206: = <2 G, 2G, 4 G + 2 \lambda> - <alpha cos(t), alpha cos(t), alpha cos(t)>
207: */
208: static PetscErrorCode linear_trig_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
209: {
210: PetscReal sum = 0.0;
211: PetscInt d;
213: for (d = 0; d < dim; ++d) sum += x[d];
214: u[0] = sum * PetscCosReal(time);
215: return PETSC_SUCCESS;
216: }
218: static PetscErrorCode linear_trig_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
219: {
220: PetscReal sum = 0.0;
221: PetscInt d;
223: for (d = 0; d < dim; ++d) sum += x[d];
224: u[0] = -sum * PetscSinReal(time);
225: return PETSC_SUCCESS;
226: }
228: 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[])
229: {
230: const PetscReal G = PetscRealPart(constants[0]);
231: const PetscReal K_u = PetscRealPart(constants[1]);
232: const PetscReal alpha = PetscRealPart(constants[2]);
233: const PetscReal M = PetscRealPart(constants[3]);
234: const PetscReal K_d = K_u - alpha * alpha * M;
235: const PetscReal lambda = K_d - (2.0 * G) / 3.0;
236: PetscInt d;
238: for (d = 0; d < dim - 1; ++d) f0[d] -= 2.0 * G - alpha * PetscCosReal(t);
239: f0[dim - 1] -= 2.0 * lambda + 4.0 * G - alpha * PetscCosReal(t);
240: }
242: 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[])
243: {
244: const PetscReal alpha = PetscRealPart(constants[2]);
245: const PetscReal M = PetscRealPart(constants[3]);
246: PetscReal sum = 0.0;
247: PetscInt d;
249: for (d = 0; d < dim; ++d) sum += x[d];
251: f0[0] += u_t ? alpha * u_t[uOff[1]] : 0.0;
252: f0[0] += u_t ? u_t[uOff[2]] / M : 0.0;
253: f0[0] += PetscSinReal(t) * sum / M;
254: }
256: /* Trigonometric space and linear time solution
258: u = sin(2 pi x)
259: v = sin(2 pi y) - 2xy
260: \varepsilon = / 2 pi cos(2 pi x) -y \
261: \ -y 2 pi cos(2 pi y) - 2x /
262: Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
263: div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
264: = \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) >
265: = \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) >
267: 2D:
268: u = sin(2 pi x)
269: v = sin(2 pi y) - 2xy
270: p = (cos(2 pi x) + cos(2 pi y)) t
271: e = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
272: 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)>
273: g = 0
274: \varepsilon = / 2 pi cos(2 pi x) -y \
275: \ -y 2 pi cos(2 pi y) - 2x /
276: Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
277: div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
278: = 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>
279: = < -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)>
280: \frac{1}{M} \frac{dp}{dt} + \alpha \frac{d\varepsilon}{dt} - \nabla \cdot \kappa \nabla p
281: = \frac{1}{M} \frac{dp}{dt} + \kappa \Delta p
282: = (cos(2 pi x) + cos(2 pi y))/M - 4 pi^2 \kappa (cos(2 pi x) + cos(2 pi y)) t
284: 3D:
285: u = sin(2 pi x)
286: v = sin(2 pi y) - 2xy
287: v = sin(2 pi y) - 2yz
288: p = (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) t
289: e = 2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2y
290: 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)>
291: g = 0
292: \varepsilon = / 2 pi cos(2 pi x) -y 0 \
293: | -y 2 pi cos(2 pi y) - 2x -z |
294: \ 0 -z 2 pi cos(2 pi z) - 2y /
295: Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2 y
296: div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
297: = 2 G < -4 pi^2 sin(2 pi x) - 1, -4 pi^2 sin(2 pi y) - 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>
298: = < -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)>
299: \frac{1}{M} \frac{dp}{dt} + \alpha \frac{d\varepsilon}{dt} - \nabla \cdot \kappa \nabla p
300: = \frac{1}{M} \frac{dp}{dt} + \kappa \Delta p
301: = (cos(2 pi x) + cos(2 pi y) + cos(2 pi z))/M - 4 pi^2 \kappa (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) t
302: */
303: static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
304: {
305: for (PetscInt d = 0; d < dim; ++d) u[d] = PetscSinReal(2. * PETSC_PI * x[d]) - (d > 0 ? 2.0 * x[d - 1] * x[d] : 0.0);
306: return PETSC_SUCCESS;
307: }
309: static PetscErrorCode trig_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
310: {
311: PetscReal sum = 0.0;
313: for (PetscInt d = 0; d < dim; ++d) sum += 2. * PETSC_PI * PetscCosReal(2. * PETSC_PI * x[d]) - (d < dim - 1 ? 2. * x[d] : 0.0);
314: u[0] = sum;
315: return PETSC_SUCCESS;
316: }
318: static PetscErrorCode trig_linear_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
319: {
320: PetscReal sum = 0.0;
322: for (PetscInt d = 0; d < dim; ++d) sum += PetscCosReal(2. * PETSC_PI * x[d]);
323: u[0] = sum * time;
324: return PETSC_SUCCESS;
325: }
327: static PetscErrorCode trig_linear_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
328: {
329: PetscReal sum = 0.0;
331: for (PetscInt d = 0; d < dim; ++d) sum += PetscCosReal(2. * PETSC_PI * x[d]);
332: u[0] = sum;
333: return PETSC_SUCCESS;
334: }
336: 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[])
337: {
338: const PetscReal G = PetscRealPart(constants[0]);
339: const PetscReal K_u = PetscRealPart(constants[1]);
340: const PetscReal alpha = PetscRealPart(constants[2]);
341: const PetscReal M = PetscRealPart(constants[3]);
342: const PetscReal K_d = K_u - alpha * alpha * M;
343: const PetscReal lambda = K_d - (2.0 * G) / 3.0;
345: for (PetscInt 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;
346: 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;
347: }
349: 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[])
350: {
351: const PetscReal alpha = PetscRealPart(constants[2]);
352: const PetscReal M = PetscRealPart(constants[3]);
353: const PetscReal kappa = PetscRealPart(constants[4]);
354: PetscReal sum = 0.0;
356: for (PetscInt d = 0; d < dim; ++d) sum += PetscCosReal(2. * PETSC_PI * x[d]);
357: f0[0] += u_t ? alpha * u_t[uOff[1]] : 0.0;
358: f0[0] += u_t ? u_t[uOff[2]] / M : 0.0;
359: f0[0] -= sum / M - 4 * PetscSqr(PETSC_PI) * kappa * sum * t;
360: }
362: /* Terzaghi Solutions */
363: /* The analytical solutions given here are drawn from chapter 7, section 3, */
364: /* "One-Dimensional Consolidation Problem," from Poroelasticity, by Cheng. */
365: static PetscErrorCode terzaghi_drainage_pressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
366: {
367: AppCtx *user = (AppCtx *)ctx;
368: Parameter *param;
370: PetscCall(PetscBagGetData(user->bag, ¶m));
371: if (time <= 0.0) {
372: PetscScalar alpha = param->alpha; /* - */
373: PetscScalar K_u = param->K_u; /* Pa */
374: PetscScalar M = param->M; /* Pa */
375: PetscScalar G = param->mu; /* Pa */
376: PetscScalar P_0 = param->P_0; /* Pa */
377: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
378: PetscScalar eta = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G); /* -, Cheng (B.11) */
379: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
381: u[0] = ((P_0 * eta) / (G * S));
382: } else {
383: u[0] = 0.0;
384: }
385: return PETSC_SUCCESS;
386: }
388: static PetscErrorCode terzaghi_initial_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
389: {
390: AppCtx *user = (AppCtx *)ctx;
391: Parameter *param;
393: PetscCall(PetscBagGetData(user->bag, ¶m));
394: {
395: PetscScalar K_u = param->K_u; /* Pa */
396: PetscScalar G = param->mu; /* Pa */
397: PetscScalar P_0 = param->P_0; /* Pa */
398: PetscReal L = user->xmax[1] - user->xmin[1]; /* m */
399: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
400: PetscReal zstar = x[1] / L; /* - */
402: u[0] = 0.0;
403: u[1] = ((P_0 * L * (1.0 - 2.0 * nu_u)) / (2.0 * G * (1.0 - nu_u))) * (1.0 - zstar);
404: }
405: return PETSC_SUCCESS;
406: }
408: static PetscErrorCode terzaghi_initial_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
409: {
410: AppCtx *user = (AppCtx *)ctx;
411: Parameter *param;
413: PetscCall(PetscBagGetData(user->bag, ¶m));
414: {
415: PetscScalar K_u = param->K_u; /* Pa */
416: PetscScalar G = param->mu; /* Pa */
417: PetscScalar P_0 = param->P_0; /* Pa */
418: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
420: u[0] = -(P_0 * (1.0 - 2.0 * nu_u)) / (2.0 * G * (1.0 - nu_u));
421: }
422: return PETSC_SUCCESS;
423: }
425: static PetscErrorCode terzaghi_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
426: {
427: AppCtx *user = (AppCtx *)ctx;
428: Parameter *param;
430: PetscCall(PetscBagGetData(user->bag, ¶m));
431: if (time < 0.0) {
432: PetscCall(terzaghi_initial_u(dim, time, x, Nc, u, ctx));
433: } else {
434: PetscScalar alpha = param->alpha; /* - */
435: PetscScalar K_u = param->K_u; /* Pa */
436: PetscScalar M = param->M; /* Pa */
437: PetscScalar G = param->mu; /* Pa */
438: PetscScalar P_0 = param->P_0; /* Pa */
439: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
440: PetscReal L = user->xmax[1] - user->xmin[1]; /* m */
441: PetscInt N = user->niter, m;
443: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
444: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
445: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
446: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
447: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
449: PetscReal zstar = x[1] / L; /* - */
450: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
451: PetscScalar F2 = 0.0;
453: for (m = 1; m < 2 * N + 1; ++m) {
454: 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));
455: }
456: u[0] = 0.0;
457: 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 */
458: }
459: return PETSC_SUCCESS;
460: }
462: static PetscErrorCode terzaghi_2d_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
463: {
464: AppCtx *user = (AppCtx *)ctx;
465: Parameter *param;
467: PetscCall(PetscBagGetData(user->bag, ¶m));
468: if (time < 0.0) {
469: PetscCall(terzaghi_initial_eps(dim, time, x, Nc, u, ctx));
470: } else {
471: PetscScalar alpha = param->alpha; /* - */
472: PetscScalar K_u = param->K_u; /* Pa */
473: PetscScalar M = param->M; /* Pa */
474: PetscScalar G = param->mu; /* Pa */
475: PetscScalar P_0 = param->P_0; /* Pa */
476: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
477: PetscReal L = user->xmax[1] - user->xmin[1]; /* m */
478: PetscInt N = user->niter, m;
480: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
481: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
482: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
483: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
484: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
486: PetscReal zstar = x[1] / L; /* - */
487: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
488: PetscScalar F2_z = 0.0;
490: for (m = 1; m < 2 * N + 1; ++m) {
491: 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));
492: }
493: 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; /* - */
494: }
495: return PETSC_SUCCESS;
496: }
498: // Pressure
499: static PetscErrorCode terzaghi_2d_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
500: {
501: AppCtx *user = (AppCtx *)ctx;
502: Parameter *param;
504: PetscCall(PetscBagGetData(user->bag, ¶m));
505: if (time <= 0.0) {
506: PetscCall(terzaghi_drainage_pressure(dim, time, x, Nc, u, ctx));
507: } else {
508: PetscScalar alpha = param->alpha; /* - */
509: PetscScalar K_u = param->K_u; /* Pa */
510: PetscScalar M = param->M; /* Pa */
511: PetscScalar G = param->mu; /* Pa */
512: PetscScalar P_0 = param->P_0; /* Pa */
513: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
514: PetscReal L = user->xmax[1] - user->xmin[1]; /* m */
515: PetscInt N = user->niter, m;
517: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
518: PetscScalar eta = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G); /* -, Cheng (B.11) */
519: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
520: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
522: PetscReal zstar = x[1] / L; /* - */
523: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
524: PetscScalar F1 = 0.0;
526: 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));
528: for (m = 1; m < 2 * N + 1; ++m) {
529: if (m % 2 == 1) F1 += (4.0 / (m * PETSC_PI)) * PetscSinReal(0.5 * m * PETSC_PI * zstar) * PetscExpReal(-PetscSqr(m * PETSC_PI) * tstar);
530: }
531: u[0] = ((P_0 * eta) / (G * S)) * F1; /* Pa */
532: }
533: return PETSC_SUCCESS;
534: }
536: static PetscErrorCode terzaghi_2d_u_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
537: {
538: AppCtx *user = (AppCtx *)ctx;
539: Parameter *param;
541: PetscCall(PetscBagGetData(user->bag, ¶m));
542: if (time <= 0.0) {
543: u[0] = 0.0;
544: u[1] = 0.0;
545: } else {
546: PetscScalar alpha = param->alpha; /* - */
547: PetscScalar K_u = param->K_u; /* Pa */
548: PetscScalar M = param->M; /* Pa */
549: PetscScalar G = param->mu; /* Pa */
550: PetscScalar P_0 = param->P_0; /* Pa */
551: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
552: PetscReal L = user->xmax[1] - user->xmin[1]; /* m */
553: PetscInt N = user->niter, m;
555: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
556: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
557: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
558: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
559: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
561: PetscReal zstar = x[1] / L; /* - */
562: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
563: PetscScalar F2_t = 0.0;
565: for (m = 1; m < 2 * N + 1; ++m) {
566: if (m % 2 == 1) F2_t += (2.0 * c / PetscSqr(L)) * PetscCosReal(0.5 * m * PETSC_PI * zstar) * PetscExpReal(-PetscSqr(m * PETSC_PI) * tstar);
567: }
568: u[0] = 0.0;
569: u[1] = ((P_0 * L * (nu_u - nu)) / (2.0 * G * (1.0 - nu_u) * (1.0 - nu))) * F2_t; /* m / s */
570: }
571: return PETSC_SUCCESS;
572: }
574: static PetscErrorCode terzaghi_2d_eps_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
575: {
576: AppCtx *user = (AppCtx *)ctx;
577: Parameter *param;
579: PetscCall(PetscBagGetData(user->bag, ¶m));
580: if (time <= 0.0) {
581: u[0] = 0.0;
582: } else {
583: PetscScalar alpha = param->alpha; /* - */
584: PetscScalar K_u = param->K_u; /* Pa */
585: PetscScalar M = param->M; /* Pa */
586: PetscScalar G = param->mu; /* Pa */
587: PetscScalar P_0 = param->P_0; /* Pa */
588: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
589: PetscReal L = user->xmax[1] - user->xmin[1]; /* m */
590: PetscInt N = user->niter, m;
592: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
593: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
594: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
595: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
596: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
598: PetscReal zstar = x[1] / L; /* - */
599: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
600: PetscScalar F2_zt = 0.0;
602: for (m = 1; m < 2 * N + 1; ++m) {
603: 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);
604: }
605: u[0] = ((P_0 * L * (nu_u - nu)) / (2.0 * G * (1.0 - nu_u) * (1.0 - nu))) * F2_zt; /* 1 / s */
606: }
607: return PETSC_SUCCESS;
608: }
610: static PetscErrorCode terzaghi_2d_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
611: {
612: AppCtx *user = (AppCtx *)ctx;
613: Parameter *param;
615: PetscCall(PetscBagGetData(user->bag, ¶m));
616: if (time <= 0.0) {
617: PetscScalar alpha = param->alpha; /* - */
618: PetscScalar K_u = param->K_u; /* Pa */
619: PetscScalar M = param->M; /* Pa */
620: PetscScalar G = param->mu; /* Pa */
621: PetscScalar P_0 = param->P_0; /* Pa */
622: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
623: PetscReal L = user->xmax[1] - user->xmin[1]; /* m */
625: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
626: PetscScalar eta = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G); /* -, Cheng (B.11) */
627: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
628: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
630: u[0] = -((P_0 * eta) / (G * S)) * PetscSqr(0 * PETSC_PI) * c / PetscSqr(2.0 * L); /* Pa / s */
631: } else {
632: PetscScalar alpha = param->alpha; /* - */
633: PetscScalar K_u = param->K_u; /* Pa */
634: PetscScalar M = param->M; /* Pa */
635: PetscScalar G = param->mu; /* Pa */
636: PetscScalar P_0 = param->P_0; /* Pa */
637: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
638: PetscReal L = user->xmax[1] - user->xmin[1]; /* m */
639: PetscInt N = user->niter, m;
641: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
642: PetscScalar eta = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G); /* -, Cheng (B.11) */
643: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
644: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
646: PetscReal zstar = x[1] / L; /* - */
647: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
648: PetscScalar F1_t = 0.0;
650: 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));
652: for (m = 1; m < 2 * N + 1; ++m) {
653: 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);
654: }
655: u[0] = ((P_0 * eta) / (G * S)) * F1_t; /* Pa / s */
656: }
657: return PETSC_SUCCESS;
658: }
660: /* Mandel Solutions */
661: static PetscErrorCode mandel_drainage_pressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
662: {
663: AppCtx *user = (AppCtx *)ctx;
664: Parameter *param;
666: PetscCall(PetscBagGetData(user->bag, ¶m));
667: if (time <= 0.0) {
668: PetscScalar alpha = param->alpha; /* - */
669: PetscScalar K_u = param->K_u; /* Pa */
670: PetscScalar M = param->M; /* Pa */
671: PetscScalar G = param->mu; /* Pa */
672: PetscScalar P_0 = param->P_0; /* Pa */
673: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
674: PetscReal a = 0.5 * (user->xmax[0] - user->xmin[0]); /* m */
675: PetscInt N = user->niter, n;
677: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
678: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
679: PetscScalar B = alpha * M / K_u; /* -, Cheng (B.12) */
680: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
681: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
683: PetscScalar A1 = 3.0 / (B * (1.0 + nu_u));
684: PetscReal aa = 0.0;
685: PetscReal p = 0.0;
686: PetscReal time = 0.0;
688: for (n = 1; n < N + 1; ++n) {
689: aa = user->zeroArray[n - 1];
690: 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));
691: }
692: u[0] = ((2.0 * P_0) / (a * A1)) * p;
693: } else {
694: u[0] = 0.0;
695: }
696: return PETSC_SUCCESS;
697: }
699: static PetscErrorCode mandel_initial_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
700: {
701: AppCtx *user = (AppCtx *)ctx;
702: Parameter *param;
704: PetscCall(PetscBagGetData(user->bag, ¶m));
705: {
706: PetscScalar alpha = param->alpha; /* - */
707: PetscScalar K_u = param->K_u; /* Pa */
708: PetscScalar M = param->M; /* Pa */
709: PetscScalar G = param->mu; /* Pa */
710: PetscScalar P_0 = param->P_0; /* Pa */
711: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
712: PetscReal a = 0.5 * (user->xmax[0] - user->xmin[0]); /* m */
713: PetscInt N = user->niter, n;
715: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
716: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
717: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
718: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
719: PetscReal c = PetscRealPart(kappa / S); /* m^2 / s, Cheng (B.16) */
721: PetscReal A_s = 0.0;
722: PetscReal B_s = 0.0;
723: for (n = 1; n < N + 1; ++n) {
724: PetscReal alpha_n = user->zeroArray[n - 1];
725: 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));
726: 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));
727: }
728: u[0] = ((P_0 * nu) / (2.0 * G * a) - (P_0 * nu_u) / (G * a) * A_s) * x[0] + P_0 / G * B_s;
729: u[1] = (-1 * (P_0 * (1.0 - nu)) / (2 * G * a) + (P_0 * (1 - nu_u)) / (G * a) * A_s) * x[1];
730: }
731: return PETSC_SUCCESS;
732: }
734: static PetscErrorCode mandel_initial_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
735: {
736: AppCtx *user = (AppCtx *)ctx;
737: Parameter *param;
739: PetscCall(PetscBagGetData(user->bag, ¶m));
740: {
741: PetscScalar alpha = param->alpha; /* - */
742: PetscScalar K_u = param->K_u; /* Pa */
743: PetscScalar M = param->M; /* Pa */
744: PetscScalar G = param->mu; /* Pa */
745: PetscScalar P_0 = param->P_0; /* Pa */
746: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
747: PetscReal a = 0.5 * (user->xmax[0] - user->xmin[0]); /* m */
748: PetscInt N = user->niter, n;
750: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
751: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
752: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
753: PetscReal c = PetscRealPart(kappa / S); /* m^2 / s, Cheng (B.16) */
755: PetscReal aa = 0.0;
756: PetscReal eps_A = 0.0;
757: PetscReal eps_B = 0.0;
758: PetscReal eps_C = 0.0;
759: PetscReal time = 0.0;
761: for (n = 1; n < N + 1; ++n) {
762: aa = user->zeroArray[n - 1];
763: 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)));
764: eps_B += (PetscExpReal((-1.0 * aa * aa * c * time) / (a * a)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa));
765: eps_C += (PetscExpReal((-1.0 * aa * aa * c * time) / (aa * aa)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa));
766: }
767: 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);
768: }
769: return PETSC_SUCCESS;
770: }
772: // Displacement
773: static PetscErrorCode mandel_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
774: {
775: Parameter *param;
777: AppCtx *user = (AppCtx *)ctx;
779: PetscCall(PetscBagGetData(user->bag, ¶m));
780: if (time <= 0.0) PetscCall(mandel_initial_u(dim, time, x, Nc, u, ctx));
781: else {
782: PetscInt NITER = user->niter;
783: PetscScalar alpha = param->alpha;
784: PetscScalar K_u = param->K_u;
785: PetscScalar M = param->M;
786: PetscScalar G = param->mu;
787: PetscScalar k = param->k;
788: PetscScalar mu_f = param->mu_f;
789: PetscScalar F = param->P_0;
791: PetscScalar K_d = K_u - alpha * alpha * M;
792: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
793: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
794: PetscScalar kappa = k / mu_f;
795: PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0;
796: PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u)));
798: // Series term
799: PetscScalar A_x = 0.0;
800: PetscScalar B_x = 0.0;
802: for (PetscInt n = 1; n < NITER + 1; n++) {
803: PetscReal alpha_n = user->zeroArray[n - 1];
805: 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));
806: 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));
807: }
808: u[0] = ((F * nu) / (2.0 * G * a) - (F * nu_u) / (G * a) * A_x) * x[0] + F / G * B_x;
809: u[1] = (-1 * (F * (1.0 - nu)) / (2 * G * a) + (F * (1 - nu_u)) / (G * a) * A_x) * x[1];
810: }
811: return PETSC_SUCCESS;
812: }
814: // Trace strain
815: static PetscErrorCode mandel_2d_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
816: {
817: Parameter *param;
819: AppCtx *user = (AppCtx *)ctx;
821: PetscCall(PetscBagGetData(user->bag, ¶m));
822: if (time <= 0.0) PetscCall(mandel_initial_eps(dim, time, x, Nc, u, ctx));
823: else {
824: PetscInt NITER = user->niter;
825: PetscScalar alpha = param->alpha;
826: PetscScalar K_u = param->K_u;
827: PetscScalar M = param->M;
828: PetscScalar G = param->mu;
829: PetscScalar k = param->k;
830: PetscScalar mu_f = param->mu_f;
831: PetscScalar F = param->P_0;
833: PetscScalar K_d = K_u - alpha * alpha * M;
834: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
835: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
836: PetscScalar kappa = k / mu_f;
837: //const PetscScalar B = (alpha*M)/(K_d + alpha*alpha * M);
839: //const PetscScalar b = (YMAX - YMIN) / 2.0;
840: PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0;
841: PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u)));
843: // Series term
844: PetscReal eps_A = 0.0;
845: PetscReal eps_B = 0.0;
846: PetscReal eps_C = 0.0;
848: for (PetscInt n = 1; n < NITER + 1; n++) {
849: PetscReal aa = user->zeroArray[n - 1];
851: 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)));
853: eps_B += (PetscExpReal((-1.0 * aa * aa * c * time) / (a * a)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa));
855: eps_C += (PetscExpReal((-1.0 * aa * aa * c * time) / (aa * aa)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa));
856: }
858: 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);
859: }
860: return PETSC_SUCCESS;
861: }
863: // Pressure
864: static PetscErrorCode mandel_2d_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
865: {
866: Parameter *param;
868: AppCtx *user = (AppCtx *)ctx;
870: PetscCall(PetscBagGetData(user->bag, ¶m));
871: if (time <= 0.0) {
872: PetscCall(mandel_drainage_pressure(dim, time, x, Nc, u, ctx));
873: } else {
874: PetscInt NITER = user->niter;
876: PetscScalar alpha = param->alpha;
877: PetscScalar K_u = param->K_u;
878: PetscScalar M = param->M;
879: PetscScalar G = param->mu;
880: PetscScalar k = param->k;
881: PetscScalar mu_f = param->mu_f;
882: PetscScalar F = param->P_0;
884: PetscScalar K_d = K_u - alpha * alpha * M;
885: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
886: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
887: PetscScalar kappa = k / mu_f;
888: PetscScalar B = (alpha * M) / (K_d + alpha * alpha * M);
890: PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0;
891: PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u)));
892: PetscScalar A1 = 3.0 / (B * (1.0 + nu_u));
893: //PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu);
895: // Series term
896: PetscReal p = 0.0;
898: for (PetscInt n = 1; n < NITER + 1; n++) {
899: PetscReal aa = user->zeroArray[n - 1];
900: p += (PetscSinReal(aa) / (aa - PetscSinReal(aa) * PetscCosReal(aa))) * (PetscCosReal((aa * x[0]) / a) - PetscCosReal(aa)) * PetscExpReal(-1.0 * (aa * aa * c * time) / (a * a));
901: }
902: u[0] = ((2.0 * F) / (a * A1)) * p;
903: }
904: return PETSC_SUCCESS;
905: }
907: // Time derivative of displacement
908: static PetscErrorCode mandel_2d_u_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
909: {
910: Parameter *param;
912: AppCtx *user = (AppCtx *)ctx;
914: PetscCall(PetscBagGetData(user->bag, ¶m));
916: PetscInt NITER = user->niter;
917: PetscScalar alpha = param->alpha;
918: PetscScalar K_u = param->K_u;
919: PetscScalar M = param->M;
920: PetscScalar G = param->mu;
921: PetscScalar F = param->P_0;
923: PetscScalar K_d = K_u - alpha * alpha * M;
924: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
925: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
926: PetscScalar kappa = param->k / param->mu_f;
927: PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0;
928: PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u)));
930: // Series term
931: PetscScalar A_s_t = 0.0;
932: PetscScalar B_s_t = 0.0;
934: for (PetscInt n = 1; n < NITER + 1; n++) {
935: PetscReal alpha_n = user->zeroArray[n - 1];
937: 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)));
938: 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)));
939: }
941: u[0] = (F / G) * A_s_t - ((F * nu_u * x[0]) / (G * a)) * B_s_t;
942: u[1] = ((F * x[1] * (1 - nu_u)) / (G * a)) * B_s_t;
944: return PETSC_SUCCESS;
945: }
947: // Time derivative of trace strain
948: static PetscErrorCode mandel_2d_eps_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
949: {
950: Parameter *param;
952: AppCtx *user = (AppCtx *)ctx;
954: PetscCall(PetscBagGetData(user->bag, ¶m));
956: PetscInt NITER = user->niter;
957: PetscScalar alpha = param->alpha;
958: PetscScalar K_u = param->K_u;
959: PetscScalar M = param->M;
960: PetscScalar G = param->mu;
961: PetscScalar k = param->k;
962: PetscScalar mu_f = param->mu_f;
963: PetscScalar F = param->P_0;
965: PetscScalar K_d = K_u - alpha * alpha * M;
966: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
967: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
968: PetscScalar kappa = k / mu_f;
969: //const PetscScalar B = (alpha*M)/(K_d + alpha*alpha * M);
971: //const PetscScalar b = (YMAX - YMIN) / 2.0;
972: PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0;
973: PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u)));
975: // Series term
976: PetscScalar eps_As = 0.0;
977: PetscScalar eps_Bs = 0.0;
978: PetscScalar eps_Cs = 0.0;
980: for (PetscInt n = 1; n < NITER + 1; n++) {
981: PetscReal alpha_n = user->zeroArray[n - 1];
983: 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)));
984: 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)));
985: 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)));
986: }
988: u[0] = (F / G) * eps_As - ((F * nu_u) / (G * a)) * eps_Bs + ((F * (1 - nu_u)) / (G * a)) * eps_Cs;
989: return PETSC_SUCCESS;
990: }
992: // Time derivative of pressure
993: static PetscErrorCode mandel_2d_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
994: {
995: Parameter *param;
997: AppCtx *user = (AppCtx *)ctx;
999: PetscCall(PetscBagGetData(user->bag, ¶m));
1001: PetscScalar alpha = param->alpha;
1002: PetscScalar K_u = param->K_u;
1003: PetscScalar M = param->M;
1004: PetscScalar G = param->mu;
1005: PetscScalar F = param->P_0;
1007: PetscScalar K_d = K_u - alpha * alpha * M;
1008: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
1009: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
1011: PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0;
1012: //PetscScalar A1 = 3.0 / (B * (1.0 + nu_u));
1013: //PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu);
1015: u[0] = ((2.0 * F * (-2.0 * nu + 3.0 * nu_u)) / (3.0 * a * alpha * (1.0 - 2.0 * nu)));
1017: return PETSC_SUCCESS;
1018: }
1020: /* Cryer Solutions */
1021: static PetscErrorCode cryer_drainage_pressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
1022: {
1023: AppCtx *user = (AppCtx *)ctx;
1024: Parameter *param;
1026: PetscCall(PetscBagGetData(user->bag, ¶m));
1027: if (time <= 0.0) {
1028: PetscScalar alpha = param->alpha; /* - */
1029: PetscScalar K_u = param->K_u; /* Pa */
1030: PetscScalar M = param->M; /* Pa */
1031: PetscScalar P_0 = param->P_0; /* Pa */
1032: PetscScalar B = alpha * M / K_u; /* -, Cheng (B.12) */
1034: u[0] = P_0 * B;
1035: } else {
1036: u[0] = 0.0;
1037: }
1038: return PETSC_SUCCESS;
1039: }
1041: static PetscErrorCode cryer_initial_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
1042: {
1043: AppCtx *user = (AppCtx *)ctx;
1044: Parameter *param;
1046: PetscCall(PetscBagGetData(user->bag, ¶m));
1047: {
1048: PetscScalar K_u = param->K_u; /* Pa */
1049: PetscScalar G = param->mu; /* Pa */
1050: PetscScalar P_0 = param->P_0; /* Pa */
1051: PetscReal R_0 = user->xmax[1]; /* m */
1052: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
1054: PetscScalar u_0 = -P_0 * R_0 * (1. - 2. * nu_u) / (2. * G * (1. + nu_u)); /* Cheng (7.407) */
1055: PetscReal u_sc = PetscRealPart(u_0) / R_0;
1057: u[0] = u_sc * x[0];
1058: u[1] = u_sc * x[1];
1059: u[2] = u_sc * x[2];
1060: }
1061: return PETSC_SUCCESS;
1062: }
1064: static PetscErrorCode cryer_initial_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
1065: {
1066: AppCtx *user = (AppCtx *)ctx;
1067: Parameter *param;
1069: PetscCall(PetscBagGetData(user->bag, ¶m));
1070: {
1071: PetscScalar K_u = param->K_u; /* Pa */
1072: PetscScalar G = param->mu; /* Pa */
1073: PetscScalar P_0 = param->P_0; /* Pa */
1074: PetscReal R_0 = user->xmax[1]; /* m */
1075: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
1077: PetscScalar u_0 = -P_0 * R_0 * (1. - 2. * nu_u) / (2. * G * (1. + nu_u)); /* Cheng (7.407) */
1078: PetscReal u_sc = PetscRealPart(u_0) / R_0;
1080: /* div R = 1/R^2 d/dR R^2 R = 3 */
1081: u[0] = 3. * u_sc;
1082: u[1] = 3. * u_sc;
1083: u[2] = 3. * u_sc;
1084: }
1085: return PETSC_SUCCESS;
1086: }
1088: // Displacement
1089: static PetscErrorCode cryer_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
1090: {
1091: AppCtx *user = (AppCtx *)ctx;
1092: Parameter *param;
1094: PetscCall(PetscBagGetData(user->bag, ¶m));
1095: if (time <= 0.0) {
1096: PetscCall(cryer_initial_u(dim, time, x, Nc, u, ctx));
1097: } else {
1098: PetscScalar alpha = param->alpha; /* - */
1099: PetscScalar K_u = param->K_u; /* Pa */
1100: PetscScalar M = param->M; /* Pa */
1101: PetscScalar G = param->mu; /* Pa */
1102: PetscScalar P_0 = param->P_0; /* Pa */
1103: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
1104: PetscReal R_0 = user->xmax[1]; /* m */
1105: PetscInt N = user->niter, n;
1107: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
1108: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
1109: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
1110: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
1111: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
1112: PetscScalar u_inf = -P_0 * R_0 * (1. - 2. * nu) / (2. * G * (1. + nu)); /* m, Cheng (7.388) */
1114: PetscReal R = PetscSqrtReal(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
1115: PetscReal R_star = R / R_0;
1116: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(R_0); /* - */
1117: PetscReal A_n = 0.0;
1118: PetscScalar u_sc;
1120: for (n = 1; n < N + 1; ++n) {
1121: const PetscReal x_n = user->zeroArray[n - 1];
1122: const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u));
1124: /* m , Cheng (7.404) */
1125: if (R_star != 0) {
1126: 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));
1127: }
1128: }
1129: if (R_star != 0) u_sc = PetscRealPart(u_inf) * (R_star - A_n) / R;
1130: else u_sc = PetscRealPart(u_inf) / R_0;
1131: u[0] = u_sc * x[0];
1132: u[1] = u_sc * x[1];
1133: u[2] = u_sc * x[2];
1134: }
1135: return PETSC_SUCCESS;
1136: }
1138: // Volumetric Strain
1139: static PetscErrorCode cryer_3d_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
1140: {
1141: AppCtx *user = (AppCtx *)ctx;
1142: Parameter *param;
1144: PetscCall(PetscBagGetData(user->bag, ¶m));
1145: if (time <= 0.0) {
1146: PetscCall(cryer_initial_eps(dim, time, x, Nc, u, ctx));
1147: } else {
1148: PetscScalar alpha = param->alpha; /* - */
1149: PetscScalar K_u = param->K_u; /* Pa */
1150: PetscScalar M = param->M; /* Pa */
1151: PetscScalar G = param->mu; /* Pa */
1152: PetscScalar P_0 = param->P_0; /* Pa */
1153: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
1154: PetscReal R_0 = user->xmax[1]; /* m */
1155: PetscInt N = user->niter, n;
1157: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
1158: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
1159: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
1160: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
1161: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
1162: PetscScalar u_inf = -P_0 * R_0 * (1. - 2. * nu) / (2. * G * (1. + nu)); /* m, Cheng (7.388) */
1164: PetscReal R = PetscSqrtReal(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
1165: PetscReal R_star = R / R_0;
1166: PetscReal tstar = PetscRealPart(c * time) / PetscSqr(R_0); /* - */
1167: PetscReal divA_n = 0.0;
1169: if (R_star < PETSC_SMALL) {
1170: for (n = 1; n < N + 1; ++n) {
1171: const PetscReal x_n = user->zeroArray[n - 1];
1172: const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u));
1174: 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));
1175: }
1176: } else {
1177: for (n = 1; n < N + 1; ++n) {
1178: const PetscReal x_n = user->zeroArray[n - 1];
1179: const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u));
1181: 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));
1182: }
1183: }
1184: 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));
1185: u[0] = PetscRealPart(u_inf) / R_0 * (3.0 - divA_n);
1186: }
1187: return PETSC_SUCCESS;
1188: }
1190: // Pressure
1191: static PetscErrorCode cryer_3d_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
1192: {
1193: AppCtx *user = (AppCtx *)ctx;
1194: Parameter *param;
1196: PetscCall(PetscBagGetData(user->bag, ¶m));
1197: if (time <= 0.0) {
1198: PetscCall(cryer_drainage_pressure(dim, time, x, Nc, u, ctx));
1199: } else {
1200: PetscScalar alpha = param->alpha; /* - */
1201: PetscScalar K_u = param->K_u; /* Pa */
1202: PetscScalar M = param->M; /* Pa */
1203: PetscScalar G = param->mu; /* Pa */
1204: PetscScalar P_0 = param->P_0; /* Pa */
1205: PetscReal R_0 = user->xmax[1]; /* m */
1206: PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
1207: PetscInt N = user->niter, n;
1209: PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
1210: PetscScalar eta = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G); /* -, Cheng (B.11) */
1211: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
1212: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
1213: PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
1214: PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */
1215: PetscReal R = PetscSqrtReal(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
1217: PetscReal R_star = R / R_0;
1218: PetscReal t_star = PetscRealPart(c * time) / PetscSqr(R_0);
1219: PetscReal A_x = 0.0;
1221: for (n = 1; n < N + 1; ++n) {
1222: const PetscReal x_n = user->zeroArray[n - 1];
1223: const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u));
1225: 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) */
1226: }
1227: u[0] = P_0 * A_x;
1228: }
1229: return PETSC_SUCCESS;
1230: }
1232: /* Boundary Kernels */
1233: 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[])
1234: {
1235: const PetscReal P = PetscRealPart(constants[5]);
1237: f0[0] = 0.0;
1238: f0[1] = P;
1239: }
1241: #if 0
1242: static void f0_mandel_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1243: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1244: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1245: PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1246: {
1247: // Uniform stress distribution
1248: /* PetscScalar xmax = 0.5;
1249: PetscScalar xmin = -0.5;
1250: PetscScalar ymax = 0.5;
1251: PetscScalar ymin = -0.5;
1252: PetscScalar P = constants[5];
1253: PetscScalar aL = (xmax - xmin) / 2.0;
1254: PetscScalar sigma_zz = -1.0*P / aL; */
1256: // Analytical (parabolic) stress distribution
1257: PetscReal a1, a2, am;
1258: PetscReal y1, y2, ym;
1260: PetscInt NITER = 500;
1261: PetscReal EPS = 0.000001;
1262: PetscReal zeroArray[500]; /* NITER */
1263: PetscReal xmax = 1.0;
1264: PetscReal xmin = 0.0;
1265: PetscReal ymax = 0.1;
1266: PetscReal ymin = 0.0;
1267: PetscReal lower[2], upper[2];
1269: lower[0] = xmin - (xmax - xmin) / 2.0;
1270: lower[1] = ymin - (ymax - ymin) / 2.0;
1271: upper[0] = xmax - (xmax - xmin) / 2.0;
1272: upper[1] = ymax - (ymax - ymin) / 2.0;
1274: xmin = lower[0];
1275: ymin = lower[1];
1276: xmax = upper[0];
1277: ymax = upper[1];
1279: PetscScalar G = constants[0];
1280: PetscScalar K_u = constants[1];
1281: PetscScalar alpha = constants[2];
1282: PetscScalar M = constants[3];
1283: PetscScalar kappa = constants[4];
1284: PetscScalar F = constants[5];
1286: PetscScalar K_d = K_u - alpha*alpha*M;
1287: PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G));
1288: PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G));
1289: PetscReal aL = (xmax - xmin) / 2.0;
1290: PetscReal c = PetscRealPart(((2.0*kappa*G) * (1.0 - nu) * (nu_u - nu)) / (alpha*alpha * (1.0 - 2.0*nu) * (1.0 - nu_u)));
1291: PetscScalar B = (3.0 * (nu_u - nu)) / ( alpha * (1.0 - 2.0*nu) * (1.0 + nu_u));
1292: PetscScalar A1 = 3.0 / (B * (1.0 + nu_u));
1293: PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu);
1295: // Generate zero values
1296: for (PetscInt i=1; i < NITER+1; i++)
1297: {
1298: a1 = ((PetscReal) i - 1.0) * PETSC_PI * PETSC_PI / 4.0 + EPS;
1299: a2 = a1 + PETSC_PI/2;
1300: for (PetscInt j=0; j<NITER; j++)
1301: {
1302: y1 = PetscTanReal(a1) - PetscRealPart(A1/A2)*a1;
1303: y2 = PetscTanReal(a2) - PetscRealPart(A1/A2)*a2;
1304: am = (a1 + a2)/2.0;
1305: ym = PetscTanReal(am) - PetscRealPart(A1/A2)*am;
1306: if ((ym*y1) > 0)
1307: {
1308: a1 = am;
1309: } else {
1310: a2 = am;
1311: }
1312: if (PetscAbsReal(y2) < EPS)
1313: {
1314: am = a2;
1315: }
1316: }
1317: zeroArray[i-1] = am;
1318: }
1320: // Solution for sigma_zz
1321: PetscScalar A_x = 0.0;
1322: PetscScalar B_x = 0.0;
1324: for (PetscInt n=1; n < NITER+1; n++)
1325: {
1326: PetscReal alpha_n = zeroArray[n-1];
1328: 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)));
1329: 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)));
1330: }
1332: PetscScalar sigma_zz = -1.0*(F/aL) - ((2.0*F)/aL) * (A2/A1) * A_x + ((2.0*F)/aL) * B_x;
1334: if (x[1] == ymax) {
1335: f0[1] += sigma_zz;
1336: } else if (x[1] == ymin) {
1337: f0[1] -= sigma_zz;
1338: }
1339: }
1340: #endif
1342: 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[])
1343: {
1344: const PetscReal P_0 = PetscRealPart(constants[5]);
1345: PetscInt d;
1347: for (d = 0; d < dim; ++d) f0[d] = -P_0 * n[d];
1348: }
1350: /* Standard Kernels - Residual */
1351: /* f0_e */
1352: 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[])
1353: {
1354: PetscInt d;
1356: for (d = 0; d < dim; ++d) f0[0] += u_x[d * dim + d];
1357: f0[0] -= u[uOff[1]];
1358: }
1360: /* f0_p */
1361: 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[])
1362: {
1363: const PetscReal alpha = PetscRealPart(constants[2]);
1364: const PetscReal M = PetscRealPart(constants[3]);
1366: f0[0] += alpha * u_t[uOff[1]];
1367: f0[0] += u_t[uOff[2]] / M;
1368: if (f0[0] != f0[0]) abort();
1369: }
1371: /* f1_u */
1372: 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[])
1373: {
1374: const PetscInt Nc = dim;
1375: const PetscReal G = PetscRealPart(constants[0]);
1376: const PetscReal K_u = PetscRealPart(constants[1]);
1377: const PetscReal alpha = PetscRealPart(constants[2]);
1378: const PetscReal M = PetscRealPart(constants[3]);
1379: const PetscReal K_d = K_u - alpha * alpha * M;
1380: const PetscReal lambda = K_d - (2.0 * G) / 3.0;
1381: PetscInt d;
1383: for (PetscInt c = 0; c < Nc; ++c) {
1384: for (d = 0; d < dim; ++d) f1[c * dim + d] -= G * (u_x[c * dim + d] + u_x[d * dim + c]);
1385: f1[c * dim + c] -= lambda * u[uOff[1]];
1386: f1[c * dim + c] += alpha * u[uOff[2]];
1387: }
1388: }
1390: /* f1_p */
1391: 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[])
1392: {
1393: const PetscReal kappa = PetscRealPart(constants[4]);
1394: PetscInt d;
1396: for (d = 0; d < dim; ++d) f1[d] += kappa * u_x[uOff_x[2] + d];
1397: }
1399: /*
1400: \partial_df \phi_fc g_{fc,gc,df,dg} \partial_dg \phi_gc
1402: \partial_df \phi_fc \lambda \delta_{fc,df} \sum_gc \partial_dg \phi_gc \delta_{gc,dg}
1403: = \partial_fc \phi_fc \sum_gc \partial_gc \phi_gc
1404: */
1406: /* Standard Kernels - Jacobian */
1407: /* g0_ee */
1408: 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[])
1409: {
1410: g0[0] = -1.0;
1411: }
1413: /* g0_pe */
1414: 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[])
1415: {
1416: const PetscReal alpha = PetscRealPart(constants[2]);
1418: g0[0] = u_tShift * alpha;
1419: }
1421: /* g0_pp */
1422: 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[])
1423: {
1424: const PetscReal M = PetscRealPart(constants[3]);
1426: g0[0] = u_tShift / M;
1427: }
1429: /* g1_eu */
1430: 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[])
1431: {
1432: PetscInt d;
1433: for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */
1434: }
1436: /* g2_ue */
1437: 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[])
1438: {
1439: const PetscReal G = PetscRealPart(constants[0]);
1440: const PetscReal K_u = PetscRealPart(constants[1]);
1441: const PetscReal alpha = PetscRealPart(constants[2]);
1442: const PetscReal M = PetscRealPart(constants[3]);
1443: const PetscReal K_d = K_u - alpha * alpha * M;
1444: const PetscReal lambda = K_d - (2.0 * G) / 3.0;
1445: PetscInt d;
1447: for (d = 0; d < dim; ++d) g2[d * dim + d] -= lambda;
1448: }
1449: /* g2_up */
1450: 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[])
1451: {
1452: const PetscReal alpha = PetscRealPart(constants[2]);
1453: PetscInt d;
1455: for (d = 0; d < dim; ++d) g2[d * dim + d] += alpha;
1456: }
1458: /* g3_uu */
1459: 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[])
1460: {
1461: const PetscInt Nc = dim;
1462: const PetscReal G = PetscRealPart(constants[0]);
1464: for (PetscInt c = 0; c < Nc; ++c) {
1465: for (PetscInt d = 0; d < dim; ++d) {
1466: g3[((c * Nc + c) * dim + d) * dim + d] -= G;
1467: g3[((c * Nc + d) * dim + d) * dim + c] -= G;
1468: }
1469: }
1470: }
1472: /* g3_pp */
1473: 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[])
1474: {
1475: const PetscReal kappa = PetscRealPart(constants[4]);
1476: for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] += kappa;
1477: }
1479: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
1480: {
1481: PetscInt sol;
1483: PetscFunctionBeginUser;
1484: options->solType = SOL_QUADRATIC_TRIG;
1485: options->niter = 500;
1486: options->eps = PETSC_SMALL;
1487: options->dtInitial = -1.0;
1488: PetscOptionsBegin(comm, "", "Biot Poroelasticity Options", "DMPLEX");
1489: PetscCall(PetscOptionsInt("-niter", "Number of series term iterations in exact solutions", "ex53.c", options->niter, &options->niter, NULL));
1490: sol = options->solType;
1491: PetscCall(PetscOptionsEList("-sol_type", "Type of exact solution", "ex53.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL));
1492: options->solType = (SolutionType)sol;
1493: PetscCall(PetscOptionsReal("-eps", "Precision value for root finding", "ex53.c", options->eps, &options->eps, NULL));
1494: PetscCall(PetscOptionsReal("-dt_initial", "Override the initial timestep", "ex53.c", options->dtInitial, &options->dtInitial, NULL));
1495: PetscOptionsEnd();
1496: PetscFunctionReturn(PETSC_SUCCESS);
1497: }
1499: static PetscErrorCode mandelZeros(MPI_Comm comm, AppCtx *ctx, Parameter *param)
1500: {
1501: //PetscBag bag;
1502: PetscReal a1, a2, am;
1503: PetscReal y1, y2, ym;
1505: PetscFunctionBeginUser;
1506: //PetscCall(PetscBagGetData(ctx->bag, ¶m));
1507: PetscInt NITER = ctx->niter;
1508: PetscReal EPS = ctx->eps;
1509: //const PetscScalar YMAX = param->ymax;
1510: //const PetscScalar YMIN = param->ymin;
1511: PetscScalar alpha = param->alpha;
1512: PetscScalar K_u = param->K_u;
1513: PetscScalar M = param->M;
1514: PetscScalar G = param->mu;
1515: //const PetscScalar k = param->k;
1516: //const PetscScalar mu_f = param->mu_f;
1517: //const PetscScalar P_0 = param->P_0;
1519: PetscScalar K_d = K_u - alpha * alpha * M;
1520: PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
1521: PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
1522: //const PetscScalar kappa = k / mu_f;
1524: // Generate zero values
1525: for (PetscInt i = 1; i < NITER + 1; i++) {
1526: a1 = ((PetscReal)i - 1.0) * PETSC_PI * PETSC_PI / 4.0 + EPS;
1527: a2 = a1 + PETSC_PI / 2;
1528: am = a1;
1529: for (PetscInt j = 0; j < NITER; j++) {
1530: y1 = PetscTanReal(a1) - PetscRealPart((1.0 - nu) / (nu_u - nu)) * a1;
1531: y2 = PetscTanReal(a2) - PetscRealPart((1.0 - nu) / (nu_u - nu)) * a2;
1532: am = (a1 + a2) / 2.0;
1533: ym = PetscTanReal(am) - PetscRealPart((1.0 - nu) / (nu_u - nu)) * am;
1534: if ((ym * y1) > 0) {
1535: a1 = am;
1536: } else {
1537: a2 = am;
1538: }
1539: if (PetscAbsReal(y2) < EPS) am = a2;
1540: }
1541: ctx->zeroArray[i - 1] = am;
1542: }
1543: PetscFunctionReturn(PETSC_SUCCESS);
1544: }
1546: static PetscReal CryerFunction(PetscReal nu_u, PetscReal nu, PetscReal x)
1547: {
1548: return PetscTanReal(PetscSqrtReal(x)) * (6.0 * (nu_u - nu) - (1.0 - nu) * (1.0 + nu_u) * x) - (6.0 * (nu_u - nu) * PetscSqrtReal(x));
1549: }
1551: static PetscErrorCode cryerZeros(MPI_Comm comm, AppCtx *ctx, Parameter *param)
1552: {
1553: PetscReal alpha = PetscRealPart(param->alpha); /* - */
1554: PetscReal K_u = PetscRealPart(param->K_u); /* Pa */
1555: PetscReal M = PetscRealPart(param->M); /* Pa */
1556: PetscReal G = PetscRealPart(param->mu); /* Pa */
1557: PetscInt N = ctx->niter, n;
1559: PetscReal K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */
1560: PetscReal nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */
1561: PetscReal nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */
1563: PetscFunctionBeginUser;
1564: for (n = 1; n < N + 1; ++n) {
1565: PetscReal tol = PetscPowReal(n, 1.5) * ctx->eps;
1566: PetscReal a1 = 0., a2 = 0., am = 0.;
1567: PetscReal y1, y2, ym;
1568: PetscInt j, k = n - 1;
1570: y1 = y2 = 1.;
1571: while (y1 * y2 > 0) {
1572: ++k;
1573: a1 = PetscSqr(n * PETSC_PI) - k * PETSC_PI;
1574: a2 = PetscSqr(n * PETSC_PI) + k * PETSC_PI;
1575: y1 = CryerFunction(nu_u, nu, a1);
1576: y2 = CryerFunction(nu_u, nu, a2);
1577: }
1578: for (j = 0; j < 50000; ++j) {
1579: y1 = CryerFunction(nu_u, nu, a1);
1580: y2 = CryerFunction(nu_u, nu, a2);
1581: 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);
1582: am = (a1 + a2) / 2.0;
1583: ym = CryerFunction(nu_u, nu, am);
1584: if ((ym * y1) < 0) a2 = am;
1585: else a1 = am;
1586: if (PetscAbsReal(ym) < tol) break;
1587: }
1588: PetscCheck(PetscAbsReal(ym) < tol, comm, PETSC_ERR_PLIB, "Root finding did not converge for root %" PetscInt_FMT " (%g)", n, (double)PetscAbsReal(ym));
1589: ctx->zeroArray[n - 1] = am;
1590: }
1591: PetscFunctionReturn(PETSC_SUCCESS);
1592: }
1594: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
1595: {
1596: PetscBag bag;
1597: Parameter *p;
1599: PetscFunctionBeginUser;
1600: /* setup PETSc parameter bag */
1601: PetscCall(PetscBagGetData(ctx->bag, &p));
1602: PetscCall(PetscBagSetName(ctx->bag, "par", "Poroelastic Parameters"));
1603: bag = ctx->bag;
1604: if (ctx->solType == SOL_TERZAGHI) {
1605: // Realistic values - Terzaghi
1606: PetscCall(PetscBagRegisterScalar(bag, &p->mu, 3.0, "mu", "Shear Modulus, Pa"));
1607: PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 9.76, "K_u", "Undrained Bulk Modulus, Pa"));
1608: PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -"));
1609: PetscCall(PetscBagRegisterScalar(bag, &p->M, 16.0, "M", "Biot Modulus, Pa"));
1610: PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2"));
1611: PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s"));
1612: PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa"));
1613: } else if (ctx->solType == SOL_MANDEL) {
1614: // Realistic values - Mandel
1615: PetscCall(PetscBagRegisterScalar(bag, &p->mu, 0.75, "mu", "Shear Modulus, Pa"));
1616: PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 2.6941176470588233, "K_u", "Undrained Bulk Modulus, Pa"));
1617: PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -"));
1618: PetscCall(PetscBagRegisterScalar(bag, &p->M, 4.705882352941176, "M", "Biot Modulus, Pa"));
1619: PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2"));
1620: PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s"));
1621: PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa"));
1622: } else if (ctx->solType == SOL_CRYER) {
1623: // Realistic values - Mandel
1624: PetscCall(PetscBagRegisterScalar(bag, &p->mu, 0.75, "mu", "Shear Modulus, Pa"));
1625: PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 2.6941176470588233, "K_u", "Undrained Bulk Modulus, Pa"));
1626: PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -"));
1627: PetscCall(PetscBagRegisterScalar(bag, &p->M, 4.705882352941176, "M", "Biot Modulus, Pa"));
1628: PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2"));
1629: PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s"));
1630: PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa"));
1631: } else {
1632: // Nonsense values
1633: PetscCall(PetscBagRegisterScalar(bag, &p->mu, 1.0, "mu", "Shear Modulus, Pa"));
1634: PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 1.0, "K_u", "Undrained Bulk Modulus, Pa"));
1635: PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 1.0, "alpha", "Biot Effective Stress Coefficient, -"));
1636: PetscCall(PetscBagRegisterScalar(bag, &p->M, 1.0, "M", "Biot Modulus, Pa"));
1637: PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.0, "k", "Isotropic Permeability, m**2"));
1638: PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s"));
1639: PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa"));
1640: }
1641: PetscCall(PetscBagSetFromOptions(bag));
1642: {
1643: PetscScalar K_d = p->K_u - p->alpha * p->alpha * p->M;
1644: PetscScalar nu_u = (3.0 * p->K_u - 2.0 * p->mu) / (2.0 * (3.0 * p->K_u + p->mu));
1645: PetscScalar nu = (3.0 * K_d - 2.0 * p->mu) / (2.0 * (3.0 * K_d + p->mu));
1646: PetscScalar S = (3.0 * p->K_u + 4.0 * p->mu) / (p->M * (3.0 * K_d + 4.0 * p->mu));
1647: PetscReal c = PetscRealPart((p->k / p->mu_f) / S);
1649: PetscViewer viewer;
1650: PetscViewerFormat format;
1651: PetscBool flg;
1653: switch (ctx->solType) {
1654: case SOL_QUADRATIC_LINEAR:
1655: case SOL_QUADRATIC_TRIG:
1656: case SOL_TRIG_LINEAR:
1657: ctx->t_r = PetscSqr(ctx->xmax[0] - ctx->xmin[0]) / c;
1658: break;
1659: case SOL_TERZAGHI:
1660: ctx->t_r = PetscSqr(2.0 * (ctx->xmax[1] - ctx->xmin[1])) / c;
1661: break;
1662: case SOL_MANDEL:
1663: ctx->t_r = PetscSqr(2.0 * (ctx->xmax[1] - ctx->xmin[1])) / c;
1664: break;
1665: case SOL_CRYER:
1666: ctx->t_r = PetscSqr(ctx->xmax[1]) / c;
1667: break;
1668: default:
1669: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%d)", solutionTypes[PetscMin(ctx->solType, NUM_SOLUTION_TYPES)], ctx->solType);
1670: }
1671: PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg));
1672: if (flg) {
1673: PetscCall(PetscViewerPushFormat(viewer, format));
1674: PetscCall(PetscBagView(bag, viewer));
1675: PetscCall(PetscViewerFlush(viewer));
1676: PetscCall(PetscViewerPopFormat(viewer));
1677: PetscCall(PetscViewerDestroy(&viewer));
1678: 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)))));
1679: PetscCall(PetscPrintf(comm, " Relaxation time: %g\n", (double)ctx->t_r));
1680: }
1681: }
1682: PetscFunctionReturn(PETSC_SUCCESS);
1683: }
1685: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
1686: {
1687: PetscFunctionBeginUser;
1688: PetscCall(DMCreate(comm, dm));
1689: PetscCall(DMSetType(*dm, DMPLEX));
1690: PetscCall(DMSetFromOptions(*dm));
1691: PetscCall(DMSetApplicationContext(*dm, user));
1692: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
1693: PetscCall(DMGetBoundingBox(*dm, user->xmin, user->xmax));
1694: PetscFunctionReturn(PETSC_SUCCESS);
1695: }
1697: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
1698: {
1699: PetscErrorCode (*exact[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
1700: PetscErrorCode (*exact_t[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
1701: PetscDS ds;
1702: DMLabel label;
1703: PetscWeakForm wf;
1704: Parameter *param;
1705: PetscInt id_mandel[2];
1706: PetscInt comp[1];
1707: PetscInt comp_mandel[2];
1708: PetscInt dim, id, bd, f;
1710: PetscFunctionBeginUser;
1711: PetscCall(DMGetLabel(dm, "marker", &label));
1712: PetscCall(DMGetDS(dm, &ds));
1713: PetscCall(PetscDSGetSpatialDimension(ds, &dim));
1714: PetscCall(PetscBagGetData(user->bag, ¶m));
1715: exact_t[0] = exact_t[1] = exact_t[2] = zero;
1717: /* Setup Problem Formulation and Boundary Conditions */
1718: switch (user->solType) {
1719: case SOL_QUADRATIC_LINEAR:
1720: PetscCall(PetscDSSetResidual(ds, 0, f0_quadratic_linear_u, f1_u));
1721: PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1722: PetscCall(PetscDSSetResidual(ds, 2, f0_quadratic_linear_p, f1_p));
1723: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1724: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1725: PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1726: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1727: PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1728: PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1729: PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));
1730: exact[0] = quadratic_u;
1731: exact[1] = linear_eps;
1732: exact[2] = linear_linear_p;
1733: exact_t[2] = linear_linear_p_t;
1735: id = 1;
1736: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)exact[0], NULL, user, NULL));
1737: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", label, 1, &id, 2, 0, NULL, (PetscVoidFn *)exact[2], (PetscVoidFn *)exact_t[2], user, NULL));
1738: break;
1739: case SOL_TRIG_LINEAR:
1740: PetscCall(PetscDSSetResidual(ds, 0, f0_trig_linear_u, f1_u));
1741: PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1742: PetscCall(PetscDSSetResidual(ds, 2, f0_trig_linear_p, f1_p));
1743: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1744: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1745: PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1746: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1747: PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1748: PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1749: PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));
1750: exact[0] = trig_u;
1751: exact[1] = trig_eps;
1752: exact[2] = trig_linear_p;
1753: exact_t[2] = trig_linear_p_t;
1755: id = 1;
1756: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)exact[0], NULL, user, NULL));
1757: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", label, 1, &id, 2, 0, NULL, (PetscVoidFn *)exact[2], (PetscVoidFn *)exact_t[2], user, NULL));
1758: break;
1759: case SOL_QUADRATIC_TRIG:
1760: PetscCall(PetscDSSetResidual(ds, 0, f0_quadratic_trig_u, f1_u));
1761: PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1762: PetscCall(PetscDSSetResidual(ds, 2, f0_quadratic_trig_p, f1_p));
1763: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1764: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1765: PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1766: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1767: PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1768: PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1769: PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));
1770: exact[0] = quadratic_u;
1771: exact[1] = linear_eps;
1772: exact[2] = linear_trig_p;
1773: exact_t[2] = linear_trig_p_t;
1775: id = 1;
1776: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)exact[0], NULL, user, NULL));
1777: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", label, 1, &id, 2, 0, NULL, (PetscVoidFn *)exact[2], (PetscVoidFn *)exact_t[2], user, NULL));
1778: break;
1779: case SOL_TERZAGHI:
1780: PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u));
1781: PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1782: PetscCall(PetscDSSetResidual(ds, 2, f0_p, f1_p));
1783: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1784: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1785: PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1786: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1787: PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1788: PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1789: PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));
1791: exact[0] = terzaghi_2d_u;
1792: exact[1] = terzaghi_2d_eps;
1793: exact[2] = terzaghi_2d_p;
1794: exact_t[0] = terzaghi_2d_u_t;
1795: exact_t[1] = terzaghi_2d_eps_t;
1796: exact_t[2] = terzaghi_2d_p_t;
1798: id = 1;
1799: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "vertical stress", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
1800: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1801: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_terzaghi_bd_u, 0, NULL));
1803: id = 3;
1804: comp[0] = 1;
1805: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed base", label, 1, &id, 0, 1, comp, (PetscVoidFn *)zero, NULL, user, NULL));
1806: id = 2;
1807: comp[0] = 0;
1808: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed side", label, 1, &id, 0, 1, comp, (PetscVoidFn *)zero, NULL, user, NULL));
1809: id = 4;
1810: comp[0] = 0;
1811: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed side", label, 1, &id, 0, 1, comp, (PetscVoidFn *)zero, NULL, user, NULL));
1812: id = 1;
1813: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", label, 1, &id, 2, 0, NULL, (PetscVoidFn *)terzaghi_drainage_pressure, NULL, user, NULL));
1814: break;
1815: case SOL_MANDEL:
1816: PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u));
1817: PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1818: PetscCall(PetscDSSetResidual(ds, 2, f0_p, f1_p));
1819: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1820: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1821: PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1822: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1823: PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1824: PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1825: PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));
1827: PetscCall(mandelZeros(PETSC_COMM_WORLD, user, param));
1829: exact[0] = mandel_2d_u;
1830: exact[1] = mandel_2d_eps;
1831: exact[2] = mandel_2d_p;
1832: exact_t[0] = mandel_2d_u_t;
1833: exact_t[1] = mandel_2d_eps_t;
1834: exact_t[2] = mandel_2d_p_t;
1836: id_mandel[0] = 3;
1837: id_mandel[1] = 1;
1838: //comp[0] = 1;
1839: comp_mandel[0] = 0;
1840: comp_mandel[1] = 1;
1841: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "vertical stress", label, 2, id_mandel, 0, 2, comp_mandel, (PetscVoidFn *)mandel_2d_u, NULL, user, NULL));
1842: //PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "vertical stress", "marker", 0, 1, comp, NULL, 2, id_mandel, user));
1843: //PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed base", "marker", 0, 1, comp, (PetscVoidFn *) zero, 2, id_mandel, user));
1844: //PetscCall(PetscDSSetBdResidual(ds, 0, f0_mandel_bd_u, NULL));
1846: id_mandel[0] = 2;
1847: id_mandel[1] = 4;
1848: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", label, 2, id_mandel, 2, 0, NULL, (PetscVoidFn *)zero, NULL, user, NULL));
1849: break;
1850: case SOL_CRYER:
1851: PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u));
1852: PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1853: PetscCall(PetscDSSetResidual(ds, 2, f0_p, f1_p));
1854: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1855: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1856: PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1857: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1858: PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1859: PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1860: PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));
1862: PetscCall(cryerZeros(PETSC_COMM_WORLD, user, param));
1864: exact[0] = cryer_3d_u;
1865: exact[1] = cryer_3d_eps;
1866: exact[2] = cryer_3d_p;
1868: id = 1;
1869: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "normal stress", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
1870: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1871: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_cryer_bd_u, 0, NULL));
1873: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", label, 1, &id, 2, 0, NULL, (PetscVoidFn *)cryer_drainage_pressure, NULL, user, NULL));
1874: break;
1875: default:
1876: SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%d)", solutionTypes[PetscMin(user->solType, NUM_SOLUTION_TYPES)], user->solType);
1877: }
1878: for (f = 0; f < 3; ++f) {
1879: PetscCall(PetscDSSetExactSolution(ds, f, exact[f], user));
1880: PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, f, exact_t[f], user));
1881: }
1883: /* Setup constants */
1884: {
1885: PetscScalar constants[6];
1886: constants[0] = param->mu; /* shear modulus, Pa */
1887: constants[1] = param->K_u; /* undrained bulk modulus, Pa */
1888: constants[2] = param->alpha; /* Biot effective stress coefficient, - */
1889: constants[3] = param->M; /* Biot modulus, Pa */
1890: constants[4] = param->k / param->mu_f; /* Darcy coefficient, m**2 / Pa*s */
1891: constants[5] = param->P_0; /* Magnitude of Vertical Stress, Pa */
1892: PetscCall(PetscDSSetConstants(ds, 6, constants));
1893: }
1894: PetscFunctionReturn(PETSC_SUCCESS);
1895: }
1897: static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace)
1898: {
1899: PetscFunctionBeginUser;
1900: PetscCall(DMPlexCreateRigidBody(dm, origField, nullspace));
1901: PetscFunctionReturn(PETSC_SUCCESS);
1902: }
1904: static PetscErrorCode SetupFE(DM dm, PetscInt Nf, PetscInt Nc[], const char *name[], PetscErrorCode (*setup)(DM, AppCtx *), PetscCtx ctx)
1905: {
1906: AppCtx *user = (AppCtx *)ctx;
1907: DM cdm = dm;
1908: PetscFE fe;
1909: PetscQuadrature q = NULL, fq = NULL;
1910: char prefix[PETSC_MAX_PATH_LEN];
1911: PetscInt dim;
1912: PetscBool simplex;
1914: PetscFunctionBeginUser;
1915: /* Create finite element */
1916: PetscCall(DMGetDimension(dm, &dim));
1917: PetscCall(DMPlexIsSimplex(dm, &simplex));
1918: for (PetscInt f = 0; f < Nf; ++f) {
1919: PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name[f]));
1920: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc[f], simplex, name[f] ? prefix : NULL, -1, &fe));
1921: PetscCall(PetscObjectSetName((PetscObject)fe, name[f]));
1922: if (!q) PetscCall(PetscFEGetQuadrature(fe, &q));
1923: if (!fq) PetscCall(PetscFEGetFaceQuadrature(fe, &fq));
1924: PetscCall(PetscFESetQuadrature(fe, q));
1925: PetscCall(PetscFESetFaceQuadrature(fe, fq));
1926: PetscCall(DMSetField(dm, f, NULL, (PetscObject)fe));
1927: PetscCall(PetscFEDestroy(&fe));
1928: }
1929: PetscCall(DMCreateDS(dm));
1930: PetscCall((*setup)(dm, user));
1931: while (cdm) {
1932: PetscCall(DMCopyDisc(dm, cdm));
1933: if (0) PetscCall(DMSetNearNullSpaceConstructor(cdm, 0, CreateElasticityNullSpace));
1934: /* TODO: Check whether the boundary of coarse meshes is marked */
1935: PetscCall(DMGetCoarseDM(cdm, &cdm));
1936: }
1937: PetscCall(PetscFEDestroy(&fe));
1938: PetscFunctionReturn(PETSC_SUCCESS);
1939: }
1941: static PetscErrorCode SetInitialConditions(TS ts, Vec u)
1942: {
1943: DM dm;
1944: PetscReal t;
1946: PetscFunctionBeginUser;
1947: PetscCall(TSGetDM(ts, &dm));
1948: PetscCall(TSGetTime(ts, &t));
1949: if (t <= 0.0) {
1950: PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
1951: void *ctxs[3];
1952: AppCtx *ctx;
1954: PetscCall(DMGetApplicationContext(dm, &ctx));
1955: switch (ctx->solType) {
1956: case SOL_TERZAGHI:
1957: funcs[0] = terzaghi_initial_u;
1958: ctxs[0] = ctx;
1959: funcs[1] = terzaghi_initial_eps;
1960: ctxs[1] = ctx;
1961: funcs[2] = terzaghi_drainage_pressure;
1962: ctxs[2] = ctx;
1963: PetscCall(DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u));
1964: break;
1965: case SOL_MANDEL:
1966: funcs[0] = mandel_initial_u;
1967: ctxs[0] = ctx;
1968: funcs[1] = mandel_initial_eps;
1969: ctxs[1] = ctx;
1970: funcs[2] = mandel_drainage_pressure;
1971: ctxs[2] = ctx;
1972: PetscCall(DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u));
1973: break;
1974: case SOL_CRYER:
1975: funcs[0] = cryer_initial_u;
1976: ctxs[0] = ctx;
1977: funcs[1] = cryer_initial_eps;
1978: ctxs[1] = ctx;
1979: funcs[2] = cryer_drainage_pressure;
1980: ctxs[2] = ctx;
1981: PetscCall(DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u));
1982: break;
1983: default:
1984: PetscCall(DMComputeExactSolution(dm, t, u, NULL));
1985: }
1986: } else {
1987: PetscCall(DMComputeExactSolution(dm, t, u, NULL));
1988: }
1989: PetscFunctionReturn(PETSC_SUCCESS);
1990: }
1992: /* Need to create Viewer each time because HDF5 can get corrupted */
1993: static PetscErrorCode SolutionMonitor(TS ts, PetscInt steps, PetscReal time, Vec u, void *mctx)
1994: {
1995: DM dm;
1996: Vec exact;
1997: PetscViewer viewer;
1998: PetscViewerFormat format;
1999: PetscOptions options;
2000: const char *prefix;
2002: PetscFunctionBeginUser;
2003: PetscCall(TSGetDM(ts, &dm));
2004: PetscCall(PetscObjectGetOptions((PetscObject)ts, &options));
2005: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, &prefix));
2006: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), options, prefix, "-monitor_solution", &viewer, &format, NULL));
2007: PetscCall(DMGetGlobalVector(dm, &exact));
2008: PetscCall(DMComputeExactSolution(dm, time, exact, NULL));
2009: PetscCall(DMSetOutputSequenceNumber(dm, steps, time));
2010: PetscCall(VecView(exact, viewer));
2011: PetscCall(VecView(u, viewer));
2012: PetscCall(DMRestoreGlobalVector(dm, &exact));
2013: {
2014: PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, PetscCtx ctx);
2015: void **ectxs;
2016: PetscReal *err;
2017: PetscInt Nf;
2019: PetscCall(DMGetNumFields(dm, &Nf));
2020: PetscCall(PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err));
2021: {
2022: PetscInt Nds;
2024: PetscCall(DMGetNumDS(dm, &Nds));
2025: for (PetscInt s = 0; s < Nds; ++s) {
2026: PetscDS ds;
2027: DMLabel label;
2028: IS fieldIS;
2029: const PetscInt *fields;
2030: PetscInt dsNf;
2032: PetscCall(DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds, NULL));
2033: PetscCall(PetscDSGetNumFields(ds, &dsNf));
2034: PetscCall(ISGetIndices(fieldIS, &fields));
2035: for (PetscInt f = 0; f < dsNf; ++f) {
2036: const PetscInt field = fields[f];
2037: PetscCall(PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]));
2038: }
2039: PetscCall(ISRestoreIndices(fieldIS, &fields));
2040: }
2041: }
2042: PetscCall(DMComputeL2FieldDiff(dm, time, exacts, ectxs, u, err));
2043: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "Time: %g L_2 Error: [", (double)time));
2044: for (PetscInt f = 0; f < Nf; ++f) {
2045: if (f) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), ", "));
2046: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%g", (double)err[f]));
2047: }
2048: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "]\n"));
2049: PetscCall(PetscFree3(exacts, ectxs, err));
2050: }
2051: PetscCall(PetscViewerDestroy(&viewer));
2052: PetscFunctionReturn(PETSC_SUCCESS);
2053: }
2055: static PetscErrorCode SetupMonitor(TS ts, AppCtx *ctx)
2056: {
2057: PetscViewer viewer;
2058: PetscViewerFormat format;
2059: PetscOptions options;
2060: const char *prefix;
2061: PetscBool flg;
2063: PetscFunctionBeginUser;
2064: PetscCall(PetscObjectGetOptions((PetscObject)ts, &options));
2065: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, &prefix));
2066: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), options, prefix, "-monitor_solution", &viewer, &format, &flg));
2067: if (flg) PetscCall(TSMonitorSet(ts, SolutionMonitor, ctx, NULL));
2068: PetscCall(PetscViewerDestroy(&viewer));
2069: PetscFunctionReturn(PETSC_SUCCESS);
2070: }
2072: static PetscErrorCode TSAdaptChoose_Terzaghi(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept, PetscReal *wlte, PetscReal *wltea, PetscReal *wlter)
2073: {
2074: static PetscReal dtTarget = -1.0;
2075: PetscReal dtInitial;
2076: DM dm;
2077: AppCtx *ctx;
2078: PetscInt step;
2080: PetscFunctionBeginUser;
2081: PetscCall(TSGetDM(ts, &dm));
2082: PetscCall(DMGetApplicationContext(dm, &ctx));
2083: PetscCall(TSGetStepNumber(ts, &step));
2084: dtInitial = ctx->dtInitial < 0.0 ? 1.0e-4 * ctx->t_r : ctx->dtInitial;
2085: if (!step) {
2086: if (PetscAbsReal(dtInitial - h) > PETSC_SMALL) {
2087: *accept = PETSC_FALSE;
2088: *next_h = dtInitial;
2089: dtTarget = h;
2090: } else {
2091: *accept = PETSC_TRUE;
2092: *next_h = dtTarget < 0.0 ? dtInitial : dtTarget;
2093: dtTarget = -1.0;
2094: }
2095: } else {
2096: *accept = PETSC_TRUE;
2097: *next_h = h;
2098: }
2099: *next_sc = 0; /* Reuse the same order scheme */
2100: *wlte = -1; /* Weighted local truncation error was not evaluated */
2101: *wltea = -1; /* Weighted absolute local truncation error was not evaluated */
2102: *wlter = -1; /* Weighted relative local truncation error was not evaluated */
2103: PetscFunctionReturn(PETSC_SUCCESS);
2104: }
2106: int main(int argc, char **argv)
2107: {
2108: AppCtx ctx; /* User-defined work context */
2109: DM dm; /* Problem specification */
2110: TS ts; /* Time Series / Nonlinear solver */
2111: Vec u; /* Solutions */
2112: const char *name[3] = {"displacement", "tracestrain", "pressure"};
2113: PetscReal t;
2114: PetscInt dim, Nc[3];
2116: PetscFunctionBeginUser;
2117: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
2118: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
2119: PetscCall(PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &ctx.bag));
2120: PetscCall(PetscMalloc1(ctx.niter, &ctx.zeroArray));
2121: PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm));
2122: PetscCall(SetupParameters(PETSC_COMM_WORLD, &ctx));
2123: /* Primal System */
2124: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
2125: PetscCall(DMSetApplicationContext(dm, &ctx));
2126: PetscCall(TSSetDM(ts, dm));
2128: PetscCall(DMGetDimension(dm, &dim));
2129: Nc[0] = dim;
2130: Nc[1] = 1;
2131: Nc[2] = 1;
2133: PetscCall(SetupFE(dm, 3, Nc, name, SetupPrimalProblem, &ctx));
2134: PetscCall(DMCreateGlobalVector(dm, &u));
2135: PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx));
2136: PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx));
2137: PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx));
2138: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
2139: PetscCall(TSSetFromOptions(ts));
2140: PetscCall(TSSetComputeInitialCondition(ts, SetInitialConditions));
2141: PetscCall(SetupMonitor(ts, &ctx));
2143: if (ctx.solType != SOL_QUADRATIC_TRIG) {
2144: TSAdapt adapt;
2146: PetscCall(TSGetAdapt(ts, &adapt));
2147: adapt->ops->choose = TSAdaptChoose_Terzaghi;
2148: }
2149: if (ctx.solType == SOL_CRYER) {
2150: Mat J;
2151: MatNullSpace sp;
2153: PetscCall(TSSetUp(ts));
2154: PetscCall(TSGetIJacobian(ts, &J, NULL, NULL, NULL));
2155: PetscCall(DMPlexCreateRigidBody(dm, 0, &sp));
2156: PetscCall(MatSetNullSpace(J, sp));
2157: PetscCall(MatNullSpaceDestroy(&sp));
2158: }
2159: PetscCall(TSGetTime(ts, &t));
2160: PetscCall(DMSetOutputSequenceNumber(dm, 0, t));
2161: PetscCall(DMTSCheckFromOptions(ts, u));
2162: PetscCall(SetInitialConditions(ts, u));
2163: PetscCall(PetscObjectSetName((PetscObject)u, "solution"));
2164: PetscCall(TSSolve(ts, u));
2165: PetscCall(DMTSCheckFromOptions(ts, u));
2166: PetscCall(TSGetSolution(ts, &u));
2167: PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));
2169: /* Cleanup */
2170: PetscCall(VecDestroy(&u));
2171: PetscCall(TSDestroy(&ts));
2172: PetscCall(DMDestroy(&dm));
2173: PetscCall(PetscBagDestroy(&ctx.bag));
2174: PetscCall(PetscFree(ctx.zeroArray));
2175: PetscCall(PetscFinalize());
2176: return 0;
2177: }
2179: /*TEST
2181: test:
2182: suffix: 2d_quad_linear
2183: requires: triangle
2184: args: -sol_type quadratic_linear -dm_refine 2 \
2185: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2186: -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme
2188: test:
2189: suffix: 3d_quad_linear
2190: requires: ctetgen
2191: args: -dm_plex_dim 3 -sol_type quadratic_linear -dm_refine 1 \
2192: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2193: -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme
2195: test:
2196: suffix: 2d_trig_linear
2197: requires: triangle
2198: args: -sol_type trig_linear -dm_refine 1 \
2199: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2200: -dmts_check .0001 -ts_max_steps 5 -ts_time_step 0.00001 -ts_monitor_extreme
2202: test:
2203: # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9, 2.1, 1.8]
2204: suffix: 2d_trig_linear_sconv
2205: requires: triangle
2206: args: -sol_type trig_linear -dm_refine 1 \
2207: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2208: -convest_num_refine 1 -ts_convergence_estimate -ts_convergence_temporal 0 -ts_max_steps 1 -ts_time_step 0.00001 -pc_type lu
2210: test:
2211: suffix: 3d_trig_linear
2212: requires: ctetgen
2213: args: -dm_plex_dim 3 -sol_type trig_linear -dm_refine 1 \
2214: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2215: -dmts_check .0001 -ts_max_steps 2 -ts_monitor_extreme
2217: test:
2218: # -dm_refine 1 -convest_num_refine 2 gets L_2 convergence rate: [2.0, 2.1, 1.9]
2219: suffix: 3d_trig_linear_sconv
2220: requires: ctetgen
2221: args: -dm_plex_dim 3 -sol_type trig_linear -dm_refine 1 \
2222: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2223: -convest_num_refine 1 -ts_convergence_estimate -ts_convergence_temporal 0 -ts_max_steps 1 -pc_type lu
2225: test:
2226: suffix: 2d_quad_trig
2227: requires: triangle
2228: args: -sol_type quadratic_trig -dm_refine 2 \
2229: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2230: -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme
2232: test:
2233: # Using -dm_refine 4 gets the convergence rates to [0.95, 0.97, 0.90]
2234: suffix: 2d_quad_trig_tconv
2235: requires: triangle
2236: args: -sol_type quadratic_trig -dm_refine 1 \
2237: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2238: -convest_num_refine 3 -ts_convergence_estimate -ts_max_steps 5 -pc_type lu
2240: test:
2241: suffix: 3d_quad_trig
2242: requires: ctetgen
2243: args: -dm_plex_dim 3 -sol_type quadratic_trig -dm_refine 1 \
2244: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2245: -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme
2247: test:
2248: # Using -dm_refine 2 -convest_num_refine 3 gets the convergence rates to [1.0, 1.0, 1.0]
2249: suffix: 3d_quad_trig_tconv
2250: requires: ctetgen
2251: args: -dm_plex_dim 3 -sol_type quadratic_trig -dm_refine 1 \
2252: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2253: -convest_num_refine 1 -ts_convergence_estimate -ts_max_steps 5 -pc_type lu
2255: testset:
2256: 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 \
2257: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 -niter 16000 \
2258: -pc_type lu
2260: test:
2261: suffix: 2d_terzaghi
2262: requires: double
2263: args: -ts_time_step 0.0028666667 -ts_max_steps 2 -ts_monitor -dmts_check .0001
2265: test:
2266: # -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]
2267: suffix: 2d_terzaghi_tconv
2268: args: -ts_time_step 0.023 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1
2270: test:
2271: # -dm_plex_box_faces 1,16 -convest_num_refine 4 gives L_2 convergence rate: [1.7, 1.2, 1.1]
2272: # if we add -displacement_petscspace_degree 3 -tracestrain_petscspace_degree 2 -pressure_petscspace_degree 2, we get [2.1, 1.6, 1.5]
2273: suffix: 2d_terzaghi_sconv
2274: args: -ts_time_step 1e-5 -dt_initial 1e-5 -ts_max_steps 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1
2276: testset:
2277: 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 \
2278: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2279: -pc_type lu
2281: test:
2282: suffix: 2d_mandel
2283: requires: double
2284: args: -ts_time_step 0.0028666667 -ts_max_steps 2 -ts_monitor -dmts_check .0001
2286: test:
2287: # -dm_refine 3 -ts_max_steps 4 -convest_num_refine 3 gives L_2 convergence rate: [1.6, 0.93, 1.2]
2288: suffix: 2d_mandel_sconv
2289: args: -ts_time_step 1e-5 -dt_initial 1e-5 -ts_max_steps 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1
2291: test:
2292: # -dm_refine 5 -ts_max_steps 4 -convest_num_refine 3 gives L_2 convergence rate: [0.26, -0.0058, 0.26]
2293: suffix: 2d_mandel_tconv
2294: args: -ts_time_step 0.023 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1
2296: testset:
2297: requires: ctetgen !complex
2298: args: -sol_type cryer -dm_plex_dim 3 -dm_plex_shape ball \
2299: -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1
2301: test:
2302: suffix: 3d_cryer
2303: args: -ts_time_step 0.0028666667 -ts_max_time 0.014333 -ts_max_steps 2 -dmts_check .0001 \
2304: -pc_type svd
2306: test:
2307: # -bd_dm_refine 3 -dm_refine_volume_limit_pre 0.004 -convest_num_refine 2 gives L_2 convergence rate: []
2308: suffix: 3d_cryer_sconv
2309: args: -bd_dm_refine 1 -dm_refine_volume_limit_pre 0.00666667 \
2310: -ts_time_step 1e-5 -dt_initial 1e-5 -ts_max_steps 2 \
2311: -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
2312: -pc_type lu -pc_factor_shift_type nonzero
2314: test:
2315: # Displacement and Pressure converge. The analytic expression for trace strain is inaccurate at the origin
2316: # -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]
2317: suffix: 3d_cryer_tconv
2318: args: -bd_dm_refine 1 -dm_refine_volume_limit_pre 0.00666667 \
2319: -ts_time_step 0.023 -ts_max_time 0.092 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1 \
2320: -pc_type lu -pc_factor_shift_type nonzero
2322: TEST*/