tImplement pressure solver - cngf-pf - continuum model for granular flows with pore-pressure dynamics (renamed from 1d_fd_simple_shear)
HTML git clone git://src.adamsgaard.dk/cngf-pf
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
DIR commit a8d3c51683cd130ccfbaa614ffef059c60d673fb
DIR parent d6d5abd9c90f013b39ad99462f4bdf419562942a
HTML Author: Anders Damsgaard <anders@adamsgaard.dk>
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(-)
---
DIR 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->nz+2);
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 BC */
set_bc_neumann(sim->p_f_ghost, sim->nz, -1);
free(dp_f_expl);
free(dp_f_impl);
DIR 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,
DIR 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.\n"
"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 [Pa]\n"
+ " -q, --fluid-pressure-freq VAL frequency of pressure variations [s^-1]\n"
+ " -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 [s]\n"
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:D:I:";
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+2);
+ puts("");
+
+ puts(".. normal stress before iterations:"); print_array(sim.sigma_n, sim.nz);
+ 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);
DIR 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
DIR 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* sim)
}
}
+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)
DIR 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,