1d_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 }