tOutline fluid timestep equation, use SL syntax rules for fluid.c - cngf-pf - c… | |
git clone git://src.adamsgaard.dk/cngf-pf | |
Log | |
Files | |
Refs | |
README | |
LICENSE | |
--- | |
commit dd88e99f42e4aa6e73ae7731af23d62d1a031c0e | |
parent b2f8d047df70db4aa5593bfafe9bd594312e1481 | |
Author: Anders Damsgaard <[email protected]> | |
Date: Fri, 28 Jun 2019 09:41:08 +0200 | |
Outline fluid timestep equation, use SL syntax rules for fluid.c | |
Diffstat: | |
D 1d_fd_simple_shear.png | 0 | |
D 1d_fd_simple_shear_rheology.png | 0 | |
D 1d_fd_simple_shear_rheology_iverso… | 0 | |
D 1d_fd_simple_shear_rheology_kamb.p… | 0 | |
D 1d_fd_simple_shear_rheology_tulacz… | 0 | |
M fluid.c | 113 ++++++++++++++++++-----------… | |
M fluid.h | 2 ++ | |
M main.c | 22 +++++++++------------- | |
8 files changed, 78 insertions(+), 59 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_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/fluid.c b/fluid.c | |
t@@ -3,39 +3,63 @@ | |
#include "simulation.h" | |
#include "arrays.h" | |
-void hydrostatic_fluid_pressure_distribution(struct simulation* sim) | |
+void | |
+hydrostatic_fluid_pressure_distribution(struct simulation* sim) | |
{ | |
- for (int i=0; i<sim->nz; ++i) | |
+ int i; | |
+ for (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) | |
+/* Determines the largest time step for the current simulation state. Note | |
+ * that the time step should be recalculated if cell sizes or diffusivities | |
+ * (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) | |
{ | |
- return amplitude*sin(2.0*PI*frequency*time + phase) + base_value; | |
+ double dx_min, diff_max; | |
+ | |
+ /* determine smallest cell size */ | |
+ dx_min = 0.0; | |
+ | |
+ /* determine largest diffusion size */ | |
+ diff_max = 0.0; | |
+ | |
+ sim->dt = safety*0.5*dx_min*dx_min/diff_max; | |
} | |
-static double darcy_pressure_change_1d(const int i, | |
- const int nz, | |
- const double* p_f_ghost_in, | |
- const double* phi, | |
- const double* k, | |
- const double dz, | |
- const double dt, | |
- const double beta_f, | |
- const double mu_f) | |
+static | |
+double sine_wave(const double time, | |
+ const double amplitude, | |
+ const double frequency, | |
+ const double phase, | |
+ const double base_value) | |
{ | |
- const double p = p_f_ghost_in[idx1g(i)]; | |
- const double p_zn = p_f_ghost_in[idx1g(i-1)]; | |
- const double p_zp = p_f_ghost_in[idx1g(i+1)]; | |
+ return amplitude*sin(2.0*PI*frequency*time + phase) + base_value; | |
+} | |
- const double k_ = k[i]; | |
+static | |
+double darcy_pressure_change_1d(const int i, | |
+ const int nz, | |
+ const double* p_f_ghost_in, | |
+ const double* phi, | |
+ const double* k, | |
+ const double dz, | |
+ const double dt, | |
+ const double beta_f, | |
+ const double mu_f) | |
+{ | |
+ double p, p_zn, p_zp, k_, div_k_grad_p; | |
double k_zn, k_zp; | |
+ | |
+ p = p_f_ghost_in[idx1g(i)]; | |
+ p_zn = p_f_ghost_in[idx1g(i-1)]; | |
+ p_zp = p_f_ghost_in[idx1g(i+1)]; | |
+ | |
+ k_ = k[i]; | |
if (i==0) k_zn = k_; else k_zn = k[i-1]; | |
if (i==nz-1) k_zp = k_; else k_zp = k[i+1]; | |
#ifdef DEBUG | |
t@@ -45,10 +69,9 @@ static double darcy_pressure_change_1d(const int i, | |
k_zn, k_, k_zp); | |
#endif | |
- const double div_k_grad_p = | |
- (2.0*k_zp*k_/(k_zp + k_) * (p_zp - p)/dz - | |
- 2.0*k_zn*k_/(k_zn + k_) * (p - p_zn)/dz | |
- )/dz; | |
+ div_k_grad_p = (2.0*k_zp*k_/(k_zp + k_) * (p_zp - p)/dz - | |
+ 2.0*k_zn*k_/(k_zn + k_) * (p - p_zn)/dz | |
+ )/dz; | |
#ifdef DEBUG | |
printf("phi[%d]=%g\tdiv_k_grad_p[%d]=%g\n", | |
i, phi[i], i, div_k_grad_p); | |
t@@ -62,10 +85,13 @@ int darcy_solver_1d(struct simulation* sim, | |
const int max_iter, | |
const double rel_tol) | |
{ | |
+ int i, iter; | |
+ double epsilon, theta, p_f, p_f_top, r_norm_max; | |
+ double* dp_f_expl, dp_f_impl, p_f_ghost_out, r_norm; | |
/* compute explicit solution to pressure change */ | |
- double* dp_f_expl = zeros(sim->nz); | |
- for (int i=0; i<sim->nz; ++i) | |
+ dp_f_expl = zeros(sim->nz); | |
+ for (i=0; i<sim->nz; ++i) | |
dp_f_expl[i] = darcy_pressure_change_1d(i, | |
sim->nz, | |
sim->p_f_ghost, | |
t@@ -80,29 +106,25 @@ int darcy_solver_1d(struct simulation* sim, | |
* epsilon = 0.0: explicit | |
* epsilon = 0.5: Crank-Nicolson | |
* epsilon = 1.0: implicit */ | |
- const double epsilon = 0.5; | |
+ epsilon = 0.5; | |
/* choose relaxation factor, parameter in ]0.0; 1.0] | |
* theta in ]0.0; 1.0]: underrelaxation | |
* theta = 1.0: Gauss-Seidel | |
* theta > 1.0: overrelaxation */ | |
- const double theta = 0.05; | |
- /* const double theta = 1.7; */ | |
- | |
- double p_f; | |
+ theta = 0.05; | |
+ /* theta = 1.7; */ | |
/* compute implicit solution to pressure change */ | |
- int iter; | |
- double* dp_f_impl = zeros(sim->nz); | |
- double* p_f_ghost_out = zeros(sim->nz+2); | |
- double* r_norm = zeros(sim->nz); | |
- double r_norm_max = NAN; | |
- 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); | |
+ dp_f_impl = zeros(sim->nz); | |
+ p_f_ghost_out = zeros(sim->nz+2); | |
+ r_norm = zeros(sim->nz); | |
+ r_norm_max = NAN; | |
+ 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); | |
for (iter=0; iter<max_iter; ++iter) { | |
t@@ -114,7 +136,7 @@ int darcy_solver_1d(struct simulation* sim, | |
#endif | |
/* for (int i=0; i<sim->nz; ++i) */ | |
- for (int i=0; i<sim->nz-1; ++i) | |
+ for (i=0; i<sim->nz-1; ++i) | |
dp_f_impl[i] = darcy_pressure_change_1d(i, | |
sim->nz, | |
sim->p_f_ghost, | |
t@@ -124,8 +146,7 @@ int darcy_solver_1d(struct simulation* sim, | |
sim->dt, | |
sim->beta_f, | |
sim->mu_f); | |
- /* for (int i=0; i<sim->nz; ++i) { */ | |
- for (int i=0; i<sim->nz-1; ++i) { | |
+ 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_expl[i], i, dp_f_impl[i]); | |
diff --git a/fluid.h b/fluid.h | |
t@@ -5,6 +5,8 @@ | |
void hydrostatic_fluid_pressure_distribution(struct simulation* sim); | |
+void set_largest_fluid_timestep(struct simulation* sim, const double safety); | |
+ | |
int darcy_solver_1d( | |
struct simulation* sim, | |
const int max_iter, | |
diff --git a/main.c b/main.c | |
t@@ -6,7 +6,7 @@ | |
#include "simulation.h" | |
#include "fluid.h" | |
-#define VERSION "0.2" | |
+#define VERSION "0.3.0" | |
#define PROGNAME "1d_fd_simple_shear" | |
#include "parameter_defaults.h" | |
t@@ -17,7 +17,9 @@ static void usage(void) | |
printf("%s: %s [OPTIONS] [NAME]\n" | |
"runs a simulation and outputs state in files prefixed with NA… | |
"If NAME is not specified, the default value '%s' is used.\n" | |
- "optional arguments:\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] (d… | |
" -P, --normal-stress VAL normal stress on top [Pa] (d… | |
t@@ -32,6 +34,7 @@ static void usage(void) | |
" -n, --resolution VAL number of cells in domain [-… | |
" -o, --origo VAL coordinate system origo [m] … | |
" -L, --length VAL domain length [m] (default %… | |
+ "\nOptional arguments only relevant with transient (fluid) sim… | |
" -F, --fluid enable pore fluid computatio… | |
" -c, --fluid-compressibility VAL fluid compressibility [Pa^-1… | |
" -i, --fluid-viscosity VAL fluid viscosity [Pa*s] (defa… | |
t@@ -43,10 +46,7 @@ static void usage(void) | |
" -H, --fluid-pressure-phase VAL fluid pressure at +z edge [P… | |
" -t, --time VAL simulation start time [s] (d… | |
" -T, --time-end VAL simulation end time [s] (def… | |
- " -D, --time-step VAL computational time step leng… | |
- " -I, --file-interval VAL interval between output file… | |
- " -v, --version show version information\n" | |
- " -h, --help show this message\n", | |
+ " -I, --file-interval VAL interval between output file… | |
__func__, PROGNAME, | |
sim.name, | |
sim.G, | |
t@@ -72,7 +72,6 @@ static void usage(void) | |
sim.p_f_mod_phase, | |
sim.t, | |
sim.t_end, | |
- sim.dt, | |
sim.file_dt); | |
free(sim.phi); | |
free(sim.k); | |
t@@ -125,7 +124,6 @@ int main(int argc, char* argv[]) | |
{"fluid-pressure-phase", required_argument, NULL, 'H'}, | |
{"time", required_argument, NULL, 't'}, | |
{"time-end", required_argument, NULL, 'T'}, | |
- {"time-step", required_argument, NULL, 'D'}, | |
{"file-interval", required_argument, NULL, 'I'}, | |
{NULL, 0, NULL, 0} | |
}; | |
t@@ -223,13 +221,9 @@ int main(int argc, char* argv[]) | |
case 'T': | |
sim.t_end = atof(optarg); | |
break; | |
- case 'D': | |
- sim.dt = atof(optarg); | |
- break; | |
case 'I': | |
sim.file_dt = atof(optarg); | |
break; | |
- | |
default: | |
fprintf(stderr, "%s: invalid option -- %c\n", … | |
fprintf(stderr, "Try `%s --help` for more info… | |
t@@ -257,8 +251,10 @@ int main(int argc, char* argv[]) | |
lithostatic_pressure_distribution(&sim); | |
- if (sim.fluid) | |
+ if (sim.fluid) { | |
hydrostatic_fluid_pressure_distribution(&sim); | |
+ set_largest_fluid_timestep(&sim, 0.5); | |
+ } | |
#ifdef DEBUG | |
puts(".. p_f_ghost before iterations:"); print_array(sim.p_f_ghost, si… |