max_depth_simple_shear.c (5830B)
1 #include <unistd.h> 2 #include <stdio.h> 3 #include <stdlib.h> 4 #include <math.h> 5 #include <getopt.h> 6 #include <time.h> 7 #include <err.h> 8 9 #include "simulation.h" 10 #include "arg.h" 11 12 #define EPS 1e-12 13 14 /* depth accuracy criteria in meter for solver */ 15 #define TOL 1e-10 16 #define MAX_ITER 100 17 18 /* uncomment to print time spent per time step to stdout */ 19 /* #define BENCHMARK_PERFORMANCE */ 20 21 #ifndef M_PI 22 #define M_PI 3.1415926535897932384626433832795028841971693993751058209749445923078164062 23 #endif 24 25 char *argv0; 26 27 static void 28 usage(void) 29 { 30 fprintf(stderr, "usage: %s " 31 "[-a fluid-pressure-ampl] " 32 "[-C fluid-compressibility] " 33 "[-D fluid-diffusivity] " 34 "[-g gravity-accel] " 35 "[-h] " 36 "[-i fluid-viscosity] " 37 "[-k fluid-permeability] " 38 "[-P grain-compressibility] " 39 "[-p grain-porosity] " 40 "[-q fluid-pressure-freq] " 41 "[-R fluid-density] " 42 "[-r grain-density] " 43 "[-v] " 44 "\n", argv0); 45 exit(1); 46 } 47 48 static double 49 skin_depth(const struct simulation *sim) 50 { 51 if (sim->D > 0.0) 52 return sqrt(sim->D / (M_PI * sim->p_f_mod_freq)); 53 else 54 return sqrt(sim->k[0] / (sim->mu_f 55 * (sim->alpha + sim->phi[0] * sim->beta_f) 56 * M_PI * sim->p_f_mod_freq)); 57 } 58 59 /* using alternate form: sin(x) + cos(x) = sqrt(2)*sin(x + pi/4) */ 60 static double 61 eff_normal_stress_gradient(struct simulation * sim, double d_s, double z_) 62 { 63 if (z_ / d_s > 10.0) { 64 fprintf(stderr, "error: %s: unrealistic depth: %g m " 65 "relative to skin depth %g m\n", __func__, z_, d_s); 66 free_arrays(sim); 67 exit(10); 68 } 69 return sqrt(2.0) * sin((3.0 * M_PI / 2.0 - z_ / d_s) + M_PI / 4.0) 70 + (sim->rho_s - sim->rho_f) * sim->G * d_s 71 / sim->p_f_mod_ampl * exp(z_ / d_s); 72 } 73 74 /* use Brent's method for finding roots (Press et al., 2007) */ 75 static double 76 zbrent(struct simulation *sim, 77 double d_s, 78 double (*f) (struct simulation *sim, double, double), 79 double x1, 80 double x2, 81 double tol) 82 { 83 int iter; 84 double a, b, c, d = NAN, e = NAN, min1, min2, 85 fa, fb, fc, p, q, r, s, tol1, xm; 86 87 a = x1; 88 b = x2; 89 c = x2; 90 fa = (*f) (sim, d_s, a); 91 fb = (*f) (sim, d_s, b); 92 93 if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) { 94 warnx("%s: no root in range %g,%g m\n", 95 __func__, x1, x2); 96 free_arrays(sim); 97 exit(11); 98 } 99 fc = fb; 100 for (iter = 0; iter < MAX_ITER; iter++) { 101 if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) { 102 c = a; 103 fc = fa; 104 d = b - a; 105 e = d; 106 } 107 if (fabs(fc) < fabs(fb)) { 108 a = b; 109 b = c; 110 c = a; 111 fa = fb; 112 fb = fc; 113 fc = fa; 114 } 115 tol1 = 2.0 * EPS * fabs(b) + 0.5 * tol; 116 xm = 0.5 * (c - b); 117 if (fabs(xm) <= tol1 || fb == 0.0) 118 return b; 119 if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) { 120 s = fb / fa; 121 if (a == c) { 122 p = 2.0 * xm * s; 123 q = 1.0 - s; 124 } else { 125 q = fa / fc; 126 r = fb / fc; 127 p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0)); 128 q = (q - 1.0) * (r - 1.0) * (s - 1.0); 129 } 130 if (p > 0.0) 131 q = -q; 132 p = fabs(p); 133 min1 = 3.0 * xm * q - fabs(tol1 * q); 134 min2 = fabs(e * q); 135 if (2.0 * p < (min1 < min2 ? min1 : min2)) { 136 e = d; 137 d = p / q; 138 } else { 139 d = xm; 140 e = d; 141 } 142 } else { 143 d = xm; 144 e = d; 145 } 146 a = b; 147 fa = fb; 148 if (fabs(d) > tol1) 149 b += d; 150 else 151 b += ((xm) >= 0.0 ? fabs(tol1) : -fabs(tol1)); 152 fb = (*f) (sim, d_s, b); 153 } 154 warnx("error: %s: exceeded maximum number of iterations", __func__); 155 free_arrays(sim); 156 exit(12); 157 158 return NAN; 159 } 160 161 int 162 main(int argc, char *argv[]) 163 { 164 int i; 165 double new_phi, new_k; 166 double d_s, depth, depth_limit1, depth_limit2; 167 struct simulation sim; 168 169 #ifdef BENCHMARK_PERFORMANCE 170 clock_t t_begin, t_end; 171 double t_elapsed; 172 #endif 173 174 #ifdef __OpenBSD__ 175 if (pledge("stdio", NULL) == -1) 176 err(2, "pledge failed"); 177 #endif 178 179 init_sim(&sim); 180 new_phi = sim.phi[0]; 181 new_k = sim.k[0]; 182 183 ARGBEGIN { 184 case 'a': 185 sim.p_f_mod_ampl = atof(EARGF(usage())); 186 break; 187 case 'C': 188 sim.beta_f = atof(EARGF(usage())); 189 break; 190 case 'D': 191 sim.D = atof(EARGF(usage())); 192 break; 193 case 'g': 194 sim.G = atof(EARGF(usage())); 195 break; 196 case 'h': 197 usage(); 198 break; 199 case 'i': 200 sim.mu_f = atof(EARGF(usage())); 201 break; 202 case 'k': 203 new_k = atof(EARGF(usage())); 204 break; 205 case 'P': 206 sim.alpha = atof(EARGF(usage())); 207 break; 208 case 'p': 209 new_phi = atof(EARGF(usage())); 210 break; 211 case 'q': 212 sim.p_f_mod_freq = atof(EARGF(usage())); 213 break; 214 case 'R': 215 sim.rho_f = atof(EARGF(usage())); 216 break; 217 case 'r': 218 sim.rho_s = atof(EARGF(usage())); 219 break; 220 case 'v': 221 printf("%s-" VERSION "\n", argv0); 222 return 0; 223 break; 224 default: 225 usage(); 226 } ARGEND; 227 228 if (argc > 0) 229 usage(); 230 231 sim.nz = 2; 232 prepare_arrays(&sim); 233 234 if (!isnan(new_phi)) 235 for (i = 0; i < sim.nz; ++i) 236 sim.phi[i] = new_phi; 237 if (!isnan(new_k)) 238 for (i = 0; i < sim.nz; ++i) 239 sim.k[i] = new_k; 240 241 check_simulation_parameters(&sim); 242 243 depth = 0.0; 244 d_s = skin_depth(&sim); 245 246 #ifdef BENCHMARK_PERFORMANCE 247 t_begin = clock(); 248 #endif 249 250 /* deep deformatin does not occur with approximately zero 251 * amplitude in water pressure forcing, or if the stress gradient 252 * is positive at zero depth */ 253 if (fabs(sim.p_f_mod_ampl) > 1e-6 && 254 eff_normal_stress_gradient(&sim, d_s, 0.0) < 0.0) { 255 256 depth_limit1 = 0.0; 257 depth_limit2 = 5.0 * d_s; 258 259 depth = zbrent(&sim, 260 d_s, 261 (double (*) (struct simulation *, double, double)) 262 eff_normal_stress_gradient, 263 depth_limit1, 264 depth_limit2, 265 TOL); 266 } 267 #ifdef BENCHMARK_PERFORMANCE 268 t_end = clock(); 269 t_elapsed = (double) (t_end - t_begin) / CLOCKS_PER_SEC; 270 printf("time spent = %g s\n", t_elapsed); 271 #endif 272 273 printf("%.17g\t%.17g\n", depth, d_s); 274 275 free_arrays(&sim); 276 277 return 0; 278 }