Actual source code: ex53.c

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

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

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

     14: /* This presentation of poroelasticity is taken from

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

     24: For visualization, use

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

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

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

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

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

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

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

     78: /* Quadratic space and linear time solution

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

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

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

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

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

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

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

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

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

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

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

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

    174: /* Quadratic space and trigonometric time solution

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

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

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

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

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

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

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

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

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

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

    257: /* Trigonometric space and linear time solution

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    954:   return PETSC_SUCCESS;
    955: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    1027:   return PETSC_SUCCESS;
    1028: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    2192: /*TEST

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    2335: TEST*/