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