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