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 #if 0 1359 static void f0_mandel_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1360 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1361 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1362 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1363 { 1364 // Uniform stress distribution 1365 /* PetscScalar xmax = 0.5; 1366 PetscScalar xmin = -0.5; 1367 PetscScalar ymax = 0.5; 1368 PetscScalar ymin = -0.5; 1369 PetscScalar P = constants[5]; 1370 PetscScalar aL = (xmax - xmin) / 2.0; 1371 PetscScalar sigma_zz = -1.0*P / aL; */ 1372 1373 // Analytical (parabolic) stress distribution 1374 PetscReal a1, a2, am; 1375 PetscReal y1, y2, ym; 1376 1377 PetscInt NITER = 500; 1378 PetscReal EPS = 0.000001; 1379 PetscReal zeroArray[500]; /* NITER */ 1380 PetscReal xmax = 1.0; 1381 PetscReal xmin = 0.0; 1382 PetscReal ymax = 0.1; 1383 PetscReal ymin = 0.0; 1384 PetscReal lower[2], upper[2]; 1385 1386 lower[0] = xmin - (xmax - xmin) / 2.0; 1387 lower[1] = ymin - (ymax - ymin) / 2.0; 1388 upper[0] = xmax - (xmax - xmin) / 2.0; 1389 upper[1] = ymax - (ymax - ymin) / 2.0; 1390 1391 xmin = lower[0]; 1392 ymin = lower[1]; 1393 xmax = upper[0]; 1394 ymax = upper[1]; 1395 1396 PetscScalar G = constants[0]; 1397 PetscScalar K_u = constants[1]; 1398 PetscScalar alpha = constants[2]; 1399 PetscScalar M = constants[3]; 1400 PetscScalar kappa = constants[4]; 1401 PetscScalar F = constants[5]; 1402 1403 PetscScalar K_d = K_u - alpha*alpha*M; 1404 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G)); 1405 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); 1406 PetscReal aL = (xmax - xmin) / 2.0; 1407 PetscReal c = PetscRealPart(((2.0*kappa*G) * (1.0 - nu) * (nu_u - nu)) / (alpha*alpha * (1.0 - 2.0*nu) * (1.0 - nu_u))); 1408 PetscScalar B = (3.0 * (nu_u - nu)) / ( alpha * (1.0 - 2.0*nu) * (1.0 + nu_u)); 1409 PetscScalar A1 = 3.0 / (B * (1.0 + nu_u)); 1410 PetscScalar A2 = (alpha * (1.0 - 2.0*nu)) / (1.0 - nu); 1411 1412 // Generate zero values 1413 for (PetscInt i=1; i < NITER+1; i++) 1414 { 1415 a1 = ((PetscReal) i - 1.0) * PETSC_PI * PETSC_PI / 4.0 + EPS; 1416 a2 = a1 + PETSC_PI/2; 1417 for (PetscInt j=0; j<NITER; j++) 1418 { 1419 y1 = PetscTanReal(a1) - PetscRealPart(A1/A2)*a1; 1420 y2 = PetscTanReal(a2) - PetscRealPart(A1/A2)*a2; 1421 am = (a1 + a2)/2.0; 1422 ym = PetscTanReal(am) - PetscRealPart(A1/A2)*am; 1423 if ((ym*y1) > 0) 1424 { 1425 a1 = am; 1426 } else { 1427 a2 = am; 1428 } 1429 if (PetscAbsReal(y2) < EPS) 1430 { 1431 am = a2; 1432 } 1433 } 1434 zeroArray[i-1] = am; 1435 } 1436 1437 // Solution for sigma_zz 1438 PetscScalar A_x = 0.0; 1439 PetscScalar B_x = 0.0; 1440 1441 for (PetscInt n=1; n < NITER+1; n++) 1442 { 1443 PetscReal alpha_n = zeroArray[n-1]; 1444 1445 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))); 1446 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))); 1447 } 1448 1449 PetscScalar sigma_zz = -1.0*(F/aL) - ((2.0*F)/aL) * (A2/A1) * A_x + ((2.0*F)/aL) * B_x; 1450 1451 1452 if (x[1] == ymax) { 1453 f0[1] += sigma_zz; 1454 } else if (x[1] == ymin) { 1455 f0[1] -= sigma_zz; 1456 } 1457 } 1458 #endif 1459 1460 static void f0_cryer_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1461 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1462 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1463 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1464 { 1465 const PetscReal P_0 = PetscRealPart(constants[5]); 1466 //const PetscReal R = PetscSqrtReal(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]); 1467 PetscInt d; 1468 1469 for (d = 0; d < dim; ++d) f0[d] = -P_0*n[d]; 1470 //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); 1471 //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); 1472 //for (d = 0; d < dim; ++d) f0[d] = -P_0*x[d]/R; 1473 } 1474 1475 /* Standard Kernels - Residual */ 1476 /* f0_e */ 1477 static void f0_epsilon(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1478 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1479 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1480 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1481 { 1482 PetscInt d; 1483 1484 for (d = 0; d < dim; ++d) { 1485 f0[0] += u_x[d*dim+d]; 1486 } 1487 f0[0] -= u[uOff[1]]; 1488 } 1489 1490 /* f0_p */ 1491 static void f0_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1492 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1493 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1494 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1495 { 1496 const PetscReal alpha = PetscRealPart(constants[2]); 1497 const PetscReal M = PetscRealPart(constants[3]); 1498 1499 f0[0] += alpha*u_t[uOff[1]]; 1500 f0[0] += u_t[uOff[2]]/M; 1501 } 1502 1503 /* f1_u */ 1504 static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1505 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1506 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1507 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 1508 { 1509 const PetscInt Nc = dim; 1510 const PetscReal G = PetscRealPart(constants[0]); 1511 const PetscReal K_u = PetscRealPart(constants[1]); 1512 const PetscReal alpha = PetscRealPart(constants[2]); 1513 const PetscReal M = PetscRealPart(constants[3]); 1514 const PetscReal K_d = K_u - alpha*alpha*M; 1515 const PetscReal lambda = K_d - (2.0 * G) / 3.0; 1516 PetscInt c, d; 1517 1518 for (c = 0; c < Nc; ++c) 1519 { 1520 for (d = 0; d < dim; ++d) 1521 { 1522 f1[c*dim+d] -= G*(u_x[c*dim+d] + u_x[d*dim+c]); 1523 } 1524 f1[c*dim+c] -= lambda*u[uOff[1]]; 1525 f1[c*dim+c] += alpha*u[uOff[2]]; 1526 } 1527 } 1528 1529 /* f1_p */ 1530 static void f1_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1531 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1532 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1533 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 1534 { 1535 const PetscReal kappa = PetscRealPart(constants[4]); 1536 PetscInt d; 1537 1538 for (d = 0; d < dim; ++d) { 1539 f1[d] += kappa*u_x[uOff_x[2]+d]; 1540 } 1541 } 1542 1543 /* 1544 \partial_df \phi_fc g_{fc,gc,df,dg} \partial_dg \phi_gc 1545 1546 \partial_df \phi_fc \lambda \delta_{fc,df} \sum_gc \partial_dg \phi_gc \delta_{gc,dg} 1547 = \partial_fc \phi_fc \sum_gc \partial_gc \phi_gc 1548 */ 1549 1550 1551 /* Standard Kernels - Jacobian */ 1552 /* g0_ee */ 1553 static void g0_ee(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1554 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1555 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1556 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 1557 { 1558 g0[0] = -1.0; 1559 } 1560 1561 /* g0_pe */ 1562 static void g0_pe(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1563 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1564 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1565 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 1566 { 1567 const PetscReal alpha = PetscRealPart(constants[2]); 1568 1569 g0[0] = u_tShift*alpha; 1570 } 1571 1572 /* g0_pp */ 1573 static void g0_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1574 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1575 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1576 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 1577 { 1578 const PetscReal M = PetscRealPart(constants[3]); 1579 1580 g0[0] = u_tShift/M; 1581 } 1582 1583 /* g1_eu */ 1584 static void g1_eu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1585 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1586 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1587 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 1588 { 1589 PetscInt d; 1590 for (d = 0; d < dim; ++d) g1[d*dim+d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */ 1591 } 1592 1593 /* g2_ue */ 1594 static void g2_ue(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1595 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1596 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1597 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]) 1598 { 1599 const PetscReal G = PetscRealPart(constants[0]); 1600 const PetscReal K_u = PetscRealPart(constants[1]); 1601 const PetscReal alpha = PetscRealPart(constants[2]); 1602 const PetscReal M = PetscRealPart(constants[3]); 1603 const PetscReal K_d = K_u - alpha*alpha*M; 1604 const PetscReal lambda = K_d - (2.0 * G) / 3.0; 1605 PetscInt d; 1606 1607 for (d = 0; d < dim; ++d) { 1608 g2[d*dim + d] -= lambda; 1609 } 1610 } 1611 /* g2_up */ 1612 static void g2_up(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1613 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1614 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1615 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]) 1616 { 1617 const PetscReal alpha = PetscRealPart(constants[2]); 1618 PetscInt d; 1619 1620 for (d = 0; d < dim; ++d) { 1621 g2[d*dim + d] += alpha; 1622 } 1623 } 1624 1625 /* g3_uu */ 1626 static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1627 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1628 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1629 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 1630 { 1631 const PetscInt Nc = dim; 1632 const PetscReal G = PetscRealPart(constants[0]); 1633 PetscInt c, d; 1634 1635 for (c = 0; c < Nc; ++c) { 1636 for (d = 0; d < dim; ++d) { 1637 g3[((c*Nc + c)*dim + d)*dim + d] -= G; 1638 g3[((c*Nc + d)*dim + d)*dim + c] -= G; 1639 } 1640 } 1641 } 1642 1643 /* g3_pp */ 1644 static void g3_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1645 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1646 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1647 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 1648 { 1649 const PetscReal kappa = PetscRealPart(constants[4]); 1650 PetscInt d; 1651 1652 for (d = 0; d < dim; ++d) g3[d*dim+d] += kappa; 1653 } 1654 1655 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 1656 { 1657 PetscInt sol; 1658 PetscErrorCode ierr; 1659 1660 PetscFunctionBeginUser; 1661 options->dim = 2; 1662 options->simplex = PETSC_TRUE; 1663 options->refLimit = -1.0; 1664 options->solType = SOL_QUADRATIC_TRIG; 1665 options->niter = 500; 1666 options->eps = PETSC_SMALL; 1667 options->dtInitial = -1.0; 1668 ierr = PetscStrncpy(options->dmType, DMPLEX, 256);CHKERRQ(ierr); 1669 1670 ierr = PetscOptionsBegin(comm, "", "Biot Poroelasticity Options", "DMPLEX");CHKERRQ(ierr); 1671 ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex53.c", options->dim, &options->dim, NULL);CHKERRQ(ierr); 1672 ierr = PetscOptionsInt("-niter", "Number of series term iterations in exact solutions", "ex53.c", options->niter, &options->niter, NULL);CHKERRQ(ierr); 1673 ierr = PetscOptionsBool("-simplex", "Simplicial (true) or tensor (false) mesh", "ex53.c", options->simplex, &options->simplex, NULL);CHKERRQ(ierr); 1674 ierr = PetscOptionsReal("-ref_limit", "Maximum cell volume for refined mesh", "ex53.c", options->refLimit, &options->refLimit, NULL);CHKERRQ(ierr); 1675 sol = options->solType; 1676 ierr = PetscOptionsEList("-sol_type", "Type of exact solution", "ex53.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL);CHKERRQ(ierr); 1677 options->solType = (SolutionType) sol; 1678 ierr = PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex53.c", DMList, options->dmType, options->dmType, 256, NULL);CHKERRQ(ierr); 1679 ierr = PetscOptionsReal("-eps", "Precision value for root finding", "ex53.c", options->eps, &options->eps, NULL);CHKERRQ(ierr); 1680 ierr = PetscOptionsReal("-dt_initial", "Override the initial timestep", "ex53.c", options->dtInitial, &options->dtInitial, NULL);CHKERRQ(ierr); 1681 1682 // Wrap up loose ends 1683 if (options->solType == SOL_CRYER) { 1684 options->dim = 3; 1685 } 1686 1687 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1688 PetscFunctionReturn(0); 1689 } 1690 1691 static PetscErrorCode mandelZeros(MPI_Comm comm, AppCtx *ctx, Parameter *param) 1692 { 1693 //PetscBag bag; 1694 PetscReal a1, a2, am; 1695 PetscReal y1, y2, ym; 1696 1697 PetscFunctionBeginUser; 1698 //ierr = PetscBagGetData(ctx->bag, (void **) ¶m);CHKERRQ(ierr); 1699 PetscInt NITER = ctx->niter; 1700 PetscReal EPS = ctx->eps; 1701 //const PetscScalar YMAX = param->ymax; 1702 //const PetscScalar YMIN = param->ymin; 1703 PetscScalar alpha = param->alpha; 1704 PetscScalar K_u = param->K_u; 1705 PetscScalar M = param->M; 1706 PetscScalar G = param->mu; 1707 //const PetscScalar k = param->k; 1708 //const PetscScalar mu_f = param->mu_f; 1709 //const PetscScalar P_0 = param->P_0; 1710 1711 PetscScalar K_d = K_u - alpha*alpha*M; 1712 PetscScalar nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G)); 1713 PetscScalar nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); 1714 //const PetscScalar kappa = k / mu_f; 1715 1716 // Generate zero values 1717 for (PetscInt i=1; i < NITER+1; i++) 1718 { 1719 a1 = ((PetscReal) i - 1.0) * PETSC_PI * PETSC_PI / 4.0 + EPS; 1720 a2 = a1 + PETSC_PI/2; 1721 am = a1; 1722 for (PetscInt j=0; j<NITER; j++) 1723 { 1724 y1 = PetscTanReal(a1) - PetscRealPart((1.0 - nu)/(nu_u - nu))*a1; 1725 y2 = PetscTanReal(a2) - PetscRealPart((1.0 - nu)/(nu_u - nu))*a2; 1726 am = (a1 + a2)/2.0; 1727 ym = PetscTanReal(am) - PetscRealPart((1.0 - nu)/(nu_u - nu))*am; 1728 if ((ym*y1) > 0) 1729 { 1730 a1 = am; 1731 } else { 1732 a2 = am; 1733 } 1734 if (PetscAbsReal(y2) < EPS) 1735 { 1736 am = a2; 1737 } 1738 } 1739 ctx->zeroArray[i-1] = am; 1740 } 1741 PetscFunctionReturn(0); 1742 } 1743 1744 static PetscReal CryerFunction(PetscReal nu_u, PetscReal nu, PetscReal x) 1745 { 1746 return PetscTanReal(PetscSqrtReal(x))*(6.0*(nu_u - nu) - (1.0 - nu)*(1.0 + nu_u)*x) - (6.0*(nu_u - nu)*PetscSqrtReal(x)); 1747 } 1748 1749 static PetscErrorCode cryerZeros(MPI_Comm comm, AppCtx *ctx, Parameter *param) 1750 { 1751 PetscReal alpha = PetscRealPart(param->alpha); /* - */ 1752 PetscReal K_u = PetscRealPart(param->K_u); /* Pa */ 1753 PetscReal M = PetscRealPart(param->M); /* Pa */ 1754 PetscReal G = PetscRealPart(param->mu); /* Pa */ 1755 PetscInt N = ctx->niter, n; 1756 1757 PetscReal K_d = K_u - alpha*alpha*M; /* Pa, Cheng (B.5) */ 1758 PetscReal nu = (3.0*K_d - 2.0*G) / (2.0*(3.0*K_d + G)); /* -, Cheng (B.8) */ 1759 PetscReal nu_u = (3.0*K_u - 2.0*G) / (2.0*(3.0*K_u + G)); /* -, Cheng (B.9) */ 1760 1761 PetscFunctionBeginUser; 1762 for (n = 1; n < N+1; ++n) { 1763 PetscReal tol = PetscPowReal(n, 1.5)*ctx->eps; 1764 PetscReal a1 = 0., a2 = 0., am = 0.; 1765 PetscReal y1, y2, ym; 1766 PetscInt j, k = n-1; 1767 1768 y1 = y2 = 1.; 1769 while (y1*y2 > 0) { 1770 ++k; 1771 a1 = PetscSqr(n*PETSC_PI) - k*PETSC_PI; 1772 a2 = PetscSqr(n*PETSC_PI) + k*PETSC_PI; 1773 y1 = CryerFunction(nu_u, nu, a1); 1774 y2 = CryerFunction(nu_u, nu, a2); 1775 } 1776 for (j = 0; j < 50000; ++j) { 1777 y1 = CryerFunction(nu_u, nu, a1); 1778 y2 = CryerFunction(nu_u, nu, a2); 1779 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); 1780 am = (a1 + a2) / 2.0; 1781 ym = CryerFunction(nu_u, nu, am); 1782 if ((ym * y1) < 0) a2 = am; 1783 else a1 = am; 1784 if (PetscAbsScalar(ym) < tol) break; 1785 } 1786 if (PetscAbsScalar(ym) >= tol) SETERRQ2(comm, PETSC_ERR_PLIB, "Root finding did not converge for root %D (%g)", n, PetscAbsScalar(ym)); 1787 ctx->zeroArray[n-1] = am; 1788 } 1789 PetscFunctionReturn(0); 1790 } 1791 1792 static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx) 1793 { 1794 PetscBag bag; 1795 Parameter *p; 1796 PetscErrorCode ierr; 1797 1798 PetscFunctionBeginUser; 1799 /* setup PETSc parameter bag */ 1800 ierr = PetscBagGetData(ctx->bag,(void**)&p);CHKERRQ(ierr); 1801 ierr = PetscBagSetName(ctx->bag,"par","Poroelastic Parameters");CHKERRQ(ierr); 1802 bag = ctx->bag; 1803 if (ctx->solType == SOL_TERZAGHI) { 1804 // Realistic values - Terzaghi 1805 ierr = PetscBagRegisterScalar(bag, &p->mu, 3.0, "mu", "Shear Modulus, Pa");CHKERRQ(ierr); 1806 ierr = PetscBagRegisterScalar(bag, &p->K_u, 9.76, "K_u", "Undrained Bulk Modulus, Pa");CHKERRQ(ierr); 1807 ierr = PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -");CHKERRQ(ierr); 1808 ierr = PetscBagRegisterScalar(bag, &p->M, 16.0, "M", "Biot Modulus, Pa");CHKERRQ(ierr); 1809 ierr = PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2");CHKERRQ(ierr); 1810 ierr = PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s");CHKERRQ(ierr); 1811 ierr = PetscBagRegisterScalar(bag, &p->zmax, 1.0, "zmax", "Depth Maximum Extent, m");CHKERRQ(ierr); 1812 ierr = PetscBagRegisterScalar(bag, &p->zmin, 0.0, "zmin", "Depth Minimum Extent, m");CHKERRQ(ierr); 1813 ierr = PetscBagRegisterScalar(bag, &p->ymax, 10.0, "ymax", "Vertical Maximum Extent, m");CHKERRQ(ierr); 1814 ierr = PetscBagRegisterScalar(bag, &p->ymin, 0.0, "ymin", "Vertical Minimum Extent, m");CHKERRQ(ierr); 1815 ierr = PetscBagRegisterScalar(bag, &p->xmax, 10.0, "xmax", "Horizontal Maximum Extent, m");CHKERRQ(ierr); 1816 ierr = PetscBagRegisterScalar(bag, &p->xmin, 0.0, "xmin", "Horizontal Minimum Extent, m");CHKERRQ(ierr); 1817 ierr = PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa");CHKERRQ(ierr); 1818 } else if (ctx->solType == SOL_MANDEL) { 1819 // Realistic values - Mandel 1820 ierr = PetscBagRegisterScalar(bag, &p->mu, 0.75, "mu", "Shear Modulus, Pa");CHKERRQ(ierr); 1821 ierr = PetscBagRegisterScalar(bag, &p->K_u, 2.6941176470588233, "K_u", "Undrained Bulk Modulus, Pa");CHKERRQ(ierr); 1822 ierr = PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -");CHKERRQ(ierr); 1823 ierr = PetscBagRegisterScalar(bag, &p->M, 4.705882352941176, "M", "Biot Modulus, Pa");CHKERRQ(ierr); 1824 ierr = PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2");CHKERRQ(ierr); 1825 ierr = PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s");CHKERRQ(ierr); 1826 ierr = PetscBagRegisterScalar(bag, &p->zmax, 1.0, "zmax", "Depth Maximum Extent, m");CHKERRQ(ierr); 1827 ierr = PetscBagRegisterScalar(bag, &p->zmin, 0.0, "zmin", "Depth Minimum Extent, m");CHKERRQ(ierr); 1828 ierr = PetscBagRegisterScalar(bag, &p->ymax, 0.25, "ymax", "Vertical Maximum Extent, m");CHKERRQ(ierr); 1829 ierr = PetscBagRegisterScalar(bag, &p->ymin, 0.0, "ymin", "Vertical Minimum Extent, m");CHKERRQ(ierr); 1830 ierr = PetscBagRegisterScalar(bag, &p->xmax, 1.0, "xmax", "Horizontal Maximum Extent, m");CHKERRQ(ierr); 1831 ierr = PetscBagRegisterScalar(bag, &p->xmin, 0.0, "xmin", "Horizontal Minimum Extent, m");CHKERRQ(ierr); 1832 ierr = PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa");CHKERRQ(ierr); 1833 } else if (ctx->solType == SOL_CRYER) { 1834 // Realistic values - Mandel 1835 ierr = PetscBagRegisterScalar(bag, &p->mu, 0.75, "mu", "Shear Modulus, Pa");CHKERRQ(ierr); 1836 ierr = PetscBagRegisterScalar(bag, &p->K_u, 2.6941176470588233, "K_u", "Undrained Bulk Modulus, Pa");CHKERRQ(ierr); 1837 ierr = PetscBagRegisterScalar(bag, &p->alpha, 0.6, "alpha", "Biot Effective Stress Coefficient, -");CHKERRQ(ierr); 1838 ierr = PetscBagRegisterScalar(bag, &p->M, 4.705882352941176, "M", "Biot Modulus, Pa");CHKERRQ(ierr); 1839 ierr = PetscBagRegisterScalar(bag, &p->k, 1.5, "k", "Isotropic Permeability, m**2");CHKERRQ(ierr); 1840 ierr = PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s");CHKERRQ(ierr); 1841 ierr = PetscBagRegisterScalar(bag, &p->zmax, 1.0, "zmax", "Depth Maximum Extent, m");CHKERRQ(ierr); 1842 ierr = PetscBagRegisterScalar(bag, &p->zmin, 0.0, "zmin", "Depth Minimum Extent, m");CHKERRQ(ierr); 1843 ierr = PetscBagRegisterScalar(bag, &p->ymax, 1.0, "ymax", "Vertical Maximum Extent, m");CHKERRQ(ierr); 1844 ierr = PetscBagRegisterScalar(bag, &p->ymin, 0.0, "ymin", "Vertical Minimum Extent, m");CHKERRQ(ierr); 1845 ierr = PetscBagRegisterScalar(bag, &p->xmax, 1.0, "xmax", "Horizontal Maximum Extent, m");CHKERRQ(ierr); 1846 ierr = PetscBagRegisterScalar(bag, &p->xmin, 0.0, "xmin", "Horizontal Minimum Extent, m");CHKERRQ(ierr); 1847 ierr = PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa");CHKERRQ(ierr); 1848 } else { 1849 // Nonsense values 1850 ierr = PetscBagRegisterScalar(bag, &p->mu, 1.0, "mu", "Shear Modulus, Pa");CHKERRQ(ierr); 1851 ierr = PetscBagRegisterScalar(bag, &p->K_u, 1.0, "K_u", "Undrained Bulk Modulus, Pa");CHKERRQ(ierr); 1852 ierr = PetscBagRegisterScalar(bag, &p->alpha, 1.0, "alpha", "Biot Effective Stress Coefficient, -");CHKERRQ(ierr); 1853 ierr = PetscBagRegisterScalar(bag, &p->M, 1.0, "M", "Biot Modulus, Pa");CHKERRQ(ierr); 1854 ierr = PetscBagRegisterScalar(bag, &p->k, 1.0, "k", "Isotropic Permeability, m**2");CHKERRQ(ierr); 1855 ierr = PetscBagRegisterScalar(bag, &p->mu_f, 1.0, "mu_f", "Fluid Dynamic Viscosity, Pa*s");CHKERRQ(ierr); 1856 ierr = PetscBagRegisterScalar(bag, &p->zmax, 1.0, "zmax", "Depth Maximum Extent, m");CHKERRQ(ierr); 1857 ierr = PetscBagRegisterScalar(bag, &p->zmin, 0.0, "zmin", "Depth Minimum Extent, m");CHKERRQ(ierr); 1858 ierr = PetscBagRegisterScalar(bag, &p->ymax, 1.0, "ymax", "Vertical Maximum Extent, m");CHKERRQ(ierr); 1859 ierr = PetscBagRegisterScalar(bag, &p->ymin, 0.0, "ymin", "Vertical Minimum Extent, m");CHKERRQ(ierr); 1860 ierr = PetscBagRegisterScalar(bag, &p->xmax, 1.0, "xmax", "Horizontal Maximum Extent, m");CHKERRQ(ierr); 1861 ierr = PetscBagRegisterScalar(bag, &p->xmin, 0.0, "xmin", "Horizontal Minimum Extent, m");CHKERRQ(ierr); 1862 ierr = PetscBagRegisterScalar(bag, &p->P_0, 1.0, "P_0", "Magnitude of Vertical Stress, Pa");CHKERRQ(ierr); 1863 } 1864 ierr = PetscBagSetFromOptions(bag);CHKERRQ(ierr); 1865 { 1866 PetscScalar K_d = p->K_u - p->alpha*p->alpha*p->M; 1867 PetscScalar nu_u = (3.0*p->K_u - 2.0*p->mu) / (2.0*(3.0*p->K_u + p->mu)); 1868 PetscScalar nu = (3.0*K_d - 2.0*p->mu) / (2.0*(3.0*K_d + p->mu)); 1869 PetscScalar S = (3.0*p->K_u + 4.0*p->mu) / (p->M*(3.0*K_d + 4.0*p->mu)); 1870 PetscReal c = PetscRealPart((p->k/p->mu_f) / S); 1871 1872 PetscViewer viewer; 1873 PetscViewerFormat format; 1874 PetscBool flg; 1875 1876 switch (ctx->solType) { 1877 case SOL_QUADRATIC_LINEAR: 1878 case SOL_QUADRATIC_TRIG: 1879 case SOL_TRIG_LINEAR: ctx->t_r = PetscSqr(p->xmax - p->xmin)/c; break; 1880 case SOL_TERZAGHI: ctx->t_r = PetscSqr(2.0*(p->ymax - p->ymin))/c; break; 1881 case SOL_MANDEL: ctx->t_r = PetscSqr(2.0*(p->ymax - p->ymin))/c; break; 1882 case SOL_CRYER: ctx->t_r = PetscSqr(p->ymax)/c; break; 1883 default: SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%D)", solutionTypes[PetscMin(ctx->solType, NUM_SOLUTION_TYPES)], ctx->solType); 1884 } 1885 ierr = PetscOptionsGetViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg);CHKERRQ(ierr); 1886 if (flg) { 1887 ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr); 1888 ierr = PetscBagView(bag, viewer);CHKERRQ(ierr); 1889 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1890 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 1891 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 1892 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))); 1893 ierr = PetscPrintf(comm, " Relaxation time: %g\n", ctx->t_r); 1894 } 1895 } 1896 PetscFunctionReturn(0); 1897 } 1898 1899 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 1900 { 1901 Parameter *param; 1902 PetscBool flg; 1903 PetscErrorCode ierr; 1904 1905 PetscFunctionBeginUser; 1906 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr); 1907 if (user->solType == SOL_CRYER) { 1908 DM rdm; 1909 1910 if (!user->simplex) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot create ball with cubic cells"); 1911 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); 1912 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); 1913 ierr = DMPlexCreateBallMesh(comm, user->dim, param->xmax, dm);CHKERRQ(ierr); 1914 1915 ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr); 1916 ierr = DMPlexSetRefinementLimit(*dm, user->refLimit);CHKERRQ(ierr); 1917 ierr = DMRefine(*dm, comm, &rdm);CHKERRQ(ierr); 1918 if (rdm) { 1919 ierr = DMDestroy(dm);CHKERRQ(ierr); 1920 *dm = rdm; 1921 } 1922 ierr = DMPlexSetRefinementUniform(*dm, PETSC_TRUE);CHKERRQ(ierr); 1923 } else if (user->solType == SOL_MANDEL) { 1924 PetscReal lower[2], upper[2]; 1925 1926 lower[0] = param->xmin - (param->xmax - param->xmin) / 2.0; 1927 lower[1] = param->ymin - (param->ymax - param->ymin) / 2.0; 1928 upper[0] = param->xmax - (param->xmax - param->xmin) / 2.0; 1929 upper[1] = param->ymax - (param->ymax - param->ymin) / 2.0; 1930 //reset min / max values for mandel 1931 param->xmin = lower[0]; 1932 param->ymin = lower[1]; 1933 param->xmax = upper[0]; 1934 param->ymax = upper[1]; 1935 ierr = DMPlexCreateBoxMesh(comm, user->dim, user->simplex, NULL, lower, upper, NULL, PETSC_TRUE, dm);CHKERRQ(ierr); 1936 } else { 1937 Parameter *param; 1938 PetscReal lower[3], upper[3]; 1939 1940 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr); 1941 lower[0] = param->xmin; 1942 lower[1] = param->ymin; 1943 lower[2] = param->zmin; 1944 upper[0] = param->xmax; 1945 upper[1] = param->ymax; 1946 upper[2] = param->zmax; 1947 ierr = DMPlexCreateBoxMesh(comm, user->dim, user->simplex, NULL, lower, upper, NULL, PETSC_TRUE, dm);CHKERRQ(ierr); 1948 } 1949 { 1950 DM pdm = NULL; 1951 PetscPartitioner part; 1952 1953 ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr); 1954 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 1955 ierr = DMPlexDistribute(*dm, 0, NULL, &pdm);CHKERRQ(ierr); 1956 if (pdm) { 1957 ierr = DMDestroy(dm);CHKERRQ(ierr); 1958 *dm = pdm; 1959 } 1960 } 1961 ierr = PetscStrcmp(user->dmType, DMPLEX, &flg);CHKERRQ(ierr); 1962 if (flg) { 1963 DM ndm; 1964 1965 ierr = DMConvert(*dm, user->dmType, &ndm);CHKERRQ(ierr); 1966 if (ndm) { 1967 ierr = DMDestroy(dm);CHKERRQ(ierr); 1968 *dm = ndm; 1969 } 1970 } 1971 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 1972 1973 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr); 1974 ierr = DMSetApplicationContext(*dm, user);CHKERRQ(ierr); 1975 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 1976 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 1977 PetscFunctionReturn(0); 1978 } 1979 1980 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 1981 { 1982 PetscErrorCode (*exact[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 1983 PetscErrorCode (*exact_t[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 1984 PetscDS ds; 1985 DMLabel label; 1986 PetscWeakForm wf; 1987 Parameter *param; 1988 PetscInt id_mandel[2]; 1989 PetscInt comp[1]; 1990 PetscInt comp_mandel[2]; 1991 PetscInt dim, id, bd, f; 1992 PetscErrorCode ierr; 1993 1994 PetscFunctionBeginUser; 1995 ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr); 1996 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 1997 ierr = PetscDSGetSpatialDimension(ds, &dim);CHKERRQ(ierr); 1998 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr); 1999 exact_t[0] = exact_t[1] = exact_t[2] = zero; 2000 2001 /* Setup Problem Formulation and Boundary Conditions */ 2002 switch (user->solType) { 2003 case SOL_QUADRATIC_LINEAR: 2004 ierr = PetscDSSetResidual(ds, 0, f0_quadratic_linear_u, f1_u);CHKERRQ(ierr); 2005 ierr = PetscDSSetResidual(ds, 1, f0_epsilon, NULL);CHKERRQ(ierr); 2006 ierr = PetscDSSetResidual(ds, 2, f0_quadratic_linear_p, f1_p);CHKERRQ(ierr); 2007 ierr = PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr); 2008 ierr = PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL);CHKERRQ(ierr); 2009 ierr = PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL);CHKERRQ(ierr); 2010 ierr = PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL);CHKERRQ(ierr); 2011 ierr = PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL);CHKERRQ(ierr); 2012 ierr = PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL);CHKERRQ(ierr); 2013 ierr = PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp);CHKERRQ(ierr); 2014 exact[0] = quadratic_u; 2015 exact[1] = linear_eps; 2016 exact[2] = linear_linear_p; 2017 exact_t[2] = linear_linear_p_t; 2018 2019 id = 1; 2020 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", label, 1, &id, 0, 0, NULL, (void (*)(void)) exact[0], NULL, user, NULL);CHKERRQ(ierr); 2021 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", label, 1, &id, 2, 0, NULL, (void (*)(void)) exact[2], (void (*)(void)) exact_t[2], user, NULL);CHKERRQ(ierr); 2022 break; 2023 case SOL_TRIG_LINEAR: 2024 ierr = PetscDSSetResidual(ds, 0, f0_trig_linear_u, f1_u);CHKERRQ(ierr); 2025 ierr = PetscDSSetResidual(ds, 1, f0_epsilon, NULL);CHKERRQ(ierr); 2026 ierr = PetscDSSetResidual(ds, 2, f0_trig_linear_p, f1_p);CHKERRQ(ierr); 2027 ierr = PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr); 2028 ierr = PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL);CHKERRQ(ierr); 2029 ierr = PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL);CHKERRQ(ierr); 2030 ierr = PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL);CHKERRQ(ierr); 2031 ierr = PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL);CHKERRQ(ierr); 2032 ierr = PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL);CHKERRQ(ierr); 2033 ierr = PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp);CHKERRQ(ierr); 2034 exact[0] = trig_u; 2035 exact[1] = trig_eps; 2036 exact[2] = trig_linear_p; 2037 exact_t[2] = trig_linear_p_t; 2038 2039 id = 1; 2040 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", label, 1, &id, 0, 0, NULL, (void (*)(void)) exact[0], NULL, user, NULL);CHKERRQ(ierr); 2041 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", label, 1, &id, 2, 0, NULL, (void (*)(void)) exact[2], (void (*)(void)) exact_t[2], user, NULL);CHKERRQ(ierr); 2042 break; 2043 case SOL_QUADRATIC_TRIG: 2044 ierr = PetscDSSetResidual(ds, 0, f0_quadratic_trig_u, f1_u);CHKERRQ(ierr); 2045 ierr = PetscDSSetResidual(ds, 1, f0_epsilon, NULL);CHKERRQ(ierr); 2046 ierr = PetscDSSetResidual(ds, 2, f0_quadratic_trig_p, f1_p);CHKERRQ(ierr); 2047 ierr = PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr); 2048 ierr = PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL);CHKERRQ(ierr); 2049 ierr = PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL);CHKERRQ(ierr); 2050 ierr = PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL);CHKERRQ(ierr); 2051 ierr = PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL);CHKERRQ(ierr); 2052 ierr = PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL);CHKERRQ(ierr); 2053 ierr = PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp);CHKERRQ(ierr); 2054 exact[0] = quadratic_u; 2055 exact[1] = linear_eps; 2056 exact[2] = linear_trig_p; 2057 exact_t[2] = linear_trig_p_t; 2058 2059 id = 1; 2060 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall displacement", label, 1, &id, 0, 0, NULL, (void (*)(void)) exact[0], NULL, user, NULL);CHKERRQ(ierr); 2061 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall pressure", label, 1, &id, 2, 0, NULL, (void (*)(void)) exact[2], (void (*)(void)) exact_t[2], user, NULL);CHKERRQ(ierr); 2062 break; 2063 case SOL_TERZAGHI: 2064 ierr = PetscDSSetResidual(ds, 0, NULL, f1_u);CHKERRQ(ierr); 2065 ierr = PetscDSSetResidual(ds, 1, f0_epsilon, NULL);CHKERRQ(ierr); 2066 ierr = PetscDSSetResidual(ds, 2, f0_p, f1_p);CHKERRQ(ierr); 2067 ierr = PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr); 2068 ierr = PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL);CHKERRQ(ierr); 2069 ierr = PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL);CHKERRQ(ierr); 2070 ierr = PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL);CHKERRQ(ierr); 2071 ierr = PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL);CHKERRQ(ierr); 2072 ierr = PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL);CHKERRQ(ierr); 2073 ierr = PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp);CHKERRQ(ierr); 2074 2075 exact[0] = terzaghi_2d_u; 2076 exact[1] = terzaghi_2d_eps; 2077 exact[2] = terzaghi_2d_p; 2078 exact_t[0] = terzaghi_2d_u_t; 2079 exact_t[1] = terzaghi_2d_eps_t; 2080 exact_t[2] = terzaghi_2d_p_t; 2081 2082 id = 1; 2083 ierr = DMAddBoundary(dm, DM_BC_NATURAL, "vertical stress", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd);CHKERRQ(ierr); 2084 ierr = PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 2085 ierr = PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, f0_terzaghi_bd_u, 0, NULL);CHKERRQ(ierr); 2086 2087 id = 3; 2088 comp[0] = 1; 2089 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed base", label, 1, &id, 0, 1, comp, (void (*)(void)) zero, NULL, user, NULL);CHKERRQ(ierr); 2090 id = 2; 2091 comp[0] = 0; 2092 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed side", label, 1, &id, 0, 1, comp, (void (*)(void)) zero, NULL, user, NULL);CHKERRQ(ierr); 2093 id = 4; 2094 comp[0] = 0; 2095 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed side", label, 1, &id, 0, 1, comp, (void (*)(void)) zero, NULL, user, NULL);CHKERRQ(ierr); 2096 id = 1; 2097 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", label, 1, &id, 2, 0, NULL, (void (*)(void)) terzaghi_drainage_pressure, NULL, user, NULL);CHKERRQ(ierr); 2098 break; 2099 case SOL_MANDEL: 2100 ierr = PetscDSSetResidual(ds, 0, NULL, f1_u);CHKERRQ(ierr); 2101 ierr = PetscDSSetResidual(ds, 1, f0_epsilon, NULL);CHKERRQ(ierr); 2102 ierr = PetscDSSetResidual(ds, 2, f0_p, f1_p);CHKERRQ(ierr); 2103 ierr = PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr); 2104 ierr = PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL);CHKERRQ(ierr); 2105 ierr = PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL);CHKERRQ(ierr); 2106 ierr = PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL);CHKERRQ(ierr); 2107 ierr = PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL);CHKERRQ(ierr); 2108 ierr = PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL);CHKERRQ(ierr); 2109 ierr = PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp);CHKERRQ(ierr); 2110 2111 ierr = mandelZeros(PETSC_COMM_WORLD, user, param);CHKERRQ(ierr); 2112 2113 exact[0] = mandel_2d_u; 2114 exact[1] = mandel_2d_eps; 2115 exact[2] = mandel_2d_p; 2116 exact_t[0] = mandel_2d_u_t; 2117 exact_t[1] = mandel_2d_eps_t; 2118 exact_t[2] = mandel_2d_p_t; 2119 2120 id_mandel[0] = 3; 2121 id_mandel[1] = 1; 2122 //comp[0] = 1; 2123 comp_mandel[0] = 0; 2124 comp_mandel[1] = 1; 2125 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "vertical stress", label, 2, id_mandel, 0, 2, comp_mandel, (void (*)(void)) mandel_2d_u, NULL, user, NULL);CHKERRQ(ierr); 2126 //ierr = DMAddBoundary(dm, DM_BC_NATURAL, "vertical stress", "marker", 0, 1, comp, NULL, 2, id_mandel, user);CHKERRQ(ierr); 2127 //ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed base", "marker", 0, 1, comp, (void (*)(void)) zero, 2, id_mandel, user);CHKERRQ(ierr); 2128 //ierr = PetscDSSetBdResidual(ds, 0, f0_mandel_bd_u, NULL);CHKERRQ(ierr); 2129 2130 id_mandel[0] = 2; 2131 id_mandel[1] = 4; 2132 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", label, 2, id_mandel, 2, 0, NULL, (void (*)(void)) zero, NULL, user, NULL);CHKERRQ(ierr); 2133 break; 2134 case SOL_CRYER: 2135 ierr = PetscDSSetResidual(ds, 0, NULL, f1_u);CHKERRQ(ierr); 2136 ierr = PetscDSSetResidual(ds, 1, f0_epsilon, NULL);CHKERRQ(ierr); 2137 ierr = PetscDSSetResidual(ds, 2, f0_p, f1_p);CHKERRQ(ierr); 2138 ierr = PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr); 2139 ierr = PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_ue, NULL);CHKERRQ(ierr); 2140 ierr = PetscDSSetJacobian(ds, 0, 2, NULL, NULL, g2_up, NULL);CHKERRQ(ierr); 2141 ierr = PetscDSSetJacobian(ds, 1, 0, NULL, g1_eu, NULL, NULL);CHKERRQ(ierr); 2142 ierr = PetscDSSetJacobian(ds, 1, 1, g0_ee, NULL, NULL, NULL);CHKERRQ(ierr); 2143 ierr = PetscDSSetJacobian(ds, 2, 1, g0_pe, NULL, NULL, NULL);CHKERRQ(ierr); 2144 ierr = PetscDSSetJacobian(ds, 2, 2, g0_pp, NULL, NULL, g3_pp);CHKERRQ(ierr); 2145 2146 ierr = cryerZeros(PETSC_COMM_WORLD, user, param);CHKERRQ(ierr); 2147 2148 exact[0] = cryer_3d_u; 2149 exact[1] = cryer_3d_eps; 2150 exact[2] = cryer_3d_p; 2151 2152 id = 1; 2153 ierr = DMAddBoundary(dm, DM_BC_NATURAL, "normal stress", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd);CHKERRQ(ierr); 2154 ierr = PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 2155 ierr = PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, f0_cryer_bd_u, 0, NULL);CHKERRQ(ierr); 2156 2157 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "drained surface", label, 1, &id, 2, 0, NULL, (void (*)(void)) cryer_drainage_pressure, NULL, user, NULL);CHKERRQ(ierr); 2158 break; 2159 default: SETERRQ2(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%D)", solutionTypes[PetscMin(user->solType, NUM_SOLUTION_TYPES)], user->solType); 2160 } 2161 for (f = 0; f < 3; ++f) { 2162 ierr = PetscDSSetExactSolution(ds, f, exact[f], user);CHKERRQ(ierr); 2163 ierr = PetscDSSetExactSolutionTimeDerivative(ds, f, exact_t[f], user);CHKERRQ(ierr); 2164 } 2165 2166 /* Setup constants */ 2167 { 2168 PetscScalar constants[6]; 2169 constants[0] = param->mu; /* shear modulus, Pa */ 2170 constants[1] = param->K_u; /* undrained bulk modulus, Pa */ 2171 constants[2] = param->alpha; /* Biot effective stress coefficient, - */ 2172 constants[3] = param->M; /* Biot modulus, Pa */ 2173 constants[4] = param->k/param->mu_f; /* Darcy coefficient, m**2 / Pa*s */ 2174 constants[5] = param->P_0; /* Magnitude of Vertical Stress, Pa */ 2175 ierr = PetscDSSetConstants(ds, 6, constants);CHKERRQ(ierr); 2176 } 2177 PetscFunctionReturn(0); 2178 } 2179 2180 static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace) 2181 { 2182 PetscErrorCode ierr; 2183 2184 PetscFunctionBegin; 2185 ierr = DMPlexCreateRigidBody(dm, origField, nullspace);CHKERRQ(ierr); 2186 PetscFunctionReturn(0); 2187 } 2188 2189 static PetscErrorCode SetupFE(DM dm, PetscBool simplex, PetscInt Nf, PetscInt Nc[], const char *name[], PetscErrorCode (*setup)(DM, AppCtx *), void *ctx) 2190 { 2191 AppCtx *user = (AppCtx *) ctx; 2192 DM cdm = dm; 2193 PetscFE fe; 2194 PetscQuadrature q = NULL; 2195 char prefix[PETSC_MAX_PATH_LEN]; 2196 PetscInt dim, f; 2197 PetscErrorCode ierr; 2198 2199 PetscFunctionBegin; 2200 /* Create finite element */ 2201 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2202 for (f = 0; f < Nf; ++f) { 2203 ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name[f]);CHKERRQ(ierr); 2204 ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc[f], simplex, name[f] ? prefix : NULL, -1, &fe);CHKERRQ(ierr); 2205 ierr = PetscObjectSetName((PetscObject) fe, name[f]);CHKERRQ(ierr); 2206 if (!q) {ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);} 2207 ierr = PetscFESetQuadrature(fe, q);CHKERRQ(ierr); 2208 ierr = DMSetField(dm, f, NULL, (PetscObject) fe);CHKERRQ(ierr); 2209 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 2210 } 2211 ierr = DMCreateDS(dm);CHKERRQ(ierr); 2212 ierr = (*setup)(dm, user);CHKERRQ(ierr); 2213 while (cdm) { 2214 ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr); 2215 if (0) {ierr = DMSetNearNullSpaceConstructor(cdm, 0, CreateElasticityNullSpace);CHKERRQ(ierr);} 2216 /* TODO: Check whether the boundary of coarse meshes is marked */ 2217 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 2218 } 2219 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 2220 PetscFunctionReturn(0); 2221 } 2222 2223 static PetscErrorCode SetInitialConditions(TS ts, Vec u) 2224 { 2225 DM dm; 2226 PetscReal t; 2227 PetscErrorCode ierr; 2228 2229 PetscFunctionBegin; 2230 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 2231 ierr = TSGetTime(ts, &t);CHKERRQ(ierr); 2232 if (t <= 0.0) { 2233 PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *); 2234 void *ctxs[3]; 2235 AppCtx *ctx; 2236 2237 ierr = DMGetApplicationContext(dm, (void **) &ctx);CHKERRQ(ierr); 2238 switch (ctx->solType) { 2239 case SOL_TERZAGHI: 2240 funcs[0] = terzaghi_initial_u; ctxs[0] = ctx; 2241 funcs[1] = terzaghi_initial_eps; ctxs[1] = ctx; 2242 funcs[2] = terzaghi_drainage_pressure; ctxs[2] = ctx; 2243 ierr = DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u);CHKERRQ(ierr); 2244 break; 2245 case SOL_MANDEL: 2246 funcs[0] = mandel_initial_u; ctxs[0] = ctx; 2247 funcs[1] = mandel_initial_eps; ctxs[1] = ctx; 2248 funcs[2] = mandel_drainage_pressure; ctxs[2] = ctx; 2249 ierr = DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u);CHKERRQ(ierr); 2250 break; 2251 case SOL_CRYER: 2252 funcs[0] = cryer_initial_u; ctxs[0] = ctx; 2253 funcs[1] = cryer_initial_eps; ctxs[1] = ctx; 2254 funcs[2] = cryer_drainage_pressure; ctxs[2] = ctx; 2255 ierr = DMProjectFunction(dm, t, funcs, ctxs, INSERT_VALUES, u);CHKERRQ(ierr); 2256 break; 2257 default: 2258 ierr = DMComputeExactSolution(dm, t, u, NULL);CHKERRQ(ierr); 2259 } 2260 } else { 2261 ierr = DMComputeExactSolution(dm, t, u, NULL);CHKERRQ(ierr); 2262 } 2263 PetscFunctionReturn(0); 2264 } 2265 2266 /* Need to create Viewer each time because HDF5 can get corrupted */ 2267 static PetscErrorCode SolutionMonitor(TS ts, PetscInt steps, PetscReal time, Vec u, void *mctx) 2268 { 2269 DM dm; 2270 Vec exact; 2271 PetscViewer viewer; 2272 PetscViewerFormat format; 2273 PetscOptions options; 2274 const char *prefix; 2275 PetscErrorCode ierr; 2276 2277 PetscFunctionBegin; 2278 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 2279 ierr = PetscObjectGetOptions((PetscObject) ts, &options);CHKERRQ(ierr); 2280 ierr = PetscObjectGetOptionsPrefix((PetscObject) ts, &prefix);CHKERRQ(ierr); 2281 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject) ts), options, prefix, "-monitor_solution", &viewer, &format, NULL);CHKERRQ(ierr); 2282 ierr = DMGetGlobalVector(dm, &exact);CHKERRQ(ierr); 2283 ierr = DMComputeExactSolution(dm, time, exact, NULL);CHKERRQ(ierr); 2284 ierr = DMSetOutputSequenceNumber(dm, steps, time);CHKERRQ(ierr); 2285 ierr = VecView(exact, viewer);CHKERRQ(ierr); 2286 ierr = VecView(u, viewer);CHKERRQ(ierr); 2287 ierr = DMRestoreGlobalVector(dm, &exact);CHKERRQ(ierr); 2288 { 2289 PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 2290 void **ectxs; 2291 PetscReal *err; 2292 PetscInt Nf, f; 2293 2294 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 2295 ierr = PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err);CHKERRQ(ierr); 2296 { 2297 PetscInt Nds, s; 2298 2299 ierr = DMGetNumDS(dm, &Nds);CHKERRQ(ierr); 2300 for (s = 0; s < Nds; ++s) { 2301 PetscDS ds; 2302 DMLabel label; 2303 IS fieldIS; 2304 const PetscInt *fields; 2305 PetscInt dsNf, f; 2306 2307 ierr = DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds);CHKERRQ(ierr); 2308 ierr = PetscDSGetNumFields(ds, &dsNf);CHKERRQ(ierr); 2309 ierr = ISGetIndices(fieldIS, &fields);CHKERRQ(ierr); 2310 for (f = 0; f < dsNf; ++f) { 2311 const PetscInt field = fields[f]; 2312 ierr = PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]);CHKERRQ(ierr); 2313 } 2314 ierr = ISRestoreIndices(fieldIS, &fields);CHKERRQ(ierr); 2315 } 2316 } 2317 ierr = DMComputeL2FieldDiff(dm, time, exacts, ectxs, u, err);CHKERRQ(ierr); 2318 ierr = PetscPrintf(PetscObjectComm((PetscObject) ts), "Time: %g L_2 Error: [", time);CHKERRQ(ierr); 2319 for (f = 0; f < Nf; ++f) { 2320 if (f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) ts), ", ");CHKERRQ(ierr);} 2321 ierr = PetscPrintf(PetscObjectComm((PetscObject) ts), "%g", (double) err[f]);CHKERRQ(ierr); 2322 } 2323 ierr = PetscPrintf(PetscObjectComm((PetscObject) ts), "]\n");CHKERRQ(ierr); 2324 ierr = PetscFree3(exacts, ectxs, err);CHKERRQ(ierr); 2325 } 2326 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2327 PetscFunctionReturn(0); 2328 } 2329 2330 static PetscErrorCode SetupMonitor(TS ts, AppCtx *ctx) 2331 { 2332 PetscViewer viewer; 2333 PetscViewerFormat format; 2334 PetscOptions options; 2335 const char *prefix; 2336 PetscBool flg; 2337 PetscErrorCode ierr; 2338 2339 PetscFunctionBegin; 2340 ierr = PetscObjectGetOptions((PetscObject) ts, &options);CHKERRQ(ierr); 2341 ierr = PetscObjectGetOptionsPrefix((PetscObject) ts, &prefix);CHKERRQ(ierr); 2342 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject) ts), options, prefix, "-monitor_solution", &viewer, &format, &flg);CHKERRQ(ierr); 2343 if (flg) {ierr = TSMonitorSet(ts, SolutionMonitor, ctx, NULL);CHKERRQ(ierr);} 2344 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2345 PetscFunctionReturn(0); 2346 } 2347 2348 static PetscErrorCode TSAdaptChoose_Terzaghi(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept, PetscReal *wlte, PetscReal *wltea, PetscReal *wlter) 2349 { 2350 static PetscReal dtTarget = -1.0; 2351 PetscReal dtInitial; 2352 DM dm; 2353 AppCtx *ctx; 2354 PetscInt step; 2355 PetscErrorCode ierr; 2356 2357 PetscFunctionBegin; 2358 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 2359 ierr = DMGetApplicationContext(dm, (void **) &ctx);CHKERRQ(ierr); 2360 ierr = TSGetStepNumber(ts, &step);CHKERRQ(ierr); 2361 dtInitial = ctx->dtInitial < 0.0 ? 1.0e-4*ctx->t_r : ctx->dtInitial; 2362 if (!step) { 2363 if (PetscAbsReal(dtInitial - h) > PETSC_SMALL) { 2364 *accept = PETSC_FALSE; 2365 *next_h = dtInitial; 2366 dtTarget = h; 2367 } else { 2368 *accept = PETSC_TRUE; 2369 *next_h = dtTarget < 0.0 ? dtInitial : dtTarget; 2370 dtTarget = -1.0; 2371 } 2372 } else { 2373 *accept = PETSC_TRUE; 2374 *next_h = h; 2375 } 2376 *next_sc = 0; /* Reuse the same order scheme */ 2377 *wlte = -1; /* Weighted local truncation error was not evaluated */ 2378 *wltea = -1; /* Weighted absolute local truncation error was not evaluated */ 2379 *wlter = -1; /* Weighted relative local truncation error was not evaluated */ 2380 PetscFunctionReturn(0); 2381 } 2382 2383 int main(int argc, char **argv) 2384 { 2385 AppCtx ctx; /* User-defined work context */ 2386 DM dm; /* Problem specification */ 2387 TS ts; /* Time Series / Nonlinear solver */ 2388 Vec u; /* Solutions */ 2389 const char *name[3] = {"displacement", "tracestrain", "pressure"}; 2390 PetscReal t; 2391 PetscInt Nc[3]; 2392 PetscErrorCode ierr; 2393 2394 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 2395 ierr = ProcessOptions(PETSC_COMM_WORLD, &ctx);CHKERRQ(ierr); 2396 ierr = PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &ctx.bag);CHKERRQ(ierr); 2397 ierr = PetscMalloc1(ctx.niter, &ctx.zeroArray);CHKERRQ(ierr); 2398 ierr = SetupParameters(PETSC_COMM_WORLD, &ctx);CHKERRQ(ierr); 2399 /* Primal System */ 2400 ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr); 2401 ierr = CreateMesh(PETSC_COMM_WORLD, &ctx, &dm);CHKERRQ(ierr); 2402 ierr = DMSetApplicationContext(dm, &ctx);CHKERRQ(ierr); 2403 ierr = TSSetDM(ts, dm);CHKERRQ(ierr); 2404 2405 Nc[0] = ctx.dim; 2406 Nc[1] = 1; 2407 Nc[2] = 1; 2408 2409 ierr = SetupFE(dm, ctx.simplex, 3, Nc, name, SetupPrimalProblem, &ctx);CHKERRQ(ierr); 2410 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr); 2411 ierr = DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);CHKERRQ(ierr); 2412 ierr = DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);CHKERRQ(ierr); 2413 ierr = DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);CHKERRQ(ierr); 2414 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 2415 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 2416 ierr = TSSetComputeInitialCondition(ts, SetInitialConditions);CHKERRQ(ierr); 2417 ierr = SetupMonitor(ts, &ctx);CHKERRQ(ierr); 2418 2419 if (ctx.solType != SOL_QUADRATIC_TRIG) { 2420 TSAdapt adapt; 2421 2422 ierr = TSGetAdapt(ts, &adapt);CHKERRQ(ierr); 2423 adapt->ops->choose = TSAdaptChoose_Terzaghi; 2424 } 2425 if (ctx.solType == SOL_CRYER) { 2426 Mat J; 2427 MatNullSpace sp; 2428 2429 ierr = TSSetUp(ts);CHKERRQ(ierr); 2430 ierr = TSGetIJacobian(ts, &J, NULL, NULL, NULL);CHKERRQ(ierr); 2431 ierr = DMPlexCreateRigidBody(dm, 0, &sp);CHKERRQ(ierr); 2432 ierr = MatSetNullSpace(J, sp);CHKERRQ(ierr); 2433 ierr = MatNullSpaceDestroy(&sp);CHKERRQ(ierr); 2434 } 2435 ierr = TSGetTime(ts, &t);CHKERRQ(ierr); 2436 ierr = DMSetOutputSequenceNumber(dm, 0, t);CHKERRQ(ierr); 2437 ierr = DMTSCheckFromOptions(ts, u);CHKERRQ(ierr); 2438 ierr = SetInitialConditions(ts, u);CHKERRQ(ierr); 2439 ierr = PetscObjectSetName((PetscObject) u, "solution");CHKERRQ(ierr); 2440 ierr = TSSolve(ts, u);CHKERRQ(ierr); 2441 ierr = DMTSCheckFromOptions(ts, u);CHKERRQ(ierr); 2442 ierr = TSGetSolution(ts, &u);CHKERRQ(ierr); 2443 ierr = VecViewFromOptions(u, NULL, "-sol_vec_view");CHKERRQ(ierr); 2444 2445 /* Cleanup */ 2446 ierr = VecDestroy(&u);CHKERRQ(ierr); 2447 ierr = TSDestroy(&ts);CHKERRQ(ierr); 2448 ierr = DMDestroy(&dm);CHKERRQ(ierr); 2449 ierr = PetscBagDestroy(&ctx.bag);CHKERRQ(ierr); 2450 ierr = PetscFree(ctx.zeroArray);CHKERRQ(ierr); 2451 ierr = PetscFinalize(); 2452 return ierr; 2453 } 2454 2455 /*TEST 2456 2457 test: 2458 suffix: 2d_quad_linear 2459 requires: triangle 2460 args: -sol_type quadratic_linear -dm_refine 2 \ 2461 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2462 -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme 2463 2464 test: 2465 suffix: 3d_quad_linear 2466 requires: ctetgen 2467 args: -dim 3 -sol_type quadratic_linear -dm_refine 1 \ 2468 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2469 -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme 2470 2471 test: 2472 suffix: 2d_trig_linear 2473 requires: triangle 2474 args: -sol_type trig_linear -dm_refine 1 \ 2475 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2476 -dmts_check .0001 -ts_max_steps 5 -ts_dt 0.00001 -ts_monitor_extreme 2477 2478 test: 2479 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9, 2.1, 1.8] 2480 suffix: 2d_trig_linear_sconv 2481 requires: triangle 2482 args: -sol_type trig_linear -dm_refine 1 \ 2483 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2484 -convest_num_refine 1 -ts_convergence_estimate -ts_convergence_temporal 0 -ts_max_steps 1 -ts_dt 0.00001 -pc_type lu 2485 2486 test: 2487 suffix: 3d_trig_linear 2488 requires: ctetgen 2489 args: -dim 3 -sol_type trig_linear -dm_refine 1 \ 2490 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2491 -dmts_check .0001 -ts_max_steps 2 -ts_monitor_extreme 2492 2493 test: 2494 # -dm_refine 1 -convest_num_refine 2 gets L_2 convergence rate: [2.0, 2.1, 1.9] 2495 suffix: 3d_trig_linear_sconv 2496 requires: ctetgen 2497 args: -dim 3 -sol_type trig_linear -dm_refine 1 \ 2498 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2499 -convest_num_refine 1 -ts_convergence_estimate -ts_convergence_temporal 0 -ts_max_steps 1 -pc_type lu 2500 2501 test: 2502 suffix: 2d_quad_trig 2503 requires: triangle 2504 args: -sol_type quadratic_trig -dm_refine 2 \ 2505 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2506 -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme 2507 2508 test: 2509 # Using -dm_refine 4 gets the convergence rates to [0.95, 0.97, 0.90] 2510 suffix: 2d_quad_trig_tconv 2511 requires: triangle 2512 args: -sol_type quadratic_trig -dm_refine 1 \ 2513 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2514 -convest_num_refine 3 -ts_convergence_estimate -ts_max_steps 5 -pc_type lu 2515 2516 test: 2517 suffix: 3d_quad_trig 2518 requires: ctetgen 2519 args: -dim 3 -sol_type quadratic_trig -dm_refine 1 \ 2520 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2521 -dmts_check .0001 -ts_max_steps 5 -ts_monitor_extreme 2522 2523 test: 2524 # Using -dm_refine 2 -convest_num_refine 3 gets the convergence rates to [1.0, 1.0, 1.0] 2525 suffix: 3d_quad_trig_tconv 2526 requires: ctetgen 2527 args: -dim 3 -sol_type quadratic_trig -dm_refine 1 \ 2528 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2529 -convest_num_refine 1 -ts_convergence_estimate -ts_max_steps 5 -pc_type lu 2530 2531 test: 2532 suffix: 2d_terzaghi 2533 requires: triangle 2534 args: -sol_type terzaghi -dm_plex_separate_marker -dm_plex_box_faces 1,8 -simplex 0 \ 2535 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 -niter 16000 \ 2536 -ts_dt 0.0028666667 -ts_max_steps 2 -ts_monitor -dmts_check .0001 -pc_type lu 2537 2538 test: 2539 # -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] 2540 suffix: 2d_terzaghi_tconv 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 0.023 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1 -pc_type lu 2545 2546 test: 2547 # -dm_plex_box_faces 1,16 -convest_num_refine 4 gives L_2 convergence rate: [1.7, 1.2, 1.1] 2548 # 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 2549 suffix: 2d_terzaghi_sconv 2550 requires: triangle 2551 args: -sol_type terzaghi -dm_plex_separate_marker -dm_plex_box_faces 1,8 -simplex 0 \ 2552 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 -niter 16000 \ 2553 -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 2554 2555 test: 2556 suffix: 2d_mandel 2557 requires: triangle 2558 args: -sol_type mandel -dm_plex_separate_marker -dm_refine 1 \ 2559 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2560 -ts_dt 0.0028666667 -ts_max_steps 2 -ts_monitor -dmts_check .0001 -pc_type lu 2561 2562 test: 2563 # -dm_refine 5 -ts_max_steps 4 -convest_num_refine 3 gives L_2 convergence rate: [0.26, -0.0058, 0.26] 2564 suffix: 2d_mandel_tconv 2565 requires: triangle 2566 args: -sol_type mandel -dm_plex_separate_marker -dm_refine 1 \ 2567 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2568 -ts_dt 0.023 -ts_max_steps 2 -ts_convergence_estimate -convest_num_refine 1 -pc_type lu 2569 2570 test: 2571 suffix: 3d_cryer 2572 requires: ctetgen !complex 2573 args: -sol_type cryer \ 2574 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2575 -ts_dt 0.0028666667 -ts_max_time 0.014333 -ts_max_steps 2 -dmts_check .0001 -pc_type svd 2576 2577 test: 2578 # Displacement and Pressure converge. The analytic expression for trace strain is inaccurate at the origin 2579 # -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] 2580 suffix: 3d_cryer_tconv 2581 requires: ctetgen !complex 2582 args: -sol_type cryer -bd_dm_refine 1 -ref_limit 0.00666667 \ 2583 -displacement_petscspace_degree 2 -tracestrain_petscspace_degree 1 -pressure_petscspace_degree 1 \ 2584 -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 2585 2586 TEST*/ 2587