commit acf33cd7db63d47c5a7ebbba9a568371ca3fe34f
parent dc46674d016883f15b65d5d0150429d1b142eebe
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Thu, 5 Mar 2020 16:33:41 +0100
Only allow short options, rework option parsing, global sim
Diffstat:
| M | 1d_fd_simple_shear.c | | | 474 | +++++++++++++++++++++++++++++-------------------------------------------------- |
| M | Makefile | | | 15 | +-------------- |
| A | arg.h | | | 38 | ++++++++++++++++++++++++++++++++++++++ |
| M | fluid.c | | | 151 | +++++++++++++++++++++++++++++++++++++++---------------------------------------- |
| M | fluid.h | | | 7 | ++++--- |
| M | max_depth_simple_shear.c | | | 213 | +++++++++++++++++++++++++++---------------------------------------------------- |
| M | parameter_defaults.h | | | 11 | ++++++----- |
| M | simulation.c | | | 432 | ++++++++++++++++++++++++++++++++++++++++--------------------------------------- |
| M | simulation.h | | | 38 | ++++++++++++++++++++++---------------- |
9 files changed, 606 insertions(+), 773 deletions(-)
diff --git a/1d_fd_simple_shear.c b/1d_fd_simple_shear.c
@@ -2,14 +2,12 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
-#include <getopt.h>
#include <string.h>
#include <time.h>
#include "simulation.h"
#include "fluid.h"
-
-#define PROGNAME "1d_fd_simple_shear"
+#include "arg.h"
#include "parameter_defaults.h"
@@ -28,119 +26,34 @@
/* uncomment to print time spent per time step to stdout */
/*#define BENCHMARK_PERFORMANCE*/
-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"
- "\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"
- " -k, --fluid-permeability VAL fluid intrinsic permeability [m^2] (default %g)\n"
- " -O, --fluid-pressure-top VAL fluid pressure at +z edge [Pa] (default %g)\n"
- " -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"
- ,
- __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);
- free(sim.phi);
- free(sim.k);
-}
+char *argv0;
+struct simulation sim;
static void
-version(void)
+usage(void)
{
- printf("%s v%s\n"
- "Licensed under the ISC License, see LICENSE for details.\n"
- "written by Anders Damsgaard, anders@adamsgaard.dk\n"
- "https://src.adamsgaard.dk/1d_fd_simple_shear\n"
- , PROGNAME, VERSION);
+ dprintf(2, "usage: %s [-FNThnv] "
+ "[-g gravity-accel] [-d grain-size] [-r grain-density] "
+ "[-m grain-friction] [-c grain-cohesion] "
+ "[-A grain-nonlocal-ampl] [-b grain-rate-dependance] "
+ "[-p grain-porosity] [-y min-porosity] [-Y max-porosity] "
+ "[-K dilatancy-constant] "
+ "[-n normal-stress] [-f applied-shear-friction] "
+ "[-s applied-shear-vel] [-l applied-shear-vel-limit] "
+ "[-o origo] [-L length] [-U resolution] "
+ "[-t curr-time] [-e end-time] [-I file-interval] "
+ "[-O fluid-pressure-top] [-a fluid-pressure-ampl] "
+ "[-q fluid-pressure-freq] [-H fluid-pressure-phase] "
+ "[-S fluid-pressure-pulse-shape] [-u fluid-pulse-time] "
+ "[-k fluid-permeability] [-R fluid-density] [-i fluid-viscosity] "
+ "[-C fluid-compressibility] [name]\n", argv0);
+ exit(1);
}
int
main(int argc, char* argv[])
{
- int norm, opt, i;
- struct simulation sim;
- const char* optstring;
+ int i, normalize;
unsigned long iter, stressiter;
double new_phi, new_k, filetimeclock, res_norm, mu_wall_orig;
#ifdef BENCHMARK_PERFORMANCE
@@ -150,193 +63,154 @@ main(int argc, char* argv[])
#ifdef __OpenBSD__
if (pledge("stdio wpath cpath", NULL) == -1) {
- fprintf(stderr, "error: pledge failed");
+ dprintf(2, "error: pledge failed");
exit(1);
}
#endif
- /* load with default values */
- sim = init_sim();
-
- norm = 0;
-
- 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:";
- 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'},
- {"fluid-permeability", required_argument, NULL, 'k'},
- {"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}
- };
+ atexit(free_arrays);
+ init_sim();
+ normalize = 0;
new_phi = sim.phi[0];
new_k = sim.k[0];
- while ((opt = getopt_long(argc, argv, optstring, longopts, NULL)) != -1) {
- switch (opt) {
- case -1: /* no more arguments */
- case 0: /* long options toggles */
- break;
-
- case 'h':
- usage();
- free(sim.phi);
- free(sim.k);
- return 0;
- case 'v':
- version();
- 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;
- case 'i':
- sim.mu_f = atof(optarg);
- break;
- case 'R':
- sim.rho_f = atof(optarg);
- break;
- case 'k':
- new_k = atof(optarg);
- break;
- case 'O':
- sim.p_f_top = atof(optarg);
- break;
- case 'a':
- sim.p_f_mod_ampl = atof(optarg);
- break;
- 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",
- argv[0]);
- return -2;
- }
- }
- for (i=optind; i<argc; ++i) {
- if (i>optind) {
- fprintf(stderr,
- "error: more than one simulation name specified\n");
+
+ ARGBEGIN {
+ case 'A':
+ sim.A = atof(EARGF(usage()));
+ break;
+ case 'C':
+ sim.beta_f = atof(EARGF(usage()));
+ break;
+ case 'F':
+ sim.fluid = 1;
+ break;
+ case 'H':
+ sim.p_f_mod_phase = atof(EARGF(usage()));
+ break;
+ case 'I':
+ sim.file_dt = atof(EARGF(usage()));
+ break;
+ case 'K':
+ sim.dilatancy_angle = atof(EARGF(usage()));
+ break;
+ case 'L':
+ sim.L_z = atof(EARGF(usage()));
+ break;
+ case 'N':
+ normalize = 1;
+ break;
+ case 'O':
+ sim.p_f_top = atof(EARGF(usage()));
+ break;
+ case 'R':
+ sim.rho_f = atof(EARGF(usage()));
+ break;
+ case 'S':
+ if (argv[1] == NULL)
+ usage();
+ else if (!strncmp(argv[1], "triangle", 8))
+ sim.p_f_mod_pulse_shape = 0;
+ else if (!strncmp(argv[1], "square", 6))
+ sim.p_f_mod_pulse_shape = 1;
+ else {
+ dprintf(2, "error: invalid pulse shape '%s'\n",
+ argv[1]);
return 1;
}
- snprintf(sim.name, sizeof(sim.name), "%s", argv[i]);
- }
+ argv++;
+ break;
+ case 'T':
+ sim.transient = 1;
+ break;
+ case 'U':
+ sim.nz = atoi(EARGF(usage()));
+ break;
+ case 'V':
+ sim.v_x_bot = atof(EARGF(usage()));
+ break;
+ case 'Y':
+ sim.phi_max = atof(EARGF(usage()));
+ break;
+ case 'a':
+ sim.p_f_mod_ampl = atof(EARGF(usage()));
+ break;
+ case 'b':
+ sim.b = atof(EARGF(usage()));
+ break;
+ case 'c':
+ sim.C = atof(EARGF(usage()));
+ break;
+ case 'd':
+ sim.d = atof(EARGF(usage()));
+ break;
+ case 'e':
+ sim.t_end = atof(EARGF(usage()));
+ break;
+ case 'f':
+ sim.mu_wall = atof(EARGF(usage()));
+ break;
+ case 'g':
+ sim.G = atof(EARGF(usage()));
+ break;
+ case 'h':
+ usage();
+ break;
+ case 'i':
+ sim.mu_f = atof(EARGF(usage()));
+ break;
+ case 'k':
+ new_k = atof(EARGF(usage()));
+ break;
+ case 'l':
+ sim.v_x_limit = atof(EARGF(usage()));
+ break;
+ case 'm':
+ sim.mu_s = atof(EARGF(usage()));
+ break;
+ case 'n':
+ sim.P_wall = atof(EARGF(usage()));
+ break;
+ case 'o':
+ sim.origo_z = atof(EARGF(usage()));
+ break;
+ case 'p':
+ new_phi = atof(EARGF(usage()));
+ break;
+ case 'q':
+ sim.p_f_mod_freq = atof(EARGF(usage()));
+ break;
+ case 'r':
+ sim.rho_s = atof(EARGF(usage()));
+ break;
+ case 's':
+ sim.v_x_fix = atof(EARGF(usage()));
+ break;
+ case 't':
+ sim.t = atof(EARGF(usage()));
+ break;
+ case 'u':
+ sim.p_f_mod_pulse_time = atof(EARGF(usage()));
+ break;
+ case 'v':
+ printf("%s-"VERSION"\n", argv0);
+ return 0;
+ case 'y':
+ sim.phi_min = atof(EARGF(usage()));
+ break;
+ default:
+ usage();
+ } ARGEND;
+
+ if (argc == 1 && argv[1])
+ snprintf(sim.name, sizeof(sim.name), "%s", argv[1]);
+ else if (argc > 1)
+ usage();
if (sim.nz < 1)
sim.nz = (int)ceil(sim.L_z/sim.d);
- prepare_arrays(&sim);
+ prepare_arrays();
if (!isnan(new_phi))
for (i=0; i<sim.nz; ++i)
@@ -345,14 +219,14 @@ main(int argc, char* argv[])
for (i=0; i<sim.nz; ++i)
sim.k[i] = new_k;
- lithostatic_pressure_distribution(&sim);
+ lithostatic_pressure_distribution();
if (sim.fluid) {
- hydrostatic_fluid_pressure_distribution(&sim);
- set_largest_fluid_timestep(&sim, 0.5);
+ hydrostatic_fluid_pressure_distribution();
+ set_largest_fluid_timestep(0.5);
}
- check_simulation_parameters(&sim);
+ check_simulation_parameters();
filetimeclock = 0.0;
iter = 0;
@@ -367,20 +241,20 @@ main(int argc, char* argv[])
stressiter = 0;
do {
if (sim.fluid) {
- if (darcy_solver_1d(&sim, MAX_ITER_DARCY, RTOL_DARCY))
+ if (darcy_solver_1d(MAX_ITER_DARCY, RTOL_DARCY))
exit(1);
}
- compute_effective_stress(&sim);
- compute_friction(&sim);
- compute_cooperativity_length(&sim);
+ compute_effective_stress();
+ compute_friction();
+ compute_cooperativity_length();
- if (implicit_1d_jacobian_poisson_solver(&sim, MAX_ITER_GRANULAR,
+ if (implicit_1d_jacobian_poisson_solver(MAX_ITER_GRANULAR,
RTOL_GRANULAR))
exit(1);
- compute_shear_strain_rate_plastic(&sim);
- compute_shear_velocity(&sim);
+ compute_shear_strain_rate_plastic();
+ compute_shear_velocity();
if (!isnan(sim.v_x_limit) || !isnan(sim.v_x_fix)) {
if (!isnan(sim.v_x_limit)) {
@@ -396,13 +270,11 @@ main(int argc, char* argv[])
}
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, "
+ dprintf(2, "error: stress solution did not converge:\n");
+ dprintf(2, "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;
}
@@ -423,18 +295,14 @@ main(int argc, char* argv[])
if ((filetimeclock - sim.dt >= sim.file_dt || iter == 1) &&
strcmp(sim.name, DEFAULT_SIMULATION_NAME) != 0) {
- write_output_file(&sim, norm);
+ write_output_file(normalize);
filetimeclock = 0.0;
if (iter == 1)
sim.t = 0.0;
}
}
- if (sim.fluid)
- print_wet_output(stdout, &sim, norm);
- else
- print_dry_output(stdout, &sim, norm);
+ print_output(stdout, normalize);
- free_arrays(&sim);
return 0;
}
diff --git a/Makefile b/Makefile
@@ -7,7 +7,7 @@ PREFIX ?= /usr/local
INSTALL ?= install
STRIP ?= strip
-VERSION = 0.4.6
+VERSION = 0.5.0
GLOBALCONST = -DVERSION=\"${VERSION}\"
default: ${BIN}
@@ -32,22 +32,9 @@ install: ${BIN}
uninstall:
rm -f ${DESTDIR}${PREFIX}/bin/${BIN}
-watch:
- echo ${SRC} ${HDR} | tr ' ' '\n' | entr -s 'make && ./1d_fd_simple_shear'
-
test: ${BIN}
make -C test/
-memtest: 1d_fd_simple_shear
- valgrind --error-exitcode=1 --leak-check=full 1d_fd_simple_shear -h
- valgrind --error-exitcode=1 --leak-check=full 1d_fd_simple_shear
- valgrind --error-exitcode=1 --leak-check=full 1d_fd_simple_shear -F
- valgrind --error-exitcode=1 --leak-check=full 1d_fd_simple_shear -F \
- --fluid-pressure-top 50e3 \
- --fluid-pressure-ampl 50e3 \
- --fluid-pressure-freq $(echo "1/3600" | bc -l) \
- --time-end 3600
-
clean:
rm -f *.txt
rm -f *.o
diff --git a/arg.h b/arg.h
@@ -0,0 +1,38 @@
+/* by 20h */
+#ifndef ARG_H
+#define ARG_H
+
+#define USED(x) ((void)(x))
+
+extern char *argv0;
+
+#define ARGBEGIN for(argv0 = *argv, argv++, argc--;\
+ argv[0] && argv[0][0] == '-'\
+ && argv[0][1];\
+ argc--, argv++) {\
+ char _argc;\
+ char **_argv;\
+ if(argv[0][1] == '-' && argv[0][2] == '\0') {\
+ argv++;\
+ argc--;\
+ break;\
+ }\
+ int i_;\
+ for(i_ = 1, _argv = argv; argv[0][i_];\
+ i_++) {\
+ if(_argv != argv)\
+ break;\
+ _argc = argv[0][i_];\
+ switch(_argc)
+
+#define ARGEND }\
+ USED(_argc);\
+ }\
+ USED(argv);\
+ USED(argc);
+
+#define EARGF(x) ((argv[1] == NULL)? ((x), abort(), (char *)0) :\
+ (argc--, argv++, argv[0]))
+
+#endif
+
diff --git a/fluid.c b/fluid.c
@@ -6,13 +6,13 @@
/* #define DEBUG true */
void
-hydrostatic_fluid_pressure_distribution(struct simulation* sim)
+hydrostatic_fluid_pressure_distribution()
{
int i;
- for (i=0; i<sim->nz; ++i)
- sim->p_f_ghost[i+1] = sim->p_f_top +
- sim->phi[i]*sim->rho_f*sim->G*
- (sim->L_z - sim->z[i]);
+ for (i=0; i<sim.nz; ++i)
+ sim.p_f_ghost[i+1] = sim.p_f_top +
+ sim.phi[i]*sim.rho_f*sim.G*
+ (sim.L_z - sim.z[i]);
}
/* Determines the largest time step for the current simulation state. Note
@@ -20,15 +20,15 @@ hydrostatic_fluid_pressure_distribution(struct simulation* sim)
* (i.e., permeabilities, porosities, viscosities, or compressibilities)
* change. The safety factor should be in ]0;1] */
void
-set_largest_fluid_timestep(struct simulation* sim, const double safety)
+set_largest_fluid_timestep(const double safety)
{
int i;
double dx_min, diff, diff_max;
double *dx;
- dx = spacing(sim->z, sim->nz);
+ dx = spacing(sim.z, sim.nz);
dx_min = INFINITY;
- for (i=0; i<sim->nz-1; ++i) {
+ for (i=0; i<sim.nz-1; ++i) {
if (dx[i] < 0.0) {
fprintf(stderr, "error: cell spacing negative (%g) in cell %d\n",
dx[i], i);
@@ -40,14 +40,14 @@ set_largest_fluid_timestep(struct simulation* sim, const double safety)
free(dx);
diff_max = -INFINITY;
- for (i=0; i<sim->nz; ++i) {
- diff = sim->k[i]/(sim->phi[i]*sim->beta_f*sim->mu_f);
+ for (i=0; i<sim.nz; ++i) {
+ diff = sim.k[i]/(sim.phi[i]*sim.beta_f*sim.mu_f);
if (diff > diff_max) diff_max = diff;
}
- sim->dt = safety*0.5*dx_min*dx_min/diff_max;
- if (sim->file_dt*0.5 < sim->dt)
- sim->dt = sim->file_dt;
+ sim.dt = safety*0.5*dx_min*dx_min/diff_max;
+ if (sim.file_dt*0.5 < sim.dt)
+ sim.dt = sim.file_dt;
}
static double
@@ -87,15 +87,15 @@ square_pulse(const double time,
}
static void
-set_fluid_bcs(struct simulation* sim, const double p_f_top)
+set_fluid_bcs(const double p_f_top)
{
- set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top);
- sim->p_f_ghost[sim->nz] = p_f_top; /* Include top node in BC */
- set_bc_neumann(sim->p_f_ghost,
- sim->nz,
+ set_bc_dirichlet(sim.p_f_ghost, sim.nz, +1, p_f_top);
+ sim.p_f_ghost[sim.nz] = p_f_top; /* Include top node in BC */
+ set_bc_neumann(sim.p_f_ghost,
+ sim.nz,
-1,
- sim->phi[0]*sim->rho_f*sim->G,
- sim->dz);
+ sim.phi[0]*sim.rho_f*sim.G,
+ sim.dz);
}
static double
@@ -169,8 +169,7 @@ darcy_pressure_change_1d(double* dp_f_dt,
}*/
int
-darcy_solver_1d(struct simulation* sim,
- const int max_iter,
+darcy_solver_1d(const int max_iter,
const double rel_tol)
{
int i, iter, solved;
@@ -193,99 +192,99 @@ darcy_solver_1d(struct simulation* sim,
/* TODO: values other than 1.0 do not work! */
theta = 1.0;
- if (isnan(sim->p_f_mod_pulse_time))
- p_f_top = sim->p_f_top + sine_wave(sim->t,
- sim->p_f_mod_ampl,
- sim->p_f_mod_freq,
- sim->p_f_mod_phase);
+ if (isnan(sim.p_f_mod_pulse_time))
+ p_f_top = sim.p_f_top + sine_wave(sim.t,
+ sim.p_f_mod_ampl,
+ sim.p_f_mod_freq,
+ sim.p_f_mod_phase);
else
- if (sim->p_f_mod_pulse_shape == 1)
- p_f_top = sim->p_f_top + square_pulse(sim->t,
- sim->p_f_mod_ampl,
- sim->p_f_mod_freq,
- sim->p_f_mod_pulse_time);
+ if (sim.p_f_mod_pulse_shape == 1)
+ p_f_top = sim.p_f_top + square_pulse(sim.t,
+ sim.p_f_mod_ampl,
+ sim.p_f_mod_freq,
+ sim.p_f_mod_pulse_time);
else
- p_f_top = sim->p_f_top + triangular_pulse(sim->t,
- sim->p_f_mod_ampl,
- sim->p_f_mod_freq,
- sim->p_f_mod_pulse_time);
+ p_f_top = sim.p_f_top + triangular_pulse(sim.t,
+ sim.p_f_mod_ampl,
+ sim.p_f_mod_freq,
+ sim.p_f_mod_pulse_time);
/* set fluid BCs (1 of 2) */
- set_fluid_bcs(sim, p_f_top);
+ set_fluid_bcs(p_f_top);
solved = 0;
if (epsilon < 1.0) {
/* compute explicit solution to pressure change */
- dp_f_dt_expl = zeros(sim->nz);
- for (i=0; i<sim->nz; ++i)
+ dp_f_dt_expl = zeros(sim.nz);
+ for (i=0; i<sim.nz; ++i)
dp_f_dt_expl[i] = darcy_pressure_change_1d(i,
- sim->nz,
- sim->p_f_ghost,
- sim->phi,
- sim->k,
- sim->dz,
- sim->beta_f,
- sim->mu_f);
+ sim.nz,
+ sim.p_f_ghost,
+ sim.phi,
+ sim.k,
+ sim.dz,
+ sim.beta_f,
+ sim.mu_f);
}
if (epsilon > 0.0) {
/* compute implicit solution with Jacobian iterations */
- dp_f_dt_impl = zeros(sim->nz);
- p_f_ghost_new = zeros(sim->nz+2);
- r_norm = zeros(sim->nz);
- p_f_ghost_old = empty(sim->nz+2);
- copy_values(sim->p_f_ghost, p_f_ghost_old, sim->nz+2);
+ dp_f_dt_impl = zeros(sim.nz);
+ p_f_ghost_new = zeros(sim.nz+2);
+ r_norm = zeros(sim.nz);
+ p_f_ghost_old = empty(sim.nz+2);
+ copy_values(sim.p_f_ghost, p_f_ghost_old, sim.nz+2);
for (iter=0; iter<max_iter; ++iter) {
/* set fluid BCs (2 of 2) */
- set_fluid_bcs(sim, p_f_top);
+ set_fluid_bcs( p_f_top);
#ifdef DEBUG
puts(".. p_f_ghost after BC:");
- print_array(sim->p_f_ghost, sim->nz+2);
+ print_array(sim.p_f_ghost, sim.nz+2);
#endif
- for (i=0; i<sim->nz-1; ++i)
+ for (i=0; i<sim.nz-1; ++i)
dp_f_dt_impl[i] = darcy_pressure_change_1d(i,
- sim->nz,
- sim->p_f_ghost,
- sim->phi,
- sim->k,
- sim->dz,
- sim->beta_f,
- sim->mu_f);
-
- for (i=0; i<sim->nz-1; ++i) {
+ sim.nz,
+ sim.p_f_ghost,
+ sim.phi,
+ sim.k,
+ sim.dz,
+ sim.beta_f,
+ sim.mu_f);
+
+ for (i=0; i<sim.nz-1; ++i) {
#ifdef DEBUG
printf("dp_f_expl[%d] = %g\ndp_f_impl[%d] = %g\n",
i, dp_f_dt_expl[i], i, dp_f_dt_impl[i]);
#endif
p_f_ghost_new[i+1] = p_f_ghost_old[i+1]
- + epsilon*dp_f_dt_impl[i]*sim->dt;
+ + epsilon*dp_f_dt_impl[i]*sim.dt;
if (epsilon < 1.0)
p_f_ghost_new[i+1] += (1.0 - epsilon)
- *dp_f_dt_expl[i]*sim->dt;
+ *dp_f_dt_expl[i]*sim.dt;
p_f_ghost_new[i+1] = p_f_ghost_old[i+1]*(1.0 - theta)
+ p_f_ghost_new[i+1]*theta;
- r_norm[i] = fabs((p_f_ghost_new[i+1] - sim->p_f_ghost[i+1])
- /(sim->p_f_ghost[i+1] + 1e-16));
+ r_norm[i] = fabs((p_f_ghost_new[i+1] - sim.p_f_ghost[i+1])
+ /(sim.p_f_ghost[i+1] + 1e-16));
}
- r_norm_max = max(r_norm, sim->nz);
+ r_norm_max = max(r_norm, sim.nz);
#ifdef DEBUG
puts(".. p_f_ghost_new:");
- print_array(p_f_ghost_new, sim->nz+2);
+ print_array(p_f_ghost_new, sim.nz+2);
#endif
- copy_values(p_f_ghost_new, sim->p_f_ghost, sim->nz+2);
+ copy_values(p_f_ghost_new, sim.p_f_ghost, sim.nz+2);
#ifdef DEBUG
puts(".. p_f_ghost after update:");
- print_array(sim->p_f_ghost, sim->nz+2);
+ print_array(sim.p_f_ghost, sim.nz+2);
#endif
if (r_norm_max <= rel_tol) {
@@ -307,17 +306,17 @@ darcy_solver_1d(struct simulation* sim,
fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max);
}
} else {
- for (i=0; i<sim->nz; ++i)
- sim->p_f_ghost[i+1] += dp_f_dt_expl[i]*sim->dt;
+ for (i=0; i<sim.nz; ++i)
+ sim.p_f_ghost[i+1] += dp_f_dt_expl[i]*sim.dt;
solved = 1;
#ifdef DEBUG
puts(".. dp_f_dt_expl:");
- print_array(dp_f_dt_expl, sim->nz);
+ print_array(dp_f_dt_expl, sim.nz);
puts(".. p_f_ghost after explicit solution:");
- print_array(sim->p_f_ghost, sim->nz+2);
+ print_array(sim.p_f_ghost, sim.nz+2);
#endif
}
- set_fluid_bcs(sim, p_f_top);
+ set_fluid_bcs(p_f_top);
if (epsilon < 1.0)
free(dp_f_dt_expl);
diff --git a/fluid.h b/fluid.h
@@ -3,12 +3,13 @@
#include "simulation.h"
-void hydrostatic_fluid_pressure_distribution(struct simulation* sim);
+extern struct simulation sim;
-void set_largest_fluid_timestep(struct simulation* sim, const double safety);
+void hydrostatic_fluid_pressure_distribution();
+
+void set_largest_fluid_timestep(const double safety);
int darcy_solver_1d(
- struct simulation* sim,
const int max_iter,
const double rel_tol);
diff --git a/max_depth_simple_shear.c b/max_depth_simple_shear.c
@@ -6,10 +6,8 @@
#include <time.h>
#include "simulation.h"
-
-#define PROGNAME "max_depth_simple_shear"
-
#include "parameter_defaults.h"
+#include "arg.h"
#define EPS 1e-12
@@ -24,64 +22,32 @@
#define M_PI 3.1415926535897932384626433832795028841971693993751058209749445923078164062
#endif
-static void
-usage(void)
-{
- struct simulation sim = init_sim();
- 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"
- "\nOptional arguments:\n"
- " -v, --version show version information\n"
- " -h, --help show this message\n"
- " -G, --gravity VAL gravity magnitude [m/s^2] (default %g)\n"
- " -p, --porosity VAL porosity fraction [-] (default %g)\n"
- " -r, --density VAL grain material density [kg/m^3] (default %g)\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"
- " -k, --fluid-permeability VAL fluid intrinsic permeability [m^2] (default %g)\n"
- " -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"
- "\n"
- "The first column of the output is the max. deformation depth from the top in meters.\n"
- "The second column is the skin depth in meters.\n"
- ,
- __func__, PROGNAME,
- sim.G,
- sim.phi[0],
- sim.rho_s,
- sim.beta_f,
- sim.mu_f,
- sim.rho_f,
- sim.k[0],
- sim.p_f_mod_ampl,
- sim.p_f_mod_freq);
- free(sim.phi);
- free(sim.k);
-}
+char *argv0;
+struct simulation sim;
static void
-version(void)
+usage(void)
{
- printf("%s v%s\n"
- "Licensed under the ISC License, see LICENSE for details.\n"
- "written by Anders Damsgaard, anders@adamsgaard.dk\n"
- "https://src.adamsgaard.dk/1d_fd_simple_shear\n"
- , PROGNAME, VERSION);
+ dprintf(2, "usage: %s [-hv] "
+ "[-g gravity-accel] [-r grain-density] "
+ "[-p grain-porosity] "
+ "[-O fluid-pressure-top] [-a fluid-pressure-ampl] "
+ "[-q fluid-pressure-freq]"
+ "[-k fluid-permeability] [-R fluid-density] [-i fluid-viscosity] "
+ "[-C fluid-compressibility]\n", argv0);
+ exit(1);
}
double
-skin_depth(const struct simulation* sim)
+skin_depth()
{
- return sqrt(sim->k[0]/
- (sim->phi[0]*sim->mu_f*sim->beta_f*M_PI*sim->p_f_mod_freq));
+ return sqrt(sim.k[0]/
+ (sim.phi[0]*sim.mu_f*sim.beta_f*M_PI*sim.p_f_mod_freq));
}
/* using alternate form: sin(x) + cos(x) = sqrt(2)*sin(x + pi/4) */
double
-eff_normal_stress_gradient(const struct simulation* sim, double d_s, double z_)
+eff_normal_stress_gradient(double d_s, double z_)
{
if (z_/d_s > 10.0) {
fprintf(stderr, "error: %s: unrealistic depth: %g m "
@@ -90,13 +56,12 @@ eff_normal_stress_gradient(const struct simulation* sim, double d_s, double z_)
}
return sqrt(2.0)*sin((3.0*M_PI/2.0 - z_/d_s) + M_PI/4.0)
- + (sim->rho_s - sim->rho_f)*sim->G*d_s/sim->p_f_mod_ampl*exp(z_/d_s);
+ + (sim.rho_s - sim.rho_f)*sim.G*d_s/sim.p_f_mod_ampl*exp(z_/d_s);
}
/* use Brent's method for finding roots (Press et al., 2007) */
-double zbrent(struct simulation* sim,
- double d_s,
- double (*f)(struct simulation* sim, double, double),
+double zbrent(double d_s,
+ double (*f)(double, double),
double x1,
double x2,
double tol)
@@ -105,8 +70,8 @@ double zbrent(struct simulation* sim,
double a, b, c, d, e, min1, min2, fa, fb, fc, p, q, r, s, tol1, xm;
a = x1; b=x2; c=x2;
- fa = (*f)(sim, d_s, a);
- fb = (*f)(sim, d_s, b);
+ fa = (*f)(d_s, a);
+ fb = (*f)(d_s, b);
if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) {
fprintf(stderr,
@@ -167,7 +132,7 @@ double zbrent(struct simulation* sim,
b += d;
else
b += ((xm) >= 0.0 ? fabs(tol1) : -fabs(tol1));
- fb = (*f)(sim, d_s, b);
+ fb = (*f)(d_s, b);
}
fprintf(stderr,
"error: %s: exceeded maximum number of iterations",
@@ -179,9 +144,7 @@ double zbrent(struct simulation* sim,
int
main(int argc, char* argv[])
{
- int opt, i;
- struct simulation sim;
- const char* optstring;
+ int i;
double new_phi, new_k;
double d_s, depth, depth_limit1, depth_limit2;
#ifdef BENCHMARK_PERFORMANCE
@@ -195,86 +158,56 @@ main(int argc, char* argv[])
exit(1);
}
#endif
+ atexit(free_arrays);
- /* load with default values */
- sim = init_sim();
-
- optstring = "hvG:p:r:c:i:R:k:a:q:";
- const struct option longopts[] = {
- {"help", no_argument, NULL, 'h'},
- {"version", no_argument, NULL, 'v'},
- {"gravity", required_argument, NULL, 'G'},
- {"porosity", required_argument, NULL, 'p'},
- {"density", required_argument, NULL, 'r'},
- {"fluid-compressibility", required_argument, NULL, 'c'},
- {"fluid-viscosity", required_argument, NULL, 'i'},
- {"fluid-density", required_argument, NULL, 'R'},
- {"fluid-permeability", required_argument, NULL, 'k'},
- {"fluid-pressure-ampl", required_argument, NULL, 'a'},
- {"fluid-pressure-freq", required_argument, NULL, 'q'},
- {NULL, 0, NULL, 0}
- };
-
+ init_sim();
new_phi = sim.phi[0];
new_k = sim.k[0];
- while ((opt = getopt_long(argc, argv, optstring, longopts, NULL)) != -1) {
- switch (opt) {
- case -1: /* no more arguments */
- case 0: /* long options toggles */
- break;
-
- case 'h':
- usage();
- free(sim.phi);
- free(sim.k);
- return 0;
- case 'v':
- version();
- free(sim.phi);
- free(sim.k);
- return 0;
- case 'G':
- sim.G = atof(optarg);
- break;
- case 'p':
- new_phi = atof(optarg);
- break;
- case 'r':
- sim.rho_s = atof(optarg);
- break;
- case 'c':
- sim.beta_f = atof(optarg);
- break;
- case 'i':
- sim.mu_f = atof(optarg);
- break;
- case 'R':
- sim.rho_f = atof(optarg);
- break;
- case 'k':
- new_k = atof(optarg);
- break;
- case 'a':
- sim.p_f_mod_ampl = atof(optarg);
- break;
- case 'q':
- sim.p_f_mod_freq = atof(optarg);
- break;
- default:
- fprintf(stderr, "%s: invalid option -- %c\n", argv[0], opt);
- fprintf(stderr, "Try `%s --help` for more information\n",
- argv[0]);
- return -2;
- }
- }
- if (optind < argc) {
- fprintf(stderr, "error: unknown argument specified\n");
- return 1;
- }
+ ARGBEGIN {
+ case 'C':
+ sim.beta_f = atof(EARGF(usage()));
+ break;
+ case 'R':
+ sim.rho_f = atof(EARGF(usage()));
+ break;
+ case 'a':
+ sim.p_f_mod_ampl = atof(EARGF(usage()));
+ break;
+ case 'g':
+ sim.G = atof(EARGF(usage()));
+ break;
+ case 'h':
+ usage();
+ break;
+ case 'i':
+ sim.mu_f = atof(EARGF(usage()));
+ break;
+ case 'k':
+ new_k = atof(EARGF(usage()));
+ break;
+ case 'p':
+ new_phi = atof(EARGF(usage()));
+ break;
+ case 'q':
+ sim.p_f_mod_freq = atof(EARGF(usage()));
+ break;
+ case 'r':
+ sim.rho_s = atof(EARGF(usage()));
+ break;
+ case 'v':
+ printf("%s-"VERSION"\n", argv0);
+ return 0;
+ break;
+ default:
+ usage();
+ } ARGEND;
+
+ if (argc > 0)
+ usage();
sim.nz = 2;
- prepare_arrays(&sim);
+ prepare_arrays();
if (!isnan(new_phi))
for (i=0; i<sim.nz; ++i)
@@ -283,10 +216,10 @@ main(int argc, char* argv[])
for (i=0; i<sim.nz; ++i)
sim.k[i] = new_k;
- check_simulation_parameters(&sim);
+ check_simulation_parameters();
depth = 0.0;
- d_s = skin_depth(&sim);
+ d_s = skin_depth();
#ifdef BENCHMARK_PERFORMANCE
t_begin = clock();
@@ -296,14 +229,13 @@ main(int argc, char* argv[])
* water pressure forcing, or if the stress gradient is positive at
* zero depth */
if (fabs(sim.p_f_mod_ampl) > 1e-6 &&
- eff_normal_stress_gradient(&sim, d_s, 0.0) < 0.0) {
+ eff_normal_stress_gradient(d_s, 0.0) < 0.0) {
depth_limit1 = 0.0;
depth_limit2 = 5.0*d_s;
- depth = zbrent(&sim,
- d_s,
- (double (*)(struct simulation*, double, double))
+ depth = zbrent(d_s,
+ (double (*)(double, double))
eff_normal_stress_gradient,
depth_limit1,
depth_limit2,
@@ -318,6 +250,5 @@ main(int argc, char* argv[])
printf("%.17g\t%.17g\n", depth, d_s);
- free_arrays(&sim);
return 0;
}
diff --git a/parameter_defaults.h b/parameter_defaults.h
@@ -9,10 +9,9 @@
#define DEFAULT_SIMULATION_NAME "unnamed_simulation"
/* Simulation settings */
-struct simulation init_sim(void)
+void
+init_sim()
{
- struct simulation sim;
-
snprintf(sim.name, sizeof(sim.name), DEFAULT_SIMULATION_NAME);
sim.G = 9.81;
@@ -43,6 +42,10 @@ struct simulation init_sim(void)
sim.phi = initval(0.25, 1);
sim.d = 0.04; /* Damsgaard et al 2013 */
+ sim.phi_min = 0.2;
+ sim.phi_max = 0.8;
+ sim.dilatancy_angle = 1.0;
+
/* Iverson et al 1997, 1998: Storglaciaren till */
/* sim.mu_s = tan(DEG2RAD(26.3)); */
/* sim.C = 5.0e3; */
@@ -107,8 +110,6 @@ struct simulation init_sim(void)
sim.p_f_mod_phase = 0.0;
sim.p_f_mod_pulse_time = NAN;
sim.p_f_mod_pulse_shape = 0;
-
- return sim;
}
#endif
diff --git a/simulation.c b/simulation.c
@@ -7,45 +7,49 @@
/* #define SHOW_PARAMETERS */
void
-prepare_arrays(struct simulation* sim)
+prepare_arrays()
{
- if (sim->nz < 2) {
+ if (sim.nz < 2) {
fprintf(stderr, "error: grid size (nz) must be at least 2 but is %d\n",
- sim->nz);
+ sim.nz);
exit(1);
}
- sim->z = linspace(sim->origo_z, /* spatial coordinates */
- sim->origo_z + sim->L_z,
- sim->nz);
- sim->dz = sim->z[1] - sim->z[0]; /* cell spacing */
- sim->mu = zeros(sim->nz); /* stress ratio */
- sim->sigma_n_eff = zeros(sim->nz); /* effective normal stress */
- sim->sigma_n = zeros(sim->nz); /* normal stess */
- sim->p_f_ghost = zeros(sim->nz+2); /* fluid pressure with ghost nodes */
- free(sim->phi);
- sim->phi = zeros(sim->nz); /* porosity */
- free(sim->k);
- sim->k = zeros(sim->nz); /* permeability */
- sim->xi = zeros(sim->nz); /* cooperativity length */
- sim->gamma_dot_p = zeros(sim->nz); /* shear velocity */
- sim->v_x = zeros(sim->nz); /* shear velocity */
- sim->g_ghost = zeros(sim->nz+2); /* fluidity with ghost nodes */
+ sim.z = linspace(sim.origo_z, /* spatial coordinates */
+ sim.origo_z + sim.L_z,
+ sim.nz);
+ sim.dz = sim.z[1] - sim.z[0]; /* cell spacing */
+ sim.mu = zeros(sim.nz); /* stress ratio */
+ sim.sigma_n_eff = zeros(sim.nz); /* effective normal stress */
+ sim.sigma_n = zeros(sim.nz); /* normal stess */
+ sim.p_f_ghost = zeros(sim.nz+2); /* fluid pressure with ghost nodes */
+ free(sim.phi);
+ sim.phi = zeros(sim.nz); /* porosity */
+ free(sim.k);
+ sim.k = zeros(sim.nz); /* permeability */
+ sim.xi = zeros(sim.nz); /* cooperativity length */
+ sim.gamma_dot_p = zeros(sim.nz); /* shear velocity */
+ sim.v_x = zeros(sim.nz); /* shear velocity */
+ sim.g_ghost = zeros(sim.nz+2); /* fluidity with ghost nodes */
+ sim.I = zeros(sim.nz); /* inertia number */
+ sim.tau = zeros(sim.nz); /* shear stress */
}
void
-free_arrays(struct simulation* sim)
+free_arrays()
{
- free(sim->z);
- free(sim->mu);
- free(sim->sigma_n_eff);
- free(sim->sigma_n);
- free(sim->p_f_ghost);
- free(sim->k);
- free(sim->phi);
- free(sim->xi);
- free(sim->gamma_dot_p);
- free(sim->v_x);
- free(sim->g_ghost);
+ free(sim.z);
+ free(sim.mu);
+ free(sim.sigma_n_eff);
+ free(sim.sigma_n);
+ free(sim.p_f_ghost);
+ free(sim.k);
+ free(sim.phi);
+ free(sim.xi);
+ free(sim.gamma_dot_p);
+ free(sim.v_x);
+ free(sim.g_ghost);
+ free(sim.I);
+ free(sim.tau);
}
static void
@@ -80,137 +84,156 @@ check_float(const char name[], const double value, int* return_status)
}
void
-check_simulation_parameters(const struct simulation* sim)
+check_simulation_parameters()
{
int return_status;
return_status = 0;
- check_float("sim.G", sim->G, &return_status);
- if (sim->G < 0.0)
- warn_parameter_value("sim.G is negative", sim->G, &return_status);
+ check_float("sim.G", sim.G, &return_status);
+ if (sim.G < 0.0)
+ warn_parameter_value("sim.G is negative", sim.G, &return_status);
- check_float("sim.P_wall", sim->P_wall, &return_status);
- if (sim->P_wall < 0.0)
- warn_parameter_value("sim.P_wall is negative", sim->P_wall,
+ check_float("sim.P_wall", sim.P_wall, &return_status);
+ if (sim.P_wall < 0.0)
+ warn_parameter_value("sim.P_wall is negative", sim.P_wall,
&return_status);
- check_float("sim.v_x_bot", sim->v_x_bot, &return_status);
+ check_float("sim.v_x_bot", sim.v_x_bot, &return_status);
- check_float("sim.mu_wall", sim->mu_wall, &return_status);
- if (sim->mu_wall < 0.0)
- warn_parameter_value("sim.mu_wall is negative", sim->mu_wall,
+ check_float("sim.mu_wall", sim.mu_wall, &return_status);
+ if (sim.mu_wall < 0.0)
+ warn_parameter_value("sim.mu_wall is negative", sim.mu_wall,
&return_status);
- check_float("sim.A", sim->A, &return_status);
- if (sim->A < 0.0)
- warn_parameter_value("sim.A is negative", sim->A, &return_status);
+ check_float("sim.A", sim.A, &return_status);
+ if (sim.A < 0.0)
+ warn_parameter_value("sim.A is negative", sim.A, &return_status);
- check_float("sim.b", sim->b, &return_status);
- if (sim->b < 0.0)
- warn_parameter_value("sim.b is negative", sim->b, &return_status);
+ check_float("sim.b", sim.b, &return_status);
+ if (sim.b < 0.0)
+ warn_parameter_value("sim.b is negative", sim.b, &return_status);
- check_float("sim.mu_s", sim->mu_s, &return_status);
- if (sim->mu_s < 0.0)
- warn_parameter_value("sim.mu_s is negative", sim->mu_s,
+ check_float("sim.mu_s", sim.mu_s, &return_status);
+ if (sim.mu_s < 0.0)
+ warn_parameter_value("sim.mu_s is negative", sim.mu_s,
&return_status);
- check_float("sim.C", sim->C, &return_status);
+ check_float("sim.C", sim.C, &return_status);
- check_float("sim.d", sim->d, &return_status);
- if (sim->d <= 0.0)
- warn_parameter_value("sim.d is not a positive number", sim->d,
+ check_float("sim.d", sim.d, &return_status);
+ if (sim.d <= 0.0)
+ warn_parameter_value("sim.d is not a positive number", sim.d,
&return_status);
- check_float("sim.rho_s", sim->rho_s, &return_status);
- if (sim->rho_s <= 0.0)
- warn_parameter_value("sim.rho_s is not a positive number", sim->rho_s,
+ check_float("sim.rho_s", sim.rho_s, &return_status);
+ if (sim.rho_s <= 0.0)
+ warn_parameter_value("sim.rho_s is not a positive number", sim.rho_s,
&return_status);
- if (sim->nz <= 0)
- warn_parameter_value("sim.nz is not a positive number", sim->nz,
+ if (sim.nz <= 0)
+ warn_parameter_value("sim.nz is not a positive number", sim.nz,
&return_status);
- check_float("sim.origo_z", sim->origo_z, &return_status);
- check_float("sim.L_z", sim->L_z, &return_status);
- if (sim->L_z <= sim->origo_z)
+ check_float("sim.origo_z", sim.origo_z, &return_status);
+ check_float("sim.L_z", sim.L_z, &return_status);
+ if (sim.L_z <= sim.origo_z)
warn_parameter_value("sim.L_z is smaller or equal to sim.origo_z",
- sim->L_z, &return_status);
+ sim.L_z, &return_status);
- if (sim->nz <= 0)
- warn_parameter_value("sim.nz is not a positive number", sim->nz,
+ if (sim.nz <= 0)
+ warn_parameter_value("sim.nz is not a positive number", sim.nz,
&return_status);
- check_float("sim.dz", sim->dz, &return_status);
- if (sim->dz <= 0.0)
- warn_parameter_value("sim.dz is not a positive number", sim->dz,
+ check_float("sim.dz", sim.dz, &return_status);
+ if (sim.dz <= 0.0)
+ warn_parameter_value("sim.dz is not a positive number", sim.dz,
&return_status);
- check_float("sim.t", sim->t, &return_status);
- if (sim->t < 0.0)
+ check_float("sim.t", sim.t, &return_status);
+ if (sim.t < 0.0)
warn_parameter_value("sim.t is a negative number",
- sim->t, &return_status);
+ sim.t, &return_status);
- check_float("sim.t_end", sim->t_end, &return_status);
- if (sim->t > sim->t_end)
+ check_float("sim.t_end", sim.t_end, &return_status);
+ if (sim.t > sim.t_end)
warn_parameter_value("sim.t_end is smaller than sim.t",
- sim->t, &return_status);
+ sim.t, &return_status);
- check_float("sim.dt", sim->dt, &return_status);
- if (sim->dt <= 0.0)
+ check_float("sim.dt", sim.dt, &return_status);
+ if (sim.dt <= 0.0)
warn_parameter_value("sim.dt is not a positive number",
- sim->dt, &return_status);
+ sim.dt, &return_status);
- check_float("sim.file_dt", sim->file_dt, &return_status);
- if (sim->file_dt < 0.0)
+ check_float("sim.file_dt", sim.file_dt, &return_status);
+ if (sim.file_dt < 0.0)
warn_parameter_value("sim.file_dt is a negative number",
- sim->file_dt, &return_status);
+ sim.file_dt, &return_status);
- check_float("sim.phi[0]", sim->phi[0], &return_status);
- if (sim->phi[0] < 0.0 || sim->phi[0] > 1.0)
+ check_float("sim.phi[0]", sim.phi[0], &return_status);
+ if (sim.phi[0] < 0.0 || sim.phi[0] > 1.0)
warn_parameter_value("sim.phi[0] is not within [0;1]",
- sim->phi[0], &return_status);
+ sim.phi[0], &return_status);
- if (sim->fluid != 0 && sim->fluid != 1)
+ check_float("sim.phi_min", sim.phi_min, &return_status);
+ if (sim.phi_min < 0.0 || sim.phi_min > 1.0)
+ warn_parameter_value("sim.phi_min is not within [0;1]",
+ sim.phi_min, &return_status);
+
+ check_float("sim.phi_max", sim.phi_max, &return_status);
+ if (sim.phi_max < 0.0 || sim.phi_max > 1.0)
+ warn_parameter_value("sim.phi_max is not within [0;1]",
+ sim.phi_max, &return_status);
+
+ check_float("sim.dilatancy_angle", sim.dilatancy_angle, &return_status);
+ if (sim.dilatancy_angle < 0.0 || sim.dilatancy_angle > 100.0)
+ warn_parameter_value("sim.dilatancy_angle is not within [0;100]",
+ sim.dilatancy_angle, &return_status);
+
+ if (sim.fluid != 0 && sim.fluid != 1)
warn_parameter_value("sim.fluid has an invalid value",
- (double)sim->fluid, &return_status);
+ (double)sim.fluid, &return_status);
- if (sim->fluid) {
+ if (sim.transient != 0 && sim.transient != 1)
+ warn_parameter_value("sim.fluid has an invalid value",
+ (double)sim.fluid, &return_status);
+
+ if (sim.fluid) {
- check_float("sim.p_f_mod_ampl", sim->p_f_mod_ampl, &return_status);
- if (sim->p_f_mod_ampl < 0.0)
+ check_float("sim.p_f_mod_ampl", sim.p_f_mod_ampl, &return_status);
+ if (sim.p_f_mod_ampl < 0.0)
warn_parameter_value("sim.p_f_mod_ampl is not a zero or positive",
- sim->p_f_mod_ampl, &return_status);
+ sim.p_f_mod_ampl, &return_status);
- if (sim->P_wall - sim->p_f_mod_ampl < 0.0)
+ if (sim.P_wall - sim.p_f_mod_ampl < 0.0)
warn_parameter_value("sim.P_wall - sim.p_f_mod_ampl is negative",
- sim->P_wall - sim->p_f_mod_ampl,
+ sim.P_wall - sim.p_f_mod_ampl,
&return_status);
- check_float("sim.p_f_mod_freq", sim->p_f_mod_freq, &return_status);
- if (sim->p_f_mod_freq < 0.0)
+ check_float("sim.p_f_mod_freq", sim.p_f_mod_freq, &return_status);
+ if (sim.p_f_mod_freq < 0.0)
warn_parameter_value("sim.p_f_mod_freq is not a zero or positive",
- sim->p_f_mod_freq, &return_status);
+ sim.p_f_mod_freq, &return_status);
- check_float("sim.beta_f", sim->beta_f, &return_status);
- if (sim->beta_f <= 0.0)
+ check_float("sim.beta_f", sim.beta_f, &return_status);
+ if (sim.beta_f <= 0.0)
warn_parameter_value("sim.beta_f is not positive",
- sim->beta_f, &return_status);
+ sim.beta_f, &return_status);
- check_float("sim.mu_f", sim->mu_f, &return_status);
- if (sim->mu_f <= 0.0)
+ check_float("sim.mu_f", sim.mu_f, &return_status);
+ if (sim.mu_f <= 0.0)
warn_parameter_value("sim.mu_f is not positive",
- sim->mu_f, &return_status);
+ sim.mu_f, &return_status);
- check_float("sim.rho_f", sim->rho_f, &return_status);
- if (sim->rho_f <= 0.0)
+ check_float("sim.rho_f", sim.rho_f, &return_status);
+ if (sim.rho_f <= 0.0)
warn_parameter_value("sim.rho_f is not positive",
- sim->rho_f, &return_status);
+ sim.rho_f, &return_status);
- check_float("sim.k[0]", sim->k[0], &return_status);
- if (sim->k[0] <= 0.0)
+ check_float("sim.k[0]", sim.k[0], &return_status);
+ if (sim.k[0] <= 0.0)
warn_parameter_value("sim.k[0] is not positive",
- sim->k[0], &return_status);
+ sim.k[0], &return_status);
}
if (return_status != 0) {
@@ -220,28 +243,28 @@ check_simulation_parameters(const struct simulation* sim)
}
void
-lithostatic_pressure_distribution(struct simulation* sim)
+lithostatic_pressure_distribution()
{
int i;
- for (i=0; i<sim->nz; ++i)
- sim->sigma_n[i] = sim->P_wall +
- (1.0 - sim->phi[i])*sim->rho_s*sim->G*
- (sim->L_z - sim->z[i]);
+ for (i=0; i<sim.nz; ++i)
+ sim.sigma_n[i] = sim.P_wall +
+ (1.0 - sim.phi[i])*sim.rho_s*sim.G*
+ (sim.L_z - sim.z[i]);
}
void
-compute_friction(struct simulation* sim)
+compute_friction()
{
int i;
- if (sim->fluid)
- for (i=0; i<sim->nz; ++i)
- sim->mu[i] = sim->mu_wall/
- (sim->sigma_n_eff[i]/(sim->P_wall - sim->p_f_top));
+ if (sim.fluid)
+ for (i=0; i<sim.nz; ++i)
+ sim.mu[i] = sim.mu_wall/
+ (sim.sigma_n_eff[i]/(sim.P_wall - sim.p_f_top));
else
- for (i=0; i<sim->nz; ++i)
- sim->mu[i] = sim->mu_wall/
- (1.0 + (1.0 - sim->phi[i])*sim->rho_s*sim->G*
- (sim->L_z - sim->z[i])/sim->P_wall);
+ for (i=0; i<sim.nz; ++i)
+ sim.mu[i] = sim.mu_wall/
+ (1.0 + (1.0 - sim.phi[i])*sim.rho_s*sim.G*
+ (sim.L_z - sim.z[i])/sim.P_wall);
}
double
@@ -251,37 +274,37 @@ shear_strain_rate_plastic(const double fluidity, const double friction)
}
void
-compute_shear_strain_rate_plastic(struct simulation* sim)
+compute_shear_strain_rate_plastic()
{
int i;
- for (i=0; i<sim->nz; ++i)
- sim->gamma_dot_p[i] = shear_strain_rate_plastic(sim->g_ghost[i+1],
- sim->mu[i]);
+ for (i=0; i<sim.nz; ++i)
+ sim.gamma_dot_p[i] = shear_strain_rate_plastic(sim.g_ghost[i+1],
+ sim.mu[i]);
}
void
-compute_shear_velocity(struct simulation* sim)
+compute_shear_velocity()
{
int i;
// TODO: implement iterative solver
// Dirichlet BC at bottom
- sim->v_x[0] = sim->v_x_bot;
+ sim.v_x[0] = sim.v_x_bot;
- for (i=1; i<sim->nz; ++i)
- sim->v_x[i] = sim->v_x[i-1] + sim->gamma_dot_p[i]*sim->dz;
+ for (i=1; i<sim.nz; ++i)
+ sim.v_x[i] = sim.v_x[i-1] + sim.gamma_dot_p[i]*sim.dz;
}
void
-compute_effective_stress(struct simulation* sim)
+compute_effective_stress()
{
int i;
- if (sim->fluid)
- for (i=0; i<sim->nz; ++i)
- sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_ghost[i+1];
+ if (sim.fluid)
+ for (i=0; i<sim.nz; ++i)
+ sim.sigma_n_eff[i] = sim.sigma_n[i] - sim.p_f_ghost[i+1];
else
- for (i=0; i<sim->nz; ++i)
- sim->sigma_n_eff[i] = sim->sigma_n[i];
+ for (i=0; i<sim.nz; ++i)
+ sim.sigma_n_eff[i] = sim.sigma_n[i];
}
static double
@@ -296,16 +319,16 @@ cooperativity_length(const double A,
}
void
-compute_cooperativity_length(struct simulation* sim)
+compute_cooperativity_length()
{
int i;
- for (i=0; i<sim->nz; ++i)
- sim->xi[i] = cooperativity_length(sim->A,
- sim->d,
- sim->mu[i],
- sim->sigma_n_eff[i],
- sim->mu_s,
- sim->C);
+ for (i=0; i<sim.nz; ++i)
+ sim.xi[i] = cooperativity_length(sim.A,
+ sim.d,
+ sim.mu[i],
+ sim.sigma_n_eff[i],
+ sim.mu_s,
+ sim.C);
}
static double
@@ -324,17 +347,17 @@ local_fluidity(const double p,
}
void
-compute_local_fluidity(struct simulation* sim)
+compute_local_fluidity()
{
int i;
- for (i=0; i<sim->nz; ++i)
- sim->g_ghost[i+1] = local_fluidity(sim->sigma_n_eff[i],
- sim->mu[i],
- sim->mu_s,
- sim->C,
- sim->b,
- sim->rho_s,
- sim->d);
+ for (i=0; i<sim.nz; ++i)
+ sim.g_ghost[i+1] = local_fluidity(sim.sigma_n_eff[i],
+ sim.mu[i],
+ sim.mu_s,
+ sim.C,
+ sim.b,
+ sim.rho_s,
+ sim.d);
}
void
@@ -406,8 +429,7 @@ poisson_solver_1d_cell_update(int i,
}
int
-implicit_1d_jacobian_poisson_solver(struct simulation* sim,
- const int max_iter,
+implicit_1d_jacobian_poisson_solver(const int max_iter,
const double rel_tol)
{
double r_norm_max;
@@ -415,41 +437,41 @@ implicit_1d_jacobian_poisson_solver(struct simulation* sim,
int iter, i;
r_norm_max = NAN;
- g_ghost_out = empty(sim->nz+2);
- r_norm = empty(sim->nz);
+ g_ghost_out = empty(sim.nz+2);
+ r_norm = empty(sim.nz);
for (iter=0; iter<max_iter; ++iter) {
#ifdef DEBUG
printf("\n@@@ ITERATION %d @@@\n", iter);
#endif
/* Dirichlet BCs resemble fixed particle velocities */
- set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0);
- set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0);
+ set_bc_dirichlet(sim.g_ghost, sim.nz, -1, 0.0);
+ set_bc_dirichlet(sim.g_ghost, sim.nz, +1, 0.0);
/* Neumann BCs resemble free surfaces */
- /* set_bc_neumann(sim->g_ghost, sim->nz, +1, 0.0, sim->dz); */
+ /* set_bc_neumann(sim.g_ghost, sim.nz, +1, 0.0, sim.dz); */
- for (i=0; i<sim->nz; ++i)
+ for (i=0; i<sim.nz; ++i)
poisson_solver_1d_cell_update(i,
- sim->g_ghost,
+ sim.g_ghost,
g_ghost_out,
r_norm,
- sim->dz,
- sim->mu,
- sim->sigma_n_eff,
- sim->xi,
- sim->mu_s,
- sim->C,
- sim->b,
- sim->rho_s,
- sim->d);
- r_norm_max = max(r_norm, sim->nz);
-
- copy_values(g_ghost_out, sim->g_ghost, sim->nz+2);
+ sim.dz,
+ sim.mu,
+ sim.sigma_n_eff,
+ sim.xi,
+ sim.mu_s,
+ sim.C,
+ sim.b,
+ sim.rho_s,
+ sim.d);
+ r_norm_max = max(r_norm, sim.nz);
+
+ copy_values(g_ghost_out, sim.g_ghost, sim.nz+2);
if (r_norm_max <= rel_tol) {
- set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0);
- set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0);
+ set_bc_dirichlet(sim.g_ghost, sim.nz, -1, 0.0);
+ set_bc_dirichlet(sim.g_ghost, sim.nz, +1, 0.0);
free(g_ghost_out);
free(r_norm);
/* printf(".. Solution converged after %d iterations\n", iter); */
@@ -466,64 +488,44 @@ implicit_1d_jacobian_poisson_solver(struct simulation* sim,
}
void
-write_output_file(struct simulation* sim, const int normalize)
+write_output_file(const int normalize)
{
char outfile[200];
FILE *fp;
snprintf(outfile, sizeof(outfile), "%s.output%05d.txt",
- sim->name, sim->n_file++);
+ sim.name, sim.n_file++);
fp = fopen(outfile, "w");
- if (sim->fluid)
- print_wet_output(fp, sim, normalize);
- else
- print_dry_output(fp, sim, normalize);
-
+ print_output(fp, normalize);
fclose(fp);
}
void
-print_dry_output(FILE* fp, struct simulation* sim, const int norm)
+print_output(FILE* fp, const int norm)
{
int i;
double *v_x_out;
if (norm)
- v_x_out = normalize(sim->v_x, sim->nz);
+ v_x_out = normalize(sim.v_x, sim.nz);
else
- v_x_out = copy(sim->v_x, sim->nz);
-
- for (i=0; i<sim->nz; ++i)
- fprintf(fp, "%.17g\t%.17g\t%.17g\t%.17g\t%.17g\n",
- sim->z[i],
- v_x_out[i],
- sim->sigma_n_eff[i],
- sim->mu[i],
- sim->gamma_dot_p[i]);
-
- free(v_x_out);
-}
-
-void
-print_wet_output(FILE* fp, struct simulation* sim, const int norm)
-{
- int i;
- double *v_x_out;
-
- if (norm)
- v_x_out = normalize(sim->v_x, sim->nz);
- else
- v_x_out = copy(sim->v_x, sim->nz);
-
- for (i=0; i<sim->nz; ++i)
- fprintf(fp, "%.17g\t%.17g\t%.17g\t%.17g\t%.17g\t%.17g\n",
- sim->z[i],
+ v_x_out = copy(sim.v_x, sim.nz);
+
+ for (i=0; i<sim.nz; ++i)
+ fprintf(fp,
+ "%.17g\t%.17g\t%.17g\t"
+ "%.17g\t%.17g\t%.17g\t"
+ "%.17g\t%.17g\t%.17g\n",
+ sim.z[i],
v_x_out[i],
- sim->sigma_n_eff[i],
- sim->p_f_ghost[i+1],
- sim->mu[i],
- sim->gamma_dot_p[i]);
+ sim.sigma_n_eff[i],
+ sim.p_f_ghost[i+1],
+ sim.mu[i],
+ sim.gamma_dot_p[i],
+ sim.phi[i],
+ sim.I[i],
+ sim.tau[i]);
free(v_x_out);
}
diff --git a/simulation.h b/simulation.h
@@ -6,6 +6,8 @@
#define PI 3.14159265358979323846
#define DEG2RAD(x) (x*PI/180.0)
+extern struct simulation sim;
+
/* Simulation settings */
struct simulation {
@@ -30,8 +32,7 @@ struct simulation {
/* stress ratio at top wall */
double mu_wall;
- /* nonlocal amplitude [-] */
- double A;
+ /* nonlocal amplitude [-] */ double A;
/* rate dependence beyond yield [-] */
double b;
@@ -78,6 +79,11 @@ struct simulation {
/* output file number */
int n_file;
+ double transient;
+ double phi_min;
+ double phi_max;
+ double dilatancy_angle;
+
/* Fluid parameters */
int fluid; /* flag to switch fluid on (1) or off (0) */
double p_f_top; /* fluid pressure at the top [Pa] */
@@ -101,15 +107,17 @@ struct simulation {
double* gamma_dot_p; /* plastic shear strain rate [1/s] */
double* v_x; /* shear velocity [m/s] */
double* g_ghost; /* fluidity with ghost nodes */
+ double* I; /* inertia number [-] */
+ double* tau; /* shear stress [Pa] */
};
-void prepare_arrays(struct simulation* sim);
-void free_arrays(struct simulation* sim);
+void prepare_arrays();
+void free_arrays();
-void check_simulation_parameters(const struct simulation* sim);
+void check_simulation_parameters();
-void lithostatic_pressure_distribution(struct simulation* sim);
+void lithostatic_pressure_distribution();
void set_bc_neumann(double* g_ghost,
const int nz,
@@ -122,18 +130,16 @@ void set_bc_dirichlet(double* g_ghost,
const int boundary,
const double value);
-void compute_cooperativity_length(struct simulation* sim);
-void compute_shear_strain_rate_plastic(struct simulation* sim);
-void compute_shear_velocity(struct simulation* sim);
-void compute_effective_stress(struct simulation* sim);
-void compute_friction(struct simulation* sim);
+void compute_cooperativity_length();
+void compute_shear_strain_rate_plastic();
+void compute_shear_velocity();
+void compute_effective_stress();
+void compute_friction();
-int implicit_1d_jacobian_poisson_solver(struct simulation* sim,
- const int max_iter,
+int implicit_1d_jacobian_poisson_solver(const int max_iter,
const double rel_tol);
-void write_output_file(struct simulation* sim, const int normalize);
-void print_dry_output(FILE* fp, struct simulation* sim, const int normalize);
-void print_wet_output(FILE* fp, struct simulation* sim, const int normalize);
+void write_output_file(const int normalize);
+void print_output(FILE* fp, const int normalize);
#endif
