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