tAdd fluid checks, change default make target - cngf-pf - continuum model for g… | |
git clone git://src.adamsgaard.dk/cngf-pf | |
Log | |
Files | |
Refs | |
README | |
LICENSE | |
--- | |
commit 750bf283e3904793069d4850c5faffdb74b37df9 | |
parent cdf51f0472dbf466e9c6da72c90a42a77cbc80b0 | |
Author: Anders Damsgaard <[email protected]> | |
Date: Mon, 15 Apr 2019 14:27:51 +0200 | |
Add fluid checks, change default make target | |
Diffstat: | |
M .gitlab-ci.yml | 2 +- | |
D 1d_fd_simple_shear_damsgaard2013.h | 78 -----------------------------… | |
M Makefile | 6 +++++- | |
M main.c | 6 +++--- | |
A parameter_defaults.h | 81 ++++++++++++++++++++++++++++++ | |
M simulation.c | 40 +++++++++++++++++++++++++++++… | |
M simulation.h | 14 +++++++++----- | |
7 files changed, 137 insertions(+), 90 deletions(-) | |
--- | |
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml | |
t@@ -4,7 +4,7 @@ build-alpine: | |
before_script: | |
- apk --no-cache add build-base gnuplot bash | |
script: | |
- - make | |
+ - make plots | |
artifacts: | |
paths: | |
- 1d_fd_simple_shear.png | |
diff --git a/1d_fd_simple_shear_damsgaard2013.h b/1d_fd_simple_shear_damsgaard2… | |
t@@ -1,78 +0,0 @@ | |
-#ifndef ONED_FD_SIMPLE_SHEAR_ | |
-#define ONED_FD_SIMPLE_SHEAR_ | |
- | |
-#include <math.h> | |
-#include <stdio.h> | |
-#include "arrays.h" | |
-#include "simulation.h" | |
- | |
-#define PI 3.14159265358979323846 | |
-#define DEG2RAD(x) (x*PI/180.0) | |
- | |
-/* Simulation settings */ | |
-struct simulation init_sim(void) | |
-{ | |
- struct simulation sim; | |
- | |
- sim.G = 9.81; | |
- | |
- sim.P_wall = 120e3; /* larger normal stress deepens the shear depth */ | |
- sim.mu_wall = 0.40; | |
- sim.v_x_bot = 0.0; | |
- | |
- sim.nz = 100; | |
- | |
- /* lower values of A mean that the velocity curve can have sharper curves, | |
- * e.g. at the transition from μ ≈ μ_s */ | |
- sim.A = 0.40; /* Loose fit to Damsgaard et al 2013 */ | |
- | |
- /* lower values of b mean larger shear velocity for a given stress ratio | |
- * above yield */ | |
- sim.b = 0.9377; /* Henann and Kamrin 2016 */ | |
- | |
- sim.mu_s = atan(DEG2RAD(22.0)); /* Damsgaard et al 2013 */ | |
- | |
- sim.phi = initval(0.25, sim.nz); /* Damsgaard et al 2013 */ | |
- | |
- /* lower values of d mean that the shear velocity curve can have sharper | |
- * curves, e.g. at the transition from μ ≈ μ_s */ | |
- sim.d = 0.04; /* Damsgaard et al 2013 */ | |
- | |
- /* grain material density [kg/m^3] */ | |
- sim.rho_s = 2.6e3; /* Damsgaard et al 2013 */ | |
- | |
- /* Spatial settings */ | |
- sim.origo_z = 0.0; | |
- sim.L_z = 0.7; /* Damsgaard et al 2013 */ | |
- | |
- /* Temporal settings */ | |
- sim.t = 0.0; | |
- sim.t_end = 10.0; | |
- sim.file_dt = 0.1; | |
- sim.n_file = 0; | |
- | |
- return sim; | |
-} | |
- | |
-void init_normal_stress(struct simulation* sim) | |
-{ | |
- for (int 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 init_friction(struct simulation* sim) | |
-{ | |
- for (int 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); | |
-} | |
- | |
-void init_water_pressure(struct simulation* sim) | |
-{ | |
- for (int i=0; i<sim->nz; ++i) | |
- sim->p_f_ghost[idx1g(i)] = 0.0; | |
-} | |
- | |
-#endif | |
diff --git a/Makefile b/Makefile | |
t@@ -4,7 +4,11 @@ SRC=$(wildcard *.c) | |
OBJ=$(patsubst %.c,%.o,$(SRC)) | |
HDR=$(wildcard *.h) | |
-default: 1d_fd_simple_shear.png \ | |
+.PHONY: default | |
+default: 1d_fd_simple_shear | |
+ | |
+.PHONY: plots | |
+plots: 1d_fd_simple_shear.png \ | |
1d_fd_simple_shear_rheology.png \ | |
1d_fd_simple_shear_rheology_kamb.png \ | |
1d_fd_simple_shear_rheology_iverson.png \ | |
diff --git a/main.c b/main.c | |
t@@ -8,8 +8,7 @@ | |
#define VERSION "0.1" | |
#define PROGNAME "1d_fd_simple_shear" | |
-/* set default simulation parameter values */ | |
-#include "1d_fd_simple_shear_damsgaard2013.h" | |
+#include "parameter_defaults.h" | |
static void usage(void) | |
{ | |
t@@ -143,7 +142,8 @@ int main(int argc, char* argv[]) | |
sim.phi[i] = new_phi; | |
init_normal_stress(&sim); | |
- init_water_pressure(&sim); | |
+ if (sim.fluid) | |
+ init_water_pressure(&sim); | |
compute_effective_stress(&sim); | |
init_friction(&sim); | |
compute_cooperativity_length(&sim); | |
diff --git a/parameter_defaults.h b/parameter_defaults.h | |
t@@ -0,0 +1,81 @@ | |
+#ifndef ONED_FD_SIMPLE_SHEAR_ | |
+#define ONED_FD_SIMPLE_SHEAR_ | |
+ | |
+#include <math.h> | |
+#include <stdio.h> | |
+#include "arrays.h" | |
+#include "simulation.h" | |
+ | |
+#define PI 3.14159265358979323846 | |
+#define DEG2RAD(x) (x*PI/180.0) | |
+ | |
+/* Simulation settings */ | |
+struct simulation init_sim(void) | |
+{ | |
+ struct simulation sim; | |
+ | |
+ sim.G = 9.81; | |
+ | |
+ sim.P_wall = 120e3; /* larger normal stress deepens the shear depth */ | |
+ sim.mu_wall = 0.40; | |
+ sim.v_x_bot = 0.0; | |
+ | |
+ sim.nz = 100; | |
+ | |
+ /* lower values of A mean that the velocity curve can have sharper curves, | |
+ * e.g. at the transition from μ ≈ μ_s */ | |
+ sim.A = 0.40; /* Loose fit to Damsgaard et al 2013 */ | |
+ | |
+ /* lower values of b mean larger shear velocity for a given stress ratio | |
+ * above yield */ | |
+ sim.b = 0.9377; /* Henann and Kamrin 2016 */ | |
+ | |
+ sim.mu_s = atan(DEG2RAD(22.0)); /* Damsgaard et al 2013 */ | |
+ | |
+ sim.phi = initval(0.25, sim.nz); /* Damsgaard et al 2013 */ | |
+ | |
+ /* lower values of d mean that the shear velocity curve can have sharper | |
+ * curves, e.g. at the transition from μ ≈ μ_s */ | |
+ sim.d = 0.04; /* Damsgaard et al 2013 */ | |
+ | |
+ /* grain material density [kg/m^3] */ | |
+ sim.rho_s = 2.6e3; /* Damsgaard et al 2013 */ | |
+ | |
+ /* spatial settings */ | |
+ sim.origo_z = 0.0; | |
+ sim.L_z = 0.7; /* Damsgaard et al 2013 */ | |
+ | |
+ /* temporal settings */ | |
+ sim.t = 0.0; | |
+ sim.t_end = 10.0; | |
+ sim.file_dt = 0.1; | |
+ sim.n_file = 0; | |
+ | |
+ /* fluid settings */ | |
+ sim.fluid = 0; | |
+ | |
+ return sim; | |
+} | |
+ | |
+void init_normal_stress(struct simulation* sim) | |
+{ | |
+ for (int 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 init_friction(struct simulation* sim) | |
+{ | |
+ for (int 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); | |
+} | |
+ | |
+void init_water_pressure(struct simulation* sim) | |
+{ | |
+ for (int i=0; i<sim->nz; ++i) | |
+ sim->p_f_ghost[idx1g(i)] = 0.0; | |
+} | |
+ | |
+#endif | |
diff --git a/simulation.c b/simulation.c | |
t@@ -143,6 +143,38 @@ void check_simulation_parameters(const struct simulation*… | |
check_float("sim.file_dt", sim->file_dt, &return_status); | |
+ if (sim->fluid != 0 && sim->fluid != 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) | |
+ warn_parameter_value("sim.p_f_mod_ampl is not a zero or positive", | |
+ 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) | |
+ warn_parameter_value("sim.p_f_mod_freq is not a zero or positive", | |
+ sim->p_f_mod_freq, &return_status); | |
+ | |
+ 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); | |
+ | |
+ 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); | |
+ | |
+ 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); | |
+ | |
+ } | |
+ | |
if (return_status != 0) { | |
fprintf(stderr, "error: aborting due to invalid parameter choices\n"); | |
exit(return_status); | |
t@@ -175,8 +207,12 @@ void compute_shear_velocity(struct simulation* sim) | |
void compute_effective_stress(struct simulation* sim) | |
{ | |
- for (int i=0; i<sim->nz; ++i) | |
- sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_ghost[idx1g(i)]; | |
+ if (sim->fluid) | |
+ for (int i=0; i<sim->nz; ++i) | |
+ sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_ghost[idx1g(i)]; | |
+ else | |
+ for (int i=0; i<sim->nz; ++i) | |
+ sim->sigma_n_eff[i] = sim->sigma_n[i]; | |
} | |
double cooperativity_length( | |
diff --git a/simulation.h b/simulation.h | |
t@@ -64,11 +64,15 @@ struct simulation { | |
/* output file number */ | |
int n_file; | |
- /* Water-pressure perturbations */ | |
- double p_f_mod_ampl; /* Amplitude of water pressure variations [Pa] */ | |
- double p_f_mod_freq; /* Frequency of water pressure variations [s^-1] */ | |
- | |
- /* other arrays */ | |
+ /* Fluid parameters */ | |
+ int fluid; /* flag to switch fluid on (1) or off (0) */ | |
+ double p_f_mod_ampl; /* amplitude of water pressure variations [Pa] */ | |
+ double p_f_mod_freq; /* frequency of water pressure variations [s^-1] */ | |
+ double beta_f; /* adiabatic fluid compressibility */ | |
+ double mu_f; /* fluid dynamic viscosity */ | |
+ double rho_f; /* fluid density */ | |
+ | |
+ /* arrays */ | |
double* mu; /* static yield friction [-] */ | |
double* sigma_n_eff; /* effective normal pressure [Pa] */ | |
double* sigma_n; /* normal stress [Pa] */ |