tImplement pressure solver - cngf-pf - continuum model for granular flows with … | |
git clone git://src.adamsgaard.dk/cngf-pf | |
Log | |
Files | |
Refs | |
README | |
LICENSE | |
--- | |
commit a8d3c51683cd130ccfbaa614ffef059c60d673fb | |
parent d6d5abd9c90f013b39ad99462f4bdf419562942a | |
Author: Anders Damsgaard <[email protected]> | |
Date: Tue, 16 Apr 2019 12:50:15 +0200 | |
Implement pressure solver | |
Diffstat: | |
M fluid.c | 47 +++++++++++++++++++++++++----… | |
M fluid.h | 2 ++ | |
M main.c | 47 ++++++++++++++++++++++++++---… | |
M parameter_defaults.h | 23 +++++------------------ | |
M simulation.c | 13 +++++++++---- | |
M simulation.h | 14 ++++++++++---- | |
6 files changed, 105 insertions(+), 41 deletions(-) | |
--- | |
diff --git a/fluid.c b/fluid.c | |
t@@ -1,7 +1,25 @@ | |
#include <stdlib.h> | |
+#include <math.h> | |
#include "simulation.h" | |
#include "arrays.h" | |
+void hydrostatic_fluid_pressure_distribution(struct simulation* sim) | |
+{ | |
+ for (int i=0; i<sim->nz; ++i) | |
+ sim->p_f_ghost[idx1g(i)] = sim->p_f_top + | |
+ sim->phi[i]*sim->rho_f*sim->G*(sim->L_z - sim->z[i]); | |
+} | |
+ | |
+static double sine_wave( | |
+ const double time, | |
+ const double amplitude, | |
+ const double frequency, | |
+ const double phase, | |
+ const double base_value) | |
+{ | |
+ return amplitude*sin(2.0*PI*frequency*time + phase) + base_value; | |
+} | |
+ | |
static double darcy_pressure_change_1d( | |
const int i, | |
const int nz, | |
t@@ -62,7 +80,8 @@ int darcy_solver_1d( | |
* theta in ]0.0; 1.0]: underrelaxation | |
* theta = 1.0: Gauss-Seidel | |
* theta > 1.0: overrrelaxation */ | |
- const double theta = 0.05; | |
+ /* const double theta = 0.05; */ | |
+ const double theta = 1.7; | |
double p_f; | |
t@@ -72,15 +91,23 @@ int darcy_solver_1d( | |
double* p_f_ghost_out = zeros(sim->nz+2); | |
double* r_norm = zeros(sim->nz); | |
double r_norm_max = NAN; | |
- for (iter=0; iter<max_iter; ++iter) { | |
+ double p_f_top = sine_wave( | |
+ sim->t, | |
+ sim->p_f_mod_ampl, | |
+ sim->p_f_mod_freq, | |
+ sim->p_f_mod_phase, | |
+ sim->p_f_top); | |
- /* Dirichlet BCs resemble fixed particle velocities */ | |
- set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, 0.0); | |
- /* Neumann BCs resemble free surfaces */ | |
+ for (iter=0; iter<max_iter; ++iter) { | |
+ | |
+ set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top); | |
+ sim->p_f_ghost[idx1g(sim->nz-1)] = p_f_top; /* Include top node in BC … | |
set_bc_neumann(sim->p_f_ghost, sim->nz, -1); | |
+ puts(".. p_f_ghost after BC:"); print_array(sim->p_f_ghost, sim->nz+2); | |
- for (int i=0; i<sim->nz; ++i) | |
+ /* for (int i=0; i<sim->nz; ++i) */ | |
+ for (int i=0; i<sim->nz-1; ++i) | |
dp_f_impl[i] = darcy_pressure_change_1d( | |
i, | |
sim->nz, | |
t@@ -91,7 +118,8 @@ int darcy_solver_1d( | |
sim->dt, | |
sim->beta_f, | |
sim->mu_f); | |
- for (int i=0; i<sim->nz; ++i) { | |
+ /* for (int i=0; i<sim->nz; ++i) { */ | |
+ for (int i=0; i<sim->nz-1; ++i) { | |
p_f = sim->p_f_ghost[idx1g(i)]; | |
p_f_ghost_out[idx1g(i)] = p_f | |
t@@ -105,11 +133,14 @@ int darcy_solver_1d( | |
} | |
r_norm_max = max(r_norm, sim->nz); | |
+ puts(".. p_f_ghost_out:"); print_array(p_f_ghost_out, sim->nz+2); | |
copy_values(p_f_ghost_out, sim->p_f_ghost, sim->nz+2); | |
+ puts(".. p_f_ghost after update:"); print_array(sim->p_f_ghost, sim->n… | |
if (r_norm_max <= rel_tol) { | |
- set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, 0.0); | |
+ set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top); | |
+ sim->p_f_ghost[idx1g(sim->nz-1)] = p_f_top; /* Include top node in… | |
set_bc_neumann(sim->p_f_ghost, sim->nz, -1); | |
free(dp_f_expl); | |
free(dp_f_impl); | |
diff --git a/fluid.h b/fluid.h | |
t@@ -3,6 +3,8 @@ | |
#include "simulation.h" | |
+void hydrostatic_fluid_pressure_distribution(struct simulation* sim); | |
+ | |
int darcy_solver_1d( | |
struct simulation* sim, | |
const int max_iter, | |
diff --git a/main.c b/main.c | |
t@@ -13,11 +13,12 @@ | |
static void usage(void) | |
{ | |
- printf("%s: %s [OPTIONS]\n" | |
+ printf("%s: %s [OPTIONS] [NAME]\n" | |
+ "runs a simulation and outputs state in files prefixed with NAME.\… | |
"optional arguments:\n" | |
" -N, --normalize normalize output velocity\n" | |
" -G, --gravity VAL gravity magnitude [m/s^2]\n" | |
- " -P, --pressure VAL normal stress on top [Pa]\n" | |
+ " -P, --normal-stress 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@@ -33,6 +34,11 @@ static void usage(void) | |
" -c, --fluid-compressibility VAL fluid compressibility [Pa^-1]\n" | |
" -i, --fluid-viscosity VAL fluid viscosity [Pa*s]\n" | |
" -R, --fluid-density VAL fluid density [kg/m^3]\n" | |
+ " -k, --fluid-permeability VAL fluid permeability [m^2]\n" | |
+ " -O, --fluid-pressure-top VAL fluid pressure at +z edge [Pa]\n" | |
+ " -a, --fluid-pressure-ampl VAL amplitude of pressure variations… | |
+ " -q, --fluid-pressure-freq VAL frequency of pressure variations… | |
+ " -H, --fluid-pressure-phase VAL fluid pressure at +z edge [Pa]\n" | |
" -t, --time VAL simulation start time [s]\n" | |
" -T, --time-end VAL simulation end time [s]\n" | |
" -D, --time-step VAL computational time step length [… | |
t@@ -60,12 +66,12 @@ int main(int argc, char* argv[]) | |
int normalize = 0; | |
int opt; | |
- const char* optstring = "hvNn:G:P:m:V:A:b:f:Fp:d:r:o:L:c:i:R:k:t:T:D:I:"; | |
+ const char* optstring = "hvNn:G:P:m:V:A:b:f:Fp:d:r:o:L:c:i:R:k:O:a:q:H:t:T… | |
const struct option longopts[] = { | |
{"help", no_argument, NULL, 'h'}, | |
{"version", no_argument, NULL, 'v'}, | |
{"gravity", required_argument, NULL, 'G'}, | |
- {"pressure", required_argument, NULL, 'P'}, | |
+ {"normal-stress", required_argument, NULL, 'P'}, | |
{"stress-ratio", required_argument, NULL, 'm'}, | |
{"velocity-bottom", required_argument, NULL, 'V'}, | |
{"nonlocal-amplitude", required_argument, NULL, 'A'}, | |
t@@ -81,7 +87,11 @@ int main(int argc, char* argv[]) | |
{"fluid-compressiblity", required_argument, NULL, 'c'}, | |
{"fluid-viscosity", required_argument, NULL, 'i'}, | |
{"fluid-density", required_argument, NULL, 'R'}, | |
- {"fluid-permeability", 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'}, | |
{"time", required_argument, NULL, 't'}, | |
{"time-end", required_argument, NULL, 'T'}, | |
{"time-step", required_argument, NULL, 'D'}, | |
t@@ -160,6 +170,18 @@ int main(int argc, char* argv[]) | |
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 't': | |
sim.t = atof(optarg); | |
break; | |
t@@ -197,14 +219,25 @@ int main(int argc, char* argv[]) | |
for (int i=0; i<sim.nz; ++i) | |
sim.k[i] = new_k; | |
+ lithostatic_pressure_distribution(&sim); | |
+ | |
+ if (sim.fluid) | |
+ hydrostatic_fluid_pressure_distribution(&sim); | |
+ | |
+ puts(".. p_f_ghost before iterations:"); print_array(sim.p_f_ghost, sim.nz… | |
+ puts(""); | |
+ | |
+ puts(".. normal stress before iterations:"); print_array(sim.sigma_n, sim.… | |
+ puts(""); | |
+ | |
double filetimeclock = 0.0; | |
while (sim.t <= sim.t_end) { | |
- init_normal_stress(&sim); | |
if (sim.fluid) { | |
- init_water_pressure(&sim); | |
if (darcy_solver_1d(&sim, 10000, 1e-5)) | |
exit(1); | |
+ puts(".. p_f_ghost:"); print_array(sim.p_f_ghost, sim.nz+2); | |
+ puts(""); | |
} | |
compute_effective_stress(&sim); | |
diff --git a/parameter_defaults.h b/parameter_defaults.h | |
t@@ -6,9 +6,6 @@ | |
#include "arrays.h" | |
#include "simulation.h" | |
-#define PI 3.14159265358979323846 | |
-#define DEG2RAD(x) (x*PI/180.0) | |
- | |
/* Simulation settings */ | |
struct simulation init_sim(void) | |
{ | |
t@@ -62,19 +59,15 @@ struct simulation init_sim(void) | |
sim.rho_f = 1e3; /* Water */ | |
sim.k = initval(1.9e-15, 1); /* Damsgaard et al 2015 */ | |
- sim.p_f_mod_ampl = 0.0; /* no fluid-pressure variations */ | |
- sim.p_f_mod_freq = 0.0; /* no fluid-pressure variations */ | |
+ /* no fluid-pressure variations */ | |
+ sim.p_f_top = 0.0; | |
+ sim.p_f_mod_ampl = 0.0; | |
+ sim.p_f_mod_freq = 1.0; | |
+ sim.p_f_mod_phase = 0.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) | |
t@@ -83,10 +76,4 @@ 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_f_ghost[idx1g(i)] = 0.0; | |
-} | |
- | |
#endif | |
diff --git a/simulation.c b/simulation.c | |
t@@ -1,5 +1,6 @@ | |
#include <stdio.h> | |
#include <stdlib.h> | |
+#include <math.h> | |
#include "arrays.h" | |
#include "simulation.h" | |
t@@ -12,7 +13,7 @@ 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_f_ghost = zeros(sim->nz+2); /* water pressure with ghost nodes */ | |
+ 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); | |
t@@ -38,7 +39,6 @@ void free_arrays(struct simulation* sim) | |
free(sim->g_ghost); | |
} | |
- | |
static void warn_parameter_value( | |
const char message[], | |
const double value, | |
t@@ -49,7 +49,6 @@ static void warn_parameter_value( | |
*return_status = 1; | |
} | |
- | |
static void check_float( | |
const char name[], | |
const double value, | |
t@@ -64,7 +63,6 @@ static void check_float( | |
} | |
} | |
- | |
void check_simulation_parameters(const struct simulation* sim) | |
{ | |
int return_status = 0; | |
t@@ -196,6 +194,13 @@ void check_simulation_parameters(const struct simulation*… | |
} | |
} | |
+void lithostatic_pressure_distribution(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]); | |
+} | |
+ | |
double shear_strain_rate_plastic( | |
const double fluidity, | |
const double friction) | |
diff --git a/simulation.h b/simulation.h | |
t@@ -1,9 +1,11 @@ | |
#ifndef SIMULATION_ | |
#define SIMULATION_ | |
-#include <math.h> | |
#include "arrays.h" | |
+#define PI 3.14159265358979323846 | |
+#define DEG2RAD(x) (x*PI/180.0) | |
+ | |
/* Simulation settings */ | |
struct simulation { | |
t@@ -69,8 +71,10 @@ struct simulation { | |
/* 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 p_f_top; /* fluid pressure at the top [Pa] */ | |
+ double p_f_mod_ampl; /* amplitude of fluid pressure variations [Pa] */ | |
+ double p_f_mod_freq; /* frequency of fluid pressure variations [s^-1] */ | |
+ double p_f_mod_phase; /* phase of fluid pressure variations [s^-1] */ | |
double beta_f; /* adiabatic fluid compressibility [Pa^-1] */ | |
double mu_f; /* fluid dynamic viscosity [Pa*s] */ | |
double rho_f; /* fluid density [kg/m^3] */ | |
t@@ -80,7 +84,7 @@ struct simulation { | |
double* sigma_n_eff; /* effective normal pressure [Pa] */ | |
double* sigma_n; /* normal stress [Pa] */ | |
double* p_f_ghost; /* fluid pressure [Pa] */ | |
- double* k; /* hydraulic permeability */ | |
+ double* k; /* hydraulic permeability [m^2] */ | |
double* phi; /* porosity [-] */ | |
double* xi; /* cooperativity length */ | |
double* gamma_dot_p; /* plastic shear strain rate [1/s] */ | |
t@@ -94,6 +98,8 @@ void free_arrays(struct simulation* sim); | |
void check_simulation_parameters(const struct simulation* sim); | |
+void lithostatic_pressure_distribution(struct simulation* sim); | |
+ | |
void set_bc_neumann( | |
double* g_ghost, | |
const int nz, |