Actual source code: ex53.c

  1: static char help[] = "Time dependent Biot Poroelasticity problem with finite elements.\n\
  2: We solve three field, quasi-static poroelasticity in a rectangular\n\
  3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\
  4: Contributed by: Robert Walker <rwalker6@buffalo.edu>\n\n\n";

  6: #include <petscdmplex.h>
  7: #include <petscsnes.h>
  8: #include <petscts.h>
  9: #include <petscds.h>
 10: #include <petscbag.h>

 12: #include <petsc/private/tsimpl.h>

 14: /* This presentation of poroelasticity is taken from

 16: @book{Cheng2016,
 17:   title={Poroelasticity},
 18:   author={Cheng, Alexander H-D},
 19:   volume={27},
 20:   year={2016},
 21:   publisher={Springer}
 22: }

 24: For visualization, use

 26:   -dm_view hdf5:${PETSC_DIR}/sol.h5 -monitor_solution hdf5:${PETSC_DIR}/sol.h5::append

 28: The weak form would then be, using test function $(v, q, \tau)$,

 30:             (q, \frac{1}{M} \frac{dp}{dt}) + (q, \alpha \frac{d\varepsilon}{dt}) + (\nabla q, \kappa \nabla p) = (q, g)
 31:  -(\nabla v, 2 G \epsilon) - (\nabla\cdot v, \frac{2 G \nu}{1 - 2\nu} \varepsilon) + \alpha (\nabla\cdot v, p) = (v, f)
 32:                                                                           (\tau, \nabla \cdot u - \varepsilon) = 0
 33: */

 35: typedef enum {
 36:   SOL_QUADRATIC_LINEAR,
 37:   SOL_QUADRATIC_TRIG,
 38:   SOL_TRIG_LINEAR,
 39:   SOL_TERZAGHI,
 40:   SOL_MANDEL,
 41:   SOL_CRYER,
 42:   NUM_SOLUTION_TYPES
 43: } SolutionType;
 44: const char *solutionTypes[NUM_SOLUTION_TYPES + 1] = {"quadratic_linear", "quadratic_trig", "trig_linear", "terzaghi", "mandel", "cryer", "unknown"};

 46: typedef struct {
 47:   PetscScalar mu;    /* shear modulus */
 48:   PetscScalar K_u;   /* undrained bulk modulus */
 49:   PetscScalar alpha; /* Biot effective stress coefficient */
 50:   PetscScalar M;     /* Biot modulus */
 51:   PetscScalar k;     /* (isotropic) permeability */
 52:   PetscScalar mu_f;  /* fluid dynamic viscosity */
 53:   PetscScalar P_0;   /* magnitude of vertical stress */
 54: } Parameter;

 56: typedef struct {
 57:   /* Domain and mesh definition */
 58:   PetscReal xmin[3]; /* Lower left bottom corner of bounding box */
 59:   PetscReal xmax[3]; /* Upper right top corner of bounding box */
 60:   /* Problem definition */
 61:   SolutionType solType;   /* Type of exact solution */
 62:   PetscBag     bag;       /* Problem parameters */
 63:   PetscReal    t_r;       /* Relaxation time: 4 L^2 / c */
 64:   PetscReal    dtInitial; /* Override the choice for first timestep */
 65:   /* Exact solution terms */
 66:   PetscInt   niter;     /* Number of series term iterations in exact solutions */
 67:   PetscReal  eps;       /* Precision value for root finding */
 68:   PetscReal *zeroArray; /* Array of root locations */
 69: } AppCtx;

 71: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 72: {
 73:   PetscInt c;
 74:   for (c = 0; c < Nc; ++c) u[c] = 0.0;
 75:   return PETSC_SUCCESS;
 76: }

 78: /* Quadratic space and linear time solution

 80:   2D:
 81:   u = x^2
 82:   v = y^2 - 2xy
 83:   p = (x + y) t
 84:   e = 2y
 85:   f = <2 G, 4 G + 2 \lambda > - <alpha t, alpha t>
 86:   g = 0
 87:   \epsilon = / 2x     -y    \
 88:              \ -y   2y - 2x /
 89:   Tr(\epsilon) = e = div u = 2y
 90:   div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
 91:     = 2 G < 2-1, 2 > + \lambda <0, 2> - alpha <t, t>
 92:     = <2 G, 4 G + 2 \lambda> - <alpha t, alpha t>
 93:   \frac{1}{M} \frac{dp}{dt} + \alpha \frac{d\varepsilon}{dt} - \nabla \cdot \kappa \nabla p
 94:     = \frac{1}{M} \frac{dp}{dt} + \kappa \Delta p
 95:     = (x + y)/M

 97:   3D:
 98:   u = x^2
 99:   v = y^2 - 2xy
100:   w = z^2 - 2yz
101:   p = (x + y + z) t
102:   e = 2z
103:   f = <2 G, 4 G + 2 \lambda > - <alpha t, alpha t, alpha t>
104:   g = 0
105:   \varepsilon = / 2x     -y       0   \
106:                 | -y   2y - 2x   -z   |
107:                 \  0     -z    2z - 2y/
108:   Tr(\varepsilon) = div u = 2z
109:   div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
110:     = 2 G < 2-1, 2-1, 2 > + \lambda <0, 0, 2> - alpha <t, t, t>
111:     = <2 G, 2G, 4 G + 2 \lambda> - <alpha t, alpha t, alpha t>
112: */
113: static PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
114: {
115:   PetscInt d;

117:   for (d = 0; d < dim; ++d) u[d] = PetscSqr(x[d]) - (d > 0 ? 2.0 * x[d - 1] * x[d] : 0.0);
118:   return PETSC_SUCCESS;
119: }

121: static PetscErrorCode linear_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
122: {
123:   u[0] = 2.0 * x[dim - 1];
124:   return PETSC_SUCCESS;
125: }

127: static PetscErrorCode linear_linear_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
128: {
129:   PetscReal sum = 0.0;
130:   PetscInt  d;

132:   for (d = 0; d < dim; ++d) sum += x[d];
133:   u[0] = sum * time;
134:   return PETSC_SUCCESS;
135: }

137: static PetscErrorCode linear_linear_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
138: {
139:   PetscReal sum = 0.0;
140:   PetscInt  d;

142:   for (d = 0; d < dim; ++d) sum += x[d];
143:   u[0] = sum;
144:   return PETSC_SUCCESS;
145: }

147: static void f0_quadratic_linear_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
148: {
149:   const PetscReal G      = PetscRealPart(constants[0]);
150:   const PetscReal K_u    = PetscRealPart(constants[1]);
151:   const PetscReal alpha  = PetscRealPart(constants[2]);
152:   const PetscReal M      = PetscRealPart(constants[3]);
153:   const PetscReal K_d    = K_u - alpha * alpha * M;
154:   const PetscReal lambda = K_d - (2.0 * G) / 3.0;
155:   PetscInt        d;

157:   for (d = 0; d < dim - 1; ++d) f0[d] -= 2.0 * G - alpha * t;
158:   f0[dim - 1] -= 2.0 * lambda + 4.0 * G - alpha * t;
159: }

161: static void f0_quadratic_linear_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
162: {
163:   const PetscReal alpha = PetscRealPart(constants[2]);
164:   const PetscReal M     = PetscRealPart(constants[3]);
165:   PetscReal       sum   = 0.0;
166:   PetscInt        d;

168:   for (d = 0; d < dim; ++d) sum += x[d];
169:   f0[0] += u_t ? alpha * u_t[uOff[1]] : 0.0;
170:   f0[0] += u_t ? u_t[uOff[2]] / M : 0.0;
171:   f0[0] -= sum / M;
172: }

174: /* Quadratic space and trigonometric time solution

176:   2D:
177:   u = x^2
178:   v = y^2 - 2xy
179:   p = (x + y) cos(t)
180:   e = 2y
181:   f = <2 G, 4 G + 2 \lambda > - <alpha cos(t), alpha cos(t)>
182:   g = 0
183:   \epsilon = / 2x     -y    \
184:              \ -y   2y - 2x /
185:   Tr(\epsilon) = e = div u = 2y
186:   div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
187:     = 2 G < 2-1, 2 > + \lambda <0, 2> - alpha <cos(t), cos(t)>
188:     = <2 G, 4 G + 2 \lambda> - <alpha cos(t), alpha cos(t)>
189:   \frac{1}{M} \frac{dp}{dt} + \alpha \frac{d\varepsilon}{dt} - \nabla \cdot \kappa \nabla p
190:     = \frac{1}{M} \frac{dp}{dt} + \kappa \Delta p
191:     = -(x + y)/M sin(t)

193:   3D:
194:   u = x^2
195:   v = y^2 - 2xy
196:   w = z^2 - 2yz
197:   p = (x + y + z) cos(t)
198:   e = 2z
199:   f = <2 G, 4 G + 2 \lambda > - <alpha cos(t), alpha cos(t), alpha cos(t)>
200:   g = 0
201:   \varepsilon = / 2x     -y       0   \
202:                 | -y   2y - 2x   -z   |
203:                 \  0     -z    2z - 2y/
204:   Tr(\varepsilon) = div u = 2z
205:   div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
206:     = 2 G < 2-1, 2-1, 2 > + \lambda <0, 0, 2> - alpha <cos(t), cos(t), cos(t)>
207:     = <2 G, 2G, 4 G + 2 \lambda> - <alpha cos(t), alpha cos(t), alpha cos(t)>
208: */
209: static PetscErrorCode linear_trig_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
210: {
211:   PetscReal sum = 0.0;
212:   PetscInt  d;

214:   for (d = 0; d < dim; ++d) sum += x[d];
215:   u[0] = sum * PetscCosReal(time);
216:   return PETSC_SUCCESS;
217: }

219: static PetscErrorCode linear_trig_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
220: {
221:   PetscReal sum = 0.0;
222:   PetscInt  d;

224:   for (d = 0; d < dim; ++d) sum += x[d];
225:   u[0] = -sum * PetscSinReal(time);
226:   return PETSC_SUCCESS;
227: }

229: static void f0_quadratic_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
230: {
231:   const PetscReal G      = PetscRealPart(constants[0]);
232:   const PetscReal K_u    = PetscRealPart(constants[1]);
233:   const PetscReal alpha  = PetscRealPart(constants[2]);
234:   const PetscReal M      = PetscRealPart(constants[3]);
235:   const PetscReal K_d    = K_u - alpha * alpha * M;
236:   const PetscReal lambda = K_d - (2.0 * G) / 3.0;
237:   PetscInt        d;

239:   for (d = 0; d < dim - 1; ++d) f0[d] -= 2.0 * G - alpha * PetscCosReal(t);
240:   f0[dim - 1] -= 2.0 * lambda + 4.0 * G - alpha * PetscCosReal(t);
241: }

243: static void f0_quadratic_trig_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
244: {
245:   const PetscReal alpha = PetscRealPart(constants[2]);
246:   const PetscReal M     = PetscRealPart(constants[3]);
247:   PetscReal       sum   = 0.0;
248:   PetscInt        d;

250:   for (d = 0; d < dim; ++d) sum += x[d];

252:   f0[0] += u_t ? alpha * u_t[uOff[1]] : 0.0;
253:   f0[0] += u_t ? u_t[uOff[2]] / M : 0.0;
254:   f0[0] += PetscSinReal(t) * sum / M;
255: }

257: /* Trigonometric space and linear time solution

259: u = sin(2 pi x)
260: v = sin(2 pi y) - 2xy
261: \varepsilon = / 2 pi cos(2 pi x)             -y        \
262:               \      -y          2 pi cos(2 pi y) - 2x /
263: Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
264: div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
265:   = \lambda \partial_j 2 pi (cos(2 pi x) + cos(2 pi y)) + 2\mu < -4 pi^2 sin(2 pi x) - 1, -4 pi^2 sin(2 pi y) >
266:   = \lambda < -4 pi^2 sin(2 pi x) - 2, -4 pi^2 sin(2 pi y) > + \mu < -8 pi^2 sin(2 pi x) - 2, -8 pi^2 sin(2 pi y) >

268:   2D:
269:   u = sin(2 pi x)
270:   v = sin(2 pi y) - 2xy
271:   p = (cos(2 pi x) + cos(2 pi y)) t
272:   e = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
273:   f = < -4 pi^2 sin(2 pi x) (2 G + lambda) - (2 G - 2 lambda), -4 pi^2 sin(2 pi y) (2G + lambda) > + 2 pi alpha t <sin(2 pi x), sin(2 pi y)>
274:   g = 0
275:   \varepsilon = / 2 pi cos(2 pi x)             -y        \
276:                 \      -y          2 pi cos(2 pi y) - 2x /
277:   Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
278:   div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
279:     = 2 G < -4 pi^2 sin(2 pi x) - 1, -4 pi^2 sin(2 pi y) > + \lambda <-4 pi^2 sin(2 pi x) - 2, -4 pi^2 sin(2 pi y)> - alpha <-2 pi sin(2 pi x) t, -2 pi sin(2 pi y) t>
280:     = < -4 pi^2 sin(2 pi x) (2 G + lambda) - (2 G + 2 lambda), -4 pi^2 sin(2 pi y) (2G + lambda) > + 2 pi alpha t <sin(2 pi x), sin(2 pi y)>
281:   \frac{1}{M} \frac{dp}{dt} + \alpha \frac{d\varepsilon}{dt} - \nabla \cdot \kappa \nabla p
282:     = \frac{1}{M} \frac{dp}{dt} + \kappa \Delta p
283:     = (cos(2 pi x) + cos(2 pi y))/M - 4 pi^2 \kappa (cos(2 pi x) + cos(2 pi y)) t

285:   3D:
286:   u = sin(2 pi x)
287:   v = sin(2 pi y) - 2xy
288:   v = sin(2 pi y) - 2yz
289:   p = (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) t
290:   e = 2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2y
291:   f = < -4 pi^2 sin(2 pi x) (2 G + lambda) - (2 G + 2 lambda),  -4 pi^2 sin(2 pi y) (2 G + lambda) - (2 G + 2 lambda), -4 pi^2 sin(2 pi z) (2G + lambda) > + 2 pi alpha t <sin(2 pi x), sin(2 pi y), , sin(2 pi z)>
292:   g = 0
293:   \varepsilon = / 2 pi cos(2 pi x)            -y                     0         \
294:                 |         -y       2 pi cos(2 pi y) - 2x            -z         |
295:                 \          0                  -z         2 pi cos(2 pi z) - 2y /
296:   Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2 y
297:   div \sigma = \partial_i 2 G \epsilon_{ij} + \partial_i \lambda \varepsilon \delta_{ij} - \partial_i \alpha p \delta_{ij}
298:     = 2 G < -4 pi^2 sin(2 pi x) - 1, -4 pi^2 sin(2 pi y) - 1, -4 pi^2 sin(2 pi z) > + \lambda <-4 pi^2 sin(2 pi x) - 2, 4 pi^2 sin(2 pi y) - 2, -4 pi^2 sin(2 pi y)> - alpha <-2 pi sin(2 pi x) t, -2 pi sin(2 pi y) t, -2 pi sin(2 pi z) t>
299:     = < -4 pi^2 sin(2 pi x) (2 G + lambda) - (2 G + 2 lambda),  -4 pi^2 sin(2 pi y) (2 G + lambda) - (2 G + 2 lambda), -4 pi^2 sin(2 pi z) (2G + lambda) > + 2 pi alpha t <sin(2 pi x), sin(2 pi y), , sin(2 pi z)>
300:   \frac{1}{M} \frac{dp}{dt} + \alpha \frac{d\varepsilon}{dt} - \nabla \cdot \kappa \nabla p
301:     = \frac{1}{M} \frac{dp}{dt} + \kappa \Delta p
302:     = (cos(2 pi x) + cos(2 pi y) + cos(2 pi z))/M - 4 pi^2 \kappa (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) t
303: */
304: static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
305: {
306:   PetscInt d;

308:   for (d = 0; d < dim; ++d) u[d] = PetscSinReal(2. * PETSC_PI * x[d]) - (d > 0 ? 2.0 * x[d - 1] * x[d] : 0.0);
309:   return PETSC_SUCCESS;
310: }

