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