Actual source code: ex53.c

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

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

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

 14: /* This presentation of poroelasticity is taken from

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

 24: For visualization, use

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

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

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

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

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

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

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

 78: /* Quadratic space and linear time solution

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

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

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

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

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

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

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

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

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

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

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

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

174: /* Quadratic space and trigonometric time solution

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

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

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

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

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

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

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

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

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

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

257: /* Trigonometric space and linear time solution

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

378:   PetscBagGetData(user->bag, (void **)&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 0;
394: }

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

401:   PetscBagGetData(user->bag, (void **)&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 0;
414: }

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

421:   PetscBagGetData(user->bag, (void **)&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 0;
431: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

549:   PetscBagGetData(user->bag, (void **)&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 0;
580: }

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

587:   PetscBagGetData(user->bag, (void **)&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 0;
616: }

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

623:   PetscBagGetData(user->bag, (void **)&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;


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

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

674:   PetscBagGetData(user->bag, (void **)&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 0;
705: }

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

712:   PetscBagGetData(user->bag, (void **)&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:     PetscScalar a     = 0.5 * (user->xmax[0] - user->xmin[0]); /* m */
721:     PetscInt    N     = user->niter, n;

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

729:     PetscScalar A_s     = 0.0;
730:     PetscScalar B_s     = 0.0;
731:     PetscScalar time    = 0.0;
732:     PetscScalar alpha_n = 0.0;

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

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

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

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

766:     PetscReal aa    = 0.0;
767:     PetscReal eps_A = 0.0;
768:     PetscReal eps_B = 0.0;
769:     PetscReal eps_C = 0.0;
770:     PetscReal time  = 0.0;

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

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

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

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

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

810:     // Series term
811:     PetscScalar A_x = 0.0;
812:     PetscScalar B_x = 0.0;

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

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

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

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

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

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

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

856:     // Series term
857:     PetscScalar eps_A = 0.0;
858:     PetscScalar eps_B = 0.0;
859:     PetscScalar eps_C = 0.0;

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

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

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

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

871:     u[0] = (F / G) * eps_A + ((F * nu) / (2.0 * G * a)) - eps_B / (G * a) - (F * (1 - nu)) / (2 * G * a) + eps_C / (G * a);
872:   }
873:   return 0;
874: }

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

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

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

889:     PetscScalar alpha = param->alpha;
890:     PetscScalar K_u   = param->K_u;
891:     PetscScalar M     = param->M;
892:     PetscScalar G     = param->mu;
893:     PetscScalar k     = param->k;
894:     PetscScalar mu_f  = param->mu_f;
895:     PetscScalar F     = param->P_0;

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

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

908:     // Series term
909:     PetscScalar aa = 0.0;
910:     PetscScalar p  = 0.0;

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

921: // Time derivative of displacement
922: static PetscErrorCode mandel_2d_u_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
923: {
924:   Parameter *param;

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

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

930:   PetscInt    NITER = user->niter;
931:   PetscScalar alpha = param->alpha;
932:   PetscScalar K_u   = param->K_u;
933:   PetscScalar M     = param->M;
934:   PetscScalar G     = param->mu;
935:   PetscScalar F     = param->P_0;

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

944:   // Series term
945:   PetscScalar A_s_t = 0.0;
946:   PetscScalar B_s_t = 0.0;

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

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

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

958:   return 0;
959: }

961: // Time derivative of trace strain
962: static PetscErrorCode mandel_2d_eps_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
963: {
964:   Parameter *param;

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

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

970:   PetscInt    NITER = user->niter;
971:   PetscScalar alpha = param->alpha;
972:   PetscScalar K_u   = param->K_u;
973:   PetscScalar M     = param->M;
974:   PetscScalar G     = param->mu;
975:   PetscScalar k     = param->k;
976:   PetscScalar mu_f  = param->mu_f;
977:   PetscScalar F     = param->P_0;

979:   PetscScalar K_d   = K_u - alpha * alpha * M;
980:   PetscScalar nu    = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
981:   PetscScalar nu_u  = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));
982:   PetscScalar kappa = k / mu_f;
983:   //const PetscScalar B = (alpha*M)/(K_d + alpha*alpha * M);

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

989:   // Series term
990:   PetscScalar eps_As = 0.0;
991:   PetscScalar eps_Bs = 0.0;
992:   PetscScalar eps_Cs = 0.0;

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

997:     eps_As += (-1.0 * alpha_n * alpha_n * alpha_n * c * PetscExpReal((-1.0 * alpha_n * alpha_n * c * time) / (a * a)) * PetscCosReal(alpha_n) * PetscCosReal((alpha_n * x[0]) / a)) / (alpha_n * alpha_n * alpha_n * (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n)));
998:     eps_Bs += (-1.0 * alpha_n * alpha_n * c * PetscExpReal((-1.0 * alpha_n * alpha_n * c * time) / (a * a)) * PetscSinReal(alpha_n) * PetscCosReal(alpha_n)) / (alpha_n * alpha_n * (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n)));
999:     eps_Cs += (-1.0 * alpha_n * alpha_n * c * PetscExpReal((-1.0 * alpha_n * alpha_n * c * time) / (a * a)) * PetscSinReal(alpha_n) * PetscCosReal(alpha_n)) / (alpha_n * alpha_n * (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n)));
1000:   }