312: static PetscErrorCode trig_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
313: {
314:   PetscReal sum = 0.0;
315:   PetscInt  d;

317:   for (d = 0; d < dim; ++d) sum += 2. * PETSC_PI * PetscCosReal(2. * PETSC_PI * x[d]) - (d < dim - 1 ? 2. * x[d] : 0.0);
318:   u[0] = sum;
319:   return PETSC_SUCCESS;
320: }

322: static PetscErrorCode trig_linear_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
323: {
324:   PetscReal sum = 0.0;
325:   PetscInt  d;

327:   for (d = 0; d < dim; ++d) sum += PetscCosReal(2. * PETSC_PI * x[d]);
328:   u[0] = sum * time;
329:   return PETSC_SUCCESS;
330: }

332: static PetscErrorCode trig_linear_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
333: {
334:   PetscReal sum = 0.0;
335:   PetscInt  d;

337:   for (d = 0; d < dim; ++d) sum += PetscCosReal(2. * PETSC_PI * x[d]);
338:   u[0] = sum;
339:   return PETSC_SUCCESS;
340: }

342: static void f0_trig_linear_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
343: {
344:   const PetscReal G      = PetscRealPart(constants[0]);
345:   const PetscReal K_u    = PetscRealPart(constants[1]);
346:   const PetscReal alpha  = PetscRealPart(constants[2]);
347:   const PetscReal M      = PetscRealPart(constants[3]);
348:   const PetscReal K_d    = K_u - alpha * alpha * M;
349:   const PetscReal lambda = K_d - (2.0 * G) / 3.0;
350:   PetscInt        d;

352:   for (d = 0; d < dim - 1; ++d) f0[d] += PetscSqr(2. * PETSC_PI) * PetscSinReal(2. * PETSC_PI * x[d]) * (2. * G + lambda) + 2.0 * (G + lambda) - 2. * PETSC_PI * alpha * PetscSinReal(2. * PETSC_PI * x[d]) * t;
353:   f0[dim - 1] += PetscSqr(2. * PETSC_PI) * PetscSinReal(2. * PETSC_PI * x[dim - 1]) * (2. * G + lambda) - 2. * PETSC_PI * alpha * PetscSinReal(2. * PETSC_PI * x[dim - 1]) * t;
354: }

356: static void f0_trig_linear_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
357: {
358:   const PetscReal alpha = PetscRealPart(constants[2]);
359:   const PetscReal M     = PetscRealPart(constants[3]);
360:   const PetscReal kappa = PetscRealPart(constants[4]);
361:   PetscReal       sum   = 0.0;
362:   PetscInt        d;

364:   for (d = 0; d < dim; ++d) sum += PetscCosReal(2. * PETSC_PI * x[d]);
365:   f0[0] += u_t ? alpha * u_t[uOff[1]] : 0.0;
366:   f0[0] += u_t ? u_t[uOff[2]] / M : 0.0;
367:   f0[0] -= sum / M - 4 * PetscSqr(PETSC_PI) * kappa * sum * t;
368: }

370: /* Terzaghi Solutions */
371: /* The analytical solutions given here are drawn from chapter 7, section 3, */
372: /* "One-Dimensional Consolidation Problem," from Poroelasticity, by Cheng.  */
373: static PetscErrorCode terzaghi_drainage_pressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
374: {
375:   AppCtx    *user = (AppCtx *)ctx;
376:   Parameter *param;

378:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
379:   if (time <= 0.0) {
380:     PetscScalar alpha = param->alpha;                                        /* -  */
381:     PetscScalar K_u   = param->K_u;                                          /* Pa */
382:     PetscScalar M     = param->M;                                            /* Pa */
383:     PetscScalar G     = param->mu;                                           /* Pa */
384:     PetscScalar P_0   = param->P_0;                                          /* Pa */
385:     PetscScalar K_d   = K_u - alpha * alpha * M;                             /* Pa,      Cheng (B.5)  */
386:     PetscScalar eta   = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G);           /* -,       Cheng (B.11) */
387:     PetscScalar S     = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */

389:     u[0] = ((P_0 * eta) / (G * S));
390:   } else {
391:     u[0] = 0.0;
392:   }
393:   return PETSC_SUCCESS;
394: }

396: static PetscErrorCode terzaghi_initial_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
397: {
398:   AppCtx    *user = (AppCtx *)ctx;
399:   Parameter *param;

401:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
402:   {
403:     PetscScalar K_u   = param->K_u;                                      /* Pa */
404:     PetscScalar G     = param->mu;                                       /* Pa */
405:     PetscScalar P_0   = param->P_0;                                      /* Pa */
406:     PetscReal   L     = user->xmax[1] - user->xmin[1];                   /* m */
407:     PetscScalar nu_u  = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -,       Cheng (B.9)  */
408:     PetscReal   zstar = x[1] / L;                                        /* - */

410:     u[0] = 0.0;
411:     u[1] = ((P_0 * L * (1.0 - 2.0 * nu_u)) / (2.0 * G * (1.0 - nu_u))) * (1.0 - zstar);
412:   }
413:   return PETSC_SUCCESS;
414: }

416: static PetscErrorCode terzaghi_initial_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
417: {
418:   AppCtx    *user = (AppCtx *)ctx;
419:   Parameter *param;

421:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
422:   {
423:     PetscScalar K_u  = param->K_u;                                      /* Pa */
424:     PetscScalar G    = param->mu;                                       /* Pa */
425:     PetscScalar P_0  = param->P_0;                                      /* Pa */
426:     PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -,       Cheng (B.9)  */

428:     u[0] = -(P_0 * (1.0 - 2.0 * nu_u)) / (2.0 * G * (1.0 - nu_u));
429:   }
430:   return PETSC_SUCCESS;
431: }

433: static PetscErrorCode terzaghi_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
434: {
435:   AppCtx    *user = (AppCtx *)ctx;
436:   Parameter *param;

438:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
439:   if (time < 0.0) {
440:     PetscCall(terzaghi_initial_u(dim, time, x, Nc, u, ctx));
441:   } else {
442:     PetscScalar alpha = param->alpha;                  /* -  */
443:     PetscScalar K_u   = param->K_u;                    /* Pa */
444:     PetscScalar M     = param->M;                      /* Pa */
445:     PetscScalar G     = param->mu;                     /* Pa */
446:     PetscScalar P_0   = param->P_0;                    /* Pa */
447:     PetscScalar kappa = param->k / param->mu_f;        /* m^2 / (Pa s) */
448:     PetscReal   L     = user->xmax[1] - user->xmin[1]; /* m */
449:     PetscInt    N     = user->niter, m;

451:     PetscScalar K_d  = K_u - alpha * alpha * M;                             /* Pa,      Cheng (B.5)  */
452:     PetscScalar nu   = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));     /* -,       Cheng (B.8)  */
453:     PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));     /* -,       Cheng (B.9)  */
454:     PetscScalar S    = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
455:     PetscScalar c    = kappa / S;                                           /* m^2 / s, Cheng (B.16) */

457:     PetscReal   zstar = x[1] / L;                                    /* - */
458:     PetscReal   tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
459:     PetscScalar F2    = 0.0;

461:     for (m = 1; m < 2 * N + 1; ++m) {
462:       if (m % 2 == 1) F2 += (8.0 / PetscSqr(m * PETSC_PI)) * PetscCosReal(0.5 * m * PETSC_PI * zstar) * (1.0 - PetscExpReal(-PetscSqr(m * PETSC_PI) * tstar));
463:     }
464:     u[0] = 0.0;
465:     u[1] = ((P_0 * L * (1.0 - 2.0 * nu_u)) / (2.0 * G * (1.0 - nu_u))) * (1.0 - zstar) + ((P_0 * L * (nu_u - nu)) / (2.0 * G * (1.0 - nu_u) * (1.0 - nu))) * F2; /* m */
466:   }
467:   return PETSC_SUCCESS;
468: }

470: static PetscErrorCode terzaghi_2d_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
471: {
472:   AppCtx    *user = (AppCtx *)ctx;
473:   Parameter *param;

475:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
476:   if (time < 0.0) {
477:     PetscCall(terzaghi_initial_eps(dim, time, x, Nc, u, ctx));
478:   } else {
479:     PetscScalar alpha = param->alpha;                  /* -  */
480:     PetscScalar K_u   = param->K_u;                    /* Pa */
481:     PetscScalar M     = param->M;                      /* Pa */
482:     PetscScalar G     = param->mu;                     /* Pa */
483:     PetscScalar P_0   = param->P_0;                    /* Pa */
484:     PetscScalar kappa = param->k / param->mu_f;        /* m^2 / (Pa s) */
485:     PetscReal   L     = user->xmax[1] - user->xmin[1]; /* m */
486:     PetscInt    N     = user->niter, m;

488:     PetscScalar K_d  = K_u - alpha * alpha * M;                             /* Pa,      Cheng (B.5)  */
489:     PetscScalar nu   = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));     /* -,       Cheng (B.8)  */
490:     PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));     /* -,       Cheng (B.9)  */
491:     PetscScalar S    = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
492:     PetscScalar c    = kappa / S;                                           /* m^2 / s, Cheng (B.16) */

494:     PetscReal   zstar = x[1] / L;                                    /* - */
495:     PetscReal   tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
496:     PetscScalar F2_z  = 0.0;

498:     for (m = 1; m < 2 * N + 1; ++m) {
499:       if (m % 2 == 1) F2_z += (-4.0 / (m * PETSC_PI * L)) * PetscSinReal(0.5 * m * PETSC_PI * zstar) * (1.0 - PetscExpReal(-PetscSqr(m * PETSC_PI) * tstar));
500:     }
501:     u[0] = -((P_0 * L * (1.0 - 2.0 * nu_u)) / (2.0 * G * (1.0 - nu_u) * L)) + ((P_0 * L * (nu_u - nu)) / (2.0 * G * (1.0 - nu_u) * (1.0 - nu))) * F2_z; /* - */
502:   }
503:   return PETSC_SUCCESS;
504: }

506: // Pressure
507: static PetscErrorCode terzaghi_2d_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
508: {
509:   AppCtx    *user = (AppCtx *)ctx;
510:   Parameter *param;

512:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
513:   if (time <= 0.0) {
514:     PetscCall(terzaghi_drainage_pressure(dim, time, x, Nc, u, ctx));
515:   } else {
516:     PetscScalar alpha = param->alpha;                  /* -  */
517:     PetscScalar K_u   = param->K_u;                    /* Pa */
518:     PetscScalar M     = param->M;                      /* Pa */
519:     PetscScalar G     = param->mu;                     /* Pa */
520:     PetscScalar P_0   = param->P_0;                    /* Pa */
521:     PetscScalar kappa = param->k / param->mu_f;        /* m^2 / (Pa s) */
522:     PetscReal   L     = user->xmax[1] - user->xmin[1]; /* m */
523:     PetscInt    N     = user->niter, m;

525:     PetscScalar K_d = K_u - alpha * alpha * M;                             /* Pa,      Cheng (B.5)  */
526:     PetscScalar eta = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G);           /* -,       Cheng (B.11) */
527:     PetscScalar S   = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
528:     PetscScalar c   = kappa / S;                                           /* m^2 / s, Cheng (B.16) */

530:     PetscReal   zstar = x[1] / L;                                    /* - */
531:     PetscReal   tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
532:     PetscScalar F1    = 0.0;

534:     PetscCheck(PetscAbsScalar((1 / M + (alpha * eta) / G) - S) <= 1.0e-10, PETSC_COMM_SELF, PETSC_ERR_PLIB, "S %g != check %g", (double)PetscAbsScalar(S), (double)PetscAbsScalar(1 / M + (alpha * eta) / G));

536:     for (m = 1; m < 2 * N + 1; ++m) {
537:       if (m % 2 == 1) F1 += (4.0 / (m * PETSC_PI)) * PetscSinReal(0.5 * m * PETSC_PI * zstar) * PetscExpReal(-PetscSqr(m * PETSC_PI) * tstar);
538:     }
539:     u[0] = ((P_0 * eta) / (G * S)) * F1; /* Pa */
540:   }
541:   return PETSC_SUCCESS;
542: }

544: static PetscErrorCode terzaghi_2d_u_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
545: {
546:   AppCtx    *user = (AppCtx *)ctx;
547:   Parameter *param;

549:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
550:   if (time <= 0.0) {
551:     u[0] = 0.0;
552:     u[1] = 0.0;
553:   } else {
554:     PetscScalar alpha = param->alpha;                  /* -  */
555:     PetscScalar K_u   = param->K_u;                    /* Pa */
556:     PetscScalar M     = param->M;                      /* Pa */
557:     PetscScalar G     = param->mu;                     /* Pa */
558:     PetscScalar P_0   = param->P_0;                    /* Pa */
559:     PetscScalar kappa = param->k / param->mu_f;        /* m^2 / (Pa s) */
560:     PetscReal   L     = user->xmax[1] - user->xmin[1]; /* m */
561:     PetscInt    N     = user->niter, m;

563:     PetscScalar K_d  = K_u - alpha * alpha * M;                             /* Pa,      Cheng (B.5)  */
564:     PetscScalar nu   = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));     /* -,       Cheng (B.8)  */
565:     PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));     /* -,       Cheng (B.9)  */
566:     PetscScalar S    = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
567:     PetscScalar c    = kappa / S;                                           /* m^2 / s, Cheng (B.16) */

