t1d_fd_simple_shear.c - cngf-pf - continuum model for granular flows with pore-… | |
git clone git://src.adamsgaard.dk/1d_fd_simple_shear | |
Log | |
Files | |
Refs | |
README | |
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'… | |
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_SHE… | |
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 } |