1002:   u[0] = (F / G) * eps_As - ((F * nu_u) / (G * a)) * eps_Bs + ((F * (1 - nu_u)) / (G * a)) * eps_Cs;
1003:   return 0;
1004: }

1006: // Time derivative of pressure
1007: static PetscErrorCode mandel_2d_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1008: {
1009:   Parameter *param;

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

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

1015:   PetscScalar alpha = param->alpha;
1016:   PetscScalar K_u   = param->K_u;
1017:   PetscScalar M     = param->M;
1018:   PetscScalar G     = param->mu;
1019:   PetscScalar F     = param->P_0;

1021:   PetscScalar K_d  = K_u - alpha * alpha * M;
1022:   PetscScalar nu   = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G));
1023:   PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G));

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

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

1031:   return 0;
1032: }

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

1040:   PetscBagGetData(user->bag, (void **)&param);
1041:   if (time <= 0.0) {
1042:     PetscScalar alpha = param->alpha;    /* -  */
1043:     PetscScalar K_u   = param->K_u;      /* Pa */
1044:     PetscScalar M     = param->M;        /* Pa */
1045:     PetscScalar P_0   = param->P_0;      /* Pa */
1046:     PetscScalar B     = alpha * M / K_u; /* -, Cheng (B.12) */

1048:     u[0] = P_0 * B;
1049:   } else {
1050:     u[0] = 0.0;
1051:   }
1052:   return 0;
1053: }

1055: static PetscErrorCode cryer_initial_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1056: {
1057:   AppCtx    *user = (AppCtx *)ctx;
1058:   Parameter *param;

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

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

1071:     u[0] = u_sc * x[0];
1072:     u[1] = u_sc * x[1];
1073:     u[2] = u_sc * x[2];
1074:   }
1075:   return 0;
1076: }

1078: static PetscErrorCode cryer_initial_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1079: {
1080:   AppCtx    *user = (AppCtx *)ctx;
1081:   Parameter *param;

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

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

1094:     /* div R = 1/R^2 d/dR R^2 R = 3 */
1095:     u[0] = 3. * u_sc;
1096:     u[1] = 3. * u_sc;
1097:     u[2] = 3. * u_sc;
1098:   }
1099:   return 0;
1100: }

1102: // Displacement
1103: static PetscErrorCode cryer_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1104: {
1105:   AppCtx    *user = (AppCtx *)ctx;
1106:   Parameter *param;

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

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

1128:     PetscReal   R      = PetscSqrtReal(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
1129:     PetscReal   R_star = R / R_0;
1130:     PetscReal   tstar  = PetscRealPart(c * time) / PetscSqr(R_0); /* - */
1131:     PetscReal   A_n    = 0.0;
1132:     PetscScalar u_sc;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1584:     y1 = y2 = 1.;
1585:     while (y1 * y2 > 0) {
1586:       ++k;
1587:       a1 = PetscSqr(n * PETSC_PI) - k * PETSC_PI;
1588:       a2 = PetscSqr(n * PETSC_PI) + k * PETSC_PI;
1589:       y1 = CryerFunction(nu_u, nu, a1);
1590:       y2 = CryerFunction(nu_u, nu, a2);
1591:     }
1592:     for (j = 0; j < 50000; ++j) {
1593:       y1 = CryerFunction(nu_u, nu, a1);
1594:       y2 = CryerFunction(nu_u, nu, a2);
1596:       am = (a1 + a2) / 2.0;
1597:       ym = CryerFunction(nu_u, nu, am);
1598:       if ((ym * y1) < 0) a2 = am;
1599:       else a1 = am;
1600:       if (PetscAbsReal(ym) < tol) break;
1601:     }
1603:     ctx->zeroArray[n - 1] = am;
1604:   }
1605:   return 0;
1606: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1841:     mandelZeros(PETSC_COMM_WORLD, user, param);

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

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

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

1876:     cryerZeros(PETSC_COMM_WORLD, user, param);

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

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

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

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

1911: static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace)
1912: {
1914:   DMPlexCreateRigidBody(dm, origField, nullspace);
1915:   return 0;
1916: }

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

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

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

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

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

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

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

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

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

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

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

2076:   PetscObjectGetOptions((PetscObject)ts, &options);
2077:   PetscObjectGetOptionsPrefix((PetscObject)ts, &prefix);
2078:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts), options, prefix, "-monitor_solution", &viewer, &format, &flg);
2079:   if (flg) TSMonitorSet(ts, SolutionMonitor, ctx, NULL);
2080:   PetscViewerDestroy(&viewer);
2081:   return 0;
2082: }

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

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

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

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

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

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

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

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

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

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

2191: /*TEST

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2334: TEST*/