569:     PetscReal   zstar = x[1] / L;                                    /* - */
570:     PetscReal   tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
571:     PetscScalar F2_t  = 0.0;

573:     for (m = 1; m < 2 * N + 1; ++m) {
574:       if (m % 2 == 1) F2_t += (2.0 * c / PetscSqr(L)) * PetscCosReal(0.5 * m * PETSC_PI * zstar) * PetscExpReal(-PetscSqr(m * PETSC_PI) * tstar);
575:     }
576:     u[0] = 0.0;
577:     u[1] = ((P_0 * L * (nu_u - nu)) / (2.0 * G * (1.0 - nu_u) * (1.0 - nu))) * F2_t; /* m / s */
578:   }
579:   return PETSC_SUCCESS;
580: }

582: static PetscErrorCode terzaghi_2d_eps_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
583: {
584:   AppCtx    *user = (AppCtx *)ctx;
585:   Parameter *param;

587:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
588:   if (time <= 0.0) {
589:     u[0] = 0.0;
590:   } else {
591:     PetscScalar alpha = param->alpha;                  /* -  */
592:     PetscScalar K_u   = param->K_u;                    /* Pa */
593:     PetscScalar M     = param->M;                      /* Pa */
594:     PetscScalar G     = param->mu;                     /* Pa */
595:     PetscScalar P_0   = param->P_0;                    /* Pa */
596:     PetscScalar kappa = param->k / param->mu_f;        /* m^2 / (Pa s) */
597:     PetscReal   L     = user->xmax[1] - user->xmin[1]; /* m */
598:     PetscInt    N     = user->niter, m;

600:     PetscScalar K_d  = K_u - alpha * alpha * M;                             /* Pa,      Cheng (B.5)  */
601:     PetscScalar nu   = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));     /* -,       Cheng (B.8)  */
602:     PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));     /* -,       Cheng (B.9)  */
603:     PetscScalar S    = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
604:     PetscScalar c    = kappa / S;                                           /* m^2 / s, Cheng (B.16) */

606:     PetscReal   zstar = x[1] / L;                                    /* - */
607:     PetscReal   tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
608:     PetscScalar F2_zt = 0.0;

610:     for (m = 1; m < 2 * N + 1; ++m) {
611:       if (m % 2 == 1) F2_zt += ((-m * PETSC_PI * c) / (L * L * L)) * PetscSinReal(0.5 * m * PETSC_PI * zstar) * PetscExpReal(-PetscSqr(m * PETSC_PI) * tstar);
612:     }
613:     u[0] = ((P_0 * L * (nu_u - nu)) / (2.0 * G * (1.0 - nu_u) * (1.0 - nu))) * F2_zt; /* 1 / s */
614:   }
615:   return PETSC_SUCCESS;
616: }

618: static PetscErrorCode terzaghi_2d_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
619: {
620:   AppCtx    *user = (AppCtx *)ctx;
621:   Parameter *param;

623:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
624:   if (time <= 0.0) {
625:     PetscScalar alpha = param->alpha;                  /* -  */
626:     PetscScalar K_u   = param->K_u;                    /* Pa */
627:     PetscScalar M     = param->M;                      /* Pa */
628:     PetscScalar G     = param->mu;                     /* Pa */
629:     PetscScalar P_0   = param->P_0;                    /* Pa */
630:     PetscScalar kappa = param->k / param->mu_f;        /* m^2 / (Pa s) */
631:     PetscReal   L     = user->xmax[1] - user->xmin[1]; /* m */

633:     PetscScalar K_d = K_u - alpha * alpha * M;                             /* Pa,      Cheng (B.5)  */
634:     PetscScalar eta = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G);           /* -,       Cheng (B.11) */
635:     PetscScalar S   = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
636:     PetscScalar c   = kappa / S;                                           /* m^2 / s, Cheng (B.16) */

638:     u[0] = -((P_0 * eta) / (G * S)) * PetscSqr(0 * PETSC_PI) * c / PetscSqr(2.0 * L); /* Pa / s */
639:   } else {
640:     PetscScalar alpha = param->alpha;                  /* -  */
641:     PetscScalar K_u   = param->K_u;                    /* Pa */
642:     PetscScalar M     = param->M;                      /* Pa */
643:     PetscScalar G     = param->mu;                     /* Pa */
644:     PetscScalar P_0   = param->P_0;                    /* Pa */
645:     PetscScalar kappa = param->k / param->mu_f;        /* m^2 / (Pa s) */
646:     PetscReal   L     = user->xmax[1] - user->xmin[1]; /* m */
647:     PetscInt    N     = user->niter, m;

649:     PetscScalar K_d = K_u - alpha * alpha * M;                             /* Pa,      Cheng (B.5)  */
650:     PetscScalar eta = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G);           /* -,       Cheng (B.11) */
651:     PetscScalar S   = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
652:     PetscScalar c   = kappa / S;                                           /* m^2 / s, Cheng (B.16) */

654:     PetscReal   zstar = x[1] / L;                                    /* - */
655:     PetscReal   tstar = PetscRealPart(c * time) / PetscSqr(2.0 * L); /* - */
656:     PetscScalar F1_t  = 0.0;

658:     PetscCheck(PetscAbsScalar((1 / M + (alpha * eta) / G) - S) <= 1.0e-10, PETSC_COMM_SELF, PETSC_ERR_PLIB, "S %g != check %g", (double)PetscAbsScalar(S), (double)PetscAbsScalar(1 / M + (alpha * eta) / G));

660:     for (m = 1; m < 2 * N + 1; ++m) {
661:       if (m % 2 == 1) F1_t += ((-m * PETSC_PI * c) / PetscSqr(L)) * PetscSinReal(0.5 * m * PETSC_PI * zstar) * PetscExpReal(-PetscSqr(m * PETSC_PI) * tstar);
662:     }
663:     u[0] = ((P_0 * eta) / (G * S)) * F1_t; /* Pa / s */
664:   }
665:   return PETSC_SUCCESS;
666: }

668: /* Mandel Solutions */
669: static PetscErrorCode mandel_drainage_pressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
670: {
671:   AppCtx    *user = (AppCtx *)ctx;
672:   Parameter *param;

674:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
675:   if (time <= 0.0) {
676:     PetscScalar alpha = param->alpha;                          /* -  */
677:     PetscScalar K_u   = param->K_u;                            /* Pa */
678:     PetscScalar M     = param->M;                              /* Pa */
679:     PetscScalar G     = param->mu;                             /* Pa */
680:     PetscScalar P_0   = param->P_0;                            /* Pa */
681:     PetscScalar kappa = param->k / param->mu_f;                /* m^2 / (Pa s) */
682:     PetscReal   a     = 0.5 * (user->xmax[0] - user->xmin[0]); /* m */
683:     PetscInt    N     = user->niter, n;

685:     PetscScalar K_d  = K_u - alpha * alpha * M;                             /* Pa,      Cheng (B.5)  */
686:     PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));     /* -,       Cheng (B.9)  */
687:     PetscScalar B    = alpha * M / K_u;                                     /* -,       Cheng (B.12) */
688:     PetscScalar S    = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
689:     PetscScalar c    = kappa / S;                                           /* m^2 / s, Cheng (B.16) */

691:     PetscScalar A1   = 3.0 / (B * (1.0 + nu_u));
692:     PetscReal   aa   = 0.0;
693:     PetscReal   p    = 0.0;
694:     PetscReal   time = 0.0;

696:     for (n = 1; n < N + 1; ++n) {
697:       aa = user->zeroArray[n - 1];
698:       p += (PetscSinReal(aa) / (aa - PetscSinReal(aa) * PetscCosReal(aa))) * (PetscCosReal((aa * x[0]) / a) - PetscCosReal(aa)) * PetscExpReal(-1.0 * (aa * aa * PetscRealPart(c) * time) / (a * a));
699:     }
700:     u[0] = ((2.0 * P_0) / (a * A1)) * p;
701:   } else {
702:     u[0] = 0.0;
703:   }
704:   return PETSC_SUCCESS;
705: }

707: static PetscErrorCode mandel_initial_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
708: {
709:   AppCtx    *user = (AppCtx *)ctx;
710:   Parameter *param;

712:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
713:   {
714:     PetscScalar alpha = param->alpha;                          /* -  */
715:     PetscScalar K_u   = param->K_u;                            /* Pa */
716:     PetscScalar M     = param->M;                              /* Pa */
717:     PetscScalar G     = param->mu;                             /* Pa */
718:     PetscScalar P_0   = param->P_0;                            /* Pa */
719:     PetscScalar kappa = param->k / param->mu_f;                /* m^2 / (Pa s) */
720:     PetscReal   a     = 0.5 * (user->xmax[0] - user->xmin[0]); /* m */
721:     PetscInt    N     = user->niter, n;

723:     PetscScalar K_d  = K_u - alpha * alpha * M;                             /* Pa,      Cheng (B.5)  */
724:     PetscScalar nu   = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));     /* -,       Cheng (B.8)  */
725:     PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));     /* -,       Cheng (B.9)  */
726:     PetscScalar S    = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
727:     PetscReal   c    = PetscRealPart(kappa / S);                            /* m^2 / s, Cheng (B.16) */

729:     PetscReal A_s = 0.0;
730:     PetscReal B_s = 0.0;
731:     for (n = 1; n < N + 1; ++n) {
732:       PetscReal alpha_n = user->zeroArray[n - 1];
733:       A_s += ((PetscSinReal(alpha_n) * PetscCosReal(alpha_n)) / (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscExpReal(-1 * (alpha_n * alpha_n * c * time) / (a * a));
734:       B_s += (PetscCosReal(alpha_n) / (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscSinReal((alpha_n * x[0]) / a) * PetscExpReal(-1 * (alpha_n * alpha_n * c * time) / (a * a));
735:     }
736:     u[0] = ((P_0 * nu) / (2.0 * G * a) - (P_0 * nu_u) / (G * a) * A_s) * x[0] + P_0 / G * B_s;
737:     u[1] = (-1 * (P_0 * (1.0 - nu)) / (2 * G * a) + (P_0 * (1 - nu_u)) / (G * a) * A_s) * x[1];
738:   }
739:   return PETSC_SUCCESS;
740: }

742: static PetscErrorCode mandel_initial_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
743: {
744:   AppCtx    *user = (AppCtx *)ctx;
745:   Parameter *param;

747:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
748:   {
749:     PetscScalar alpha = param->alpha;                          /* -  */
750:     PetscScalar K_u   = param->K_u;                            /* Pa */
751:     PetscScalar M     = param->M;                              /* Pa */
752:     PetscScalar G     = param->mu;                             /* Pa */
753:     PetscScalar P_0   = param->P_0;                            /* Pa */
754:     PetscScalar kappa = param->k / param->mu_f;                /* m^2 / (Pa s) */
755:     PetscReal   a     = 0.5 * (user->xmax[0] - user->xmin[0]); /* m */
756:     PetscInt    N     = user->niter, n;

758:     PetscScalar K_d = K_u - alpha * alpha * M;                             /* Pa,      Cheng (B.5)  */
759:     PetscScalar nu  = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));     /* -,       Cheng (B.8)  */
760:     PetscScalar S   = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
761:     PetscReal   c   = PetscRealPart(kappa / S);                            /* m^2 / s, Cheng (B.16) */

763:     PetscReal aa    = 0.0;
764:     PetscReal eps_A = 0.0;
765:     PetscReal eps_B = 0.0;
766:     PetscReal eps_C = 0.0;
767:     PetscReal time  = 0.0;

769:     for (n = 1; n < N + 1; ++n) {
770:       aa = user->zeroArray[n - 1];
771:       eps_A += (aa * PetscExpReal((-1.0 * aa * aa * c * time) / (a * a)) * PetscCosReal(aa) * PetscCosReal((aa * x[0]) / a)) / (a * (aa - PetscSinReal(aa) * PetscCosReal(aa)));
772:       eps_B += (PetscExpReal((-1.0 * aa * aa * c * time) / (a * a)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa));
773:       eps_C += (PetscExpReal((-1.0 * aa * aa * c * time) / (aa * aa)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa));
774:     }
775:     u[0] = (P_0 / G) * eps_A + ((P_0 * nu) / (2.0 * G * a)) - eps_B / (G * a) - (P_0 * (1 - nu)) / (2 * G * a) + eps_C / (G * a);
776:   }
777:   return PETSC_SUCCESS;
778: }

780: // Displacement
781: static PetscErrorCode mandel_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
782: {
783:   Parameter *param;

785:   AppCtx *user = (AppCtx *)ctx;

787:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
788:   if (time <= 0.0) {
789:     PetscCall(mandel_initial_u(dim, time, x, Nc, u, ctx));
790:   } else {
791:     PetscInt    NITER = user->niter;
792:     PetscScalar alpha = param->alpha;
793:     PetscScalar K_u   = param->K_u;
794:     PetscScalar M     = param->M;
795:     PetscScalar G     = param->mu;
796:     PetscScalar k     = param->k;
797:     PetscScalar mu_f  = param->mu_f;
798:     PetscScalar F     = param->P_0;

800:     PetscScalar K_d   = K_u - alpha * alpha * M;
801:     PetscScalar nu    = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
802:     PetscScalar nu_u  = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
803:     PetscScalar kappa = k / mu_f;
804:     PetscReal   a     = (user->xmax[0] - user->xmin[0]) / 2.0;
805:     PetscReal   c     = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u)));

807:     // Series term
808:     PetscScalar A_x = 0.0;
809:     PetscScalar B_x = 0.0;

