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