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