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