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