tAdd temporal parameters, parameter value check, and zeros for water pressure -… | |
git clone git://src.adamsgaard.dk/cngf-pf | |
Log | |
Files | |
Refs | |
README | |
LICENSE | |
--- | |
commit bfa41cc2e556c1fb56b794cac4505371890dc86b | |
parent be5a38658a041750de71b977e28948a03f5cc510 | |
Author: Anders Damsgaard <[email protected]> | |
Date: Mon, 15 Apr 2019 11:11:05 +0200 | |
Add temporal parameters, parameter value check, and zeros for water pressure | |
Diffstat: | |
M 1d_fd_simple_shear_damsgaard2013.h | 22 +++++++++++++++++----- | |
M main.c | 12 ++++++++---- | |
M simulation.c | 138 +++++++++++++++++++++++++++++… | |
M simulation.h | 21 ++++++++++++++++++++- | |
4 files changed, 175 insertions(+), 18 deletions(-) | |
--- | |
diff --git a/1d_fd_simple_shear_damsgaard2013.h b/1d_fd_simple_shear_damsgaard2… | |
t@@ -32,26 +32,32 @@ struct simulation init_sim(void) | |
sim.mu_s = atan(DEG2RAD(22.0)); /* Damsgaard et al 2013 */ | |
- sim.phi = 0.25; /* Damsgaard et al 2013 */ | |
+ sim.phi = 0.25; /* 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 */ | |
+ sim.d = 0.04; /* Damsgaard et al 2013 */ | |
/* grain material density [kg/m^3] */ | |
- sim.rho_s = 2.6e3; /* Damsgaard et al 2013 */ | |
+ 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_pressure(struct simulation* sim) | |
+void init_normal_stress(struct simulation* sim) | |
{ | |
for (int i=0; i<sim->nz; ++i) | |
- sim->p[i] = sim->P_wall + | |
+ sim->sigma_n[i] = sim->P_wall + | |
(1.0 - sim->phi)*sim->rho_s*sim->G*(sim->L_z - sim->z[i]); | |
} | |
t@@ -63,4 +69,10 @@ void init_friction(struct simulation* sim) | |
sim->P_wall); | |
} | |
+void init_water_pressure(struct simulation* sim) | |
+{ | |
+ for (int i=0; i<sim->nz; ++i) | |
+ sim->p_w_ghost[idx1g(i)] = 0.0; | |
+} | |
+ | |
#endif | |
diff --git a/main.c b/main.c | |
t@@ -17,7 +17,7 @@ static void usage(void) | |
"optional arguments:\n" | |
" -N, --normalize normalize output velocity\n" | |
" -G, --gravity VAL gravity magnitude [m/s^2]\n" | |
- " -P, --pressure VAL normal stress [Pa]\n" | |
+ " -P, --pressure VAL normal stress on top [Pa]\n" | |
" -m, --stress-ratio VAL applied stress ratio [-]\n" | |
" -V, --velocity-bottom VAL base velocity at bottom [m/s]\n" | |
" -A, --nonlocal-amplitude VAL amplitude of nonlocality [-]\n" | |
t@@ -137,13 +137,17 @@ int main(int argc, char* argv[]) | |
prepare_arrays(&sim); | |
- init_pressure(&sim); | |
+ init_normal_stress(&sim); | |
+ init_water_pressure(&sim); | |
+ compute_effective_stress(&sim); | |
init_friction(&sim); | |
compute_cooperativity_length(&sim); | |
+ check_simulation_parameters(&sim); | |
+ | |
#ifdef DEBUG | |
puts("\n## Before solver"); | |
- puts(".. p:"); print_array(sim.p, sim.nz); | |
+ puts(".. sigma_n_eff:"); print_array(sim.sigma_n_eff, sim.nz); | |
puts(".. mu:"); print_array(sim.mu, sim.nz); | |
#endif | |
t@@ -156,7 +160,7 @@ int main(int argc, char* argv[]) | |
if (normalize) | |
print_arrays_2nd_normalized(sim.z, sim.v_x, sim.nz); | |
else | |
- print_arrays(sim.z, sim.v_x, sim.nz); | |
+ print_three_arrays(sim.z, sim.v_x, sim.sigma_n_eff, sim.nz); | |
free_arrays(&sim); | |
return 0; | |
diff --git a/simulation.c b/simulation.c | |
t@@ -9,20 +9,22 @@ void prepare_arrays(struct simulation* sim) | |
sim->origo_z + sim->L_z, | |
sim->nz); | |
sim->dz = sim->z[1] - sim->z[0]; /* cell spacing */ | |
- sim->mu = zeros(sim->nz); /* local stress ratio */ | |
- sim->p = zeros(sim->nz); /* local pressure */ | |
+ 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_w_ghost = zeros(sim->nz+2); /* water pressure with ghost nodes */ | |
sim->xi = zeros(sim->nz); /* cooperativity length */ | |
- sim->gamma_dot_p = zeros(sim->nz); /* local shear velocity */ | |
- sim->v_x = zeros(sim->nz); /* local shear velocity */ | |
- sim->g_ghost = zeros(sim->nz+2); /* local fluidity with ghost nodes */ | |
+ 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 */ | |
} | |
void free_arrays(struct simulation* sim) | |
{ | |
free(sim->z); | |
free(sim->mu); | |
- free(sim->p); | |
+ free(sim->sigma_n_eff); | |
+ free(sim->sigma_n); | |
free(sim->p_w_ghost); | |
free(sim->xi); | |
free(sim->gamma_dot_p); | |
t@@ -30,6 +32,119 @@ void free_arrays(struct simulation* sim) | |
free(sim->g_ghost); | |
} | |
+ | |
+static void warn_parameter_value( | |
+ const char message[], | |
+ const double value, | |
+ int* return_status) | |
+{ | |
+ fprintf(stderr, "check_simulation_parameters: %s (%.17g)\n", | |
+ message, value); | |
+ *return_status = 1; | |
+} | |
+ | |
+ | |
+static void check_float( | |
+ const char name[], | |
+ const double value, | |
+ int* return_status) | |
+{ | |
+ if (isnan(value)) { | |
+ warn_parameter_value("%s is NaN", value, return_status); | |
+ *return_status = 1; | |
+ } else if (isinf(value)) { | |
+ warn_parameter_value("%s is infinite", value, return_status); | |
+ *return_status = 1; | |
+ } | |
+} | |
+ | |
+ | |
+void check_simulation_parameters(const struct simulation* sim) | |
+{ | |
+ int 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.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.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.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, | |
+ &return_status); | |
+ | |
+ check_float("sim.phi", sim->phi, &return_status); | |
+ if (sim->phi <= 0.0 || sim->phi > 1.0) | |
+ warn_parameter_value("sim.phi is not in ]0;1]", sim->phi, | |
+ &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, | |
+ &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, | |
+ &return_status); | |
+ | |
+ 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) | |
+ warn_parameter_value("sim.L_z is smaller or equal to sim.origo_z", | |
+ sim->L_z, &return_status); | |
+ | |
+ 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) | |
+ warn_parameter_value("sim.dz is not a positive number", sim->dz, | |
+ &return_status); | |
+ | |
+ check_float("sim.t", sim->t, &return_status); | |
+ 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); | |
+ | |
+ check_float("sim.dt", sim->t_end, &return_status); | |
+ if (sim->dt <= 0.0) | |
+ warn_parameter_value("sim.dt is not a positive number", | |
+ sim->dt, &return_status); | |
+ | |
+ check_float("sim.file_dt", sim->file_dt, &return_status); | |
+ | |
+ if (return_status != 0) { | |
+ fprintf(stderr, "error: aborting due to invalid parameter choices\n"); | |
+ exit(return_status); | |
+ } | |
+} | |
+ | |
double shear_strain_rate_plastic( | |
const double fluidity, | |
const double friction) | |
t@@ -54,6 +169,12 @@ void compute_shear_velocity(struct simulation* sim) | |
sim->v_x[i] = sim->v_x[i-1] + sim->gamma_dot_p[i]*sim->dz; | |
} | |
+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_w_ghost[idx1g(i)]; | |
+} | |
+ | |
double cooperativity_length( | |
const double A, | |
const double d, | |
t@@ -88,7 +209,8 @@ void compute_local_fluidity(struct simulation* sim) | |
{ | |
for (int i=0; i<sim->nz; ++i) | |
sim->g_ghost[idx1g(i)] = local_fluidity( | |
- sim->p[i], sim->mu[i], sim->mu_s, sim->b, sim->rho_s, sim->d); | |
+ sim->sigma_n_eff[i], sim->mu[i], sim->mu_s, sim->b, sim->rho_s, | |
+ sim->d); | |
} | |
void set_bc_neumann(double* g_ghost, const int nz, const int boundary) | |
t@@ -182,7 +304,7 @@ int implicit_1d_jacobian_poisson_solver( | |
r_norm, | |
sim->dz, | |
sim->mu, | |
- sim->p, | |
+ sim->sigma_n_eff, | |
sim->xi, | |
sim->mu_s, | |
sim->b, | |
diff --git a/simulation.h b/simulation.h | |
t@@ -52,9 +52,25 @@ struct simulation { | |
/* cell spacing [m] */ | |
double dz; | |
+ /* current time [s] */ | |
+ double t; | |
+ | |
+ /* end time [s] */ | |
+ double t_end; | |
+ | |
+ /* time step length [s] */ | |
+ double dt; | |
+ | |
+ /* interval between output files [s] */ | |
+ double file_dt; | |
+ | |
+ /* output file number */ | |
+ int n_file; | |
+ | |
/* other arrays */ | |
double* mu; /* static yield friction [-] */ | |
- double* p; /* effective normal pressure [Pa] */ | |
+ double* sigma_n_eff; /* effective normal pressure [Pa] */ | |
+ double* sigma_n; /* normal stress [Pa] */ | |
double* p_w_ghost; /* water pressure [Pa] */ | |
double* xi; /* cooperativity length */ | |
double* gamma_dot_p; /* plastic shear strain rate [1/s] */ | |
t@@ -66,6 +82,8 @@ struct simulation { | |
void prepare_arrays(struct simulation* sim); | |
void free_arrays(struct simulation* sim); | |
+void check_simulation_parameters(const struct simulation* sim); | |
+ | |
void set_bc_neumann( | |
double* g_ghost, | |
const int nz, | |
t@@ -80,6 +98,7 @@ void set_bc_dirichlet( | |
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); | |
int implicit_1d_jacobian_poisson_solver( | |
struct simulation* sim, |