tcngf-pf.c - cngf-pf - continuum model for granular flows with pore-pressure dynamics (renamed from 1d_fd_simple_shear)
HTML git clone git://src.adamsgaard.dk/cngf-pf
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
tcngf-pf.c (6450B)
---
1 #include <stdlib.h>
2 #include <math.h>
3 #include <string.h>
4 #include <time.h>
5 #include <unistd.h>
6 #include <err.h>
7
8 #include "simulation.h"
9 #include "fluid.h"
10
11 #include "arg.h"
12
13 /* uncomment to print time spent per time step to stdout */
14 /* #define BENCHMARK_PERFORMANCE */
15
16 char *argv0;
17
18 static void
19 usage(void)
20 {
21 fprintf(stderr, "usage: %s "
22 "[-A grain-nonlocal-ampl] "
23 "[-a fluid-pressure-ampl] "
24 "[-b grain-rate-dependence] "
25 "[-C fluid-compressibility] "
26 "[-c grain-cohesion] "
27 "[-d grain-size] "
28 "[-D fluid-diffusivity] "
29 "[-e end-time] "
30 "[-F] "
31 "[-f applied-shear-friction] "
32 "[-g gravity-accel] "
33 "[-H fluid-pressure-phase] "
34 "[-h] "
35 "[-I file-interval] "
36 "[-i fluid-viscosity] "
37 "[-j time-step] "
38 "[-K dilatancy-constant] "
39 "[-k fluid-permeability] "
40 "[-L length] "
41 "[-l applied-shear-vel-limit] "
42 "[-m grain-friction] "
43 "[-N] "
44 "[-n normal-stress] "
45 "[-O fluid-pressure-top] "
46 "[-o origo] "
47 "[-P grain-compressibility] "
48 "[-p grain-porosity] "
49 "[-q fluid-pressure-freq] "
50 "[-R fluid-density] "
51 "[-r grain-density] "
52 "[-S fluid-pressure-pulse-shape] "
53 "[-s applied-shear-vel] "
54 "[-T] "
55 "[-t curr-time] "
56 "[-U resolution] "
57 "[-u fluid-pulse-time] "
58 "[-v] "
59 "[-Y max-porosity] "
60 "[-y min-porosity] "
61 "[-X relative-tolerance] "
62 "[-x max-iter] "
63 "[name]\n", argv0);
64 exit(1);
65 }
66
67 int
68 main(int argc, char *argv[])
69 {
70 int i, normalize = 0, dt_override = 0, ret, iter, max_iter = 100000;
71 double new_phi, new_k, filetimeclock;
72 struct simulation sim;
73 double rtol = 1e-5;
74
75 #ifdef BENCHMARK_PERFORMANCE
76 clock_t t_begin, t_end;
77 double t_elapsed;
78 #endif
79
80 #ifdef __OpenBSD__
81 if (pledge("stdio wpath cpath", NULL) == -1)
82 err(2, "pledge failed");
83 #endif
84
85 init_sim(&sim);
86 new_phi = sim.phi[0];
87 new_k = sim.k[0];
88
89 ARGBEGIN {
90 case 'A':
91 sim.A = atof(EARGF(usage()));
92 break;
93 case 'a':
94 sim.p_f_mod_ampl = atof(EARGF(usage()));
95 break;
96 case 'b':
97 sim.b = atof(EARGF(usage()));
98 break;
99 case 'C':
100 sim.beta_f = atof(EARGF(usage()));
101 break;
102 case 'c':
103 sim.C = atof(EARGF(usage()));
104 break;
105 case 'd':
106 sim.d = atof(EARGF(usage()));
107 break;
108 case 'D':
109 sim.D = atof(EARGF(usage()));
110 break;
111 case 'e':
112 sim.t_end = atof(EARGF(usage()));
113 break;
114 case 'F':
115 sim.fluid = 1;
116 break;
117 case 'f':
118 sim.mu_wall = atof(EARGF(usage()));
119 break;
120 case 'g':
121 sim.G = atof(EARGF(usage()));
122 break;
123 case 'H':
124 sim.p_f_mod_phase = atof(EARGF(usage()));
125 break;
126 case 'h':
127 usage();
128 break;
129 case 'I':
130 sim.file_dt = atof(EARGF(usage()));
131 break;
132 case 'i':
133 sim.mu_f = atof(EARGF(usage()));
134 break;
135 case 'j':
136 sim.dt = atof(EARGF(usage()));
137 dt_override = 1;
138 break;
139 case 'K':
140 sim.dilatancy_constant = 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 case 'X':
224 rtol = atoi(EARGF(usage()));
225 break;
226 case 'x':
227 max_iter = atoi(EARGF(usage()));
228 break;
229 default:
230 usage();
231 } ARGEND;
232
233 if (argc == 1 && argv[0]) {
234 ret = snprintf(sim.name, sizeof(sim.name), "%s", argv[0]);
235 if (ret < 0 || (size_t)ret >= sizeof(sim.name))
236 errx(1, "%s: could not write sim.name", __func__);
237 } else if (argc > 1)
238 usage();
239
240 if (sim.nz < 1)
241 sim.nz = (int) ceil(sim.L_z / sim.d);
242
243 prepare_arrays(&sim);
244
245 if (!isnan(new_phi))
246 for (i = 0; i < sim.nz; ++i)
247 sim.phi[i] = new_phi;
248 if (!isnan(new_k))
249 for (i = 0; i < sim.nz; ++i)
250 sim.k[i] = new_k;
251
252 lithostatic_pressure_distribution(&sim);
253
254 if (sim.fluid) {
255 hydrostatic_fluid_pressure_distribution(&sim);
256 if (!dt_override) {
257 if (set_largest_fluid_timestep(&sim, 0.5)) {
258 free_arrays(&sim);
259 return 20;
260 }
261 if (sim.transient)
262 set_coupled_fluid_transient_timestep(&sim, 0.5);
263 }
264 }
265 #ifdef DEBUG
266 fprintf(stderr, "t_val = %g\n", sim.dt);
267 #endif
268
269 if (sim.dt > sim.file_dt)
270 sim.dt = sim.file_dt;
271 if (sim.dt > sim.t_end)
272 sim.dt = sim.t_end;
273 compute_effective_stress(&sim);
274
275 check_simulation_parameters(&sim);
276
277 filetimeclock = 0.0;
278 iter = 0;
279 do {
280
281 #ifdef BENCHMARK_PERFORMANCE
282 t_begin = clock();
283 #endif
284
285 if (coupled_shear_solver(&sim, max_iter, rtol)) {
286 free_arrays(&sim);
287 exit(10);
288 }
289 #ifdef BENCHMARK_PERFORMANCE
290 t_end = clock();
291 t_elapsed = (double) (t_end - t_begin) / CLOCKS_PER_SEC;
292 printf("time spent per time step = %g s\n", t_elapsed);
293 #endif
294
295 if ((filetimeclock >= sim.file_dt || iter == 1) &&
296 strncmp(sim.name, DEFAULT_SIMULATION_NAME,
297 sizeof(DEFAULT_SIMULATION_NAME)) != 0) {
298 write_output_file(&sim, normalize);
299 filetimeclock = 0.0;
300 }
301 filetimeclock += sim.dt;
302 iter++;
303
304 } while (sim.t < sim.t_end);
305
306 print_output(&sim, stdout, normalize);
307
308 free_arrays(&sim);
309 return 0;
310 }