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