tmax_depth_simple_shear.c - cngf-pf - continuum model for granular flows with pore-pressure dynamics (renamed from 1d_fd_simple_shear)
(HTM) git clone git://src.adamsgaard.dk/cngf-pf
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tmax_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 }