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