811:     for (PetscInt n = 1; n < NITER + 1; n++) {
812:       PetscReal alpha_n = user->zeroArray[n - 1];

814:       A_x += ((PetscSinReal(alpha_n) * PetscCosReal(alpha_n)) / (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscExpReal(-1 * (alpha_n * alpha_n * c * time) / (a * a));
815:       B_x += (PetscCosReal(alpha_n) / (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscSinReal((alpha_n * x[0]) / a) * PetscExpReal(-1 * (alpha_n * alpha_n * c * time) / (a * a));
816:     }
817:     u[0] = ((F * nu) / (2.0 * G * a) - (F * nu_u) / (G * a) * A_x) * x[0] + F / G * B_x;
818:     u[1] = (-1 * (F * (1.0 - nu)) / (2 * G * a) + (F * (1 - nu_u)) / (G * a) * A_x) * x[1];
819:   }
820:   return PETSC_SUCCESS;
821: }

823: // Trace strain
824: static PetscErrorCode mandel_2d_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
825: {
826:   Parameter *param;

828:   AppCtx *user = (AppCtx *)ctx;

830:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
831:   if (time <= 0.0) {
832:     PetscCall(mandel_initial_eps(dim, time, x, Nc, u, ctx));
833:   } else {
834:     PetscInt    NITER = user->niter;
835:     PetscScalar alpha = param->alpha;
836:     PetscScalar K_u   = param->K_u;
837:     PetscScalar M     = param->M;
838:     PetscScalar G     = param->mu;
839:     PetscScalar k     = param->k;
840:     PetscScalar mu_f  = param->mu_f;
841:     PetscScalar F     = param->P_0;

843:     PetscScalar K_d   = K_u - alpha * alpha * M;
844:     PetscScalar nu    = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
845:     PetscScalar nu_u  = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
846:     PetscScalar kappa = k / mu_f;
847:     //const PetscScalar B = (alpha*M)/(K_d + alpha*alpha * M);

849:     //const PetscScalar b = (YMAX - YMIN) / 2.0;
850:     PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0;
851:     PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u)));

853:     // Series term
854:     PetscReal eps_A = 0.0;
855:     PetscReal eps_B = 0.0;
856:     PetscReal eps_C = 0.0;

858:     for (PetscInt n = 1; n < NITER + 1; n++) {
859:       PetscReal aa = user->zeroArray[n - 1];

861:       eps_A += (aa * PetscExpReal((-1.0 * aa * aa * c * time) / (a * a)) * PetscCosReal(aa) * PetscCosReal((aa * x[0]) / a)) / (a * (aa - PetscSinReal(aa) * PetscCosReal(aa)));

863:       eps_B += (PetscExpReal((-1.0 * aa * aa * c * time) / (a * a)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa));

865:       eps_C += (PetscExpReal((-1.0 * aa * aa * c * time) / (aa * aa)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa));
866:     }

868:     u[0] = (F / G) * eps_A + ((F * nu) / (2.0 * G * a)) - eps_B / (G * a) - (F * (1 - nu)) / (2 * G * a) + eps_C / (G * a);
869:   }
870:   return PETSC_SUCCESS;
871: }

873: // Pressure
874: static PetscErrorCode mandel_2d_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
875: {
876:   Parameter *param;

878:   AppCtx *user = (AppCtx *)ctx;

880:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
881:   if (time <= 0.0) {
882:     PetscCall(mandel_drainage_pressure(dim, time, x, Nc, u, ctx));
883:   } else {
884:     PetscInt NITER = user->niter;

886:     PetscScalar alpha = param->alpha;
887:     PetscScalar K_u   = param->K_u;
888:     PetscScalar M     = param->M;
889:     PetscScalar G     = param->mu;
890:     PetscScalar k     = param->k;
891:     PetscScalar mu_f  = param->mu_f;
892:     PetscScalar F     = param->P_0;

894:     PetscScalar K_d   = K_u - alpha * alpha * M;
895:     PetscScalar nu    = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
896:     PetscScalar nu_u  = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
897:     PetscScalar kappa = k / mu_f;
898:     PetscScalar B     = (alpha * M) / (K_d + alpha * alpha * M);

900:     PetscReal   a  = (user->xmax[0] - user->xmin[0]) / 2.0;
901:     PetscReal   c  = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u)));
902:     PetscScalar A1 = 3.0 / (B * (1.0 + nu_u));
903:     //PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu);

905:     // Series term
906:     PetscReal p = 0.0;

908:     for (PetscInt n = 1; n < NITER + 1; n++) {
909:       PetscReal aa = user->zeroArray[n - 1];
910:       p += (PetscSinReal(aa) / (aa - PetscSinReal(aa) * PetscCosReal(aa))) * (PetscCosReal((aa * x[0]) / a) - PetscCosReal(aa)) * PetscExpReal(-1.0 * (aa * aa * c * time) / (a * a));
911:     }
912:     u[0] = ((2.0 * F) / (a * A1)) * p;
913:   }
914:   return PETSC_SUCCESS;
915: }

917: // Time derivative of displacement
918: static PetscErrorCode mandel_2d_u_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
919: {
920:   Parameter *param;

922:   AppCtx *user = (AppCtx *)ctx;

924:   PetscCall(PetscBagGetData(user->bag, (void **)&param));

926:   PetscInt    NITER = user->niter;
927:   PetscScalar alpha = param->alpha;
928:   PetscScalar K_u   = param->K_u;
929:   PetscScalar M     = param->M;
930:   PetscScalar G     = param->mu;
931:   PetscScalar F     = param->P_0;

933:   PetscScalar K_d   = K_u - alpha * alpha * M;
934:   PetscScalar nu    = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
935:   PetscScalar nu_u  = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
936:   PetscScalar kappa = param->k / param->mu_f;
937:   PetscReal   a     = (user->xmax[0] - user->xmin[0]) / 2.0;
938:   PetscReal   c     = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u)));

940:   // Series term
941:   PetscScalar A_s_t = 0.0;
942:   PetscScalar B_s_t = 0.0;

944:   for (PetscInt n = 1; n < NITER + 1; n++) {
945:     PetscReal alpha_n = user->zeroArray[n - 1];

947:     A_s_t += (-1.0 * alpha_n * alpha_n * c * PetscExpReal((-1.0 * alpha_n * alpha_n * time) / (a * a)) * PetscSinReal((alpha_n * x[0]) / a) * PetscCosReal(alpha_n)) / (a * a * (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n)));
948:     B_s_t += (-1.0 * alpha_n * alpha_n * c * PetscExpReal((-1.0 * alpha_n * alpha_n * time) / (a * a)) * PetscSinReal(alpha_n) * PetscCosReal(alpha_n)) / (a * a * (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n)));
949:   }

951:   u[0] = (F / G) * A_s_t - ((F * nu_u * x[0]) / (G * a)) * B_s_t;
952:   u[1] = ((F * x[1] * (1 - nu_u)) / (G * a)) * B_s_t;

954:   return PETSC_SUCCESS;
955: }

957: // Time derivative of trace strain
958: static PetscErrorCode mandel_2d_eps_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
959: {
960:   Parameter *param;

962:   AppCtx *user = (AppCtx *)ctx;

964:   PetscCall(PetscBagGetData(user->bag, (void **)&param));

966:   PetscInt    NITER = user->niter;
967:   PetscScalar alpha = param->alpha;
968:   PetscScalar K_u   = param->K_u;
969:   PetscScalar M     = param->M;
970:   PetscScalar G     = param->mu;
971:   PetscScalar k     = param->k;
972:   PetscScalar mu_f  = param->mu_f;
973:   PetscScalar F     = param->P_0;

975:   PetscScalar K_d   = K_u - alpha * alpha * M;
976:   PetscScalar nu    = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
977:   PetscScalar nu_u  = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
978:   PetscScalar kappa = k / mu_f;
979:   //const PetscScalar B = (alpha*M)/(K_d + alpha*alpha * M);

981:   //const PetscScalar b = (YMAX - YMIN) / 2.0;
982:   PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0;
983:   PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u)));

985:   // Series term
986:   PetscScalar eps_As = 0.0;
987:   PetscScalar eps_Bs = 0.0;
988:   PetscScalar eps_Cs = 0.0;

990:   for (PetscInt n = 1; n < NITER + 1; n++) {
991:     PetscReal alpha_n = user->zeroArray[n - 1];

993:     eps_As += (-1.0 * alpha_n * alpha_n * alpha_n * c * PetscExpReal((-1.0 * alpha_n * alpha_n * c * time) / (a * a)) * PetscCosReal(alpha_n) * PetscCosReal((alpha_n * x[0]) / a)) / (alpha_n * alpha_n * alpha_n * (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n)));
994:     eps_Bs += (-1.0 * alpha_n * alpha_n * c * PetscExpReal((-1.0 * alpha_n * alpha_n * c * time) / (a * a)) * PetscSinReal(alpha_n) * PetscCosReal(alpha_n)) / (alpha_n * alpha_n * (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n)));
995:     eps_Cs += (-1.0 * alpha_n * alpha_n * c * PetscExpReal((-1.0 * alpha_n * alpha_n * c * time) / (a * a)) * PetscSinReal(alpha_n) * PetscCosReal(alpha_n)) / (alpha_n * alpha_n * (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n)));
996:   }

998:   u[0] = (F / G) * eps_As - ((F * nu_u) / (G * a)) * eps_Bs + ((F * (1 - nu_u)) / (G * a)) * eps_Cs;
999:   return PETSC_SUCCESS;
1000: }

1002: // Time derivative of pressure
1003: static PetscErrorCode mandel_2d_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1004: {
1005:   Parameter *param;

1007:   AppCtx *user = (AppCtx *)ctx;

1009:   PetscCall(PetscBagGetData(user->bag, (void **)&param));

1011:   PetscScalar alpha = param->alpha;
1012:   PetscScalar K_u   = param->K_u;
1013:   PetscScalar M     = param->M;
1014:   PetscScalar G     = param->mu;
1015:   PetscScalar F     = param->P_0;

1017:   PetscScalar K_d  = K_u - alpha * alpha * M;
1018:   PetscScalar nu   = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
1019:   PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));

1021:   PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0;
1022:   //PetscScalar A1 = 3.0 / (B * (1.0 + nu_u));
1023:   //PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu);

1025:   u[0] = ((2.0 * F * (-2.0 * nu + 3.0 * nu_u)) / (3.0 * a * alpha * (1.0 - 2.0 * nu)));

1027:   return PETSC_SUCCESS;
1028: }

1030: /* Cryer Solutions */
1031: static PetscErrorCode cryer_drainage_pressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1032: {
1033:   AppCtx    *user = (AppCtx *)ctx;
1034:   Parameter *param;

1036:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
1037:   if (time <= 0.0) {
1038:     PetscScalar alpha = param->alpha;    /* -  */
1039:     PetscScalar K_u   = param->K_u;      /* Pa */
1040:     PetscScalar M     = param->M;        /* Pa */
1041:     PetscScalar P_0   = param->P_0;      /* Pa */
1042:     PetscScalar B     = alpha * M / K_u; /* -, Cheng (B.12) */

1044:     u[0] = P_0 * B;
1045:   } else {
1046:     u[0] = 0.0;
1047:   }
1048:   return PETSC_SUCCESS;
1049: }

1051: static PetscErrorCode cryer_initial_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1052: {
1053:   AppCtx    *user = (AppCtx *)ctx;
1054:   Parameter *param;

1056:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
1057:   {
1058:     PetscScalar K_u  = param->K_u;                                      /* Pa */
1059:     PetscScalar G    = param->mu;                                       /* Pa */
1060:     PetscScalar P_0  = param->P_0;                                      /* Pa */
1061:     PetscReal   R_0  = user->xmax[1];                                   /* m */
1062:     PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -,       Cheng (B.9)  */

1064:     PetscScalar u_0  = -P_0 * R_0 * (1. - 2. * nu_u) / (2. * G * (1. + nu_u)); /* Cheng (7.407) */
1065:     PetscReal   u_sc = PetscRealPart(u_0) / R_0;

1067:     u[0] = u_sc * x[0];
1068:     u[1] = u_sc * x[1];
1069:     u[2] = u_sc * x[2];
1070:   }
1071:   return PETSC_SUCCESS;
1072: }

1074: static PetscErrorCode cryer_initial_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1075: {
1076:   AppCtx    *user = (AppCtx *)ctx;
1077:   Parameter *param;

1079:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
1080:   {
1081:     PetscScalar K_u  = param->K_u;                                      /* Pa */
1082:     PetscScalar G    = param->mu;                                       /* Pa */
1083:     PetscScalar P_0  = param->P_0;                                      /* Pa */
1084:     PetscReal   R_0  = user->xmax[1];                                   /* m */
1085:     PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -,       Cheng (B.9)  */

1087:     PetscScalar u_0  = -P_0 * R_0 * (1. - 2. * nu_u) / (2. * G * (1. + nu_u)); /* Cheng (7.407) */
1088:     PetscReal   u_sc = PetscRealPart(u_0) / R_0;

1090:     /* div R = 1/R^2 d/dR R^2 R = 3 */
1091:     u[0] = 3. * u_sc;
1092:     u[1] = 3. * u_sc;
1093:     u[2] = 3. * u_sc;
1094:   }
1095:   return PETSC_SUCCESS;
1096: }

1098: // Displacement
1099: static PetscErrorCode cryer_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1100: {
1101:   AppCtx    *user = (AppCtx *)ctx;
1102:   Parameter *param;

1104:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
1105:   if (time <= 0.0) {
1106:     PetscCall(cryer_initial_u(dim, time, x, Nc, u, ctx));
1107:   } else {
1108:     PetscScalar alpha = param->alpha;           /* -  */
1109:     PetscScalar K_u   = param->K_u;             /* Pa */
1110:     PetscScalar M     = param->M;               /* Pa */
1111:     PetscScalar G     = param->mu;              /* Pa */
1112:     PetscScalar P_0   = param->P_0;             /* Pa */
1113:     PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
1114:     PetscReal   R_0   = user->xmax[1];          /* m */
1115:     PetscInt    N     = user->niter, n;

1117:     PetscScalar K_d   = K_u - alpha * alpha * M;                             /* Pa,      Cheng (B.5)  */
1118:     PetscScalar nu    = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));     /* -,       Cheng (B.8)  */
1119:     PetscScalar nu_u  = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));     /* -,       Cheng (B.9)  */
1120:     PetscScalar S     = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
1121:     PetscScalar c     = kappa / S;                                           /* m^2 / s, Cheng (B.16) */
1122:     PetscScalar u_inf = -P_0 * R_0 * (1. - 2. * nu) / (2. * G * (1. + nu));  /* m,       Cheng (7.388) */

1124:     PetscReal   R      = PetscSqrtReal(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
1125:     PetscReal   R_star = R / R_0;
1126:     PetscReal   tstar  = PetscRealPart(c * time) / PetscSqr(R_0); /* - */
1127:     PetscReal   A_n    = 0.0;
1128:     PetscScalar u_sc;

1130:     for (n = 1; n < N + 1; ++n) {
1131:       const PetscReal x_n = user->zeroArray[n - 1];
1132:       const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u));

