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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PetscReal 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 PetscReal c = PetscRealPart(kappa / S); /* m^2 / s, Cheng (B.16) */ 728 729 PetscReal A_s = 0.0; 730 PetscReal B_s = 0.0; 731 for (n = 1; n < N + 1; ++n) { 732 PetscReal alpha_n = user->zeroArray[n - 1]; 733 A_s += ((PetscSinReal(alpha_n) * PetscCosReal(alpha_n)) / (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscExpReal(-1 * (alpha_n * alpha_n * c * time) / (a * a)); 734 B_s += (PetscCosReal(alpha_n) / (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscSinReal((alpha_n * x[0]) / a) * PetscExpReal(-1 * (alpha_n * alpha_n * c * time) / (a * a)); 735 } 736 u[0] = ((P_0 * nu) / (2.0 * G * a) - (P_0 * nu_u) / (G * a) * A_s) * x[0] + P_0 / G * B_s; 737 u[1] = (-1 * (P_0 * (1.0 - nu)) / (2 * G * a) + (P_0 * (1 - nu_u)) / (G * a) * A_s) * x[1]; 738 } 739 return PETSC_SUCCESS; 740 } 741 742 static PetscErrorCode mandel_initial_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 743 { 744 AppCtx *user = (AppCtx *)ctx; 745 Parameter *param; 746 747 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 748 { 749 PetscScalar alpha = param->alpha; /* - */ 750 PetscScalar K_u = param->K_u; /* Pa */ 751 PetscScalar M = param->M; /* Pa */ 752 PetscScalar G = param->mu; /* Pa */ 753 PetscScalar P_0 = param->P_0; /* Pa */ 754 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */ 755 PetscReal a = 0.5 * (user->xmax[0] - user->xmin[0]); /* m */ 756 PetscInt N = user->niter, n; 757 758 PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */ 759 PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */ 760 PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */ 761 PetscReal c = PetscRealPart(kappa / S); /* m^2 / s, Cheng (B.16) */ 762 763 PetscReal aa = 0.0; 764 PetscReal eps_A = 0.0; 765 PetscReal eps_B = 0.0; 766 PetscReal eps_C = 0.0; 767 PetscReal time = 0.0; 768 769 for (n = 1; n < N + 1; ++n) { 770 aa = user->zeroArray[n - 1]; 771 eps_A += (aa * PetscExpReal((-1.0 * aa * aa * c * time) / (a * a)) * PetscCosReal(aa) * PetscCosReal((aa * x[0]) / a)) / (a * (aa - PetscSinReal(aa) * PetscCosReal(aa))); 772 eps_B += (PetscExpReal((-1.0 * aa * aa * c * time) / (a * a)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa)); 773 eps_C += (PetscExpReal((-1.0 * aa * aa * c * time) / (aa * aa)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa)); 774 } 775 u[0] = (P_0 / G) * eps_A + ((P_0 * nu) / (2.0 * G * a)) - eps_B / (G * a) - (P_0 * (1 - nu)) / (2 * G * a) + eps_C / (G * a); 776 } 777 return PETSC_SUCCESS; 778 } 779 780 // Displacement 781 static PetscErrorCode mandel_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 782 { 783 Parameter *param; 784 785 AppCtx *user = (AppCtx *)ctx; 786 787 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 788 if (time <= 0.0) { 789 PetscCall(mandel_initial_u(dim, time, x, Nc, u, ctx)); 790 } else { 791 PetscInt NITER = user->niter; 792 PetscScalar alpha = param->alpha; 793 PetscScalar K_u = param->K_u; 794 PetscScalar M = param->M; 795 PetscScalar G = param->mu; 796 PetscScalar k = param->k; 797 PetscScalar mu_f = param->mu_f; 798 PetscScalar F = param->P_0; 799 800 PetscScalar K_d = K_u - alpha * alpha * M; 801 PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); 802 PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); 803 PetscScalar kappa = k / mu_f; 804 PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0; 805 PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u))); 806 807 // Series term 808 PetscScalar A_x = 0.0; 809 PetscScalar B_x = 0.0; 810 811 for (PetscInt n = 1; n < NITER + 1; n++) { 812 PetscReal alpha_n = user->zeroArray[n - 1]; 813 814 A_x += ((PetscSinReal(alpha_n) * PetscCosReal(alpha_n)) / (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscExpReal(-1 * (alpha_n * alpha_n * c * time) / (a * a)); 815 B_x += (PetscCosReal(alpha_n) / (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscSinReal((alpha_n * x[0]) / a) * PetscExpReal(-1 * (alpha_n * alpha_n * c * time) / (a * a)); 816 } 817 u[0] = ((F * nu) / (2.0 * G * a) - (F * nu_u) / (G * a) * A_x) * x[0] + F / G * B_x; 818 u[1] = (-1 * (F * (1.0 - nu)) / (2 * G * a) + (F * (1 - nu_u)) / (G * a) * A_x) * x[1]; 819 } 820 return PETSC_SUCCESS; 821 } 822 823 // Trace strain 824 static PetscErrorCode mandel_2d_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 825 { 826 Parameter *param; 827 828 AppCtx *user = (AppCtx *)ctx; 829 830 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 831 if (time <= 0.0) { 832 PetscCall(mandel_initial_eps(dim, time, x, Nc, u, ctx)); 833 } else { 834 PetscInt NITER = user->niter; 835 PetscScalar alpha = param->alpha; 836 PetscScalar K_u = param->K_u; 837 PetscScalar M = param->M; 838 PetscScalar G = param->mu; 839 PetscScalar k = param->k; 840 PetscScalar mu_f = param->mu_f; 841 PetscScalar F = param->P_0; 842 843 PetscScalar K_d = K_u - alpha * alpha * M; 844 PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); 845 PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); 846 PetscScalar kappa = k / mu_f; 847 //const PetscScalar B = (alpha*M)/(K_d + alpha*alpha * M); 848 849 //const PetscScalar b = (YMAX - YMIN) / 2.0; 850 PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0; 851 PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u))); 852 853 // Series term 854 PetscReal eps_A = 0.0; 855 PetscReal eps_B = 0.0; 856 PetscReal eps_C = 0.0; 857 858 for (PetscInt n = 1; n < NITER + 1; n++) { 859 PetscReal aa = user->zeroArray[n - 1]; 860 861 eps_A += (aa * PetscExpReal((-1.0 * aa * aa * c * time) / (a * a)) * PetscCosReal(aa) * PetscCosReal((aa * x[0]) / a)) / (a * (aa - PetscSinReal(aa) * PetscCosReal(aa))); 862 863 eps_B += (PetscExpReal((-1.0 * aa * aa * c * time) / (a * a)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa)); 864 865 eps_C += (PetscExpReal((-1.0 * aa * aa * c * time) / (aa * aa)) * PetscSinReal(aa) * PetscCosReal(aa)) / (aa - PetscSinReal(aa) * PetscCosReal(aa)); 866 } 867 868 u[0] = (F / G) * eps_A + ((F * nu) / (2.0 * G * a)) - eps_B / (G * a) - (F * (1 - nu)) / (2 * G * a) + eps_C / (G * a); 869 } 870 return PETSC_SUCCESS; 871 } 872 873 // Pressure 874 static PetscErrorCode mandel_2d_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 875 { 876 Parameter *param; 877 878 AppCtx *user = (AppCtx *)ctx; 879 880 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 881 if (time <= 0.0) { 882 PetscCall(mandel_drainage_pressure(dim, time, x, Nc, u, ctx)); 883 } else { 884 PetscInt NITER = user->niter; 885 886 PetscScalar alpha = param->alpha; 887 PetscScalar K_u = param->K_u; 888 PetscScalar M = param->M; 889 PetscScalar G = param->mu; 890 PetscScalar k = param->k; 891 PetscScalar mu_f = param->mu_f; 892 PetscScalar F = param->P_0; 893 894 PetscScalar K_d = K_u - alpha * alpha * M; 895 PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); 896 PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); 897 PetscScalar kappa = k / mu_f; 898 PetscScalar B = (alpha * M) / (K_d + alpha * alpha * M); 899 900 PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0; 901 PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u))); 902 PetscScalar A1 = 3.0 / (B * (1.0 + nu_u)); 903 //PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu); 904 905 // Series term 906 PetscReal p = 0.0; 907 908 for (PetscInt n = 1; n < NITER + 1; n++) { 909 PetscReal aa = user->zeroArray[n - 1]; 910 p += (PetscSinReal(aa) / (aa - PetscSinReal(aa) * PetscCosReal(aa))) * (PetscCosReal((aa * x[0]) / a) - PetscCosReal(aa)) * PetscExpReal(-1.0 * (aa * aa * c * time) / (a * a)); 911 } 912 u[0] = ((2.0 * F) / (a * A1)) * p; 913 } 914 return PETSC_SUCCESS; 915 } 916 917 // Time derivative of displacement 918 static PetscErrorCode mandel_2d_u_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 919 { 920 Parameter *param; 921 922 AppCtx *user = (AppCtx *)ctx; 923 924 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 925 926 PetscInt NITER = user->niter; 927 PetscScalar alpha = param->alpha; 928 PetscScalar K_u = param->K_u; 929 PetscScalar M = param->M; 930 PetscScalar G = param->mu; 931 PetscScalar F = param->P_0; 932 933 PetscScalar K_d = K_u - alpha * alpha * M; 934 PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); 935 PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); 936 PetscScalar kappa = param->k / param->mu_f; 937 PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0; 938 PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u))); 939 940 // Series term 941 PetscScalar A_s_t = 0.0; 942 PetscScalar B_s_t = 0.0; 943 944 for (PetscInt n = 1; n < NITER + 1; n++) { 945 PetscReal alpha_n = user->zeroArray[n - 1]; 946 947 A_s_t += (-1.0 * alpha_n * alpha_n * c * PetscExpReal((-1.0 * alpha_n * alpha_n * time) / (a * a)) * PetscSinReal((alpha_n * x[0]) / a) * PetscCosReal(alpha_n)) / (a * a * (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))); 948 B_s_t += (-1.0 * alpha_n * alpha_n * c * PetscExpReal((-1.0 * alpha_n * alpha_n * time) / (a * a)) * PetscSinReal(alpha_n) * PetscCosReal(alpha_n)) / (a * a * (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))); 949 } 950 951 u[0] = (F / G) * A_s_t - ((F * nu_u * x[0]) / (G * a)) * B_s_t; 952 u[1] = ((F * x[1] * (1 - nu_u)) / (G * a)) * B_s_t; 953 954 return PETSC_SUCCESS; 955 } 956 957 // Time derivative of trace strain 958 static PetscErrorCode mandel_2d_eps_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 959 { 960 Parameter *param; 961 962 AppCtx *user = (AppCtx *)ctx; 963 964 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 965 966 PetscInt NITER = user->niter; 967 PetscScalar alpha = param->alpha; 968 PetscScalar K_u = param->K_u; 969 PetscScalar M = param->M; 970 PetscScalar G = param->mu; 971 PetscScalar k = param->k; 972 PetscScalar mu_f = param->mu_f; 973 PetscScalar F = param->P_0; 974 975 PetscScalar K_d = K_u - alpha * alpha * M; 976 PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); 977 PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); 978 PetscScalar kappa = k / mu_f; 979 //const PetscScalar B = (alpha*M)/(K_d + alpha*alpha * M); 980 981 //const PetscScalar b = (YMAX - YMIN) / 2.0; 982 PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0; 983 PetscReal c = PetscRealPart(((2.0 * kappa * G) * (1.0 - nu) * (nu_u - nu)) / (alpha * alpha * (1.0 - 2.0 * nu) * (1.0 - nu_u))); 984 985 // Series term 986 PetscScalar eps_As = 0.0; 987 PetscScalar eps_Bs = 0.0; 988 PetscScalar eps_Cs = 0.0; 989 990 for (PetscInt n = 1; n < NITER + 1; n++) { 991 PetscReal alpha_n = user->zeroArray[n - 1]; 992 993 eps_As += (-1.0 * alpha_n * alpha_n * alpha_n * c * PetscExpReal((-1.0 * alpha_n * alpha_n * c * time) / (a * a)) * PetscCosReal(alpha_n) * PetscCosReal((alpha_n * x[0]) / a)) / (alpha_n * alpha_n * alpha_n * (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))); 994 eps_Bs += (-1.0 * alpha_n * alpha_n * c * PetscExpReal((-1.0 * alpha_n * alpha_n * c * time) / (a * a)) * PetscSinReal(alpha_n) * PetscCosReal(alpha_n)) / (alpha_n * alpha_n * (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))); 995 eps_Cs += (-1.0 * alpha_n * alpha_n * c * PetscExpReal((-1.0 * alpha_n * alpha_n * c * time) / (a * a)) * PetscSinReal(alpha_n) * PetscCosReal(alpha_n)) / (alpha_n * alpha_n * (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))); 996 } 997 998 u[0] = (F / G) * eps_As - ((F * nu_u) / (G * a)) * eps_Bs + ((F * (1 - nu_u)) / (G * a)) * eps_Cs; 999 return PETSC_SUCCESS; 1000 } 1001 1002 // Time derivative of pressure 1003 static PetscErrorCode mandel_2d_p_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 1004 { 1005 Parameter *param; 1006 1007 AppCtx *user = (AppCtx *)ctx; 1008 1009 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 1010 1011 PetscScalar alpha = param->alpha; 1012 PetscScalar K_u = param->K_u; 1013 PetscScalar M = param->M; 1014 PetscScalar G = param->mu; 1015 PetscScalar F = param->P_0; 1016 1017 PetscScalar K_d = K_u - alpha * alpha * M; 1018 PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); 1019 PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); 1020 1021 PetscReal a = (user->xmax[0] - user->xmin[0]) / 2.0; 1022 //PetscScalar A1 = 3.0 / (B * (1.0 + nu_u)); 1023 //PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu); 1024 1025 u[0] = ((2.0 * F * (-2.0 * nu + 3.0 * nu_u)) / (3.0 * a * alpha * (1.0 - 2.0 * nu))); 1026 1027 return PETSC_SUCCESS; 1028 } 1029 1030 /* Cryer Solutions */ 1031 static PetscErrorCode cryer_drainage_pressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 1032 { 1033 AppCtx *user = (AppCtx *)ctx; 1034 Parameter *param; 1035 1036 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 1037 if (time <= 0.0) { 1038 PetscScalar alpha = param->alpha; /* - */ 1039 PetscScalar K_u = param->K_u; /* Pa */ 1040 PetscScalar M = param->M; /* Pa */ 1041 PetscScalar P_0 = param->P_0; /* Pa */ 1042 PetscScalar B = alpha * M / K_u; /* -, Cheng (B.12) */ 1043 1044 u[0] = P_0 * B; 1045 } else { 1046 u[0] = 0.0; 1047 } 1048 return PETSC_SUCCESS; 1049 } 1050 1051 static PetscErrorCode cryer_initial_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 1052 { 1053 AppCtx *user = (AppCtx *)ctx; 1054 Parameter *param; 1055 1056 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 1057 { 1058 PetscScalar K_u = param->K_u; /* Pa */ 1059 PetscScalar G = param->mu; /* Pa */ 1060 PetscScalar P_0 = param->P_0; /* Pa */ 1061 PetscReal R_0 = user->xmax[1]; /* m */ 1062 PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */ 1063 1064 PetscScalar u_0 = -P_0 * R_0 * (1. - 2. * nu_u) / (2. * G * (1. + nu_u)); /* Cheng (7.407) */ 1065 PetscReal u_sc = PetscRealPart(u_0) / R_0; 1066 1067 u[0] = u_sc * x[0]; 1068 u[1] = u_sc * x[1]; 1069 u[2] = u_sc * x[2]; 1070 } 1071 return PETSC_SUCCESS; 1072 } 1073 1074 static PetscErrorCode cryer_initial_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 1075 { 1076 AppCtx *user = (AppCtx *)ctx; 1077 Parameter *param; 1078 1079 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 1080 { 1081 PetscScalar K_u = param->K_u; /* Pa */ 1082 PetscScalar G = param->mu; /* Pa */ 1083 PetscScalar P_0 = param->P_0; /* Pa */ 1084 PetscReal R_0 = user->xmax[1]; /* m */ 1085 PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */ 1086 1087 PetscScalar u_0 = -P_0 * R_0 * (1. - 2. * nu_u) / (2. * G * (1. + nu_u)); /* Cheng (7.407) */ 1088 PetscReal u_sc = PetscRealPart(u_0) / R_0; 1089 1090 /* div R = 1/R^2 d/dR R^2 R = 3 */ 1091 u[0] = 3. * u_sc; 1092 u[1] = 3. * u_sc; 1093 u[2] = 3. * u_sc; 1094 } 1095 return PETSC_SUCCESS; 1096 } 1097 1098 // Displacement 1099 static PetscErrorCode cryer_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 1100 { 1101 AppCtx *user = (AppCtx *)ctx; 1102 Parameter *param; 1103 1104 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 1105 if (time <= 0.0) { 1106 PetscCall(cryer_initial_u(dim, time, x, Nc, u, ctx)); 1107 } else { 1108 PetscScalar alpha = param->alpha; /* - */ 1109 PetscScalar K_u = param->K_u; /* Pa */ 1110 PetscScalar M = param->M; /* Pa */ 1111 PetscScalar G = param->mu; /* Pa */ 1112 PetscScalar P_0 = param->P_0; /* Pa */ 1113 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */ 1114 PetscReal R_0 = user->xmax[1]; /* m */ 1115 PetscInt N = user->niter, n; 1116 1117 PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */ 1118 PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */ 1119 PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */ 1120 PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */ 1121 PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */ 1122 PetscScalar u_inf = -P_0 * R_0 * (1. - 2. * nu) / (2. * G * (1. + nu)); /* m, Cheng (7.388) */ 1123 1124 PetscReal R = PetscSqrtReal(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]); 1125 PetscReal R_star = R / R_0; 1126 PetscReal tstar = PetscRealPart(c * time) / PetscSqr(R_0); /* - */ 1127 PetscReal A_n = 0.0; 1128 PetscScalar u_sc; 1129 1130 for (n = 1; n < N + 1; ++n) { 1131 const PetscReal x_n = user->zeroArray[n - 1]; 1132 const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u)); 1133 1134 /* m , Cheng (7.404) */ 1135 if (R_star != 0) { 1136 A_n += PetscRealPart((12.0 * (1.0 + nu) * (nu_u - nu)) / ((1.0 - 2.0 * nu) * E_n * PetscSqr(R_star) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * (3.0 * (nu_u - nu) * (PetscSinReal(R_star * PetscSqrtReal(x_n)) - R_star * PetscSqrtReal(x_n) * PetscCosReal(R_star * PetscSqrtReal(x_n))) + (1.0 - nu) * (1.0 - 2.0 * nu) * PetscPowRealInt(R_star, 3) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * PetscExpReal(-x_n * tstar)); 1137 } 1138 } 1139 if (R_star != 0) u_sc = PetscRealPart(u_inf) * (R_star - A_n) / R; 1140 else u_sc = PetscRealPart(u_inf) / R_0; 1141 u[0] = u_sc * x[0]; 1142 u[1] = u_sc * x[1]; 1143 u[2] = u_sc * x[2]; 1144 } 1145 return PETSC_SUCCESS; 1146 } 1147 1148 // Volumetric Strain 1149 static PetscErrorCode cryer_3d_eps(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 1150 { 1151 AppCtx *user = (AppCtx *)ctx; 1152 Parameter *param; 1153 1154 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 1155 if (time <= 0.0) { 1156 PetscCall(cryer_initial_eps(dim, time, x, Nc, u, ctx)); 1157 } else { 1158 PetscScalar alpha = param->alpha; /* - */ 1159 PetscScalar K_u = param->K_u; /* Pa */ 1160 PetscScalar M = param->M; /* Pa */ 1161 PetscScalar G = param->mu; /* Pa */ 1162 PetscScalar P_0 = param->P_0; /* Pa */ 1163 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */ 1164 PetscReal R_0 = user->xmax[1]; /* m */ 1165 PetscInt N = user->niter, n; 1166 1167 PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */ 1168 PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */ 1169 PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */ 1170 PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */ 1171 PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */ 1172 PetscScalar u_inf = -P_0 * R_0 * (1. - 2. * nu) / (2. * G * (1. + nu)); /* m, Cheng (7.388) */ 1173 1174 PetscReal R = PetscSqrtReal(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]); 1175 PetscReal R_star = R / R_0; 1176 PetscReal tstar = PetscRealPart(c * time) / PetscSqr(R_0); /* - */ 1177 PetscReal divA_n = 0.0; 1178 1179 if (R_star < PETSC_SMALL) { 1180 for (n = 1; n < N + 1; ++n) { 1181 const PetscReal x_n = user->zeroArray[n - 1]; 1182 const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u)); 1183 1184 divA_n += PetscRealPart((12.0 * (1.0 + nu) * (nu_u - nu)) / ((1.0 - 2.0 * nu) * E_n * PetscSqr(R_star) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * (3.0 * (nu_u - nu) * PetscSqrtReal(x_n) * ((2.0 + PetscSqr(R_star * PetscSqrtReal(x_n))) - 2.0 * PetscCosReal(R_star * PetscSqrtReal(x_n))) + 5.0 * (1.0 - nu) * (1.0 - 2.0 * nu) * PetscPowRealInt(R_star, 2) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * PetscExpReal(-x_n * tstar)); 1185 } 1186 } else { 1187 for (n = 1; n < N + 1; ++n) { 1188 const PetscReal x_n = user->zeroArray[n - 1]; 1189 const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u)); 1190 1191 divA_n += PetscRealPart((12.0 * (1.0 + nu) * (nu_u - nu)) / ((1.0 - 2.0 * nu) * E_n * PetscSqr(R_star) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * (3.0 * (nu_u - nu) * PetscSqrtReal(x_n) * ((2.0 / (R_star * PetscSqrtReal(x_n)) + R_star * PetscSqrtReal(x_n)) * PetscSinReal(R_star * PetscSqrtReal(x_n)) - 2.0 * PetscCosReal(R_star * PetscSqrtReal(x_n))) + 5.0 * (1.0 - nu) * (1.0 - 2.0 * nu) * PetscPowRealInt(R_star, 2) * x_n * PetscSinReal(PetscSqrtReal(x_n))) * PetscExpReal(-x_n * tstar)); 1192 } 1193 } 1194 if (PetscAbsReal(divA_n) > 1e3) PetscCall(PetscPrintf(PETSC_COMM_SELF, "(%g, %g, %g) divA_n: %g\n", (double)x[0], (double)x[1], (double)x[2], (double)divA_n)); 1195 u[0] = PetscRealPart(u_inf) / R_0 * (3.0 - divA_n); 1196 } 1197 return PETSC_SUCCESS; 1198 } 1199 1200 // Pressure 1201 static PetscErrorCode cryer_3d_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 1202 { 1203 AppCtx *user = (AppCtx *)ctx; 1204 Parameter *param; 1205 1206 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 1207 if (time <= 0.0) { 1208 PetscCall(cryer_drainage_pressure(dim, time, x, Nc, u, ctx)); 1209 } else { 1210 PetscScalar alpha = param->alpha; /* - */ 1211 PetscScalar K_u = param->K_u; /* Pa */ 1212 PetscScalar M = param->M; /* Pa */ 1213 PetscScalar G = param->mu; /* Pa */ 1214 PetscScalar P_0 = param->P_0; /* Pa */ 1215 PetscReal R_0 = user->xmax[1]; /* m */ 1216 PetscScalar kappa = param->k / param->mu_f; /* m^2 / (Pa s) */ 1217 PetscInt N = user->niter, n; 1218 1219 PetscScalar K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */ 1220 PetscScalar eta = (3.0 * alpha * G) / (3.0 * K_d + 4.0 * G); /* -, Cheng (B.11) */ 1221 PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */ 1222 PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */ 1223 PetscScalar S = (3.0 * K_u + 4.0 * G) / (M * (3.0 * K_d + 4.0 * G)); /* Pa^{-1}, Cheng (B.14) */ 1224 PetscScalar c = kappa / S; /* m^2 / s, Cheng (B.16) */ 1225 PetscReal R = PetscSqrtReal(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]); 1226 1227 PetscReal R_star = R / R_0; 1228 PetscReal t_star = PetscRealPart(c * time) / PetscSqr(R_0); 1229 PetscReal A_x = 0.0; 1230 1231 for (n = 1; n < N + 1; ++n) { 1232 const PetscReal x_n = user->zeroArray[n - 1]; 1233 const PetscReal E_n = PetscRealPart(PetscSqr(1 - nu) * PetscSqr(1 + nu_u) * x_n - 18.0 * (1 + nu) * (nu_u - nu) * (1 - nu_u)); 1234 1235 A_x += PetscRealPart(((18.0 * PetscSqr(nu_u - nu)) / (eta * E_n)) * (PetscSinReal(R_star * PetscSqrtReal(x_n)) / (R_star * PetscSinReal(PetscSqrtReal(x_n))) - 1.0) * PetscExpReal(-x_n * t_star)); /* Cheng (7.395) */ 1236 } 1237 u[0] = P_0 * A_x; 1238 } 1239 return PETSC_SUCCESS; 1240 } 1241 1242 /* Boundary Kernels */ 1243 static void f0_terzaghi_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1244 { 1245 const PetscReal P = PetscRealPart(constants[5]); 1246 1247 f0[0] = 0.0; 1248 f0[1] = P; 1249 } 1250 1251 #if 0 1252 static void f0_mandel_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1253 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1254 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1255 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1256 { 1257 // Uniform stress distribution 1258 /* PetscScalar xmax = 0.5; 1259 PetscScalar xmin = -0.5; 1260 PetscScalar ymax = 0.5; 1261 PetscScalar ymin = -0.5; 1262 PetscScalar P = constants[5]; 1263 PetscScalar aL = (xmax - xmin) / 2.0; 1264 PetscScalar sigma_zz = -1.0*P / aL; */ 1265 1266 // Analytical (parabolic) stress distribution 1267 PetscReal a1, a2, am; 1268 PetscReal y1, y2, ym; 1269 1270 PetscInt NITER = 500; 1271 PetscReal EPS = 0.000001; 1272 PetscReal zeroArray[500]; /* NITER */ 1273 PetscReal xmax = 1.0; 1274 PetscReal xmin = 0.0; 1275 PetscReal ymax = 0.1; 1276 PetscReal ymin = 0.0; 1277 PetscReal lower[2], upper[2]; 1278 1279 lower[0] = xmin - (xmax - xmin) / 2.0; 1280 lower[1] = ymin - (ymax - ymin) / 2.0; 1281 upper[0] = xmax - (xmax - xmin) / 2.0; 1282 upper[1] = ymax - (ymax - ymin) / 2.0; 1283 1284 xmin = lower[0]; 1285 ymin = lower[1]; 1286 xmax = upper[0]; 1287 ymax = upper[1]; 1288 1289 PetscScalar G = constants[0]; 1290 PetscScalar K_u = constants[1]; 1291 PetscScalar alpha = constants[2]; 1292 PetscScalar M = constants[3]; 1293 PetscScalar kappa = constants[4]; 1294 PetscScalar F = constants[5]; 1295 1296 PetscScalar K_d = K_u - alpha*alpha*M; 1297 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G)); 1298 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); 1299 PetscReal aL = (xmax - xmin) / 2.0; 1300 PetscReal c = PetscRealPart(((2.0*kappa*G) * (1.0 - nu) * (nu_u - nu)) / (alpha*alpha * (1.0 - 2.0*nu) * (1.0 - nu_u))); 1301 PetscScalar B = (3.0 * (nu_u - nu)) / ( alpha * (1.0 - 2.0*nu) * (1.0 + nu_u)); 1302 PetscScalar A1 = 3.0 / (B * (1.0 + nu_u)); 1303 PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu); 1304 1305 // Generate zero values 1306 for (PetscInt i=1; i < NITER+1; i++) 1307 { 1308 a1 = ((PetscReal) i - 1.0) * PETSC_PI * PETSC_PI / 4.0 + EPS; 1309 a2 = a1 + PETSC_PI/2; 1310 for (PetscInt j=0; j<NITER; j++) 1311 { 1312 y1 = PetscTanReal(a1) - PetscRealPart(A1/A2)*a1; 1313 y2 = PetscTanReal(a2) - PetscRealPart(A1/A2)*a2; 1314 am = (a1 + a2)/2.0; 1315 ym = PetscTanReal(am) - PetscRealPart(A1/A2)*am; 1316 if ((ym*y1) > 0) 1317 { 1318 a1 = am; 1319 } else { 1320 a2 = am; 1321 } 1322 if (PetscAbsReal(y2) < EPS) 1323 { 1324 am = a2; 1325 } 1326 } 1327 zeroArray[i-1] = am; 1328 } 1329 1330 // Solution for sigma_zz 1331 PetscScalar A_x = 0.0; 1332 PetscScalar B_x = 0.0; 1333 1334 for (PetscInt n=1; n < NITER+1; n++) 1335 { 1336 PetscReal alpha_n = zeroArray[n-1]; 1337 1338 A_x += ( PetscSinReal(alpha_n) / (alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscCosReal( (alpha_n * x[0]) / aL) * PetscExpReal( -1.0*( (alpha_n*alpha_n*c*t)/(aL*aL))); 1339 B_x += ( (PetscSinReal(alpha_n) * PetscCosReal(alpha_n))/(alpha_n - PetscSinReal(alpha_n) * PetscCosReal(alpha_n))) * PetscExpReal( -1.0*( (alpha_n*alpha_n*c*t)/(aL*aL))); 1340 } 1341 1342 PetscScalar sigma_zz = -1.0*(F/aL) - ((2.0*F)/aL) * (A2/A1) * A_x + ((2.0*F)/aL) * B_x; 1343 1344 if (x[1] == ymax) { 1345 f0[1] += sigma_zz; 1346 } else if (x[1] == ymin) { 1347 f0[1] -= sigma_zz; 1348 } 1349 } 1350 #endif 1351 1352 static void f0_cryer_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1353 { 1354 const PetscReal P_0 = PetscRealPart(constants[5]); 1355 PetscInt d; 1356 1357 for (d = 0; d < dim; ++d) f0[d] = -P_0 * n[d]; 1358 } 1359 1360 /* Standard Kernels - Residual */ 1361 /* f0_e */ 1362 static void f0_epsilon(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1363 { 1364 PetscInt d; 1365 1366 for (d = 0; d < dim; ++d) f0[0] += u_x[d * dim + d]; 1367 f0[0] -= u[uOff[1]]; 1368 } 1369 1370 /* f0_p */ 1371 static void f0_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1372 { 1373 const PetscReal alpha = PetscRealPart(constants[2]); 1374 const PetscReal M = PetscRealPart(constants[3]); 1375 1376 f0[0] += alpha * u_t[uOff[1]]; 1377 f0[0] += u_t[uOff[2]] / M; 1378 if (f0[0] != f0[0]) abort(); 1379 } 1380 1381 /* f1_u */ 1382 static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 1383 { 1384 const PetscInt Nc = dim; 1385 const PetscReal G = PetscRealPart(constants[0]); 1386 const PetscReal K_u = PetscRealPart(constants[1]); 1387 const PetscReal alpha = PetscRealPart(constants[2]); 1388 const PetscReal M = PetscRealPart(constants[3]); 1389 const PetscReal K_d = K_u - alpha * alpha * M; 1390 const PetscReal lambda = K_d - (2.0 * G) / 3.0; 1391 PetscInt c, d; 1392 1393 for (c = 0; c < Nc; ++c) { 1394 for (d = 0; d < dim; ++d) f1[c * dim + d] -= G * (u_x[c * dim + d] + u_x[d * dim + c]); 1395 f1[c * dim + c] -= lambda * u[uOff[1]]; 1396 f1[c * dim + c] += alpha * u[uOff[2]]; 1397 } 1398 } 1399 1400 /* f1_p */ 1401 static void f1_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 1402 { 1403 const PetscReal kappa = PetscRealPart(constants[4]); 1404 PetscInt d; 1405 1406 for (d = 0; d < dim; ++d) f1[d] += kappa * u_x[uOff_x[2] + d]; 1407 } 1408 1409 /* 1410 \partial_df \phi_fc g_{fc,gc,df,dg} \partial_dg \phi_gc 1411 1412 \partial_df \phi_fc \lambda \delta_{fc,df} \sum_gc \partial_dg \phi_gc \delta_{gc,dg} 1413 = \partial_fc \phi_fc \sum_gc \partial_gc \phi_gc 1414 */ 1415 1416 /* Standard Kernels - Jacobian */ 1417 /* g0_ee */ 1418 static void g0_ee(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 1419 { 1420 g0[0] = -1.0; 1421 } 1422 1423 /* g0_pe */ 1424 static void g0_pe(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 1425 { 1426 const PetscReal alpha = PetscRealPart(constants[2]); 1427 1428 g0[0] = u_tShift * alpha; 1429 } 1430 1431 /* g0_pp */ 1432 static void g0_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 1433 { 1434 const PetscReal M = PetscRealPart(constants[3]); 1435 1436 g0[0] = u_tShift / M; 1437 } 1438 1439 /* g1_eu */ 1440 static void g1_eu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 1441 { 1442 PetscInt d; 1443 for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */ 1444 } 1445 1446 /* g2_ue */ 1447 static void g2_ue(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]) 1448 { 1449 const PetscReal G = PetscRealPart(constants[0]); 1450 const PetscReal K_u = PetscRealPart(constants[1]); 1451 const PetscReal alpha = PetscRealPart(constants[2]); 1452 const PetscReal M = PetscRealPart(constants[3]); 1453 const PetscReal K_d = K_u - alpha * alpha * M; 1454 const PetscReal lambda = K_d - (2.0 * G) / 3.0; 1455 PetscInt d; 1456 1457 for (d = 0; d < dim; ++d) g2[d * dim + d] -= lambda; 1458 } 1459 /* g2_up */ 1460 static void g2_up(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]) 1461 { 1462 const PetscReal alpha = PetscRealPart(constants[2]); 1463 PetscInt d; 1464 1465 for (d = 0; d < dim; ++d) g2[d * dim + d] += alpha; 1466 } 1467 1468 /* g3_uu */ 1469 static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 1470 { 1471 const PetscInt Nc = dim; 1472 const PetscReal G = PetscRealPart(constants[0]); 1473 PetscInt c, d; 1474 1475 for (c = 0; c < Nc; ++c) { 1476 for (d = 0; d < dim; ++d) { 1477 g3[((c * Nc + c) * dim + d) * dim + d] -= G; 1478 g3[((c * Nc + d) * dim + d) * dim + c] -= G; 1479 } 1480 } 1481 } 1482 1483 /* g3_pp */ 1484 static void g3_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 1485 { 1486 const PetscReal kappa = PetscRealPart(constants[4]); 1487 PetscInt d; 1488 1489 for (d = 0; d < dim; ++d) g3[d * dim + d] += kappa; 1490 } 1491 1492 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 1493 { 1494 PetscInt sol; 1495 1496 PetscFunctionBeginUser; 1497 options->solType = SOL_QUADRATIC_TRIG; 1498 options->niter = 500; 1499 options->eps = PETSC_SMALL; 1500 options->dtInitial = -1.0; 1501 PetscOptionsBegin(comm, "", "Biot Poroelasticity Options", "DMPLEX"); 1502 PetscCall(PetscOptionsInt("-niter", "Number of series term iterations in exact solutions", "ex53.c", options->niter, &options->niter, NULL)); 1503 sol = options->solType; 1504 PetscCall(PetscOptionsEList("-sol_type", "Type of exact solution", "ex53.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL)); 1505 options->solType = (SolutionType)sol; 1506 PetscCall(PetscOptionsReal("-eps", "Precision value for root finding", "ex53.c", options->eps, &options->eps, NULL)); 1507 PetscCall(PetscOptionsReal("-dt_initial", "Override the initial timestep", "ex53.c", options->dtInitial, &options->dtInitial, NULL)); 1508 PetscOptionsEnd(); 1509 PetscFunctionReturn(PETSC_SUCCESS); 1510 } 1511 1512 static PetscErrorCode mandelZeros(MPI_Comm comm, AppCtx *ctx, Parameter *param) 1513 { 1514 //PetscBag bag; 1515 PetscReal a1, a2, am; 1516 PetscReal y1, y2, ym; 1517 1518 PetscFunctionBeginUser; 1519 //PetscCall(PetscBagGetData(ctx->bag, (void **) ¶m)); 1520 PetscInt NITER = ctx->niter; 1521 PetscReal EPS = ctx->eps; 1522 //const PetscScalar YMAX = param->ymax; 1523 //const PetscScalar YMIN = param->ymin; 1524 PetscScalar alpha = param->alpha; 1525 PetscScalar K_u = param->K_u; 1526 PetscScalar M = param->M; 1527 PetscScalar G = param->mu; 1528 //const PetscScalar k = param->k; 1529 //const PetscScalar mu_f = param->mu_f; 1530 //const PetscScalar P_0 = param->P_0; 1531 1532 PetscScalar K_d = K_u - alpha * alpha * M; 1533 PetscScalar nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); 1534 PetscScalar nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); 1535 //const PetscScalar kappa = k / mu_f; 1536 1537 // Generate zero values 1538 for (PetscInt i = 1; i < NITER + 1; i++) { 1539 a1 = ((PetscReal)i - 1.0) * PETSC_PI * PETSC_PI / 4.0 + EPS; 1540 a2 = a1 + PETSC_PI / 2; 1541 am = a1; 1542 for (PetscInt j = 0; j < NITER; j++) { 1543 y1 = PetscTanReal(a1) - PetscRealPart((1.0 - nu) / (nu_u - nu)) * a1; 1544 y2 = PetscTanReal(a2) - PetscRealPart((1.0 - nu) / (nu_u - nu)) * a2; 1545 am = (a1 + a2) / 2.0; 1546 ym = PetscTanReal(am) - PetscRealPart((1.0 - nu) / (nu_u - nu)) * am; 1547 if ((ym * y1) > 0) { 1548 a1 = am; 1549 } else { 1550 a2 = am; 1551 } 1552 if (PetscAbsReal(y2) < EPS) am = a2; 1553 } 1554 ctx->zeroArray[i - 1] = am; 1555 } 1556 PetscFunctionReturn(PETSC_SUCCESS); 1557 } 1558 1559 static PetscReal CryerFunction(PetscReal nu_u, PetscReal nu, PetscReal x) 1560 { 1561 return PetscTanReal(PetscSqrtReal(x)) * (6.0 * (nu_u - nu) - (1.0 - nu) * (1.0 + nu_u) * x) - (6.0 * (nu_u - nu) * PetscSqrtReal(x)); 1562 } 1563 1564 static PetscErrorCode cryerZeros(MPI_Comm comm, AppCtx *ctx, Parameter *param) 1565 { 1566 PetscReal alpha = PetscRealPart(param->alpha); /* - */ 1567 PetscReal K_u = PetscRealPart(param->K_u); /* Pa */ 1568 PetscReal M = PetscRealPart(param->M); /* Pa */ 1569 PetscReal G = PetscRealPart(param->mu); /* Pa */ 1570 PetscInt N = ctx->niter, n; 1571 1572 PetscReal K_d = K_u - alpha * alpha * M; /* Pa, Cheng (B.5) */ 1573 PetscReal nu = (3.0 * K_d - 2.0 * G) / (2.0 * (3.0 * K_d + G)); /* -, Cheng (B.8) */ 1574 PetscReal nu_u = (3.0 * K_u - 2.0 * G) / (2.0 * (3.0 * K_u + G)); /* -, Cheng (B.9) */ 1575 1576 PetscFunctionBeginUser; 1577 for (n = 1; n < N + 1; ++n) { 1578 PetscReal tol = PetscPowReal(n, 1.5) * ctx->eps; 1579 PetscReal a1 = 0., a2 = 0., am = 0.; 1580 PetscReal y1, y2, ym; 1581 PetscInt j, k = n - 1; 1582 1583 y1 = y2 = 1.; 1584 while (y1 * y2 > 0) { 1585 ++k; 1586 a1 = PetscSqr(n * PETSC_PI) - k * PETSC_PI; 1587 a2 = PetscSqr(n * PETSC_PI) + k * PETSC_PI; 1588 y1 = CryerFunction(nu_u, nu, a1); 1589 y2 = CryerFunction(nu_u, nu, a2); 1590 } 1591 for (j = 0; j < 50000; ++j) { 1592 y1 = CryerFunction(nu_u, nu, a1); 1593 y2 = CryerFunction(nu_u, nu, a2); 1594 PetscCheck(y1 * y2 <= 0, comm, PETSC_ERR_PLIB, "Invalid root finding initialization for root %" PetscInt_FMT ", (%g, %g)--(%g, %g)", n, (double)a1, (double)y1, (double)a2, (double)y2); 1595 am = (a1 + a2) / 2.0; 1596 ym = CryerFunction(nu_u, nu, am); 1597 if ((ym * y1) < 0) a2 = am; 1598 else a1 = am; 1599 if (PetscAbsReal(ym) < tol) break; 1600 } 1601 PetscCheck(PetscAbsReal(ym) < tol, comm, PETSC_ERR_PLIB, "Root finding did not converge for root %" PetscInt_FMT " (%g)", n, (double)PetscAbsReal(ym)); 1602 ctx->zeroArray[n - 1] = am; 1603 } 1604 PetscFunctionReturn(PETSC_SUCCESS); 1605 } 1606 1607 static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx) 1608 { 1609 PetscBag bag; 1610 Parameter *p; 1611 1612 PetscFunctionBeginUser; 1613 /* setup PETSc parameter bag */ 1614 PetscCall(PetscBagGetData(ctx->bag, (void **)&p)); 1615 PetscCall(PetscBagSetName(ctx->bag, "par", "Poroelastic Parameters")); 1616 bag = ctx->bag; 1617 if (ctx->solType == SOL_TERZAGHI) { 1618 // Realistic values - Terzaghi 1619 PetscCall(PetscBagRegisterScalar(bag, &p->mu, 3.0, "mu", "Shear Modulus, Pa")); 1620 PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 9.76, "K_u", "Undrained Bulk Modulus, Pa")); 1621 PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -")); 1622 PetscCall(PetscBagRegisterScalar(bag, &p->M, 16.0, "M", "Biot Modulus, Pa")); 1623 PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2")); 1624 PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s")); 1625 PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa")); 1626 } else if (ctx->solType == SOL_MANDEL) { 1627 // Realistic values - Mandel 1628 PetscCall(PetscBagRegisterScalar(bag, &p->mu, 0.75, "mu", "Shear Modulus, Pa")); 1629 PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 2.6941176470588233, "K_u", "Undrained Bulk Modulus, Pa")); 1630 PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -")); 1631 PetscCall(PetscBagRegisterScalar(bag, &p->M, 4.705882352941176, "M", "Biot Modulus, Pa")); 1632 PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2")); 1633 PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s")); 1634 PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa")); 1635 } else if (ctx->solType == SOL_CRYER) { 1636 // Realistic values - Mandel 1637 PetscCall(PetscBagRegisterScalar(bag, &p->mu, 0.75, "mu", "Shear Modulus, Pa")); 1638 PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 2.6941176470588233, "K_u", "Undrained Bulk Modulus, Pa")); 1639 PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -")); 1640 PetscCall(PetscBagRegisterScalar(bag, &p->M, 4.705882352941176, "M", "Biot Modulus, Pa")); 1641 PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2")); 1642 PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s")); 1643 PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa")); 1644 } else { 1645 // Nonsense values 1646 PetscCall(PetscBagRegisterScalar(bag, &p->mu, 1.0, "mu", "Shear Modulus, Pa")); 1647 PetscCall(PetscBagRegisterScalar(bag, &p->K_u, 1.0, "K_u", "Undrained Bulk Modulus, Pa")); 1648 PetscCall(PetscBagRegisterScalar(bag, &p->alpha, 1.0, "alpha", "Biot Effective Stress Coefficient, -")); 1649 PetscCall(PetscBagRegisterScalar(bag, &p->M, 1.0, "M", "Biot Modulus, Pa")); 1650 PetscCall(PetscBagRegisterScalar(bag, &p->k, 1.0, "k", "Isotropic Permeability, m**2")); 1651 PetscCall(PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s")); 1652 PetscCall(PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa")); 1653 } 1654 PetscCall(PetscBagSetFromOptions(bag)); 1655 { 1656 PetscScalar K_d = p->K_u - p->alpha * p->alpha * p->M; 1657 PetscScalar nu_u = (3.0 * p->K_u - 2.0 * p->mu) / (2.0 * (3.0 * p->K_u + p->mu)); 1658 PetscScalar nu = (3.0 * K_d - 2.0 * p->mu) / (2.0 * (3.0 * K_d + p->mu)); 1659 PetscScalar S = (3.0 * p->K_u + 4.0 * p->mu) / (p->M * (3.0 * K_d + 4.0 * p->mu)); 1660 PetscReal c = PetscRealPart((p->k / p->mu_f) / S); 1661 1662 PetscViewer viewer; 1663 PetscViewerFormat format; 1664 PetscBool flg; 1665 1666 switch (ctx->solType) { 1667 case SOL_QUADRATIC_LINEAR: 1668 case SOL_QUADRATIC_TRIG: 1669 case SOL_TRIG_LINEAR: 1670 ctx->t_r = PetscSqr(ctx->xmax[0] - ctx->xmin[0]) / c; 1671 break; 1672 case SOL_TERZAGHI: 1673 ctx->t_r = PetscSqr(2.0 * (ctx->xmax[1] - ctx->xmin[1])) / c; 1674 break; 1675 case SOL_MANDEL: 1676 ctx->t_r = PetscSqr(2.0 * (ctx->xmax[1] - ctx->xmin[1])) / c; 1677 break; 1678 case SOL_CRYER: 1679 ctx->t_r = PetscSqr(ctx->xmax[1]) / c; 1680 break; 1681 default: 1682 SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%d)", solutionTypes[PetscMin(ctx->solType, NUM_SOLUTION_TYPES)], ctx->solType); 1683 } 1684 PetscCall(PetscOptionsGetViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg)); 1685 if (flg) { 1686 PetscCall(PetscViewerPushFormat(viewer, format)); 1687 PetscCall(PetscBagView(bag, viewer)); 1688 PetscCall(PetscViewerFlush(viewer)); 1689 PetscCall(PetscViewerPopFormat(viewer)); 1690 PetscCall(PetscOptionsRestoreViewer(&viewer)); 1691 PetscCall(PetscPrintf(comm, " Max displacement: %g %g\n", (double)PetscRealPart(p->P_0 * (ctx->xmax[1] - ctx->xmin[1]) * (1. - 2. * nu_u) / (2. * p->mu * (1. - nu_u))), (double)PetscRealPart(p->P_0 * (ctx->xmax[1] - ctx->xmin[1]) * (1. - 2. * nu) / (2. * p->mu * (1. - nu))))); 1692 PetscCall(PetscPrintf(comm, " Relaxation time: %g\n", (double)ctx->t_r)); 1693 } 1694 } 1695 PetscFunctionReturn(PETSC_SUCCESS); 1696 } 1697 1698 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 1699 { 1700 PetscFunctionBeginUser; 1701 PetscCall(DMCreate(comm, dm)); 1702 PetscCall(DMSetType(*dm, DMPLEX)); 1703 PetscCall(DMSetFromOptions(*dm)); 1704 PetscCall(DMSetApplicationContext(*dm, user)); 1705 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 1706 PetscCall(DMGetBoundingBox(*dm, user->xmin, user->xmax)); 1707 PetscFunctionReturn(PETSC_SUCCESS); 1708 } 1709 1710 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 1711 { 1712 PetscErrorCode (*exact[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 1713 PetscErrorCode (*exact_t[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 1714 PetscDS ds; 1715 DMLabel label; 1716 PetscWeakForm wf; 1717 Parameter *param; 1718 PetscInt id_mandel[2]; 1719 PetscInt comp[1]; 1720 PetscInt comp_mandel[2]; 1721 PetscInt dim, id, bd, f; 1722 1723 PetscFunctionBeginUser; 1724 PetscCall(DMGetLabel(dm, "marker", &label)); 1725 PetscCall(DMGetDS(dm, &ds)); 1726 PetscCall(PetscDSGetSpatialDimension(ds, &dim)); 1727 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 1728 exact_t[0] = exact_t[1] = exact_t[2] = zero; 1729 1730 /* Setup Problem Formulation and Boundary Conditions */ 1731 switch (user->solType) { 1732 case SOL_QUADRATIC_LINEAR: 1733 PetscCall(PetscDSSetResidual(ds, 0, f0_quadratic_linear_u, f1_u)); 1734 PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL)); 1735 PetscCall(PetscDSSetResidual(ds, 2, f0_quadratic_linear_p, f1_p)); 1736 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu)); 1737 PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL)); 1738 PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL)); 1739 PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL)); 1740 PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL)); 1741 PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL)); 1742 PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp)); 1743 exact[0] = quadratic_u; 1744 exact[1] = linear_eps; 1745 exact[2] = linear_linear_p; 1746 exact_t[2] = linear_linear_p_t; 1747 1748 id = 1; 1749 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", label, 1, &id, 0, 0, NULL, (void (*)(void))exact[0], NULL, user, NULL)); 1750 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", label, 1, &id, 2, 0, NULL, (void (*)(void))exact[2], (void (*)(void))exact_t[2], user, NULL)); 1751 break; 1752 case SOL_TRIG_LINEAR: 1753 PetscCall(PetscDSSetResidual(ds, 0, f0_trig_linear_u, f1_u)); 1754 PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL)); 1755 PetscCall(PetscDSSetResidual(ds, 2, f0_trig_linear_p, f1_p)); 1756 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu)); 1757 PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL)); 1758 PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL)); 1759 PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL)); 1760 PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL)); 1761 PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL)); 1762 PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp)); 1763 exact[0] = trig_u; 1764 exact[1] = trig_eps; 1765 exact[2] = trig_linear_p; 1766 exact_t[2] = trig_linear_p_t; 1767 1768 id = 1; 1769 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", label, 1, &id, 0, 0, NULL, (void (*)(void))exact[0], NULL, user, NULL)); 1770 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", label, 1, &id, 2, 0, NULL, (void (*)(void))exact[2], (void (*)(void))exact_t[2], user, NULL)); 1771 break; 1772 case SOL_QUADRATIC_TRIG: 1773 PetscCall(PetscDSSetResidual(ds, 0, f0_quadratic_trig_u, f1_u)); 1774 PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL)); 1775 PetscCall(PetscDSSetResidual(ds, 2, f0_quadratic_trig_p, f1_p)); 1776 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu)); 1777 PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL)); 1778 PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL)); 1779 PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL)); 1780 PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL)); 1781 PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL)); 1782 PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp)); 1783 exact[0] = quadratic_u; 1784 exact[1] = linear_eps; 1785 exact[2] = linear_trig_p; 1786 exact_t[2] = linear_trig_p_t; 1787 1788 id = 1; 1789 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", label, 1, &id, 0, 0, NULL, (void (*)(void))exact[0], NULL, user, NULL)); 1790 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", label, 1, &id, 2, 0, NULL, (void (*)(void))exact[2], (void (*)(void))exact_t[2], user, NULL)); 1791 break; 1792 case SOL_TERZAGHI: 1793 PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u)); 1794 PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL)); 1795 PetscCall(PetscDSSetResidual(ds, 2, f0_p, f1_p)); 1796 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu)); 1797 PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL)); 1798 PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL)); 1799 PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL)); 1800 PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL)); 1801 PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL)); 1802 PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp)); 1803 1804 exact[0] = terzaghi_2d_u; 1805 exact[1] = terzaghi_2d_eps; 1806 exact[2] = terzaghi_2d_p; 1807 exact_t[0] = terzaghi_2d_u_t; 1808 exact_t[1] = terzaghi_2d_eps_t; 1809 exact_t[2] = terzaghi_2d_p_t; 1810 1811 id = 1; 1812 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "vertical stress", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 1813 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 1814 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_terzaghi_bd_u, 0, NULL)); 1815 1816 id = 3; 1817 comp[0] = 1; 1818 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed base", label, 1, &id, 0, 1, comp, (void (*)(void))zero, NULL, user, NULL)); 1819 id = 2; 1820 comp[0] = 0; 1821 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed side", label, 1, &id, 0, 1, comp, (void (*)(void))zero, NULL, user, NULL)); 1822 id = 4; 1823 comp[0] = 0; 1824 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed side", label, 1, &id, 0, 1, comp, (void (*)(void))zero, NULL, user, NULL)); 1825 id = 1; 1826 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", label, 1, &id, 2, 0, NULL, (void (*)(void))terzaghi_drainage_pressure, NULL, user, NULL)); 1827 break; 1828 case SOL_MANDEL: 1829 PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u)); 1830 PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL)); 1831 PetscCall(PetscDSSetResidual(ds, 2, f0_p, f1_p)); 1832 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu)); 1833 PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL)); 1834 PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL)); 1835 PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL)); 1836 PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL)); 1837 PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL)); 1838 PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp)); 1839 1840 PetscCall(mandelZeros(PETSC_COMM_WORLD, user, param)); 1841 1842 exact[0] = mandel_2d_u; 1843 exact[1] = mandel_2d_eps; 1844 exact[2] = mandel_2d_p; 1845 exact_t[0] = mandel_2d_u_t; 1846 exact_t[1] = mandel_2d_eps_t; 1847 exact_t[2] = mandel_2d_p_t; 1848 1849 id_mandel[0] = 3; 1850 id_mandel[1] = 1; 1851 //comp[0] = 1; 1852 comp_mandel[0] = 0; 1853 comp_mandel[1] = 1; 1854 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "vertical stress", label, 2, id_mandel, 0, 2, comp_mandel, (void (*)(void))mandel_2d_u, NULL, user, NULL)); 1855 //PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "vertical stress", "marker", 0, 1, comp, NULL, 2, id_mandel, user)); 1856 //PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed base", "marker", 0, 1, comp, (void (*)(void)) zero, 2, id_mandel, user)); 1857 //PetscCall(PetscDSSetBdResidual(ds, 0, f0_mandel_bd_u, NULL)); 1858 1859 id_mandel[0] = 2; 1860 id_mandel[1] = 4; 1861 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", label, 2, id_mandel, 2, 0, NULL, (void (*)(void))zero, NULL, user, NULL)); 1862 break; 1863 case SOL_CRYER: 1864 PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u)); 1865 PetscCall(PetscDSSetResidual(ds, 1, f0_epsilon, NULL)); 1866 PetscCall(PetscDSSetResidual(ds, 2, f0_p, f1_p)); 1867 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu)); 1868 PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL)); 1869 PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL)); 1870 PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL)); 1871 PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL)); 1872 PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL)); 1873 PetscCall(PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp)); 1874 1875 PetscCall(cryerZeros(PETSC_COMM_WORLD, user, param)); 1876 1877 exact[0] = cryer_3d_u; 1878 exact[1] = cryer_3d_eps; 1879 exact[2] = cryer_3d_p; 1880 1881 id = 1; 1882 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "normal stress", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 1883 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 1884 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_cryer_bd_u, 0, NULL)); 1885 1886 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", label, 1, &id, 2, 0, NULL, (void (*)(void))cryer_drainage_pressure, NULL, user, NULL)); 1887 break; 1888 default: 1889 SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%d)", solutionTypes[PetscMin(user->solType, NUM_SOLUTION_TYPES)], user->solType); 1890 } 1891 for (f = 0; f < 3; ++f) { 1892 PetscCall(PetscDSSetExactSolution(ds, f, exact[f], user)); 1893 PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, f, exact_t[f], user)); 1894 } 1895 1896 /* Setup constants */ 1897 { 1898 PetscScalar constants[6]; 1899 constants[0] = param->mu; /* shear modulus, Pa */ 1900 constants[1] = param->K_u; /* undrained bulk modulus, Pa */ 1901 constants[2] = param->alpha; /* Biot effective stress coefficient, - */ 1902 constants[3] = param->M; /* Biot modulus, Pa */ 1903 constants[4] = param->k / param->mu_f; /* Darcy coefficient, m**2 / Pa*s */ 1904 constants[5] = param->P_0; /* Magnitude of Vertical Stress, Pa */ 1905 PetscCall(PetscDSSetConstants(ds, 6, constants)); 1906 } 1907 PetscFunctionReturn(PETSC_SUCCESS); 1908 } 1909 1910 static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace) 1911 { 1912 PetscFunctionBeginUser; 1913 PetscCall(DMPlexCreateRigidBody(dm, origField, nullspace)); 1914 PetscFunctionReturn(PETSC_SUCCESS); 1915 } 1916 1917 static PetscErrorCode SetupFE(DM dm, PetscInt Nf, PetscInt Nc[], const char *name[], PetscErrorCode (*setup)(DM, AppCtx *), void *ctx) 1918 { 1919 AppCtx *user = (AppCtx *)ctx; 1920 DM cdm = dm; 1921 PetscFE fe; 1922 PetscQuadrature q = NULL, fq = NULL; 1923 char prefix[PETSC_MAX_PATH_LEN]; 1924 PetscInt dim, f; 1925 PetscBool simplex; 1926 1927 PetscFunctionBeginUser; 1928 /* Create finite element */ 1929 PetscCall(DMGetDimension(dm, &dim)); 1930 PetscCall(DMPlexIsSimplex(dm, &simplex)); 1931 for (f = 0; f < Nf; ++f) { 1932 PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name[f])); 1933 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc[f], simplex, name[f] ? prefix : NULL, -1, &fe)); 1934 PetscCall(PetscObjectSetName((PetscObject)fe, name[f])); 1935 if (!q) PetscCall(PetscFEGetQuadrature(fe, &q)); 1936 if (!fq) PetscCall(PetscFEGetFaceQuadrature(fe, &fq)); 1937 PetscCall(PetscFESetQuadrature(fe, q)); 1938 PetscCall(PetscFESetFaceQuadrature(fe, fq)); 1939 PetscCall(DMSetField(dm, f, NULL, (PetscObject)fe)); 1940 PetscCall(PetscFEDestroy(&fe)); 1941 } 1942 PetscCall(DMCreateDS(dm)); 1943 PetscCall((*setup)(dm, user)); 1944 while (cdm) { 1945 PetscCall(DMCopyDisc(dm, cdm)); 1946 if (0) PetscCall(DMSetNearNullSpaceConstructor(cdm, 0, CreateElasticityNullSpace)); 1947 /* TODO: Check whether the boundary of coarse meshes is marked */ 1948 PetscCall(DMGetCoarseDM(cdm, &cdm)); 1949 } 1950 PetscCall(PetscFEDestroy(&fe)); 1951 PetscFunctionReturn(PETSC_SUCCESS); 1952 } 1953 1954 static PetscErrorCode SetInitialConditions(TS ts, Vec u) 1955 { 1956 DM dm; 1957 PetscReal t; 1958 1959 PetscFunctionBeginUser; 1960 PetscCall(TSGetDM(ts, &dm)); 1961 PetscCall(TSGetTime(ts, &t)); 1962 if (t <= 0.0) { 1963 PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 1964 void *ctxs[3]; 1965 AppCtx *ctx; 1966 1967 PetscCall(DMGetApplicationContext(dm, &ctx)); 1968 switch (ctx->solType) { 1969 case SOL_TERZAGHI: 1970 funcs[0] = terzaghi_initial_u; 1971 ctxs[0] = ctx; 1972 funcs[1] = terzaghi_initial_eps; 1973 ctxs[1] = ctx; 1974 funcs[2] = terzaghi_drainage_pressure; 1975 ctxs[2] = ctx; 1976 PetscCall(DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u)); 1977 break; 1978 case SOL_MANDEL: 1979 funcs[0] = mandel_initial_u; 1980 ctxs[0] = ctx; 1981 funcs[1] = mandel_initial_eps; 1982 ctxs[1] = ctx; 1983 funcs[2] = mandel_drainage_pressure; 1984 ctxs[2] = ctx; 1985 PetscCall(DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u)); 1986 break; 1987 case SOL_CRYER: 1988 funcs[0] = cryer_initial_u; 1989 ctxs[0] = ctx; 1990 funcs[1] = cryer_initial_eps; 1991 ctxs[1] = ctx; 1992 funcs[2] = cryer_drainage_pressure; 1993 ctxs[2] = ctx; 1994 PetscCall(DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u)); 1995 break; 1996 default: 1997 PetscCall(DMComputeExactSolution(dm, t, u, NULL)); 1998 } 1999 } else { 2000 PetscCall(DMComputeExactSolution(dm, t, u, NULL)); 2001 } 2002 PetscFunctionReturn(PETSC_SUCCESS); 2003 } 2004 2005 /* Need to create Viewer each time because HDF5 can get corrupted */ 2006 static PetscErrorCode SolutionMonitor(TS ts, PetscInt steps, PetscReal time, Vec u, void *mctx) 2007 { 2008 DM dm; 2009 Vec exact; 2010 PetscViewer viewer; 2011 PetscViewerFormat format; 2012 PetscOptions options; 2013 const char *prefix; 2014 2015 PetscFunctionBeginUser; 2016 PetscCall(TSGetDM(ts, &dm)); 2017 PetscCall(PetscObjectGetOptions((PetscObject)ts, &options)); 2018 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, &prefix)); 2019 PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts), options, prefix, "-monitor_solution", &viewer, &format, NULL)); 2020 PetscCall(DMGetGlobalVector(dm, &exact)); 2021 PetscCall(DMComputeExactSolution(dm, time, exact, NULL)); 2022 PetscCall(DMSetOutputSequenceNumber(dm, steps, time)); 2023 PetscCall(VecView(exact, viewer)); 2024 PetscCall(VecView(u, viewer)); 2025 PetscCall(DMRestoreGlobalVector(dm, &exact)); 2026 { 2027 PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 2028 void **ectxs; 2029 PetscReal *err; 2030 PetscInt Nf, f; 2031 2032 PetscCall(DMGetNumFields(dm, &Nf)); 2033 PetscCall(PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err)); 2034 { 2035 PetscInt Nds, s; 2036 2037 PetscCall(DMGetNumDS(dm, &Nds)); 2038 for (s = 0; s < Nds; ++s) { 2039 PetscDS ds; 2040 DMLabel label; 2041 IS fieldIS; 2042 const PetscInt *fields; 2043 PetscInt dsNf, f; 2044 2045 PetscCall(DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds, NULL)); 2046 PetscCall(PetscDSGetNumFields(ds, &dsNf)); 2047 PetscCall(ISGetIndices(fieldIS, &fields)); 2048 for (f = 0; f < dsNf; ++f) { 2049 const PetscInt field = fields[f]; 2050 PetscCall(PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field])); 2051 } 2052 PetscCall(ISRestoreIndices(fieldIS, &fields)); 2053 } 2054 } 2055 PetscCall(DMComputeL2FieldDiff(dm, time, exacts, ectxs, u, err)); 2056 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "Time: %g L_2 Error: [", (double)time)); 2057 for (f = 0; f < Nf; ++f) { 2058 if (f) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), ", ")); 2059 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%g", (double)err[f])); 2060 } 2061 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "]\n")); 2062 PetscCall(PetscFree3(exacts, ectxs, err)); 2063 } 2064 PetscCall(PetscOptionsRestoreViewer(&viewer)); 2065 PetscFunctionReturn(PETSC_SUCCESS); 2066 } 2067 2068 static PetscErrorCode SetupMonitor(TS ts, AppCtx *ctx) 2069 { 2070 PetscViewer viewer; 2071 PetscViewerFormat format; 2072 PetscOptions options; 2073 const char *prefix; 2074 PetscBool flg; 2075 2076 PetscFunctionBeginUser; 2077 PetscCall(PetscObjectGetOptions((PetscObject)ts, &options)); 2078 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, &prefix)); 2079 PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts), options, prefix, "-monitor_solution", &viewer, &format, &flg)); 2080 if (flg) PetscCall(TSMonitorSet(ts, SolutionMonitor, ctx, NULL)); 2081 PetscCall(PetscOptionsRestoreViewer(&viewer)); 2082 PetscFunctionReturn(PETSC_SUCCESS); 2083 } 2084 2085 static PetscErrorCode TSAdaptChoose_Terzaghi(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept, PetscReal *wlte, PetscReal *wltea, PetscReal *wlter) 2086 { 2087 static PetscReal dtTarget = -1.0; 2088 PetscReal dtInitial; 2089 DM dm; 2090 AppCtx *ctx; 2091 PetscInt step; 2092 2093 PetscFunctionBeginUser; 2094 PetscCall(TSGetDM(ts, &dm)); 2095 PetscCall(DMGetApplicationContext(dm, &ctx)); 2096 PetscCall(TSGetStepNumber(ts, &step)); 2097 dtInitial = ctx->dtInitial < 0.0 ? 1.0e-4 * ctx->t_r : ctx->dtInitial; 2098 if (!step) { 2099 if (PetscAbsReal(dtInitial - h) > PETSC_SMALL) { 2100 *accept = PETSC_FALSE; 2101 *next_h = dtInitial; 2102 dtTarget = h; 2103 } else { 2104 *accept = PETSC_TRUE; 2105 *next_h = dtTarget < 0.0 ? dtInitial : dtTarget; 2106 dtTarget = -1.0; 2107 } 2108 } else { 2109 *accept = PETSC_TRUE; 2110 *next_h = h; 2111 } 2112 *next_sc = 0; /* Reuse the same order scheme */ 2113 *wlte = -1; /* Weighted local truncation error was not evaluated */ 2114 *wltea = -1; /* Weighted absolute local truncation error was not evaluated */ 2115 *wlter = -1; /* Weighted relative local truncation error was not evaluated */ 2116 PetscFunctionReturn(PETSC_SUCCESS); 2117 } 2118 2119 int main(int argc, char **argv) 2120 { 2121 AppCtx ctx; /* User-defined work context */ 2122 DM dm; /* Problem specification */ 2123 TS ts; /* Time Series / Nonlinear solver */ 2124 Vec u; /* Solutions */ 2125 const char *name[3] = {"displacement", "tracestrain", "pressure"}; 2126 PetscReal t; 2127 PetscInt dim, Nc[3]; 2128 2129 PetscFunctionBeginUser; 2130 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 2131 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx)); 2132 PetscCall(PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &ctx.bag)); 2133 PetscCall(PetscMalloc1(ctx.niter, &ctx.zeroArray)); 2134 PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm)); 2135 PetscCall(SetupParameters(PETSC_COMM_WORLD, &ctx)); 2136 /* Primal System */ 2137 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 2138 PetscCall(DMSetApplicationContext(dm, &ctx)); 2139 PetscCall(TSSetDM(ts, dm)); 2140 2141 PetscCall(DMGetDimension(dm, &dim)); 2142 Nc[0] = dim; 2143 Nc[1] = 1; 2144 Nc[2] = 1; 2145 2146 PetscCall(SetupFE(dm, 3, Nc, name, SetupPrimalProblem, &ctx)); 2147 PetscCall(DMCreateGlobalVector(dm, &u)); 2148 PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx)); 2149 PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx)); 2150 PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx)); 2151 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 2152 PetscCall(TSSetFromOptions(ts)); 2153 PetscCall(TSSetComputeInitialCondition(ts, SetInitialConditions)); 2154 PetscCall(SetupMonitor(ts, &ctx)); 2155 2156 if (ctx.solType != SOL_QUADRATIC_TRIG) { 2157 TSAdapt adapt; 2158 2159 PetscCall(TSGetAdapt(ts, &adapt)); 2160 adapt->ops->choose = TSAdaptChoose_Terzaghi; 2161 } 2162 if (ctx.solType == SOL_CRYER) { 2163 Mat J; 2164 MatNullSpace sp; 2165 2166 PetscCall(TSSetUp(ts)); 2167 PetscCall(TSGetIJacobian(ts, &J, NULL, NULL, NULL)); 2168 PetscCall(DMPlexCreateRigidBody(dm, 0, &sp)); 2169 PetscCall(MatSetNullSpace(J, sp)); 2170 PetscCall(MatNullSpaceDestroy(&sp)); 2171 } 2172 PetscCall(TSGetTime(ts, &t)); 2173 PetscCall(DMSetOutputSequenceNumber(dm, 0, t)); 2174 PetscCall(DMTSCheckFromOptions(ts, u)); 2175 PetscCall(SetInitialConditions(ts, u)); 2176 PetscCall(PetscObjectSetName((PetscObject)u, "solution")); 2177 PetscCall(TSSolve(ts, u)); 2178 PetscCall(DMTSCheckFromOptions(ts, u)); 2179 PetscCall(TSGetSolution(ts, &u)); 2180 PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view")); 2181 2182 /* Cleanup */ 2183 PetscCall(VecDestroy(&u)); 2184 PetscCall(TSDestroy(&ts)); 2185 PetscCall(DMDestroy(&dm)); 2186 PetscCall(PetscBagDestroy(&ctx.bag)); 2187 PetscCall(PetscFree(ctx.zeroArray)); 2188 PetscCall(PetscFinalize()); 2189 return 0; 2190 } 2191 2192 /*TEST 2193 2194 test: 2195 suffix: 2d_quad_linear 2196 requires: triangle 2197 args: -sol_type quadratic_linear -dm_refine 2 \ 2198 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2199 -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme 2200 2201 test: 2202 suffix: 3d_quad_linear 2203 requires: ctetgen 2204 args: -dm_plex_dim 3 -sol_type quadratic_linear -dm_refine 1 \ 2205 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2206 -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme 2207 2208 test: 2209 suffix: 2d_trig_linear 2210 requires: triangle 2211 args: -sol_type trig_linear -dm_refine 1 \ 2212 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2213 -dmts_check .0001 -ts_max_steps 5 -ts_dt 0.00001 -ts_monitor_extreme 2214 2215 test: 2216 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9, 2.1, 1.8] 2217 suffix: 2d_trig_linear_sconv 2218 requires: triangle 2219 args: -sol_type trig_linear -dm_refine 1 \ 2220 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2221 -convest_num_refine 1 -ts_convergence_estimate -ts_convergence_temporal 0 -ts_max_steps 1 -ts_dt 0.00001 -pc_type lu 2222 2223 test: 2224 suffix: 3d_trig_linear 2225 requires: ctetgen 2226 args: -dm_plex_dim 3 -sol_type trig_linear -dm_refine 1 \ 2227 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2228 -dmts_check .0001 -ts_max_steps 2 -ts_monitor_extreme 2229 2230 test: 2231 # -dm_refine 1 -convest_num_refine 2 gets L_2 convergence rate: [2.0, 2.1, 1.9] 2232 suffix: 3d_trig_linear_sconv 2233 requires: ctetgen 2234 args: -dm_plex_dim 3 -sol_type trig_linear -dm_refine 1 \ 2235 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2236 -convest_num_refine 1 -ts_convergence_estimate -ts_convergence_temporal 0 -ts_max_steps 1 -pc_type lu 2237 2238 test: 2239 suffix: 2d_quad_trig 2240 requires: triangle 2241 args: -sol_type quadratic_trig -dm_refine 2 \ 2242 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2243 -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme 2244 2245 test: 2246 # Using -dm_refine 4 gets the convergence rates to [0.95, 0.97, 0.90] 2247 suffix: 2d_quad_trig_tconv 2248 requires: triangle 2249 args: -sol_type quadratic_trig -dm_refine 1 \ 2250 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2251 -convest_num_refine 3 -ts_convergence_estimate -ts_max_steps 5 -pc_type lu 2252 2253 test: 2254 suffix: 3d_quad_trig 2255 requires: ctetgen 2256 args: -dm_plex_dim 3 -sol_type quadratic_trig -dm_refine 1 \ 2257 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2258 -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme 2259 2260 test: 2261 # Using -dm_refine 2 -convest_num_refine 3 gets the convergence rates to [1.0, 1.0, 1.0] 2262 suffix: 3d_quad_trig_tconv 2263 requires: ctetgen 2264 args: -dm_plex_dim 3 -sol_type quadratic_trig -dm_refine 1 \ 2265 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2266 -convest_num_refine 1 -ts_convergence_estimate -ts_max_steps 5 -pc_type lu 2267 2268 testset: 2269 args: -sol_type terzaghi -dm_plex_simplex 0 -dm_plex_box_faces 1,8 -dm_plex_box_lower 0,0 -dm_plex_box_upper 10,10 -dm_plex_separate_marker \ 2270 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 -niter 16000 \ 2271 -pc_type lu 2272 2273 test: 2274 suffix: 2d_terzaghi 2275 requires: double 2276 args: -ts_dt 0.0028666667 -ts_max_steps 2 -ts_monitor -dmts_check .0001 2277 2278 test: 2279 # -dm_plex_box_faces 1,64 -ts_max_steps 4 -convest_num_refine 3 gives L_2 convergence rate: [1.1, 1.1, 1.1] 2280 suffix: 2d_terzaghi_tconv 2281 args: -ts_dt 0.023 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1 2282 2283 test: 2284 # -dm_plex_box_faces 1,16 -convest_num_refine 4 gives L_2 convergence rate: [1.7, 1.2, 1.1] 2285 # if we add -displacement_petscspace_degree 3 -tracestrain_petscspace_degree 2 -pressure_petscspace_degree 2, we get [2.1, 1.6, 1.5] 2286 suffix: 2d_terzaghi_sconv 2287 args: -ts_dt 1e-5 -dt_initial 1e-5 -ts_max_steps 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 2288 2289 testset: 2290 args: -sol_type mandel -dm_plex_simplex 0 -dm_plex_box_lower -0.5,-0.125 -dm_plex_box_upper 0.5,0.125 -dm_plex_separate_marker -dm_refine 1 \ 2291 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2292 -pc_type lu 2293 2294 test: 2295 suffix: 2d_mandel 2296 requires: double 2297 args: -ts_dt 0.0028666667 -ts_max_steps 2 -ts_monitor -dmts_check .0001 2298 2299 test: 2300 # -dm_refine 3 -ts_max_steps 4 -convest_num_refine 3 gives L_2 convergence rate: [1.6, 0.93, 1.2] 2301 suffix: 2d_mandel_sconv 2302 args: -ts_dt 1e-5 -dt_initial 1e-5 -ts_max_steps 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 2303 2304 test: 2305 # -dm_refine 5 -ts_max_steps 4 -convest_num_refine 3 gives L_2 convergence rate: [0.26, -0.0058, 0.26] 2306 suffix: 2d_mandel_tconv 2307 args: -ts_dt 0.023 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1 2308 2309 testset: 2310 requires: ctetgen !complex 2311 args: -sol_type cryer -dm_plex_dim 3 -dm_plex_shape ball \ 2312 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 2313 2314 test: 2315 suffix: 3d_cryer 2316 args: -ts_dt 0.0028666667 -ts_max_time 0.014333 -ts_max_steps 2 -dmts_check .0001 \ 2317 -pc_type svd 2318 2319 test: 2320 # -bd_dm_refine 3 -dm_refine_volume_limit_pre 0.004 -convest_num_refine 2 gives L_2 convergence rate: [] 2321 suffix: 3d_cryer_sconv 2322 args: -bd_dm_refine 1 -dm_refine_volume_limit_pre 0.00666667 \ 2323 -ts_dt 1e-5 -dt_initial 1e-5 -ts_max_steps 2 \ 2324 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 2325 -pc_type lu -pc_factor_shift_type nonzero 2326 2327 test: 2328 # Displacement and Pressure converge. The analytic expression for trace strain is inaccurate at the origin 2329 # -bd_dm_refine 3 -ref_limit 0.00666667 -ts_max_steps 5 -convest_num_refine 2 gives L_2 convergence rate: [0.47, -0.43, 1.5] 2330 suffix: 3d_cryer_tconv 2331 args: -bd_dm_refine 1 -dm_refine_volume_limit_pre 0.00666667 \ 2332 -ts_dt 0.023 -ts_max_time 0.092 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1 \ 2333 -pc_type lu -pc_factor_shift_type nonzero 2334 2335 TEST*/ 2336