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