1134:       /* m , Cheng (7.404) */
1135:       if (R_star != 0) {
1136:         A_n += PetscRealPart((12.0 * (1.0 + nu) * (nu_u - nu)) / ((1.0 - 2.0 * nu) * E_n * PetscSqr(R_star) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * (3.0 * (nu_u - nu) * (PetscSinReal(R_star * PetscSqrtReal(x_n)) - R_star * PetscSqrtReal(x_n) * PetscCosReal(R_star * PetscSqrtReal(x_n))) + (1.0 - nu) * (1.0 - 2.0 * nu) * PetscPowRealInt(R_star, 3) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * PetscExpReal(-x_n * tstar));
1137:       }
1138:     }
1139:     if (R_star != 0) u_sc = PetscRealPart(u_inf) * (R_star - A_n) / R;
1140:     else u_sc = PetscRealPart(u_inf) / R_0;
1141:     u[0] = u_sc * x[0];
1142:     u[1] = u_sc * x[1];
1143:     u[2] = u_sc * x[2];
1144:   }
1145:   return PETSC_SUCCESS;
1146: }

1148: // Volumetric Strain
1149: static PetscErrorCode cryer_3d_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1150: {
1151:   AppCtx    *user = (AppCtx *)ctx;
1152:   Parameter *param;

1154:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
1155:   if (time <= 0.0) {
1156:     PetscCall(cryer_initial_eps(dim, time, x, Nc, u, ctx));
1157:   } else {
1158:     PetscScalar alpha = param->alpha;           /* -  */
1159:     PetscScalar K_u   = param->K_u;             /* Pa */
1160:     PetscScalar M     = param->M;               /* Pa */
1161:     PetscScalar G     = param->mu;              /* Pa */
1162:     PetscScalar P_0   = param->P_0;             /* Pa */
1163:     PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
1164:     PetscReal   R_0   = user->xmax[1];          /* m */
1165:     PetscInt    N     = user->niter, n;

1167:     PetscScalar K_d   = K_u - alpha * alpha * M;                             /* Pa,      Cheng (B.5)  */
1168:     PetscScalar nu    = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));     /* -,       Cheng (B.8)  */
1169:     PetscScalar nu_u  = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));     /* -,       Cheng (B.9)  */
1170:     PetscScalar S     = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
1171:     PetscScalar c     = kappa / S;                                           /* m^2 / s, Cheng (B.16) */
1172:     PetscScalar u_inf = -P_0 * R_0 * (1. - 2. * nu) / (2. * G * (1. + nu));  /* m,       Cheng (7.388) */

1174:     PetscReal R      = PetscSqrtReal(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
1175:     PetscReal R_star = R / R_0;
1176:     PetscReal tstar  = PetscRealPart(c * time) / PetscSqr(R_0); /* - */
1177:     PetscReal divA_n = 0.0;

1179:     if (R_star < PETSC_SMALL) {
1180:       for (n = 1; n < N + 1; ++n) {
1181:         const PetscReal x_n = user->zeroArray[n - 1];
1182:         const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u));

1184:         divA_n += PetscRealPart((12.0 * (1.0 + nu) * (nu_u - nu)) / ((1.0 - 2.0 * nu) * E_n * PetscSqr(R_star) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * (3.0 * (nu_u - nu) * PetscSqrtReal(x_n) * ((2.0 + PetscSqr(R_star * PetscSqrtReal(x_n))) - 2.0 * PetscCosReal(R_star * PetscSqrtReal(x_n))) + 5.0 * (1.0 - nu) * (1.0 - 2.0 * nu) * PetscPowRealInt(R_star, 2) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * PetscExpReal(-x_n * tstar));
1185:       }
1186:     } else {
1187:       for (n = 1; n < N + 1; ++n) {
1188:         const PetscReal x_n = user->zeroArray[n - 1];
1189:         const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u));

1191:         divA_n += PetscRealPart((12.0 * (1.0 + nu) * (nu_u - nu)) / ((1.0 - 2.0 * nu) * E_n * PetscSqr(R_star) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * (3.0 * (nu_u - nu) * PetscSqrtReal(x_n) * ((2.0 / (R_star * PetscSqrtReal(x_n)) + R_star * PetscSqrtReal(x_n)) * PetscSinReal(R_star * PetscSqrtReal(x_n)) - 2.0 * PetscCosReal(R_star * PetscSqrtReal(x_n))) + 5.0 * (1.0 - nu) * (1.0 - 2.0 * nu) * PetscPowRealInt(R_star, 2) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * PetscExpReal(-x_n * tstar));
1192:       }
1193:     }
1194:     if (PetscAbsReal(divA_n) > 1e3) PetscCall(PetscPrintf(PETSC_COMM_SELF, "(%g, %g, %g) divA_n: %g\n", (double)x[0], (double)x[1], (double)x[2], (double)divA_n));
1195:     u[0] = PetscRealPart(u_inf) / R_0 * (3.0 - divA_n);
1196:   }
1197:   return PETSC_SUCCESS;
1198: }

1200: // Pressure
1201: static PetscErrorCode cryer_3d_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1202: {
1203:   AppCtx    *user = (AppCtx *)ctx;
1204:   Parameter *param;

1206:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
1207:   if (time <= 0.0) {
1208:     PetscCall(cryer_drainage_pressure(dim, time, x, Nc, u, ctx));
1209:   } else {
1210:     PetscScalar alpha = param->alpha;           /* -  */
1211:     PetscScalar K_u   = param->K_u;             /* Pa */
1212:     PetscScalar M     = param->M;               /* Pa */
1213:     PetscScalar G     = param->mu;              /* Pa */
1214:     PetscScalar P_0   = param->P_0;             /* Pa */
1215:     PetscReal   R_0   = user->xmax[1];          /* m */
1216:     PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */
1217:     PetscInt    N     = user->niter, n;

1219:     PetscScalar K_d  = K_u - alpha * alpha * M;                             /* Pa,      Cheng (B.5)  */
1220:     PetscScalar eta  = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G);           /* -,       Cheng (B.11) */
1221:     PetscScalar nu   = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));     /* -,       Cheng (B.8)  */
1222:     PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));     /* -,       Cheng (B.9)  */
1223:     PetscScalar S    = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */
1224:     PetscScalar c    = kappa / S;                                           /* m^2 / s, Cheng (B.16) */
1225:     PetscReal   R    = PetscSqrtReal(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);

1227:     PetscReal R_star = R / R_0;
1228:     PetscReal t_star = PetscRealPart(c * time) / PetscSqr(R_0);
1229:     PetscReal A_x    = 0.0;

1231:     for (n = 1; n < N + 1; ++n) {
1232:       const PetscReal x_n = user->zeroArray[n - 1];
1233:       const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u));

1235:       A_x += PetscRealPart(((18.0 * PetscSqr(nu_u - nu)) / (eta * E_n)) * (PetscSinReal(R_star * PetscSqrtReal(x_n)) / (R_star * PetscSinReal(PetscSqrtReal(x_n))) - 1.0) * PetscExpReal(-x_n * t_star)); /* Cheng (7.395) */
1236:     }
1237:     u[0] = P_0 * A_x;
1238:   }
1239:   return PETSC_SUCCESS;
1240: }

1242: /* Boundary Kernels */
1243: static void f0_terzaghi_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1244: {
1245:   const PetscReal P = PetscRealPart(constants[5]);

1247:   f0[0] = 0.0;
1248:   f0[1] = P;
1249: }

1251: #if 0
1252: static void f0_mandel_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1253:                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1254:                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1255:                                     PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1256: {
1257:   // Uniform stress distribution
1258:   /* PetscScalar xmax =  0.5;
1259:   PetscScalar xmin = -0.5;
1260:   PetscScalar ymax =  0.5;
1261:   PetscScalar ymin = -0.5;
1262:   PetscScalar P = constants[5];
1263:   PetscScalar aL = (xmax - xmin) / 2.0;
1264:   PetscScalar sigma_zz = -1.0*P / aL; */

1266:   // Analytical (parabolic) stress distribution
1267:   PetscReal a1, a2, am;
1268:   PetscReal y1, y2, ym;

1270:   PetscInt NITER = 500;
1271:   PetscReal EPS = 0.000001;
1272:   PetscReal zeroArray[500]; /* NITER */
1273:   PetscReal xmax =  1.0;
1274:   PetscReal xmin =  0.0;
1275:   PetscReal ymax =  0.1;
1276:   PetscReal ymin =  0.0;
1277:   PetscReal lower[2], upper[2];

1279:   lower[0] = xmin - (xmax - xmin) / 2.0;
1280:   lower[1] = ymin - (ymax - ymin) / 2.0;
1281:   upper[0] = xmax - (xmax - xmin) / 2.0;
1282:   upper[1] = ymax - (ymax - ymin) / 2.0;

1284:   xmin = lower[0];
1285:   ymin = lower[1];
1286:   xmax = upper[0];
1287:   ymax = upper[1];

1289:   PetscScalar G     = constants[0];
1290:   PetscScalar K_u   = constants[1];
1291:   PetscScalar alpha = constants[2];
1292:   PetscScalar M     = constants[3];
1293:   PetscScalar kappa = constants[4];
1294:   PetscScalar F     = constants[5];

1296:   PetscScalar K_d = K_u - alpha*alpha*M;
1297:   PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G));
1298:   PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G));
1299:   PetscReal   aL = (xmax - xmin) / 2.0;
1300:   PetscReal   c = PetscRealPart(((2.0*kappa*G) * (1.0 - nu) * (nu_u - nu)) / (alpha*alpha * (1.0 - 2.0*nu) * (1.0 - nu_u)));
1301:   PetscScalar B = (3.0 * (nu_u - nu)) / ( alpha * (1.0 - 2.0*nu) * (1.0 + nu_u));
1302:   PetscScalar A1 = 3.0 / (B * (1.0 + nu_u));
1303:   PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu);

1305:   // Generate zero values
1306:   for (PetscInt i=1; i < NITER+1; i++)
1307:   {
1308:     a1 = ((PetscReal) i - 1.0) * PETSC_PI * PETSC_PI / 4.0 + EPS;
1309:     a2 = a1 + PETSC_PI/2;
1310:     for (PetscInt j=0; j<NITER; j++)
1311:     {
1312:       y1 = PetscTanReal(a1) - PetscRealPart(A1/A2)*a1;
1313:       y2 = PetscTanReal(a2) - PetscRealPart(A1/A2)*a2;
1314:       am = (a1 + a2)/2.0;
1315:       ym = PetscTanReal(am) - PetscRealPart(A1/A2)*am;
1316:       if ((ym*y1) > 0)
1317:       {
1318:         a1 = am;
1319:       } else {
1320:         a2 = am;
1321:       }
1322:       if (PetscAbsReal(y2) < EPS)
1323:       {
1324:         am = a2;
1325:       }
1326:     }
1327:     zeroArray[i-1] = am;
1328:   }

1330:   // Solution for sigma_zz
1331:   PetscScalar A_x = 0.0;
1332:   PetscScalar B_x = 0.0;

