commit 84ff27ebe390512189687c4fc28adddef22a3d0a
parent 812fb3de573d8e4412dbf66aafd9f6b54828f635
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Wed, 11 Dec 2019 11:25:27 +0100
Clean up new main() and fix version info
Diffstat:
2 files changed, 21 insertions(+), 234 deletions(-)
diff --git a/1d_fd_simple_shear.c b/1d_fd_simple_shear.c
@@ -129,9 +129,9 @@ static void
version(void)
{
printf("%s v%s\n"
- "Licensed under the GNU Public License, v3+\n"
+ "Licensed under the ISC License, see LICENSE for details.\n"
"written by Anders Damsgaard, anders@adamsgaard.dk\n"
- "https://gitlab.com/admesg/1d_fd_simple_shear\n"
+ "https://src.adamsgaard.dk/1d_fd_simple_shear\n"
, PROGNAME, VERSION);
}
diff --git a/max_depth_simple_shear.c b/max_depth_simple_shear.c
@@ -9,21 +9,13 @@
#include "simulation.h"
#include "fluid.h"
-#define PROGNAME "1d_fd_simple_shear"
+#define PROGNAME "max_depth_simple_shear"
#include "parameter_defaults.h"
-/* relative tolerance criteria for granular fluidity solver */
-#define RTOL_GRANULAR 1e-5
-#define MAX_ITER_GRANULAR 10000
-
-/* relative tolerance criteria for fluid-pressure solver */
-#define RTOL_DARCY 1e-5
-#define MAX_ITER_DARCY 10000
-
-/* relative tolerance criteria when shear velocity is restricted */
-#define RTOL_STRESS 1e-3
-#define MAX_ITER_STRESS 20000
+/* relative tolerance criteria for solver */
+#define RTOL 1e-5
+#define MAX_ITER 10000
/* uncomment to print time spent per time step to stdout */
/*#define BENCHMARK_PERFORMANCE*/
@@ -32,35 +24,16 @@ static void
usage(void)
{
struct simulation sim = init_sim();
- printf("%s: %s [OPTIONS] [NAME]\n"
- "runs a simulation and outputs state in files prefixed with NAME.\n"
- "If NAME is not specified, intermediate output files are not written.\n"
+ printf("%s: %s [OPTIONS]\n"
+ "outputs the maximum deformation depth in a shear system with sinusoidal\n"
+ "fluid-pressure perturbations from the top, assuming infinite layer thickness.\n"
+ "Assumes infinite layer thickness.\n
"\nOptional arguments:\n"
" -v, --version show version information\n"
" -h, --help show this message\n"
- " -N, --normalize normalize output velocity\n"
" -G, --gravity VAL gravity magnitude [m/s^2] (default %g)\n"
" -P, --normal-stress VAL normal stress on top [Pa] (default %g)\n"
- " -m, --stress-ratio VAL applied stress ratio [-] (default %g)\n"
- " -s, --set-shear-velocity VAL set top shear velocity to this value [m/s]\n"
- " (default %g), overrides --stress-ratio\n"
- " -l, --limit-shear-velocity VAL limit top shear velocity to this value [m/s]\n"
- " (default %g), overrides --stress-ratio and\n"
- " --set-shear-velocity\n"
- " -V, --velocity-bottom VAL base velocity at bottom [m/s] (default %g)\n"
- " -A, --nonlocal-amplitude VAL amplitude of nonlocality [-] (default %g)\n"
- " -b, --rate-dependence VAL rate dependence beyond yield [-] (default %g)\n"
- " -f, --friction-coefficient VAL grain friction coefficient [-] (default %g)\n"
- " -C, --cohesion VAL grain cohesion [Pa] (default %g)\n"
- " -p, --porosity VAL porosity fraction [-] (default %g)\n"
- " -d, --grain-size VAL representative grain size [m] (default %g)\n"
" -r, --density VAL grain material density [kg/m^3] (default %g)\n"
- " -n, --resolution VAL number of cells in domain [-]\n"
- " (default cell size equals grain size)\n"
- " -o, --origo VAL coordinate system origo [m] (default %g)\n"
- " -L, --length VAL domain length [m] (default %g)\n"
- "\nOptional arguments only relevant with transient (fluid) simulation:\n"
- " -F, --fluid enable pore fluid computations\n"
" -c, --fluid-compressibility VAL fluid compressibility [Pa^-1] (default %g)\n"
" -i, --fluid-viscosity VAL fluid viscosity [Pa*s] (default %g)\n"
" -R, --fluid-density VAL fluid density [kg/m^3] (default %g)\n"
@@ -69,58 +42,20 @@ usage(void)
" -a, --fluid-pressure-ampl VAL amplitude of pressure variations [Pa] (default %g)\n"
" -q, --fluid-pressure-freq VAL frequency of sinusoidal pressure variations [s^-1]\n"
" (default %g)\n"
- " -H, --fluid-pressure-phase VAL fluid pressure at +z edge [Pa] (default %g)\n"
- " -u, --fluid-pressure-pulse-time VAL fluid pressure pulse peak time [s]\n"
- " (default %g)\n"
- " -S, --fluid-pressure-pulse-shape VAL\n"
- " where VAL is triangular (default) or square\n"
- " -t, --time VAL simulation start time [s] (default %g)\n"
- " -T, --time-end VAL simulation end time [s] (default %g)\n"
- " -I, --file-interval VAL interval between output files [s] (default %g)\n"
"\n"
- "The output (stdout or output files) consists of the following "
- "tab-delimited\nfields, with one row per cell:\n"
- " 1. z_position [m]\n"
- " 2. x_velocity [m/s]\n"
- " 3. normal_stress [Pa]\n"
- " 4. friction [-]\n"
- " 5. strain_rate [1/s]\n\n"
- "With --fluid enabled:\n"
- " 1. z_position [m]\n"
- " 2. x_velocity [m/s]\n"
- " 3. effective_normal_stress [Pa]\n"
- " 4. fluid_pressure [Pa]\n"
- " 5. friction [-]\n"
- " 6. strain_rate [1/s]\n"
+ "The output value is the max. deformation depth from the top in meters.\n"
,
__func__, PROGNAME,
sim.G,
sim.P_wall,
- sim.mu_wall,
- sim.v_x_fix,
- sim.v_x_limit,
- sim.v_x_bot,
- sim.A,
- sim.b,
- sim.mu_s,
- sim.C,
- sim.phi[0],
- sim.d,
sim.rho_s,
- sim.origo_z,
- sim.L_z,
sim.beta_f,
sim.mu_f,
sim.rho_f,
sim.k[0],
sim.p_f_top,
sim.p_f_mod_ampl,
- sim.p_f_mod_freq,
- sim.p_f_mod_phase,
- sim.p_f_mod_pulse_time,
- sim.t,
- sim.t_end,
- sim.file_dt);
+ sim.p_f_mod_freq);
free(sim.phi);
free(sim.k);
}
@@ -129,9 +64,9 @@ static void
version(void)
{
printf("%s v%s\n"
- "Licensed under the GNU Public License, v3+\n"
+ "Licensed under the ISC License, see LICENSE for details.\n"
"written by Anders Damsgaard, anders@adamsgaard.dk\n"
- "https://gitlab.com/admesg/1d_fd_simple_shear\n"
+ "https://src.adamsgaard.dk/1d_fd_simple_shear\n"
, PROGNAME, VERSION);
}
@@ -141,15 +76,14 @@ main(int argc, char* argv[])
int norm, opt, i;
struct simulation sim;
const char* optstring;
- unsigned long iter, stressiter;
- double new_phi, new_k, filetimeclock, res_norm, mu_wall_orig;
+ unsigned long iter;
#ifdef BENCHMARK_PERFORMANCE
clock_t t_begin, t_end;
double t_elapsed;
#endif
#ifdef __OpenBSD__
- if (pledge("stdio wpath cpath", NULL) == -1) {
+ if (pledge("stdio", NULL) == -1) {
fprintf(stderr, "error: pledge failed");
exit(1);
}
@@ -164,24 +98,9 @@ main(int argc, char* argv[])
const struct option longopts[] = {
{"help", no_argument, NULL, 'h'},
{"version", no_argument, NULL, 'v'},
- {"normalize", no_argument, NULL, 'N'},
{"gravity", required_argument, NULL, 'G'},
{"normal-stress", required_argument, NULL, 'P'},
- {"stress-ratio", required_argument, NULL, 'm'},
- {"set-shear-velocity", required_argument, NULL, 's'},
- {"limit-shear-velocity", required_argument, NULL, 'l'},
- {"velocity-bottom", required_argument, NULL, 'V'},
- {"nonlocal-amplitude", required_argument, NULL, 'A'},
- {"rate-dependence", required_argument, NULL, 'b'},
- {"friction-coefficient", required_argument, NULL, 'f'},
- {"cohesion", required_argument, NULL, 'C'},
- {"porosity", required_argument, NULL, 'p'},
- {"grain-size", required_argument, NULL, 'd'},
{"density", required_argument, NULL, 'r'},
- {"resolution", required_argument, NULL, 'n'},
- {"origo", required_argument, NULL, 'o'},
- {"length", required_argument, NULL, 'L'},
- {"fluid", no_argument, NULL, 'F'},
{"fluid-compressiblity", required_argument, NULL, 'c'},
{"fluid-viscosity", required_argument, NULL, 'i'},
{"fluid-density", required_argument, NULL, 'R'},
@@ -189,12 +108,6 @@ main(int argc, char* argv[])
{"fluid-pressure-top", required_argument, NULL, 'O'},
{"fluid-pressure-ampl", required_argument, NULL, 'a'},
{"fluid-pressure-freq", required_argument, NULL, 'q'},
- {"fluid-pressure-phase", required_argument, NULL, 'H'},
- {"fluid-pressure-pulse-time", required_argument, NULL, 'u'},
- {"fluid-pressure-pulse-shape",required_argument, NULL, 'S'},
- {"time", required_argument, NULL, 't'},
- {"time-end", required_argument, NULL, 'T'},
- {"file-interval", required_argument, NULL, 'I'},
{NULL, 0, NULL, 0}
};
@@ -216,60 +129,15 @@ main(int argc, char* argv[])
free(sim.phi);
free(sim.k);
return 0;
- case 'n':
- sim.nz = atoi(optarg);
- break;
- case 'N':
- norm = 1;
- break;
case 'G':
sim.G = atof(optarg);
break;
case 'P':
sim.P_wall = atof(optarg);
break;
- case 'm':
- sim.mu_wall = atof(optarg);
- break;
- case 's':
- sim.v_x_fix = atof(optarg);
- break;
- case 'l':
- sim.v_x_limit = atof(optarg);
- break;
- case 'V':
- sim.v_x_bot = atof(optarg);
- break;
- case 'A':
- sim.A = atof(optarg);
- break;
- case 'b':
- sim.b = atof(optarg);
- break;
- case 'f':
- sim.mu_s = atof(optarg);
- break;
- case 'C':
- sim.C = atof(optarg);
- break;
- case 'p':
- new_phi = atof(optarg);
- break;
- case 'd':
- sim.d = atof(optarg);
- break;
case 'r':
sim.rho_s = atof(optarg);
break;
- case 'o':
- sim.origo_z = atof(optarg);
- break;
- case 'L':
- sim.L_z = atof(optarg);
- break;
- case 'F':
- sim.fluid = 1;
- break;
case 'c':
sim.beta_f = atof(optarg);
break;
@@ -291,32 +159,6 @@ main(int argc, char* argv[])
case 'q':
sim.p_f_mod_freq = atof(optarg);
break;
- case 'H':
- sim.p_f_mod_phase = atof(optarg);
- break;
- case 'u':
- sim.p_f_mod_pulse_time = atof(optarg);
- break;
- case 'S':
- if (strcmp(optarg, "triangle") == 0)
- sim.p_f_mod_pulse_shape = 0;
- else if (strcmp(optarg, "square") == 0)
- sim.p_f_mod_pulse_shape = 1;
- else {
- fprintf(stderr, "error: invalid pulse shape '%s'\n",
- optarg);
- return 1;
- }
- break;
- case 't':
- sim.t = atof(optarg);
- break;
- case 'T':
- sim.t_end = atof(optarg);
- break;
- case 'I':
- sim.file_dt = atof(optarg);
- break;
default:
fprintf(stderr, "%s: invalid option -- %c\n", argv[0], opt);
fprintf(stderr, "Try `%s --help` for more information\n",
@@ -324,17 +166,11 @@ main(int argc, char* argv[])
return -2;
}
}
- for (i=optind; i<argc; ++i) {
- if (i>optind) {
- fprintf(stderr,
- "error: more than one simulation name specified\n");
- return 1;
- }
- snprintf(sim.name, sizeof(sim.name), "%s", argv[i]);
- }
- if (sim.nz < 1)
- sim.nz = (int)ceil(sim.L_z/sim.d);
+ if (optind>=argc) {
+ fprintf(stderr, "error: unknown argument specified\n");
+ return 1;
+ }
prepare_arrays(&sim);
@@ -354,60 +190,14 @@ main(int argc, char* argv[])
check_simulation_parameters(&sim);
- filetimeclock = 0.0;
- iter = 0;
- mu_wall_orig = sim.mu_wall;
- res_norm = NAN;
while (sim.t <= sim.t_end) {
#ifdef BENCHMARK_PERFORMANCE
t_begin = clock();
#endif
- stressiter = 0;
- do {
- if (sim.fluid) {
- if (darcy_solver_1d(&sim, MAX_ITER_DARCY, RTOL_DARCY))
- exit(1);
- }
-
- compute_effective_stress(&sim);
- compute_friction(&sim);
- compute_cooperativity_length(&sim);
-
- if (implicit_1d_jacobian_poisson_solver(&sim, MAX_ITER_GRANULAR,
- RTOL_GRANULAR))
- exit(1);
-
- compute_shear_strain_rate_plastic(&sim);
- compute_shear_velocity(&sim);
-
- if (!isnan(sim.v_x_limit) || !isnan(sim.v_x_fix)) {
- if (!isnan(sim.v_x_limit)) {
- res_norm = (sim.v_x_limit - sim.v_x[sim.nz-1])
- /(sim.v_x[sim.nz-1] + 1e-12);
- if (res_norm > 0.0)
- res_norm = 0.0;
- } else {
- res_norm = (sim.v_x_fix - sim.v_x[sim.nz-1])
- /(sim.v_x[sim.nz-1] + 1e-12);
- }
- sim.mu_wall *= 1.0 + (res_norm*1e-2);
- }
- if (++stressiter > MAX_ITER_STRESS) {
- fprintf(stderr, "error: stress solution did not converge:\n");
- fprintf(stderr,
- "v_x=%g, v_x_fix=%g, v_x_limit=%g, "
- "res_norm=%g, mu_wall=%g\n",
- sim.v_x[sim.nz-1], sim.v_x_fix, sim.v_x_limit,
- res_norm, sim.mu_wall);
- free_arrays(&sim);
- return 10;
- }
- } while ((!isnan(sim.v_x_fix) || !isnan(sim.v_x_limit))
- && fabs(res_norm) > RTOL_STRESS);
#ifdef BENCHMARK_PERFORMANCE
t_end = clock();
@@ -430,10 +220,7 @@ main(int argc, char* argv[])
}
}
- if (sim.fluid)
- print_wet_output(stdout, &sim, norm);
- else
- print_dry_output(stdout, &sim, norm);
+ printf("%s\n", max_depth);
free_arrays(&sim);
return 0;
