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