Introduction
Introduction Statistics Contact Development Disclaimer Help
tmax_depth_simple_shear.c - cngf-pf - continuum model for granular flows with p…
git clone git://src.adamsgaard.dk/cngf-pf
Log
Files
Refs
README
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.141592653589793238462643383279502884197169399375105820974…
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] * s…
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…
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.…
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…
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,…
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 }
You are viewing proxied material from mx1.adamsgaard.dk. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.