1334:   for (PetscInt n=1; n < NITER+1; n++)
1335:   {
1336:     PetscReal alpha_n = zeroArray[n-1];

1338:     A_x += ( PetscSinReal(alpha_n) / (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscCosReal( (alpha_n * x[0]) / aL) * PetscExpReal( -1.0*( (alpha_n*alpha_n*c*t)/(aL*aL)));
1339:     B_x += ( (PetscSinReal(alpha_n) * PetscCosReal(alpha_n))/(alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscExpReal( -1.0*( (alpha_n*alpha_n*c*t)/(aL*aL)));
1340:   }

1342:   PetscScalar sigma_zz = -1.0*(F/aL) - ((2.0*F)/aL) * (A2/A1) * A_x + ((2.0*F)/aL) * B_x;

1344:   if (x[1] == ymax) {
1345:     f0[1] += sigma_zz;
1346:   } else if (x[1] == ymin) {
1347:     f0[1] -= sigma_zz;
1348:   }
1349: }
1350: #endif

1352: static void f0_cryer_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1353: {
1354:   const PetscReal P_0 = PetscRealPart(constants[5]);
1355:   PetscInt        d;

1357:   for (d = 0; d < dim; ++d) f0[d] = -P_0 * n[d];
1358: }

1360: /* Standard Kernels - Residual */
1361: /* f0_e */
1362: static void f0_epsilon(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1363: {
1364:   PetscInt d;

1366:   for (d = 0; d < dim; ++d) f0[0] += u_x[d * dim + d];
1367:   f0[0] -= u[uOff[1]];
1368: }

1370: /* f0_p */
1371: static void f0_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1372: {
1373:   const PetscReal alpha = PetscRealPart(constants[2]);
1374:   const PetscReal M     = PetscRealPart(constants[3]);

1376:   f0[0] += alpha * u_t[uOff[1]];
1377:   f0[0] += u_t[uOff[2]] / M;
1378:   if (f0[0] != f0[0]) abort();
1379: }

1381: /* f1_u */
1382: static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
1383: {
1384:   const PetscInt  Nc     = dim;
1385:   const PetscReal G      = PetscRealPart(constants[0]);
1386:   const PetscReal K_u    = PetscRealPart(constants[1]);
1387:   const PetscReal alpha  = PetscRealPart(constants[2]);
1388:   const PetscReal M      = PetscRealPart(constants[3]);
1389:   const PetscReal K_d    = K_u - alpha * alpha * M;
1390:   const PetscReal lambda = K_d - (2.0 * G) / 3.0;
1391:   PetscInt        c, d;

1393:   for (c = 0; c < Nc; ++c) {
1394:     for (d = 0; d < dim; ++d) f1[c * dim + d] -= G * (u_x[c * dim + d] + u_x[d * dim + c]);
1395:     f1[c * dim + c] -= lambda * u[uOff[1]];
1396:     f1[c * dim + c] += alpha * u[uOff[2]];
1397:   }
1398: }

1400: /* f1_p */
1401: static void f1_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
1402: {
1403:   const PetscReal kappa = PetscRealPart(constants[4]);
1404:   PetscInt        d;

1406:   for (d = 0; d < dim; ++d) f1[d] += kappa * u_x[uOff_x[2] + d];
1407: }

1409: /*
1410:   \partial_df \phi_fc g_{fc,gc,df,dg} \partial_dg \phi_gc

1412:   \partial_df \phi_fc \lambda \delta_{fc,df} \sum_gc \partial_dg \phi_gc \delta_{gc,dg}
1413:   = \partial_fc \phi_fc \sum_gc \partial_gc \phi_gc
1414: */

1416: /* Standard Kernels - Jacobian */
1417: /* g0_ee */
1418: static void g0_ee(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1419: {
1420:   g0[0] = -1.0;
1421: }

1423: /* g0_pe */
1424: static void g0_pe(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1425: {
1426:   const PetscReal alpha = PetscRealPart(constants[2]);

1428:   g0[0] = u_tShift * alpha;
1429: }

1431: /* g0_pp */
1432: static void g0_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1433: {
1434:   const PetscReal M = PetscRealPart(constants[3]);

1436:   g0[0] = u_tShift / M;
1437: }

1439: /* g1_eu */
1440: static void g1_eu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
1441: {
1442:   PetscInt d;
1443:   for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */
1444: }

1446: /* g2_ue */
1447: static void g2_ue(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
1448: {
1449:   const PetscReal G      = PetscRealPart(constants[0]);
1450:   const PetscReal K_u    = PetscRealPart(constants[1]);
1451:   const PetscReal alpha  = PetscRealPart(constants[2]);
1452:   const PetscReal M      = PetscRealPart(constants[3]);
1453:   const PetscReal K_d    = K_u - alpha * alpha * M;
1454:   const PetscReal lambda = K_d - (2.0 * G) / 3.0;
1455:   PetscInt        d;

1457:   for (d = 0; d < dim; ++d) g2[d * dim + d] -= lambda;
1458: }
1459: /* g2_up */
1460: static void g2_up(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
1461: {
1462:   const PetscReal alpha = PetscRealPart(constants[2]);
1463:   PetscInt        d;

1465:   for (d = 0; d < dim; ++d) g2[d * dim + d] += alpha;
1466: }

1468: /* g3_uu */
1469: static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
1470: {
1471:   const PetscInt  Nc = dim;
1472:   const PetscReal G  = PetscRealPart(constants[0]);
1473:   PetscInt        c, d;

1475:   for (c = 0; c < Nc; ++c) {
1476:     for (d = 0; d < dim; ++d) {
1477:       g3[((c * Nc + c) * dim + d) * dim + d] -= G;
1478:       g3[((c * Nc + d) * dim + d) * dim + c] -= G;
1479:     }
1480:   }
1481: }

1483: /* g3_pp */
1484: static void g3_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
1485: {
1486:   const PetscReal kappa = PetscRealPart(constants[4]);
1487:   PetscInt        d;

1489:   for (d = 0; d < dim; ++d) g3[d * dim + d] += kappa;
1490: }

1492: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
1493: {
1494:   PetscInt sol;

1496:   PetscFunctionBeginUser;
1497:   options->solType   = SOL_QUADRATIC_TRIG;
1498:   options->niter     = 500;
1499:   options->eps       = PETSC_SMALL;
1500:   options->dtInitial = -1.0;
1501:   PetscOptionsBegin(comm, "", "Biot Poroelasticity Options", "DMPLEX");
1502:   PetscCall(PetscOptionsInt("-niter", "Number of series term iterations in exact solutions", "ex53.c", options->niter, &options->niter, NULL));
1503:   sol = options->solType;
1504:   PetscCall(PetscOptionsEList("-sol_type", "Type of exact solution", "ex53.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL));
1505:   options->solType = (SolutionType)sol;
1506:   PetscCall(PetscOptionsReal("-eps", "Precision value for root finding", "ex53.c", options->eps, &options->eps, NULL));
1507:   PetscCall(PetscOptionsReal("-dt_initial", "Override the initial timestep", "ex53.c", options->dtInitial, &options->dtInitial, NULL));
1508:   PetscOptionsEnd();
1509:   PetscFunctionReturn(PETSC_SUCCESS);
1510: }

1512: static PetscErrorCode mandelZeros(MPI_Comm comm, AppCtx *ctx, Parameter *param)
1513: {
1514:   //PetscBag       bag;
1515:   PetscReal a1, a2, am;
1516:   PetscReal y1, y2, ym;

1518:   PetscFunctionBeginUser;
1519:   //PetscCall(PetscBagGetData(ctx->bag, (void **) ¶m));
1520:   PetscInt  NITER = ctx->niter;
1521:   PetscReal EPS   = ctx->eps;
1522:   //const PetscScalar YMAX = param->ymax;
1523:   //const PetscScalar YMIN = param->ymin;
1524:   PetscScalar alpha = param->alpha;
1525:   PetscScalar K_u   = param->K_u;
1526:   PetscScalar M     = param->M;
1527:   PetscScalar G     = param->mu;
1528:   //const PetscScalar k = param->k;
1529:   //const PetscScalar mu_f = param->mu_f;
1530:   //const PetscScalar P_0 = param->P_0;

1532:   PetscScalar K_d  = K_u - alpha * alpha * M;
1533:   PetscScalar nu   = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
1534:   PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
1535:   //const PetscScalar kappa = k / mu_f;

1537:   // Generate zero values
1538:   for (PetscInt i = 1; i < NITER + 1; i++) {
1539:     a1 = ((PetscReal)i - 1.0) * PETSC_PI * PETSC_PI / 4.0 + EPS;
1540:     a2 = a1 + PETSC_PI / 2;
1541:     am = a1;
1542:     for (PetscInt j = 0; j < NITER; j++) {
1543:       y1 = PetscTanReal(a1) - PetscRealPart((1.0 - nu) / (nu_u - nu)) * a1;
1544:       y2 = PetscTanReal(a2) - PetscRealPart((1.0 - nu) / (nu_u - nu)) * a2;
1545:       am = (a1 + a2) / 2.0;
1546:       ym = PetscTanReal(am) - PetscRealPart((1.0 - nu) / (nu_u - nu)) * am;
1547:       if ((ym * y1) > 0) {
1548:         a1 = am;
1549:       } else {
1550:         a2 = am;
1551:       }
1552:       if (PetscAbsReal(y2) < EPS) am = a2;
1553:     }
1554:     ctx->zeroArray[i - 1] = am;
1555:   }
1556:   PetscFunctionReturn(PETSC_SUCCESS);
1557: }

1559: static PetscReal CryerFunction(PetscReal nu_u, PetscReal nu, PetscReal x)
1560: {
1561:   return PetscTanReal(PetscSqrtReal(x)) * (6.0 * (nu_u - nu) - (1.0 - nu) * (1.0 + nu_u) * x) - (6.0 * (nu_u - nu) * PetscSqrtReal(x));
1562: }

1564: static PetscErrorCode cryerZeros(MPI_Comm comm, AppCtx *ctx, Parameter *param)
1565: {
1566:   PetscReal alpha = PetscRealPart(param->alpha); /* -  */
1567:   PetscReal K_u   = PetscRealPart(param->K_u);   /* Pa */
1568:   PetscReal M     = PetscRealPart(param->M);     /* Pa */
1569:   PetscReal G     = PetscRealPart(param->mu);    /* Pa */
1570:   PetscInt  N     = ctx->niter, n;

1572:   PetscReal K_d  = K_u - alpha * alpha * M;                         /* Pa,      Cheng (B.5)  */
1573:   PetscReal nu   = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -,       Cheng (B.8)  */
1574:   PetscReal nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -,       Cheng (B.9)  */

1576:   PetscFunctionBeginUser;
1577:   for (n = 1; n < N + 1; ++n) {
1578:     PetscReal tol = PetscPowReal(n, 1.5) * ctx->eps;
1579:     PetscReal a1 = 0., a2 = 0., am = 0.;
1580:     PetscReal y1, y2, ym;
1581:     PetscInt  j, k = n - 1;

1583:     y1 = y2 = 1.;
1584:     while (y1 * y2 > 0) {
1585:       ++k;
1586:       a1 = PetscSqr(n * PETSC_PI) - k * PETSC_PI;
1587:       a2 = PetscSqr(n * PETSC_PI) + k * PETSC_PI;
1588:       y1 = CryerFunction(nu_u, nu, a1);
1589:       y2 = CryerFunction(nu_u, nu, a2);
1590:     }
1591:     for (j = 0; j < 50000; ++j) {
1592:       y1 = CryerFunction(nu_u, nu, a1);
1593:       y2 = CryerFunction(nu_u, nu, a2);
1594:       PetscCheck(y1 * y2 <= 0, comm, PETSC_ERR_PLIB, "Invalid root finding initialization for root %" PetscInt_FMT ", (%g, %g)--(%g, %g)", n, (double)a1, (double)y1, (double)a2, (double)y2);
1595:       am = (a1 + a2) / 2.0;
1596:       ym = CryerFunction(nu_u, nu, am);
1597:       if ((ym * y1) < 0) a2 = am;
1598:       else a1 = am;
1599:       if (PetscAbsReal(ym) < tol) break;
1600:     }
1601:     PetscCheck(PetscAbsReal(ym) < tol, comm, PETSC_ERR_PLIB, "Root finding did not converge for root %" PetscInt_FMT " (%g)", n, (double)PetscAbsReal(ym));
1602:     ctx->zeroArray[n - 1] = am;
1603:   }
1604:   PetscFunctionReturn(PETSC_SUCCESS);
1605: }

1607: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
1608: {
1609:   PetscBag   bag;
1610:   Parameter *p;

1612:   PetscFunctionBeginUser;
1613:   /* setup PETSc parameter bag */
1614:   PetscCall(PetscBagGetData(ctx->bag, (void **)&p));
1615:   PetscCall(PetscBagSetName(ctx->bag, "par", "Poroelastic Parameters"));
1616:   bag = ctx->bag;
1617:   if (ctx->solType == SOL_TERZAGHI) {
1618:     // Realistic values - Terzaghi
1619:     PetscCall(PetscBagRegisterScalar(bag, &p->mu, 3.0, "mu", "Shear Modulus, Pa"));
1620:     PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 9.76, "K_u", "Undrained Bulk Modulus, Pa"));
1621:     PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -"));
1622:     PetscCall(PetscBagRegisterScalar(bag, &p->M, 16.0, "M", "Biot Modulus, Pa"));
1623:     PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2"));
1624:     PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s"));
1625:     PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa"));
1626:   } else if (ctx->solType == SOL_MANDEL) {
1627:     // Realistic values - Mandel
1628:     PetscCall(PetscBagRegisterScalar(bag, &p->mu, 0.75, "mu", "Shear Modulus, Pa"));
1629:     PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 2.6941176470588233, "K_u", "Undrained Bulk Modulus, Pa"));
1630:     PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -"));
1631:     PetscCall(PetscBagRegisterScalar(bag, &p->M, 4.705882352941176, "M", "Biot Modulus, Pa"));
1632:     PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2"));
1633:     PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s"));
1634:     PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa"));
1635:   } else if (ctx->solType == SOL_CRYER) {
1636:     // Realistic values - Mandel
1637:     PetscCall(PetscBagRegisterScalar(bag, &p->mu, 0.75, "mu", "Shear Modulus, Pa"));
1638:     PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 2.6941176470588233, "K_u", "Undrained Bulk Modulus, Pa"));
1639:     PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -"));
1640:     PetscCall(PetscBagRegisterScalar(bag, &p->M, 4.705882352941176, "M", "Biot Modulus, Pa"));
1641:     PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2"));
1642:     PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s"));
1643:     PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa"));
1644:   } else {
1645:     // Nonsense values
1646:     PetscCall(PetscBagRegisterScalar(bag, &p->mu, 1.0, "mu", "Shear Modulus, Pa"));
1647:     PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 1.0, "K_u", "Undrained Bulk Modulus, Pa"));
1648:     PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 1.0, "alpha", "Biot Effective Stress Coefficient, -"));
1649:     PetscCall(PetscBagRegisterScalar(bag, &p->M, 1.0, "M", "Biot Modulus, Pa"));
1650:     PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.0, "k", "Isotropic Permeability, m**2"));
1651:     PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s"));
1652:     PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa"));
1653:   }
1654:   PetscCall(PetscBagSetFromOptions(bag));
1655:   {
1656:     PetscScalar K_d  = p->K_u - p->alpha * p->alpha * p->M;
1657:     PetscScalar nu_u = (3.0 * p->K_u - 2.0 * p->mu) / (2.0 * (3.0 * p->K_u + p->mu));
1658:     PetscScalar nu   = (3.0 * K_d - 2.0 * p->mu) / (2.0 * (3.0 * K_d + p->mu));
1659:     PetscScalar S    = (3.0 * p->K_u + 4.0 * p->mu) / (p->M * (3.0 * K_d + 4.0 * p->mu));
1660:     PetscReal   c    = PetscRealPart((p->k / p->mu_f) / S);

1662:     PetscViewer       viewer;
1663:     PetscViewerFormat format;
1664:     PetscBool         flg;

1666:     switch (ctx->solType) {
1667:     case SOL_QUADRATIC_LINEAR:
1668:     case SOL_QUADRATIC_TRIG:
1669:     case SOL_TRIG_LINEAR:
1670:       ctx->t_r = PetscSqr(ctx->xmax[0] - ctx->xmin[0]) / c;
1671:       break;
1672:     case SOL_TERZAGHI:
1673:       ctx->t_r = PetscSqr(2.0 * (ctx->xmax[1] - ctx->xmin[1])) / c;
1674:       break;
1675:     case SOL_MANDEL:
1676:       ctx->t_r = PetscSqr(2.0 * (ctx->xmax[1] - ctx->xmin[1])) / c;
1677:       break;
1678:     case SOL_CRYER:
1679:       ctx->t_r = PetscSqr(ctx->xmax[1]) / c;
1680:       break;
1681:     default:
1682:       SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%d)", solutionTypes[PetscMin(ctx->solType, NUM_SOLUTION_TYPES)], ctx->solType);
1683:     }
1684:     PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg));
1685:     if (flg) {
1686:       PetscCall(PetscViewerPushFormat(viewer, format));
1687:       PetscCall(PetscBagView(bag, viewer));
1688:       PetscCall(PetscViewerFlush(viewer));
1689:       PetscCall(PetscViewerPopFormat(viewer));
1690:       PetscCall(PetscViewerDestroy(&viewer));
1691:       PetscCall(PetscPrintf(comm, "  Max displacement: %g %g\n", (double)PetscRealPart(p->P_0 * (ctx->xmax[1] - ctx->xmin[1]) * (1 - 2 * nu_u) / (2 * p->mu * (1 - nu_u))), (double)PetscRealPart(p->P_0 * (ctx->xmax[1] - ctx->xmin[1]) * (1 - 2 * nu) / (2 * p->mu * (1 - nu)))));
1692:       PetscCall(PetscPrintf(comm, "  Relaxation time: %g\n", (double)ctx->t_r));
1693:     }
1694:   }
1695:   PetscFunctionReturn(PETSC_SUCCESS);
1696: }

