t1d_fd_simple_shear.c - cngf-pf - continuum model for granular flows with pore-pressure dynamics
(HTM) git clone git://src.adamsgaard.dk/1d_fd_simple_shear
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
t1d_fd_simple_shear.c (5723B)
---
1 #include <stdlib.h>
2 #include <math.h>
3 #include <string.h>
4 #include <time.h>
5 #include <unistd.h>
6
7 #include "simulation.h"
8 #include "fluid.h"
9
10 #include "arg.h"
11
12 /* relative tolerance criteria for granular fluidity solver */
13 #define RTOL 1e-5
14 #define MAX_ITER_1D_FD_SIMPLE_SHEAR 100000
15
16 /* uncomment to print time spent per time step to stdout */
17 /* #define BENCHMARK_PERFORMANCE */
18
19 char *argv0;
20
21 static void
22 usage(void)
23 {
24 fprintf(stderr, "usage: %s "
25 "[-A grain-nonlocal-ampl] "
26 "[-a fluid-pressure-ampl] "
27 "[-b grain-rate-dependence] "
28 "[-C fluid-compressibility] "
29 "[-c grain-cohesion] "
30 "[-d grain-size] "
31 "[-D fluid-diffusivity] "
32 "[-e end-time] "
33 "[-F] "
34 "[-f applied-shear-friction] "
35 "[-g gravity-accel] "
36 "[-H fluid-pressure-phase] "
37 "[-h] "
38 "[-I file-interval] "
39 "[-i fluid-viscosity] "
40 "[-K dilatancy-constant] "
41 "[-k fluid-permeability] "
42 "[-L length] "
43 "[-l applied-shear-vel-limit] "
44 "[-m grain-friction] "
45 "[-N] "
46 "[-n normal-stress] "
47 "[-O fluid-pressure-top] "
48 "[-o origo] "
49 "[-P grain-compressibility] "
50 "[-p grain-porosity] "
51 "[-q fluid-pressure-freq] "
52 "[-R fluid-density] "
53 "[-r grain-density] "
54 "[-S fluid-pressure-pulse-shape] "
55 "[-s applied-shear-vel] "
56 "[-T] "
57 "[-t curr-time] "
58 "[-U resolution] "
59 "[-u fluid-pulse-time] "
60 "[-v] "
61 "[-Y max-porosity] "
62 "[-y min-porosity] "
63 "[name]\n", argv0);
64 exit(1);
65 }
66
67 int
68 main(int argc, char *argv[])
69 {
70 int i, normalize;
71 unsigned long iter;
72 double new_phi, new_k, filetimeclock;
73 struct simulation sim;
74
75 #ifdef BENCHMARK_PERFORMANCE
76 clock_t t_begin, t_end;
77 double t_elapsed;
78
79 #endif
80
81 #ifdef __OpenBSD__
82 if (pledge("stdio wpath cpath", NULL) == -1) {
83 fprintf(stderr, "error: pledge failed");
84 exit(2);
85 }
86 #endif
87
88 init_sim(&sim);
89 normalize = 0;
90 new_phi = sim.phi[0];
91 new_k = sim.k[0];
92
93 ARGBEGIN {
94 case 'A':
95 sim.A = atof(EARGF(usage()));
96 break;
97 case 'a':
98 sim.p_f_mod_ampl = atof(EARGF(usage()));
99 break;
100 case 'b':
101 sim.b = atof(EARGF(usage()));
102 break;
103 case 'C':
104 sim.beta_f = atof(EARGF(usage()));
105 break;
106 case 'c':
107 sim.C = atof(EARGF(usage()));
108 break;
109 case 'd':
110 sim.d = atof(EARGF(usage()));
111 break;
112 case 'D':
113 sim.D = atof(EARGF(usage()));
114 break;
115 case 'e':
116 sim.t_end = atof(EARGF(usage()));
117 break;
118 case 'F':
119 sim.fluid = 1;
120 break;
121 case 'f':
122 sim.mu_wall = atof(EARGF(usage()));
123 break;
124 case 'g':
125 sim.G = atof(EARGF(usage()));
126 break;
127 case 'H':
128 sim.p_f_mod_phase = atof(EARGF(usage()));
129 break;
130 case 'h':
131 usage();
132 break;
133 case 'I':
134 sim.file_dt = atof(EARGF(usage()));
135 break;
136 case 'i':
137 sim.mu_f = atof(EARGF(usage()));
138 break;
139 case 'K':
140 sim.dilatancy_angle = atof(EARGF(usage()));
141 break;
142 case 'k':
143 new_k = atof(EARGF(usage()));
144 break;
145 case 'L':
146 sim.L_z = atof(EARGF(usage()));
147 break;
148 case 'l':
149 sim.v_x_limit = atof(EARGF(usage()));
150 break;
151 case 'm':
152 sim.mu_s = atof(EARGF(usage()));
153 break;
154 case 'N':
155 normalize = 1;
156 break;
157 case 'n':
158 sim.P_wall = atof(EARGF(usage()));
159 break;
160 case 'O':
161 sim.p_f_top = atof(EARGF(usage()));
162 break;
163 case 'o':
164 sim.origo_z = atof(EARGF(usage()));
165 break;
166 case 'P':
167 sim.alpha = atof(EARGF(usage()));
168 break;
169 case 'p':
170 new_phi = atof(EARGF(usage()));
171 break;
172 case 'q':
173 sim.p_f_mod_freq = atof(EARGF(usage()));
174 break;
175 case 'R':
176 sim.rho_f = atof(EARGF(usage()));
177 break;
178 case 'r':
179 sim.rho_s = atof(EARGF(usage()));
180 break;
181 case 'S':
182 if (argv[1] == NULL)
183 usage();
184 else if (!strncmp(argv[1], "triangle", 8))
185 sim.p_f_mod_pulse_shape = 0;
186 else if (!strncmp(argv[1], "square", 6))
187 sim.p_f_mod_pulse_shape = 1;
188 else {
189 fprintf(stderr, "error: invalid pulse shape '%s'\n",
190 argv[1]);
191 return 1;
192 }
193 argc--;
194 argv++;
195 break;
196 case 's':
197 sim.v_x_fix = atof(EARGF(usage()));
198 break;
199 case 'T':
200 sim.transient = 1;
201 break;
202 case 't':
203 sim.t = atof(EARGF(usage()));
204 break;
205 case 'U':
206 sim.nz = atoi(EARGF(usage()));
207 break;
208 case 'u':
209 sim.p_f_mod_pulse_time = atof(EARGF(usage()));
210 break;
211 case 'V':
212 sim.v_x_bot = atof(EARGF(usage()));
213 break;
214 case 'v':
215 printf("%s-" VERSION "\n", argv0);
216 return 0;
217 case 'Y':
218 sim.phi_max = atof(EARGF(usage()));
219 break;
220 case 'y':
221 sim.phi_min = atof(EARGF(usage()));
222 break;
223 default:
224 usage();
225 } ARGEND;
226
227 if (argc == 1 && argv[0])
228 snprintf(sim.name, sizeof(sim.name), "%s", argv[0]);
229 else if (argc > 1)
230 usage();
231
232 if (sim.nz < 1)
233 sim.nz = (int) ceil(sim.L_z / sim.d);
234
235 prepare_arrays(&sim);
236
237 if (!isnan(new_phi))
238 for (i = 0; i < sim.nz; ++i)
239 sim.phi[i] = new_phi;
240 if (!isnan(new_k))
241 for (i = 0; i < sim.nz; ++i)
242 sim.k[i] = new_k;
243
244 lithostatic_pressure_distribution(&sim);
245
246 if (sim.fluid) {
247 hydrostatic_fluid_pressure_distribution(&sim);
248 if (set_largest_fluid_timestep(&sim, 0.5)) {
249 free_arrays(&sim);
250 return 20;
251 }
252 }
253 if (sim.dt > sim.file_dt)
254 sim.dt = sim.file_dt;
255 compute_effective_stress(&sim);
256
257 check_simulation_parameters(&sim);
258
259 filetimeclock = 0.0;
260 iter = 0;
261 do {
262
263 #ifdef BENCHMARK_PERFORMANCE
264 t_begin = clock();
265 #endif
266
267 if (coupled_shear_solver(&sim, MAX_ITER_1D_FD_SIMPLE_SHEAR, RTOL)) {
268 free_arrays(&sim);
269 exit(10);
270 }
271 #ifdef BENCHMARK_PERFORMANCE
272 t_end = clock();
273 t_elapsed = (double) (t_end - t_begin) / CLOCKS_PER_SEC;
274 printf("time spent per time step = %g s\n", t_elapsed);
275 #endif
276
277 if ((filetimeclock >= sim.file_dt || iter == 1) &&
278 strncmp(sim.name, DEFAULT_SIMULATION_NAME,
279 sizeof(DEFAULT_SIMULATION_NAME)) != 0) {
280 write_output_file(&sim, normalize);
281 filetimeclock = 0.0;
282 }
283 sim.t += sim.dt;
284 filetimeclock += sim.dt;
285 iter++;
286
287 } while (sim.t - sim.dt < sim.t_end);
288
289 print_output(&sim, stdout, normalize);
290
291 free_arrays(&sim);
292 return 0;
293 }