tAdd new fluid-related parameters, save porosity as array - cngf-pf - continuum… | |
git clone git://src.adamsgaard.dk/cngf-pf | |
Log | |
Files | |
Refs | |
README | |
LICENSE | |
--- | |
commit cdf51f0472dbf466e9c6da72c90a42a77cbc80b0 | |
parent edc6affeb0cb2691c0db6ea8ef5043f3a43b1e3e | |
Author: Anders Damsgaard <[email protected]> | |
Date: Mon, 15 Apr 2019 13:06:23 +0200 | |
Add new fluid-related parameters, save porosity as array | |
Diffstat: | |
M 1d_fd_simple_shear.png | 0 | |
M 1d_fd_simple_shear_damsgaard2013.h | 10 +++++----- | |
M 1d_fd_simple_shear_henann_kamrin20… | 2 +- | |
M 1d_fd_simple_shear_rheology.png | 0 | |
M 1d_fd_simple_shear_rheology_iverso… | 0 | |
M 1d_fd_simple_shear_rheology_kamb.p… | 0 | |
M 1d_fd_simple_shear_rheology_tulacz… | 0 | |
M arrays.c | 15 ++++++++++++--- | |
M arrays.h | 7 ++++--- | |
M main.c | 15 ++++++++++----- | |
M simulation.c | 18 +++++++++++------- | |
M simulation.h | 11 +++++++---- | |
12 files changed, 50 insertions(+), 28 deletions(-) | |
--- | |
diff --git a/1d_fd_simple_shear.png b/1d_fd_simple_shear.png | |
Binary files differ. | |
diff --git a/1d_fd_simple_shear_damsgaard2013.h b/1d_fd_simple_shear_damsgaard2… | |
t@@ -30,9 +30,9 @@ struct simulation init_sim(void) | |
* above yield */ | |
sim.b = 0.9377; /* Henann and Kamrin 2016 */ | |
- sim.mu_s = atan(DEG2RAD(22.0)); /* Damsgaard et al 2013 */ | |
+ sim.mu_s = atan(DEG2RAD(22.0)); /* Damsgaard et al 2013 */ | |
- sim.phi = 0.25; /* 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 */ | |
t@@ -58,21 +58,21 @@ 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)*sim->rho_s*sim->G*(sim->L_z - sim->z[i]); | |
+ (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)*sim->rho_s*sim->G*(sim->L_z - sim->z[i])/ | |
+ (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_w_ghost[idx1g(i)] = 0.0; | |
+ sim->p_f_ghost[idx1g(i)] = 0.0; | |
} | |
#endif | |
diff --git a/1d_fd_simple_shear_henann_kamrin2016.h b/1d_fd_simple_shear_henann… | |
t@@ -23,7 +23,7 @@ struct simulation init_sim(void) | |
sim.A = 0.48; | |
sim.b = 0.9377; | |
sim.mu_s = 0.3819; | |
- sim.phi = 0.38; | |
+ sim.phi = initval(0.38, sim.nz); | |
sim.d = 1e-3; | |
sim.rho_s = 2.485e3; | |
diff --git a/1d_fd_simple_shear_rheology.png b/1d_fd_simple_shear_rheology.png | |
Binary files differ. | |
diff --git a/1d_fd_simple_shear_rheology_iverson.png b/1d_fd_simple_shear_rheol… | |
Binary files differ. | |
diff --git a/1d_fd_simple_shear_rheology_kamb.png b/1d_fd_simple_shear_rheology… | |
Binary files differ. | |
diff --git a/1d_fd_simple_shear_rheology_tulaczyk.png b/1d_fd_simple_shear_rheo… | |
Binary files differ. | |
diff --git a/arrays.c b/arrays.c | |
t@@ -54,7 +54,7 @@ double* linspace(const double lower, const double upper, con… | |
} | |
/* Return an array of `n` values with the value 0.0 */ | |
-double* zeros(const double n) | |
+double* zeros(const int n) | |
{ | |
double *x = malloc(n*sizeof(double)); | |
for (int i=0; i<n; ++i) | |
t@@ -63,7 +63,7 @@ double* zeros(const double n) | |
} | |
/* Return an array of `n` values with the value 1.0 */ | |
-double* ones(const double n) | |
+double* ones(const int n) | |
{ | |
double *x = malloc(n*sizeof(double)); | |
for (int i=0; i<n; ++i) | |
t@@ -71,8 +71,17 @@ double* ones(const double n) | |
return x; | |
} | |
+/* Return an array of `n` values with a specified value */ | |
+double* initval(const double value, const int n) | |
+{ | |
+ double *x = malloc(n*sizeof(double)); | |
+ for (int i=0; i<n; ++i) | |
+ x[i] = value; | |
+ return x; | |
+} | |
+ | |
/* Return an array of `n` uninitialized values */ | |
-double* empty(const double n) | |
+double* empty(const int n) | |
{ | |
return malloc(n*sizeof(double)); | |
} | |
diff --git a/arrays.h b/arrays.h | |
t@@ -18,9 +18,10 @@ unsigned int idx2g( | |
unsigned int idx1g(const unsigned int i); | |
double* linspace(const double lower, const double upper, const int n); | |
-double* zeros(const double n); | |
-double* ones(const double n); | |
-double* empty(const double n); | |
+double* zeros(const int n); | |
+double* ones(const int n); | |
+double* initval(const double value, const int n); | |
+double* empty(const int n); | |
double max(const double* a, const int n); | |
double min(const double* a, const int n); | |
diff --git a/main.c b/main.c | |
t@@ -52,7 +52,7 @@ int main(int argc, char* argv[]) | |
int normalize = 0; | |
int opt; | |
- const char* optstring = "hvNG:P:m:V:A:b:f:p:d:r:n:o:L:"; | |
+ const char* optstring = "hvNn:G:P:m:V:A:b:f:p:d:r:o:L:"; | |
const struct option longopts[] = { | |
{"help", no_argument, NULL, 'h'}, | |
{"version", no_argument, NULL, 'v'}, | |
t@@ -72,6 +72,7 @@ int main(int argc, char* argv[]) | |
{NULL, 0, NULL, 0} | |
}; | |
+ double new_phi = NAN; | |
while ((opt = getopt_long(argc, argv, optstring, longopts, NULL)) != -1) { | |
switch (opt) { | |
case -1: /* no more arguments */ | |
t@@ -84,6 +85,9 @@ int main(int argc, char* argv[]) | |
case 'v': | |
version(); | |
return 0; | |
+ case 'n': | |
+ sim.nz = atoi(optarg); | |
+ break; | |
case 'N': | |
normalize = 1; | |
break; | |
t@@ -109,7 +113,7 @@ int main(int argc, char* argv[]) | |
sim.mu_s = atof(optarg); | |
break; | |
case 'p': | |
- sim.phi = atof(optarg); | |
+ new_phi = atof(optarg); | |
break; | |
case 'd': | |
sim.d = atof(optarg); | |
t@@ -117,9 +121,6 @@ int main(int argc, char* argv[]) | |
case 'r': | |
sim.rho_s = atof(optarg); | |
break; | |
- case 'n': | |
- sim.nz = atoi(optarg); | |
- break; | |
case 'o': | |
sim.origo_z = atof(optarg); | |
break; | |
t@@ -137,6 +138,10 @@ int main(int argc, char* argv[]) | |
prepare_arrays(&sim); | |
+ if (!isnan(new_phi)) | |
+ for (int i=0; i<sim.nz; ++i) | |
+ sim.phi[i] = new_phi; | |
+ | |
init_normal_stress(&sim); | |
init_water_pressure(&sim); | |
compute_effective_stress(&sim); | |
diff --git a/simulation.c b/simulation.c | |
t@@ -12,7 +12,9 @@ void prepare_arrays(struct simulation* sim) | |
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->p_f_ghost = zeros(sim->nz+2); /* water pressure with ghost nodes */ | |
+ sim->k = zeros(sim->nz); /* hydraulic permeability */ | |
+ sim->phi = zeros(sim->nz); /* porosity */ | |
sim->xi = zeros(sim->nz); /* cooperativity length */ | |
sim->gamma_dot_p = zeros(sim->nz); /* shear velocity */ | |
sim->v_x = zeros(sim->nz); /* shear velocity */ | |
t@@ -25,7 +27,9 @@ void free_arrays(struct simulation* sim) | |
free(sim->mu); | |
free(sim->sigma_n_eff); | |
free(sim->sigma_n); | |
- free(sim->p_w_ghost); | |
+ free(sim->p_f_ghost); | |
+ free(sim->k); | |
+ free(sim->phi); | |
free(sim->xi); | |
free(sim->gamma_dot_p); | |
free(sim->v_x); | |
t@@ -92,10 +96,10 @@ void check_simulation_parameters(const struct simulation* … | |
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.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) | |
t@@ -172,7 +176,7 @@ 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_w_ghost[idx1g(i)]; | |
+ sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_ghost[idx1g(i)]; | |
} | |
double cooperativity_length( | |
diff --git a/simulation.h b/simulation.h | |
t@@ -28,9 +28,6 @@ struct simulation { | |
/* bulk and critical state static yield friction coefficient [-] */ | |
double mu_s; | |
- /* porosity [-] */ | |
- double phi; | |
- | |
/* representative grain size [m] */ | |
double d; | |
t@@ -67,11 +64,17 @@ 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 */ | |
double* mu; /* static yield friction [-] */ | |
double* sigma_n_eff; /* effective normal pressure [Pa] */ | |
double* sigma_n; /* normal stress [Pa] */ | |
- double* p_w_ghost; /* water pressure [Pa] */ | |
+ double* p_f_ghost; /* fluid pressure [Pa] */ | |
+ double* k; /* hydraulic permeability */ | |
+ double* phi; /* porosity [-] */ | |
double* xi; /* cooperativity length */ | |
double* gamma_dot_p; /* plastic shear strain rate [1/s] */ | |
double* v_x; /* shear velocity [m/s] */ |