cngf-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 }