main.c (13800B)
1 #include <unistd.h> 2 #include <stdio.h> 3 #include <stdlib.h> 4 #include <math.h> 5 #include <getopt.h> 6 #include <string.h> 7 #include <time.h> 8 9 #include "simulation.h" 10 #include "fluid.h" 11 12 #define VERSION "0.4.6" 13 #define PROGNAME "1d_fd_simple_shear" 14 15 #include "parameter_defaults.h" 16 17 /* relative tolerance criteria for granular fluidity solver */ 18 #define RTOL_GRANULAR 1e-5 19 #define MAX_ITER_GRANULAR 10000 20 21 /* relative tolerance criteria for fluid-pressure solver */ 22 #define RTOL_DARCY 1e-5 23 #define MAX_ITER_DARCY 10000 24 25 /* relative tolerance criteria when shear velocity is restricted */ 26 #define RTOL_STRESS 1e-3 27 #define MAX_ITER_STRESS 20000 28 29 /* uncomment to print time spent per time step to stdout */ 30 /*#define BENCHMARK_PERFORMANCE*/ 31 32 static void 33 usage(void) 34 { 35 struct simulation sim = init_sim(); 36 printf("%s: %s [OPTIONS] [NAME]\n" 37 "runs a simulation and outputs state in files prefixed with NAME.\n" 38 "If NAME is not specified, intermediate output files are not written.\n" 39 "\nOptional arguments:\n" 40 " -v, --version show version information\n" 41 " -h, --help show this message\n" 42 " -N, --normalize normalize output velocity\n" 43 " -G, --gravity VAL gravity magnitude [m/s^2] (default %g)\n" 44 " -P, --normal-stress VAL normal stress on top [Pa] (default %g)\n" 45 " -m, --stress-ratio VAL applied stress ratio [-] (default %g)\n" 46 " -s, --set-shear-velocity VAL set top shear velocity to this value [m/s]\n" 47 " (default %g), overrides --stress-ratio\n" 48 " -l, --limit-shear-velocity VAL limit top shear velocity to this value [m/s]\n" 49 " (default %g), overrides --stress-ratio and\n" 50 " --set-shear-velocity\n" 51 " -V, --velocity-bottom VAL base velocity at bottom [m/s] (default %g)\n" 52 " -A, --nonlocal-amplitude VAL amplitude of nonlocality [-] (default %g)\n" 53 " -b, --rate-dependence VAL rate dependence beyond yield [-] (default %g)\n" 54 " -f, --friction-coefficient VAL grain friction coefficient [-] (default %g)\n" 55 " -C, --cohesion VAL grain cohesion [Pa] (default %g)\n" 56 " -p, --porosity VAL porosity fraction [-] (default %g)\n" 57 " -d, --grain-size VAL representative grain size [m] (default %g)\n" 58 " -r, --density VAL grain material density [kg/m^3] (default %g)\n" 59 " -n, --resolution VAL number of cells in domain [-]\n" 60 " (default cell size equals grain size)\n" 61 " -o, --origo VAL coordinate system origo [m] (default %g)\n" 62 " -L, --length VAL domain length [m] (default %g)\n" 63 "\nOptional arguments only relevant with transient (fluid) simulation:\n" 64 " -F, --fluid enable pore fluid computations\n" 65 " -c, --fluid-compressibility VAL fluid compressibility [Pa^-1] (default %g)\n" 66 " -i, --fluid-viscosity VAL fluid viscosity [Pa*s] (default %g)\n" 67 " -R, --fluid-density VAL fluid density [kg/m^3] (default %g)\n" 68 " -k, --fluid-permeability VAL fluid intrinsic permeability [m^2] (default %g)\n" 69 " -O, --fluid-pressure-top VAL fluid pressure at +z edge [Pa] (default %g)\n" 70 " -a, --fluid-pressure-ampl VAL amplitude of pressure variations [Pa] (default %g)\n" 71 " -q, --fluid-pressure-freq VAL frequency of sinusoidal pressure variations [s^-1]\n" 72 " (default %g)\n" 73 " -H, --fluid-pressure-phase VAL fluid pressure at +z edge [Pa] (default %g)\n" 74 " -u, --fluid-pressure-pulse-time VAL fluid pressure pulse peak time [s]\n" 75 " (default %g)\n" 76 " -S, --fluid-pressure-pulse-shape VAL\n" 77 " where VAL is triangular (default) or square\n" 78 " -t, --time VAL simulation start time [s] (default %g)\n" 79 " -T, --time-end VAL simulation end time [s] (default %g)\n" 80 " -I, --file-interval VAL interval between output files [s] (default %g)\n" 81 "\n" 82 "The output (stdout or output files) consists of the following " 83 "tab-delimited\nfields, with one row per cell:\n" 84 " 1. z_position [m]\n" 85 " 2. x_velocity [m/s]\n" 86 " 3. normal_stress [Pa]\n" 87 " 4. friction [-]\n" 88 " 5. strain_rate [1/s]\n\n" 89 "With --fluid enabled:\n" 90 " 1. z_position [m]\n" 91 " 2. x_velocity [m/s]\n" 92 " 3. effective_normal_stress [Pa]\n" 93 " 4. fluid_pressure [Pa]\n" 94 " 5. friction [-]\n" 95 " 6. strain_rate [1/s]\n" 96 , 97 __func__, PROGNAME, 98 sim.G, 99 sim.P_wall, 100 sim.mu_wall, 101 sim.v_x_fix, 102 sim.v_x_limit, 103 sim.v_x_bot, 104 sim.A, 105 sim.b, 106 sim.mu_s, 107 sim.C, 108 sim.phi[0], 109 sim.d, 110 sim.rho_s, 111 sim.origo_z, 112 sim.L_z, 113 sim.beta_f, 114 sim.mu_f, 115 sim.rho_f, 116 sim.k[0], 117 sim.p_f_top, 118 sim.p_f_mod_ampl, 119 sim.p_f_mod_freq, 120 sim.p_f_mod_phase, 121 sim.p_f_mod_pulse_time, 122 sim.t, 123 sim.t_end, 124 sim.file_dt); 125 free(sim.phi); 126 free(sim.k); 127 } 128 129 static void 130 version(void) 131 { 132 printf("%s v%s\n" 133 "Licensed under the GNU Public License, v3+\n" 134 "written by Anders Damsgaard, anders@adamsgaard.dk\n" 135 "https://gitlab.com/admesg/1d_fd_simple_shear\n" 136 , PROGNAME, VERSION); 137 } 138 139 int 140 main(int argc, char* argv[]) 141 { 142 int norm, opt, i; 143 struct simulation sim; 144 const char* optstring; 145 unsigned long iter, stressiter; 146 double new_phi, new_k, filetimeclock, res_norm, mu_wall_orig; 147 #ifdef BENCHMARK_PERFORMANCE 148 clock_t t_begin, t_end; 149 double t_elapsed; 150 #endif 151 152 #ifdef __OpenBSD__ 153 if (pledge("stdio wpath cpath", NULL) == -1) { 154 fprintf(stderr, "error: pledge failed"); 155 exit(1); 156 } 157 #endif 158 159 /* load with default values */ 160 sim = init_sim(); 161 162 norm = 0; 163 164 optstring = "hvNn:G:P:m:s:l:V:A:b:f:C:Fp:d:r:o:L:c:i:R:k:O:a:q:H:T:S:t:T:D:I:"; 165 const struct option longopts[] = { 166 {"help", no_argument, NULL, 'h'}, 167 {"version", no_argument, NULL, 'v'}, 168 {"normalize", no_argument, NULL, 'N'}, 169 {"gravity", required_argument, NULL, 'G'}, 170 {"normal-stress", required_argument, NULL, 'P'}, 171 {"stress-ratio", required_argument, NULL, 'm'}, 172 {"set-shear-velocity", required_argument, NULL, 's'}, 173 {"limit-shear-velocity", required_argument, NULL, 'l'}, 174 {"velocity-bottom", required_argument, NULL, 'V'}, 175 {"nonlocal-amplitude", required_argument, NULL, 'A'}, 176 {"rate-dependence", required_argument, NULL, 'b'}, 177 {"friction-coefficient", required_argument, NULL, 'f'}, 178 {"cohesion", required_argument, NULL, 'C'}, 179 {"porosity", required_argument, NULL, 'p'}, 180 {"grain-size", required_argument, NULL, 'd'}, 181 {"density", required_argument, NULL, 'r'}, 182 {"resolution", required_argument, NULL, 'n'}, 183 {"origo", required_argument, NULL, 'o'}, 184 {"length", required_argument, NULL, 'L'}, 185 {"fluid", no_argument, NULL, 'F'}, 186 {"fluid-compressiblity", required_argument, NULL, 'c'}, 187 {"fluid-viscosity", required_argument, NULL, 'i'}, 188 {"fluid-density", required_argument, NULL, 'R'}, 189 {"fluid-permeability", required_argument, NULL, 'k'}, 190 {"fluid-pressure-top", required_argument, NULL, 'O'}, 191 {"fluid-pressure-ampl", required_argument, NULL, 'a'}, 192 {"fluid-pressure-freq", required_argument, NULL, 'q'}, 193 {"fluid-pressure-phase", required_argument, NULL, 'H'}, 194 {"fluid-pressure-pulse-time", required_argument, NULL, 'u'}, 195 {"fluid-pressure-pulse-shape",required_argument, NULL, 'S'}, 196 {"time", required_argument, NULL, 't'}, 197 {"time-end", required_argument, NULL, 'T'}, 198 {"file-interval", required_argument, NULL, 'I'}, 199 {NULL, 0, NULL, 0} 200 }; 201 202 new_phi = sim.phi[0]; 203 new_k = sim.k[0]; 204 while ((opt = getopt_long(argc, argv, optstring, longopts, NULL)) != -1) { 205 switch (opt) { 206 case -1: /* no more arguments */ 207 case 0: /* long options toggles */ 208 break; 209 210 case 'h': 211 usage(); 212 free(sim.phi); 213 free(sim.k); 214 return 0; 215 case 'v': 216 version(); 217 free(sim.phi); 218 free(sim.k); 219 return 0; 220 case 'n': 221 sim.nz = atoi(optarg); 222 break; 223 case 'N': 224 norm = 1; 225 break; 226 case 'G': 227 sim.G = atof(optarg); 228 break; 229 case 'P': 230 sim.P_wall = atof(optarg); 231 break; 232 case 'm': 233 sim.mu_wall = atof(optarg); 234 break; 235 case 's': 236 sim.v_x_fix = atof(optarg); 237 break; 238 case 'l': 239 sim.v_x_limit = atof(optarg); 240 break; 241 case 'V': 242 sim.v_x_bot = atof(optarg); 243 break; 244 case 'A': 245 sim.A = atof(optarg); 246 break; 247 case 'b': 248 sim.b = atof(optarg); 249 break; 250 case 'f': 251 sim.mu_s = atof(optarg); 252 break; 253 case 'C': 254 sim.C = atof(optarg); 255 break; 256 case 'p': 257 new_phi = atof(optarg); 258 break; 259 case 'd': 260 sim.d = atof(optarg); 261 break; 262 case 'r': 263 sim.rho_s = atof(optarg); 264 break; 265 case 'o': 266 sim.origo_z = atof(optarg); 267 break; 268 case 'L': 269 sim.L_z = atof(optarg); 270 break; 271 case 'F': 272 sim.fluid = 1; 273 break; 274 case 'c': 275 sim.beta_f = atof(optarg); 276 break; 277 case 'i': 278 sim.mu_f = atof(optarg); 279 break; 280 case 'R': 281 sim.rho_f = atof(optarg); 282 break; 283 case 'k': 284 new_k = atof(optarg); 285 break; 286 case 'O': 287 sim.p_f_top = atof(optarg); 288 break; 289 case 'a': 290 sim.p_f_mod_ampl = atof(optarg); 291 break; 292 case 'q': 293 sim.p_f_mod_freq = atof(optarg); 294 break; 295 case 'H': 296 sim.p_f_mod_phase = atof(optarg); 297 break; 298 case 'u': 299 sim.p_f_mod_pulse_time = atof(optarg); 300 break; 301 case 'S': 302 if (strcmp(optarg, "triangle") == 0) 303 sim.p_f_mod_pulse_shape = 0; 304 else if (strcmp(optarg, "square") == 0) 305 sim.p_f_mod_pulse_shape = 1; 306 else { 307 fprintf(stderr, "error: invalid pulse shape '%s'\n", 308 optarg); 309 return 1; 310 } 311 break; 312 case 't': 313 sim.t = atof(optarg); 314 break; 315 case 'T': 316 sim.t_end = atof(optarg); 317 break; 318 case 'I': 319 sim.file_dt = atof(optarg); 320 break; 321 default: 322 fprintf(stderr, "%s: invalid option -- %c\n", argv[0], opt); 323 fprintf(stderr, "Try `%s --help` for more information\n", 324 argv[0]); 325 return -2; 326 } 327 } 328 for (i=optind; i<argc; ++i) { 329 if (i>optind) { 330 fprintf(stderr, 331 "error: more than one simulation name specified\n"); 332 return 1; 333 } 334 snprintf(sim.name, sizeof(sim.name), "%s", argv[i]); 335 } 336 337 if (sim.nz < 1) 338 sim.nz = (int)ceil(sim.L_z/sim.d); 339 340 prepare_arrays(&sim); 341 342 if (!isnan(new_phi)) 343 for (i=0; i<sim.nz; ++i) 344 sim.phi[i] = new_phi; 345 if (!isnan(new_k)) 346 for (i=0; i<sim.nz; ++i) 347 sim.k[i] = new_k; 348 349 lithostatic_pressure_distribution(&sim); 350 351 if (sim.fluid) { 352 hydrostatic_fluid_pressure_distribution(&sim); 353 set_largest_fluid_timestep(&sim, 0.5); 354 } 355 356 check_simulation_parameters(&sim); 357 358 filetimeclock = 0.0; 359 iter = 0; 360 mu_wall_orig = sim.mu_wall; 361 res_norm = NAN; 362 while (sim.t <= sim.t_end) { 363 364 #ifdef BENCHMARK_PERFORMANCE 365 t_begin = clock(); 366 #endif 367 368 stressiter = 0; 369 do { 370 if (sim.fluid) { 371 if (darcy_solver_1d(&sim, MAX_ITER_DARCY, RTOL_DARCY)) 372 exit(1); 373 } 374 375 compute_effective_stress(&sim); 376 compute_friction(&sim); 377 compute_cooperativity_length(&sim); 378 379 if (implicit_1d_jacobian_poisson_solver(&sim, MAX_ITER_GRANULAR, 380 RTOL_GRANULAR)) 381 exit(1); 382 383 compute_shear_strain_rate_plastic(&sim); 384 compute_shear_velocity(&sim); 385 386 if (!isnan(sim.v_x_limit) || !isnan(sim.v_x_fix)) { 387 if (!isnan(sim.v_x_limit)) { 388 res_norm = (sim.v_x_limit - sim.v_x[sim.nz-1]) 389 /(sim.v_x[sim.nz-1] + 1e-12); 390 if (res_norm > 0.0) 391 res_norm = 0.0; 392 } else { 393 res_norm = (sim.v_x_fix - sim.v_x[sim.nz-1]) 394 /(sim.v_x[sim.nz-1] + 1e-12); 395 } 396 sim.mu_wall *= 1.0 + (res_norm*1e-2); 397 } 398 399 if (++stressiter > MAX_ITER_STRESS) { 400 fprintf(stderr, "error: stress solution did not converge:\n"); 401 fprintf(stderr, 402 "v_x=%g, v_x_fix=%g, v_x_limit=%g, " 403 "res_norm=%g, mu_wall=%g\n", 404 sim.v_x[sim.nz-1], sim.v_x_fix, sim.v_x_limit, 405 res_norm, sim.mu_wall); 406 free_arrays(&sim); 407 return 10; 408 } 409 410 } while ((!isnan(sim.v_x_fix) || !isnan(sim.v_x_limit)) 411 && fabs(res_norm) > RTOL_STRESS); 412 413 #ifdef BENCHMARK_PERFORMANCE 414 t_end = clock(); 415 t_elapsed = (double)(t_end - t_begin)/CLOCKS_PER_SEC; 416 printf("time spent per time step = %g s\n", t_elapsed); 417 #endif 418 419 if (!isnan(sim.v_x_limit)) 420 sim.mu_wall = mu_wall_orig; 421 sim.t += sim.dt; 422 filetimeclock += sim.dt; 423 iter++; 424 425 if ((filetimeclock - sim.dt >= sim.file_dt || iter == 1) && 426 strcmp(sim.name, DEFAULT_SIMULATION_NAME) != 0) { 427 write_output_file(&sim, norm); 428 filetimeclock = 0.0; 429 if (iter == 1) 430 sim.t = 0.0; 431 } 432 } 433 434 if (sim.fluid) 435 print_wet_output(stdout, &sim, norm); 436 else 437 print_dry_output(stdout, &sim, norm); 438 439 free_arrays(&sim); 440 return 0; 441 }