cngf-pf.c (6606B)
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 if (prepare_arrays(&sim)) 251 errx(1, "prepare_arrays failed"); 252 253 if (!isnan(new_phi)) 254 for (i = 0; i < sim.nz; ++i) 255 sim.phi[i] = new_phi; 256 if (!isnan(new_k)) 257 for (i = 0; i < sim.nz; ++i) 258 sim.k[i] = new_k; 259 260 lithostatic_pressure_distribution(&sim); 261 262 if (sim.fluid) { 263 hydrostatic_fluid_pressure_distribution(&sim); 264 if (!dt_override) { 265 if (set_largest_fluid_timestep(&sim, 0.5)) { 266 free_arrays(&sim); 267 return 20; 268 } 269 if (sim.transient) 270 set_coupled_fluid_transient_timestep(&sim, 0.5); 271 } 272 } 273 #ifdef DEBUG 274 fprintf(stderr, "t_val = %g\n", sim.dt); 275 #endif 276 277 if (sim.dt > sim.file_dt) 278 sim.dt = sim.file_dt; 279 if (sim.dt > sim.t_end) 280 sim.dt = sim.t_end; 281 compute_effective_stress(&sim); 282 283 if (check_simulation_parameters(&sim)) { 284 free_arrays(&sim); 285 exit(1); 286 } 287 288 filetimeclock = 0.0; 289 iter = 0; 290 t_begin = clock(); 291 do { 292 if ((ret = coupled_shear_solver(&sim, max_iter, rtol))) { 293 free_arrays(&sim); 294 exit(ret); 295 } 296 297 if ((filetimeclock >= sim.file_dt || iter == 1) && 298 strncmp(sim.name, DEFAULT_SIMULATION_NAME, 299 sizeof(DEFAULT_SIMULATION_NAME)) != 0) { 300 write_output_file(&sim, normalize); 301 filetimeclock = 0.0; 302 } 303 filetimeclock += sim.dt; 304 iter++; 305 306 } while (sim.t < sim.t_end); 307 t_end = clock(); 308 t_elapsed = (double) (t_end - t_begin) / CLOCKS_PER_SEC; 309 310 if (!benchmark) 311 print_output(&sim, stdout, normalize); 312 else { 313 fprintf(stderr, "benchmark: elapsed_time=%.6f\n", t_elapsed); 314 print_solver_stats(stderr); 315 } 316 317 free_arrays(&sim); 318 return 0; 319 }