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