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