1698: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
1699: {
1700:   PetscFunctionBeginUser;
1701:   PetscCall(DMCreate(comm, dm));
1702:   PetscCall(DMSetType(*dm, DMPLEX));
1703:   PetscCall(DMSetFromOptions(*dm));
1704:   PetscCall(DMSetApplicationContext(*dm, user));
1705:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
1706:   PetscCall(DMGetBoundingBox(*dm, user->xmin, user->xmax));
1707:   PetscFunctionReturn(PETSC_SUCCESS);
1708: }

1710: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
1711: {
1712:   PetscErrorCode (*exact[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
1713:   PetscErrorCode (*exact_t[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
1714:   PetscDS       ds;
1715:   DMLabel       label;
1716:   PetscWeakForm wf;
1717:   Parameter    *param;
1718:   PetscInt      id_mandel[2];
1719:   PetscInt      comp[1];
1720:   PetscInt      comp_mandel[2];
1721:   PetscInt      dim, id, bd, f;

1723:   PetscFunctionBeginUser;
1724:   PetscCall(DMGetLabel(dm, "marker", &label));
1725:   PetscCall(DMGetDS(dm, &ds));
1726:   PetscCall(PetscDSGetSpatialDimension(ds, &dim));
1727:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
1728:   exact_t[0] = exact_t[1] = exact_t[2] = zero;

1730:   /* Setup Problem Formulation and Boundary Conditions */
1731:   switch (user->solType) {
1732:   case SOL_QUADRATIC_LINEAR:
1733:     PetscCall(PetscDSSetResidual(ds, 0, f0_quadratic_linear_u, f1_u));
1734:     PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1735:     PetscCall(PetscDSSetResidual(ds, 2, f0_quadratic_linear_p, f1_p));
1736:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1737:     PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1738:     PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1739:     PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1740:     PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1741:     PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1742:     PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));
1743:     exact[0]   = quadratic_u;
1744:     exact[1]   = linear_eps;
1745:     exact[2]   = linear_linear_p;
1746:     exact_t[2] = linear_linear_p_t;

1748:     id = 1;
1749:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", label, 1, &id, 0, 0, NULL, (void (*)(void))exact[0], NULL, user, NULL));
1750:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", label, 1, &id, 2, 0, NULL, (void (*)(void))exact[2], (void (*)(void))exact_t[2], user, NULL));
1751:     break;
1752:   case SOL_TRIG_LINEAR:
1753:     PetscCall(PetscDSSetResidual(ds, 0, f0_trig_linear_u, f1_u));
1754:     PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1755:     PetscCall(PetscDSSetResidual(ds, 2, f0_trig_linear_p, f1_p));
1756:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1757:     PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1758:     PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1759:     PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1760:     PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1761:     PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1762:     PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));
1763:     exact[0]   = trig_u;
1764:     exact[1]   = trig_eps;
1765:     exact[2]   = trig_linear_p;
1766:     exact_t[2] = trig_linear_p_t;

1768:     id = 1;
1769:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", label, 1, &id, 0, 0, NULL, (void (*)(void))exact[0], NULL, user, NULL));
1770:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", label, 1, &id, 2, 0, NULL, (void (*)(void))exact[2], (void (*)(void))exact_t[2], user, NULL));
1771:     break;
1772:   case SOL_QUADRATIC_TRIG:
1773:     PetscCall(PetscDSSetResidual(ds, 0, f0_quadratic_trig_u, f1_u));
1774:     PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1775:     PetscCall(PetscDSSetResidual(ds, 2, f0_quadratic_trig_p, f1_p));
1776:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1777:     PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1778:     PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1779:     PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1780:     PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1781:     PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1782:     PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));
1783:     exact[0]   = quadratic_u;
1784:     exact[1]   = linear_eps;
1785:     exact[2]   = linear_trig_p;
1786:     exact_t[2] = linear_trig_p_t;

1788:     id = 1;
1789:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", label, 1, &id, 0, 0, NULL, (void (*)(void))exact[0], NULL, user, NULL));
1790:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", label, 1, &id, 2, 0, NULL, (void (*)(void))exact[2], (void (*)(void))exact_t[2], user, NULL));
1791:     break;
1792:   case SOL_TERZAGHI:
1793:     PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u));
1794:     PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1795:     PetscCall(PetscDSSetResidual(ds, 2, f0_p, f1_p));
1796:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1797:     PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1798:     PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1799:     PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1800:     PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1801:     PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1802:     PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));

1804:     exact[0]   = terzaghi_2d_u;
1805:     exact[1]   = terzaghi_2d_eps;
1806:     exact[2]   = terzaghi_2d_p;
1807:     exact_t[0] = terzaghi_2d_u_t;
1808:     exact_t[1] = terzaghi_2d_eps_t;
1809:     exact_t[2] = terzaghi_2d_p_t;

1811:     id = 1;
1812:     PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "vertical stress", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
1813:     PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1814:     PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_terzaghi_bd_u, 0, NULL));

1816:     id      = 3;
1817:     comp[0] = 1;
1818:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed base", label, 1, &id, 0, 1, comp, (void (*)(void))zero, NULL, user, NULL));
1819:     id      = 2;
1820:     comp[0] = 0;
1821:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed side", label, 1, &id, 0, 1, comp, (void (*)(void))zero, NULL, user, NULL));
1822:     id      = 4;
1823:     comp[0] = 0;
1824:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed side", label, 1, &id, 0, 1, comp, (void (*)(void))zero, NULL, user, NULL));
1825:     id = 1;
1826:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", label, 1, &id, 2, 0, NULL, (void (*)(void))terzaghi_drainage_pressure, NULL, user, NULL));
1827:     break;
1828:   case SOL_MANDEL:
1829:     PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u));
1830:     PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1831:     PetscCall(PetscDSSetResidual(ds, 2, f0_p, f1_p));
1832:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1833:     PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1834:     PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1835:     PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1836:     PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1837:     PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1838:     PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));

1840:     PetscCall(mandelZeros(PETSC_COMM_WORLD, user, param));

1842:     exact[0]   = mandel_2d_u;
1843:     exact[1]   = mandel_2d_eps;
1844:     exact[2]   = mandel_2d_p;
1845:     exact_t[0] = mandel_2d_u_t;
1846:     exact_t[1] = mandel_2d_eps_t;
1847:     exact_t[2] = mandel_2d_p_t;

1849:     id_mandel[0] = 3;
1850:     id_mandel[1] = 1;
1851:     //comp[0] = 1;
1852:     comp_mandel[0] = 0;
1853:     comp_mandel[1] = 1;
1854:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "vertical stress", label, 2, id_mandel, 0, 2, comp_mandel, (void (*)(void))mandel_2d_u, NULL, user, NULL));
1855:     //PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "vertical stress", "marker", 0, 1, comp, NULL, 2, id_mandel, user));
1856:     //PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed base", "marker", 0, 1, comp, (void (*)(void)) zero, 2, id_mandel, user));
1857:     //PetscCall(PetscDSSetBdResidual(ds, 0, f0_mandel_bd_u, NULL));

1859:     id_mandel[0] = 2;
1860:     id_mandel[1] = 4;
1861:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", label, 2, id_mandel, 2, 0, NULL, (void (*)(void))zero, NULL, user, NULL));
1862:     break;
1863:   case SOL_CRYER:
1864:     PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u));
1865:     PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL));
1866:     PetscCall(PetscDSSetResidual(ds, 2, f0_p, f1_p));
1867:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
1868:     PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL));
1869:     PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL));
1870:     PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL));
1871:     PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL));
1872:     PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL));
1873:     PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp));

1875:     PetscCall(cryerZeros(PETSC_COMM_WORLD, user, param));

1877:     exact[0] = cryer_3d_u;
1878:     exact[1] = cryer_3d_eps;
1879:     exact[2] = cryer_3d_p;

1881:     id = 1;
1882:     PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "normal stress", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
1883:     PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1884:     PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_cryer_bd_u, 0, NULL));

1886:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", label, 1, &id, 2, 0, NULL, (void (*)(void))cryer_drainage_pressure, NULL, user, NULL));
1887:     break;
1888:   default:
1889:     SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%d)", solutionTypes[PetscMin(user->solType, NUM_SOLUTION_TYPES)], user->solType);
1890:   }
1891:   for (f = 0; f < 3; ++f) {
1892:     PetscCall(PetscDSSetExactSolution(ds, f, exact[f], user));
1893:     PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, f, exact_t[f], user));
1894:   }

1896:   /* Setup constants */
1897:   {
1898:     PetscScalar constants[6];
1899:     constants[0] = param->mu;              /* shear modulus, Pa */
1900:     constants[1] = param->K_u;             /* undrained bulk modulus, Pa */
1901:     constants[2] = param->alpha;           /* Biot effective stress coefficient, - */
1902:     constants[3] = param->M;               /* Biot modulus, Pa */
1903:     constants[4] = param->k / param->mu_f; /* Darcy coefficient, m**2 / Pa*s */
1904:     constants[5] = param->P_0;             /* Magnitude of Vertical Stress, Pa */
1905:     PetscCall(PetscDSSetConstants(ds, 6, constants));
1906:   }
1907:   PetscFunctionReturn(PETSC_SUCCESS);
1908: }

1910: static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace)
1911: {
1912:   PetscFunctionBeginUser;
1913:   PetscCall(DMPlexCreateRigidBody(dm, origField, nullspace));
1914:   PetscFunctionReturn(PETSC_SUCCESS);
1915: }

1917: static PetscErrorCode SetupFE(DM dm, PetscInt Nf, PetscInt Nc[], const char *name[], PetscErrorCode (*setup)(DM, AppCtx *), void *ctx)
1918: {
1919:   AppCtx         *user = (AppCtx *)ctx;
1920:   DM              cdm  = dm;
1921:   PetscFE         fe;
1922:   PetscQuadrature q = NULL, fq = NULL;
1923:   char            prefix[PETSC_MAX_PATH_LEN];
1924:   PetscInt        dim, f;
1925:   PetscBool       simplex;

1927:   PetscFunctionBeginUser;
1928:   /* Create finite element */
1929:   PetscCall(DMGetDimension(dm, &dim));
1930:   PetscCall(DMPlexIsSimplex(dm, &simplex));
1931:   for (f = 0; f < Nf; ++f) {
1932:     PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name[f]));
1933:     PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc[f], simplex, name[f] ? prefix : NULL, -1, &fe));
1934:     PetscCall(PetscObjectSetName((PetscObject)fe, name[f]));
1935:     if (!q) PetscCall(PetscFEGetQuadrature(fe, &q));
1936:     if (!fq) PetscCall(PetscFEGetFaceQuadrature(fe, &fq));
1937:     PetscCall(PetscFESetQuadrature(fe, q));
1938:     PetscCall(PetscFESetFaceQuadrature(fe, fq));
1939:     PetscCall(DMSetField(dm, f, NULL, (PetscObject)fe));
1940:     PetscCall(PetscFEDestroy(&fe));
1941:   }
1942:   PetscCall(DMCreateDS(dm));
1943:   PetscCall((*setup)(dm, user));
1944:   while (cdm) {
1945:     PetscCall(DMCopyDisc(dm, cdm));
1946:     if (0) PetscCall(DMSetNearNullSpaceConstructor(cdm, 0, CreateElasticityNullSpace));
1947:     /* TODO: Check whether the boundary of coarse meshes is marked */
1948:     PetscCall(DMGetCoarseDM(cdm, &cdm));
1949:   }
1950:   PetscCall(PetscFEDestroy(&fe));
1951:   PetscFunctionReturn(PETSC_SUCCESS);
1952: }

1954: static PetscErrorCode SetInitialConditions(TS ts, Vec u)
1955: {
1956:   DM        dm;
1957:   PetscReal t;

1959:   PetscFunctionBeginUser;
1960:   PetscCall(TSGetDM(ts, &dm));
1961:   PetscCall(TSGetTime(ts, &t));
1962:   if (t <= 0.0) {
1963:     PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
1964:     void   *ctxs[3];
1965:     AppCtx *ctx;

1967:     PetscCall(DMGetApplicationContext(dm, &ctx));
1968:     switch (ctx->solType) {
1969:     case SOL_TERZAGHI:
1970:       funcs[0] = terzaghi_initial_u;
1971:       ctxs[0]  = ctx;
1972:       funcs[1] = terzaghi_initial_eps;
1973:       ctxs[1]  = ctx;
1974:       funcs[2] = terzaghi_drainage_pressure;
1975:       ctxs[2]  = ctx;
1976:       PetscCall(DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u));
1977:       break;
1978:     case SOL_MANDEL:
1979:       funcs[0] = mandel_initial_u;
1980:       ctxs[0]  = ctx;
1981:       funcs[1] = mandel_initial_eps;
1982:       ctxs[1]  = ctx;
1983:       funcs[2] = mandel_drainage_pressure;
1984:       ctxs[2]  = ctx;
1985:       PetscCall(DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u));
1986:       break;
1987:     case SOL_CRYER:
1988:       funcs[0] = cryer_initial_u;
1989:       ctxs[0]  = ctx;
1990:       funcs[1] = cryer_initial_eps;
1991:       ctxs[1]  = ctx;
1992:       funcs[2] = cryer_drainage_pressure;
1993:       ctxs[2]  = ctx;
1994:       PetscCall(DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u));
1995:       break;
1996:     default:
1997:       PetscCall(DMComputeExactSolution(dm, t, u, NULL));
1998:     }
1999:   } else {
2000:     PetscCall(DMComputeExactSolution(dm, t, u, NULL));
2001:   }
2002:   PetscFunctionReturn(PETSC_SUCCESS);
2003: }

2005: /* Need to create Viewer each time because HDF5 can get corrupted */
2006: static PetscErrorCode SolutionMonitor(TS ts, PetscInt steps, PetscReal time, Vec u, void *mctx)
2007: {
2008:   DM                dm;
2009:   Vec               exact;
2010:   PetscViewer       viewer;
2011:   PetscViewerFormat format;
2012:   PetscOptions      options;
2013:   const char       *prefix;

2015:   PetscFunctionBeginUser;
2016:   PetscCall(TSGetDM(ts, &dm));
2017:   PetscCall(PetscObjectGetOptions((PetscObject)ts, &options));
2018:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, &prefix));
2019:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), options, prefix, "-monitor_solution", &viewer, &format, NULL));
2020:   PetscCall(DMGetGlobalVector(dm, &exact));
2021:   PetscCall(DMComputeExactSolution(dm, time, exact, NULL));
2022:   PetscCall(DMSetOutputSequenceNumber(dm, steps, time));
2023:   PetscCall(VecView(exact, viewer));
2024:   PetscCall(VecView(u, viewer));
2025:   PetscCall(DMRestoreGlobalVector(dm, &exact));
2026:   {
2027:     PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
2028:     void     **ectxs;
2029:     PetscReal *err;
2030:     PetscInt   Nf, f;

2032:     PetscCall(DMGetNumFields(dm, &Nf));
2033:     PetscCall(PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err));
2034:     {
2035:       PetscInt Nds, s;

2037:       PetscCall(DMGetNumDS(dm, &Nds));
2038:       for (s = 0; s < Nds; ++s) {
2039:         PetscDS         ds;
2040:         DMLabel         label;
2041:         IS              fieldIS;
2042:         const PetscInt *fields;
2043:         PetscInt        dsNf, f;

2045:         PetscCall(DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds, NULL));
2046:         PetscCall(PetscDSGetNumFields(ds, &dsNf));
2047:         PetscCall(ISGetIndices(fieldIS, &fields));
2048:         for (f = 0; f < dsNf; ++f) {
2049:           const PetscInt field = fields[f];
2050:           PetscCall(PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]));
2051:         }
2052:         PetscCall(ISRestoreIndices(fieldIS, &fields));
2053:       }
2054:     }
2055:     PetscCall(DMComputeL2FieldDiff(dm, time, exacts, ectxs, u, err));
2056:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "Time: %g L_2 Error: [", (double)time));
2057:     for (f = 0; f < Nf; ++f) {
2058:       if (f) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), ", "));
2059:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%g", (double)err[f]));
2060:     }
2061:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "]\n"));
2062:     PetscCall(PetscFree3(exacts, ectxs, err));
2063:   }
2064:   PetscCall(PetscViewerDestroy(&viewer));
2065:   PetscFunctionReturn(PETSC_SUCCESS);
2066: }

2068: static PetscErrorCode SetupMonitor(TS ts, AppCtx *ctx)
2069: {
2070:   PetscViewer       viewer;
2071:   PetscViewerFormat format;
2072:   PetscOptions      options;
2073:   const char       *prefix;
2074:   PetscBool         flg;

2076:   PetscFunctionBeginUser;
2077:   PetscCall(PetscObjectGetOptions((PetscObject)ts, &options));
2078:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, &prefix));
2079:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), options, prefix, "-monitor_solution", &viewer, &format, &flg));
2080:   if (flg) PetscCall(TSMonitorSet(ts, SolutionMonitor, ctx, NULL));
2081:   PetscCall(PetscViewerDestroy(&viewer));
2082:   PetscFunctionReturn(PETSC_SUCCESS);
2083: }

2085: static PetscErrorCode TSAdaptChoose_Terzaghi(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept, PetscReal *wlte, PetscReal *wltea, PetscReal *wlter)
2086: {
2087:   static PetscReal dtTarget = -1.0;
2088:   PetscReal        dtInitial;
2089:   DM               dm;
2090:   AppCtx          *ctx;
2091:   PetscInt         step;

2093:   PetscFunctionBeginUser;
2094:   PetscCall(TSGetDM(ts, &dm));
2095:   PetscCall(DMGetApplicationContext(dm, &ctx));
2096:   PetscCall(TSGetStepNumber(ts, &step));
2097:   dtInitial = ctx->dtInitial < 0.0 ? 1.0e-4 * ctx->t_r : ctx->dtInitial;
2098:   if (!step) {
2099:     if (PetscAbsReal(dtInitial - h) > PETSC_SMALL) {
2100:       *accept  = PETSC_FALSE;
2101:       *next_h  = dtInitial;
2102:       dtTarget = h;
2103:     } else {
2104:       *accept  = PETSC_TRUE;
2105:       *next_h  = dtTarget < 0.0 ? dtInitial : dtTarget;
2106:       dtTarget = -1.0;
2107:     }
2108:   } else {
2109:     *accept = PETSC_TRUE;
2110:     *next_h = h;
2111:   }
2112:   *next_sc = 0;  /* Reuse the same order scheme */
2113:   *wlte    = -1; /* Weighted local truncation error was not evaluated */
2114:   *wltea   = -1; /* Weighted absolute local truncation error was not evaluated */
2115:   *wlter   = -1; /* Weighted relative local truncation error was not evaluated */
2116:   PetscFunctionReturn(PETSC_SUCCESS);
2117: }

2119: int main(int argc, char **argv)
2120: {
2121:   AppCtx      ctx; /* User-defined work context */
2122:   DM          dm;  /* Problem specification */
2123:   TS          ts;  /* Time Series / Nonlinear solver */
2124:   Vec         u;   /* Solutions */
2125:   const char *name[3] = {"displacement", "tracestrain", "pressure"};
2126:   PetscReal   t;
2127:   PetscInt    dim, Nc[3];

2129:   PetscFunctionBeginUser;
2130:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
2131:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
2132:   PetscCall(PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &ctx.bag));
2133:   PetscCall(PetscMalloc1(ctx.niter, &ctx.zeroArray));
2134:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm));
2135:   PetscCall(SetupParameters(PETSC_COMM_WORLD, &ctx));
2136:   /* Primal System */
2137:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
2138:   PetscCall(DMSetApplicationContext(dm, &ctx));
2139:   PetscCall(TSSetDM(ts, dm));

2141:   PetscCall(DMGetDimension(dm, &dim));
2142:   Nc[0] = dim;
2143:   Nc[1] = 1;
2144:   Nc[2] = 1;

2146:   PetscCall(SetupFE(dm, 3, Nc, name, SetupPrimalProblem, &ctx));
2147:   PetscCall(DMCreateGlobalVector(dm, &u));
2148:   PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx));
2149:   PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx));
2150:   PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx));
2151:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
2152:   PetscCall(TSSetFromOptions(ts));
2153:   PetscCall(TSSetComputeInitialCondition(ts, SetInitialConditions));
2154:   PetscCall(SetupMonitor(ts, &ctx));

2156:   if (ctx.solType != SOL_QUADRATIC_TRIG) {
2157:     TSAdapt adapt;

2159:     PetscCall(TSGetAdapt(ts, &adapt));
2160:     adapt->ops->choose = TSAdaptChoose_Terzaghi;
2161:   }
2162:   if (ctx.solType == SOL_CRYER) {
2163:     Mat          J;
2164:     MatNullSpace sp;

2166:     PetscCall(TSSetUp(ts));
2167:     PetscCall(TSGetIJacobian(ts, &J, NULL, NULL, NULL));
2168:     PetscCall(DMPlexCreateRigidBody(dm, 0, &sp));
2169:     PetscCall(MatSetNullSpace(J, sp));
2170:     PetscCall(MatNullSpaceDestroy(&sp));
2171:   }
2172:   PetscCall(TSGetTime(ts, &t));
2173:   PetscCall(DMSetOutputSequenceNumber(dm, 0, t));
2174:   PetscCall(DMTSCheckFromOptions(ts, u));
2175:   PetscCall(SetInitialConditions(ts, u));
2176:   PetscCall(PetscObjectSetName((PetscObject)u, "solution"));
2177:   PetscCall(TSSolve(ts, u));
2178:   PetscCall(DMTSCheckFromOptions(ts, u));
2179:   PetscCall(TSGetSolution(ts, &u));
2180:   PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));

2182:   /* Cleanup */
2183:   PetscCall(VecDestroy(&u));
2184:   PetscCall(TSDestroy(&ts));
2185:   PetscCall(DMDestroy(&dm));
2186:   PetscCall(PetscBagDestroy(&ctx.bag));
2187:   PetscCall(PetscFree(ctx.zeroArray));
2188:   PetscCall(PetscFinalize());
2189:   return 0;
2190: }

2192: /*TEST

2194:   test:
2195:     suffix: 2d_quad_linear
2196:     requires: triangle
2197:     args: -sol_type quadratic_linear -dm_refine 2 \
2198:       -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2199:       -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme

2201:   test:
2202:     suffix: 3d_quad_linear
2203:     requires: ctetgen
2204:     args: -dm_plex_dim 3 -sol_type quadratic_linear -dm_refine 1 \
2205:       -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2206:       -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme

2208:   test:
2209:     suffix: 2d_trig_linear
2210:     requires: triangle
2211:     args: -sol_type trig_linear -dm_refine 1 \
2212:       -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2213:       -dmts_check .0001 -ts_max_steps 5 -ts_dt 0.00001 -ts_monitor_extreme

2215:   test:
2216:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9, 2.1, 1.8]
2217:     suffix: 2d_trig_linear_sconv
2218:     requires: triangle
2219:     args: -sol_type trig_linear -dm_refine 1 \
2220:       -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2221:       -convest_num_refine 1 -ts_convergence_estimate -ts_convergence_temporal 0 -ts_max_steps 1 -ts_dt 0.00001 -pc_type lu

2223:   test:
2224:     suffix: 3d_trig_linear
2225:     requires: ctetgen
2226:     args: -dm_plex_dim 3 -sol_type trig_linear -dm_refine 1 \
2227:       -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2228:       -dmts_check .0001 -ts_max_steps 2 -ts_monitor_extreme

2230:   test:
2231:     # -dm_refine 1 -convest_num_refine 2 gets L_2 convergence rate: [2.0, 2.1, 1.9]
2232:     suffix: 3d_trig_linear_sconv
2233:     requires: ctetgen
2234:     args: -dm_plex_dim 3 -sol_type trig_linear -dm_refine 1 \
2235:       -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2236:       -convest_num_refine 1 -ts_convergence_estimate -ts_convergence_temporal 0 -ts_max_steps 1 -pc_type lu

2238:   test:
2239:     suffix: 2d_quad_trig
2240:     requires: triangle
2241:     args: -sol_type quadratic_trig -dm_refine 2 \
2242:       -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2243:       -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme

2245:   test:
2246:     # Using -dm_refine 4 gets the convergence rates to [0.95, 0.97, 0.90]
2247:     suffix: 2d_quad_trig_tconv
2248:     requires: triangle
2249:     args: -sol_type quadratic_trig -dm_refine 1 \
2250:       -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2251:       -convest_num_refine 3 -ts_convergence_estimate -ts_max_steps 5 -pc_type lu

2253:   test:
2254:     suffix: 3d_quad_trig
2255:     requires: ctetgen
2256:     args: -dm_plex_dim 3 -sol_type quadratic_trig -dm_refine 1 \
2257:       -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2258:       -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme

2260:   test:
2261:     # Using -dm_refine 2 -convest_num_refine 3 gets the convergence rates to [1.0, 1.0, 1.0]
2262:     suffix: 3d_quad_trig_tconv
2263:     requires: ctetgen
2264:     args: -dm_plex_dim 3 -sol_type quadratic_trig -dm_refine 1 \
2265:       -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2266:       -convest_num_refine 1 -ts_convergence_estimate -ts_max_steps 5 -pc_type lu

2268:   testset:
2269:     args: -sol_type terzaghi -dm_plex_simplex 0 -dm_plex_box_faces 1,8 -dm_plex_box_lower 0,0 -dm_plex_box_upper 10,10 -dm_plex_separate_marker \
2270:           -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 -niter 16000 \
2271:           -pc_type lu

2273:     test:
2274:       suffix: 2d_terzaghi
2275:       requires: double
2276:       args: -ts_dt 0.0028666667 -ts_max_steps 2 -ts_monitor -dmts_check .0001

2278:     test:
2279:       # -dm_plex_box_faces 1,64 -ts_max_steps 4 -convest_num_refine 3 gives L_2 convergence rate: [1.1, 1.1, 1.1]
2280:       suffix: 2d_terzaghi_tconv
2281:       args: -ts_dt 0.023 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1

2283:     test:
2284:       # -dm_plex_box_faces 1,16 -convest_num_refine 4 gives L_2 convergence rate: [1.7, 1.2, 1.1]
2285:       # if we add -displacement_petscspace_degree 3 -tracestrain_petscspace_degree 2 -pressure_petscspace_degree 2, we get [2.1, 1.6, 1.5]
2286:       suffix: 2d_terzaghi_sconv
2287:       args: -ts_dt 1e-5 -dt_initial 1e-5 -ts_max_steps 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1

2289:   testset:
2290:     args: -sol_type mandel -dm_plex_simplex 0 -dm_plex_box_lower -0.5,-0.125 -dm_plex_box_upper 0.5,0.125 -dm_plex_separate_marker -dm_refine 1 \
2291:           -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \
2292:           -pc_type lu

2294:     test:
2295:       suffix: 2d_mandel
2296:       requires: double
2297:       args: -ts_dt 0.0028666667 -ts_max_steps 2 -ts_monitor -dmts_check .0001

2299:     test:
2300:       # -dm_refine 3 -ts_max_steps 4 -convest_num_refine 3 gives L_2 convergence rate: [1.6, 0.93, 1.2]
2301:       suffix: 2d_mandel_sconv
2302:       args: -ts_dt 1e-5 -dt_initial 1e-5 -ts_max_steps 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1

2304:     test:
2305:       # -dm_refine 5 -ts_max_steps 4 -convest_num_refine 3 gives L_2 convergence rate: [0.26, -0.0058, 0.26]
2306:       suffix: 2d_mandel_tconv
2307:       args: -ts_dt 0.023 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1

2309:   testset:
2310:     requires: ctetgen !complex
2311:     args: -sol_type cryer -dm_plex_dim 3 -dm_plex_shape ball \
2312:           -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1

2314:     test:
2315:       suffix: 3d_cryer
2316:       args: -ts_dt 0.0028666667 -ts_max_time 0.014333 -ts_max_steps 2 -dmts_check .0001 \
2317:             -pc_type svd

2319:     test:
2320:       # -bd_dm_refine 3 -dm_refine_volume_limit_pre 0.004 -convest_num_refine 2 gives L_2 convergence rate: []
2321:       suffix: 3d_cryer_sconv
2322:       args: -bd_dm_refine 1 -dm_refine_volume_limit_pre 0.00666667 \
2323:             -ts_dt 1e-5 -dt_initial 1e-5 -ts_max_steps 2 \
2324:             -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
2325:             -pc_type lu -pc_factor_shift_type nonzero

2327:     test:
2328:       # Displacement and Pressure converge. The analytic expression for trace strain is inaccurate at the origin
2329:       # -bd_dm_refine 3 -ref_limit 0.00666667 -ts_max_steps 5 -convest_num_refine 2 gives L_2 convergence rate: [0.47, -0.43, 1.5]
2330:       suffix: 3d_cryer_tconv
2331:       args: -bd_dm_refine 1 -dm_refine_volume_limit_pre 0.00666667 \
2332:             -ts_dt 0.023 -ts_max_time 0.092 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1 \
2333:             -pc_type lu -pc_factor_shift_type nonzero

2335: TEST*/