commit 76f099629e26cfd4d8338b5862856ebc974d4e74
parent 9e079cfadd9105a19d6d651546fd09fe069612e3
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 29 Apr 2019 15:02:51 +0000
Merge branch 'darcy' into 'master'
Darcian diffusion of porewater
See merge request admesg/1d_fd_simple_shear!1
Diffstat:
23 files changed, 878 insertions(+), 412 deletions(-)
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
@@ -4,7 +4,7 @@ build-alpine:
before_script:
- apk --no-cache add build-base gnuplot bash
script:
- - make
+ - make plots
artifacts:
paths:
- 1d_fd_simple_shear.png
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_damsgaard2013.h
@@ -1,66 +0,0 @@
-#ifndef ONED_FD_SIMPLE_SHEAR_
-#define ONED_FD_SIMPLE_SHEAR_
-
-#include <math.h>
-#include <stdio.h>
-#include "arrays.h"
-#include "simulation.h"
-
-#define PI 3.14159265358979323846
-#define DEG2RAD(x) (x*PI/180.0)
-
-/* Simulation settings */
-struct simulation init_sim(void)
-{
- struct simulation sim;
-
- sim.G = 9.81;
-
- sim.P_wall = 120e3; /* larger normal stress deepens the shear depth */
- sim.mu_wall = 0.40;
- sim.v_x_bot = 0.0;
-
- sim.nz = 100;
-
- /* lower values of A mean that the velocity curve can have sharper curves,
- * e.g. at the transition from μ ≈ μ_s */
- sim.A = 0.40; /* Loose fit to Damsgaard et al 2013 */
-
- /* lower values of b mean larger shear velocity for a given stress ratio
- * above yield */
- sim.b = 0.9377; /* Henann and Kamrin 2016 */
-
- sim.mu_s = atan(DEG2RAD(22.0)); /* 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 */
-
- /* grain material density [kg/m^3] */
- 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 */
-
- return sim;
-}
-
-void init_pressure(struct simulation* sim)
-{
- for (int i=0; i<sim->nz; ++i)
- sim->p[i] = sim->P_wall +
- (1.0 - sim->phi)*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])/
- sim->P_wall);
-}
-
-#endif
diff --git a/1d_fd_simple_shear_fluid.gp b/1d_fd_simple_shear_fluid.gp
@@ -0,0 +1,35 @@
+#!/usr/bin/env gnuplot
+
+# specify input file with:
+# gnuplot -e "filename='file.txt'" 1d_fd_simple_shear_fluid.gp
+
+set terminal pngcairo color size 60 cm, 17.6 cm
+set multiplot layout 1,3
+# set multiplot layout 1,4
+
+
+set yrange [0.0:2.0]
+
+set key bottom right #samplen 0.9
+
+set xlabel "Water pressure, p_f [Pa]"
+set ylabel "Vertical position, z [m]" offset 2
+set xrange [p_min:p_max]
+plot filename u 4:1 w l lw 2 lc "blue" t ""
+
+set xlabel "Effective normal stress [Pa]"
+set ylabel ""
+set xrange [0:200e3]
+plot filename u 3:1 w l lw 2 lc "black" t ""
+
+# set xlabel "Friction, mu [-]"
+# set xrange [-0.05:1.5]
+# plot filename u 5:1 w l lw 2 lc "gray" t ""
+
+set xlabel "Horizontal velocity, v_x [m/s]"
+set ylabel ""
+set xrange [0.0:0.025]
+plot filename u 2:1 w l lw 2 lc "red" t ""
+
+
+unset multiplot
diff --git a/1d_fd_simple_shear_henann_kamrin2016.h b/1d_fd_simple_shear_henann_kamrin2016.h
@@ -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_rheology_iverson.png
Binary files differ.
diff --git a/1d_fd_simple_shear_rheology_kamb.png b/1d_fd_simple_shear_rheology_kamb.png
Binary files differ.
diff --git a/1d_fd_simple_shear_rheology_tulaczyk.png b/1d_fd_simple_shear_rheology_tulaczyk.png
Binary files differ.
diff --git a/Makefile b/Makefile
@@ -4,12 +4,29 @@ SRC=$(wildcard *.c)
OBJ=$(patsubst %.c,%.o,$(SRC))
HDR=$(wildcard *.h)
-default: 1d_fd_simple_shear.png \
+.PHONY: default
+default: 1d_fd_simple_shear
+
+.PHONY: plots
+plots: 1d_fd_simple_shear.png \
1d_fd_simple_shear_rheology.png \
1d_fd_simple_shear_rheology_kamb.png \
1d_fd_simple_shear_rheology_iverson.png \
1d_fd_simple_shear_rheology_tulaczyk.png
+diurnal.mp4: 1d_fd_simple_shear 1d_fd_simple_shear_fluid.gp
+ /usr/bin/env zsh -c '\
+ ./$< --resolution 50 --length 2.0 --normal-stress 150e3 --fluid --fluid-permeability 1e-17 --fluid-pressure-top 50e3 --fluid-pressure-ampl 50e3 --fluid-pressure-freq $$(( 1.0/(3600*24) )) --time-step 1e-1 --file-interval $$(( 60*10 )) --time-end $$(( 3600*24*4 )) diurnal'
+ /bin/bash -c '\
+ for f in diurnal.output*.txt; do \
+ gnuplot -e "filename=\"$$f\"; p_min=\"0\"; p_max=\"100e3\"" $<_fluid.gp > $$f.png; \
+ done'
+ ffmpeg -i diurnal.output%05d.txt.png -y diurnal.mp4
+
+diurnal.gif: diurnal.mp4
+ convert diurnal.output*.txt.png \
+ -delay 1 -loop 0 -fuzz 10% -layers Optimize $@
+
1d_fd_simple_shear: $(OBJ) $(HDR)
$(CC) $(LDFLAGS) $(OBJ) -o $@
diff --git a/README.md b/README.md
@@ -12,6 +12,44 @@ directory and an output PNG figure is generated.
Alternatively, an implementation in [Julia](https://julialang.org) resides in
**julia/**.
+### Advanced usage
+The majority of simulation parameters can be adjusted from the command line:
+
+```
+usage: 1d_fd_simple_shear [OPTIONS] [NAME]
+runs a simulation and outputs state in files prefixed with NAME.
+optional arguments:
+ -N, --normalize normalize output velocity
+ -G, --gravity VAL gravity magnitude [m/s^2]
+ -P, --normal-stress VAL normal stress on top [Pa]
+ -m, --stress-ratio VAL applied stress ratio [-]
+ -V, --velocity-bottom VAL base velocity at bottom [m/s]
+ -A, --nonlocal-amplitude VAL amplitude of nonlocality [-]
+ -b, --rate-dependence VAL rate dependence beyond yield [-]
+ -f, --friction-coefficient VAL grain friction coefficient [-]
+ -p, --porosity VAL porosity fraction [-]
+ -d, --grain-size VAL representative grain size [m]
+ -r, --density VAL grain material density [kg/m^3]
+ -n, --resolution VAL number of cells in domain [-]
+ -o, --origo VAL coordinate system origo [m]
+ -L, --length VAL domain length [m]
+ -F, --fluid enable pore fluid computations
+ -c, --fluid-compressibility VAL fluid compressibility [Pa^-1]
+ -i, --fluid-viscosity VAL fluid viscosity [Pa*s]
+ -R, --fluid-density VAL fluid density [kg/m^3]
+ -k, --fluid-permeability VAL fluid intrinsic permeability [m^2]
+ -O, --fluid-pressure-top VAL fluid pressure at +z edge [Pa]
+ -a, --fluid-pressure-ampl VAL amplitude of pressure variations [Pa]
+ -q, --fluid-pressure-freq VAL frequency of pressure variations [s^-1]
+ -H, --fluid-pressure-phase VAL fluid pressure at +z edge [Pa]
+ -t, --time VAL simulation start time [s]
+ -T, --time-end VAL simulation end time [s]
+ -D, --time-step VAL computational time step length [s]
+ -I, --file-interval VAL interval between output files [s]
+ -v, --version show version information
+ -h, --help show this message
+```
+
## Results
### Strain distribution
@@ -26,7 +64,6 @@ probably be improved further.
|  |  |
### Stress and strain rate
-
The rheology is of Bingham type, where no deformation occurs beneath the
Mohr-Coulomb yield limit. Above it, deformation is highly non-linearly viscous.
The model has a parameter *b* for rate dependence beyond yield. Glass beads
@@ -41,3 +78,10 @@ have *b* = 0.94 ([Forterre and Pouliquen
|  |  |
| Whillans Ice Plain ([Tulaczyk 2006](https://doi.org/10.3189/172756506781828601)): | |
|  |  |
+
+### Variable water pressure
+The model is expanded from the Henann and Kamrin 2013 model by including a
+diffusive porewater pressure parameterization. Below is an example of diurnal
+water-pressure variations that gradually propagate into the bed.
+
+
diff --git a/arrays.c b/arrays.c
@@ -54,7 +54,7 @@ double* linspace(const double lower, const double upper, const int n)
}
/* 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)
@@ -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)
@@ -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));
}
@@ -116,6 +125,37 @@ void print_arrays_2nd_normalized(const double* a, const double* b, const int n)
printf("%.17g\t%.17g\n", a[i], b[i]/max_b);
}
+void print_three_arrays(
+ const double* a,
+ const double* b,
+ const double* c,
+ const int n)
+{
+ for (int i=0; i<n; ++i)
+ printf("%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]);
+}
+
+void fprint_arrays(
+ FILE* fp,
+ const double* a,
+ const double* b,
+ const int n)
+{
+ for (int i=0; i<n; ++i)
+ fprintf(fp, "%.17g\t%.17g\n", a[i], b[i]);
+}
+
+void fprint_three_arrays(
+ FILE* fp,
+ const double* a,
+ const double* b,
+ const double* c,
+ const int n)
+{
+ for (int i=0; i<n; ++i)
+ fprintf(fp, "%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]);
+}
+
void copy_values(const double* in, double* out, const int n)
{
for (int i=0; i<n; ++i)
diff --git a/arrays.h b/arrays.h
@@ -1,3 +1,5 @@
+#include <stdio.h>
+
#ifndef ARRAYS_
#define ARRAYS_
@@ -18,9 +20,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);
@@ -28,6 +31,20 @@ double min(const double* a, const int n);
void print_array(const double* a, const int n);
void print_arrays(const double* a, const double* b, const int n);
void print_arrays_2nd_normalized(const double* a, const double* b, const int n);
+void print_three_arrays(
+ const double* a,
+ const double* b,
+ const double* c,
+ const int n);
+
+void fprint_arrays(FILE* fp, const double* a, const double* b, const int n);
+
+void fprint_three_arrays(
+ FILE* fp,
+ const double* a,
+ const double* b,
+ const double* c,
+ const int n);
void copy_values(const double* in, double* out, const int n);
diff --git a/diurnal.gif b/diurnal.gif
Binary files differ.
diff --git a/fluid.c b/fluid.c
@@ -0,0 +1,185 @@
+#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,
+ 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)
+{
+ 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)];
+
+ const double k_ = k[i];
+ double k_zn, k_zp;
+ 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
+ printf("%d->%d: p=[%g, %g, %g]\tk=[%g, %g, %g]\n",
+ i, idx1g(i),
+ p_zn, p, p_zp,
+ 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;
+#ifdef DEBUG
+ printf("phi[%d]=%g\tdiv_k_grad_p[%d]=%g\n",
+ i, phi[i], i, div_k_grad_p);
+#endif
+
+ /* return delta p */
+ return dt/(beta_f*phi[i]*mu_f)*div_k_grad_p;
+}
+
+int darcy_solver_1d(
+ struct simulation* sim,
+ const int max_iter,
+ const double rel_tol)
+{
+
+ /* compute explicit solution to pressure change */
+ double* dp_f_expl = zeros(sim->nz);
+ for (int i=0; i<sim->nz; ++i)
+ dp_f_expl[i] = darcy_pressure_change_1d(
+ i,
+ sim->nz,
+ sim->p_f_ghost,
+ sim->phi,
+ sim->k,
+ sim->dz,
+ sim->dt,
+ sim->beta_f,
+ sim->mu_f);
+
+ /* choose integration method, parameter in [0.0; 1.0]
+ * epsilon = 0.0: explicit
+ * epsilon = 0.5: Crank-Nicolson
+ * epsilon = 1.0: implicit */
+ const double 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: overrrelaxation */
+ const double theta = 0.05;
+ /* const double theta = 1.7; */
+
+ double p_f;
+
+ /* 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);
+
+
+ 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);
+#ifdef DEBUG
+ puts(".. p_f_ghost after BC:"); print_array(sim->p_f_ghost, sim->nz+2);
+#endif
+
+ /* 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,
+ sim->p_f_ghost,
+ sim->phi,
+ sim->k,
+ sim->dz,
+ 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) {
+#ifdef DEBUG
+ printf("dp_f_expl[%d] = %g\ndp_f_impl[%d] = %g\n",
+ i, dp_f_expl[i], i, dp_f_impl[i]);
+#endif
+
+ p_f = sim->p_f_ghost[idx1g(i)];
+
+ p_f_ghost_out[idx1g(i)] = p_f
+ + (1.0 - epsilon)*dp_f_expl[i] + epsilon*dp_f_impl[i];
+
+ /* apply relaxation */
+ p_f_ghost_out[idx1g(i)] = p_f*(1.0 - theta)
+ + p_f_ghost_out[idx1g(i)]*theta;
+
+ r_norm[i] = (p_f_ghost_out[idx1g(i)] - p_f)/(p_f + 1e-16);
+ }
+
+ r_norm_max = max(r_norm, sim->nz);
+#ifdef DEBUG
+ puts(".. p_f_ghost_out:"); print_array(p_f_ghost_out, sim->nz+2);
+#endif
+
+ copy_values(p_f_ghost_out, sim->p_f_ghost, sim->nz+2);
+#ifdef DEBUG
+ puts(".. p_f_ghost after update:");
+ print_array(sim->p_f_ghost, sim->nz+2);
+#endif
+
+ if (r_norm_max <= rel_tol) {
+ set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top);
+ sim->p_f_ghost[idx1g(sim->nz-1)] = p_f_top; /* top node in BC */
+ set_bc_neumann(sim->p_f_ghost, sim->nz, -1);
+ free(dp_f_expl);
+ free(dp_f_impl);
+ free(p_f_ghost_out);
+ free(r_norm);
+#ifdef DEBUG
+ printf(".. Solution converged after %d iterations\n", iter);
+#endif
+ return 0;
+ }
+ }
+
+ free(dp_f_expl);
+ free(dp_f_impl);
+ free(p_f_ghost_out);
+ free(r_norm);
+ fprintf(stderr, "darcy_solver_1d: ");
+ fprintf(stderr, "Solution did not converge after %d iterations\n", iter);
+ fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max);
+ return 1;
+}
diff --git a/fluid.h b/fluid.h
@@ -0,0 +1,13 @@
+#ifndef FLUID_H_
+#define FLUID_H_
+
+#include "simulation.h"
+
+void hydrostatic_fluid_pressure_distribution(struct simulation* sim);
+
+int darcy_solver_1d(
+ struct simulation* sim,
+ const int max_iter,
+ const double rel_tol);
+
+#endif
diff --git a/julia/1d_fd_simple_shear.jl b/julia/1d_fd_simple_shear.jl
@@ -1,277 +0,0 @@
-#!/usr/bin/env julia
-
-ENV["MPLBACKEND"] = "Agg"
-import PyPlot
-
-let
-
-# Simulation settings
-
-## Gravitational acceleration magnitude
-G = 9.81
-
-## Wall parameters
-
-### Effective normal stress exerted by top wall
-# A larger normal stress deepens the deformational depth
-P_wall_ = [10, 20, 40, 60, 80, 120] .* 1e3 # normal stress [Pa]
-
-### bottom velocity along x [m/s]
-v_x_bot = 0.0
-
-# stress ratio at top wall
-μ_wall = 0.40
-
-### Nodes along z
-nz = 100
-
-## Material properties
-
-### nonlocal amplitude [-]
-# lower values of A mean that the velocity curve can have sharper curves, e.g.
-# at the transition from μ ≈ μ_s
-#A = 0.48 # Henann and Kamrin 2016
-A = 0.40 # Loose fit to Damsgaard et al 2013
-
-### rate dependence beyond yield [-]
-# lower values of b mean larger shear velocity for a given stress ratio above
-# yield
-b = 0.9377 # Henann and Kamrin 2016
-
-### bulk and critical state static yield friction coefficient [-]
-#μ_s = 0.3819 # Henann and Kamrin 2016
-μ_s = atan(deg2rad(22.0)) # Damsgaard et al 2013
-
-### porosity [-]
-#ϕ = 0.38 # Henann and Kamrin 2016
-ϕ = 0.25 # Damsgaard et al 2013
-
-# representative grain size [m]
-# lower values of d mean that the shear velocity curve can have sharper curves,
-# e.g. at the transition from μ ≈ μ_s
-#d = 1e-3 # Henann and Kamrin 2016
-d = 0.04 # Damsgaard et al 2013
-
-### grain material density [kg/m^3]
-#ρ_s = 2.485e3 # Henann and Kamrin 2016
-ρ_s = 2.6e3 # Damsgaard et al 2013
-
-## Spatial settings
-origo_z = 0.0
-#L_z = 20.0*d # Henann and Kamrin 2016
-L_z = 0.7 # Damsgaard et al 2013
-z = collect(range(origo_z, L_z, length=nz))
-Δz = z[2] - z[1]
-
-## Allocate other arrays
-μ = zero(z) # local stress ratio
-p = zero(z) # local pressure
-v_x = zero(z) # local shear velocity
-g_ghost = zeros(size(z)[1]+2) # local fluidity with ghost nodes
-
-
-# Function definitions
-
-## Shear plastic strain rate (eq 2)
-γ_dot_p(g, μ) = g.*μ
-
-## Normal stress
-p_lithostatic(P_wall, z) = P_wall .+ (1 - ϕ).*ρ_s.*G.*(L_z .- z)
-
-## local cooperativity length
-ξ(μ) = A*d./sqrt.(abs.(μ .- μ_s))
-
-## Local fluidity
-function g_local(p, μ)
- if μ <= μ_s
- return 0.0
- else
- return sqrt(p./ρ_s.*d.^2.0) .* (μ .- μ_s)./(b.*μ)
- end
-end
-
-## Update ghost nodes for g from current values
-## BC: Neumann (dg/dx = 0)
-function set_bc_neumann(g_ghost, boundary)
- if boundary == "-z"
- g_ghost[1] = g_ghost[2]
- elseif boundary == "+z"
- g_ghost[end] = g_ghost[end-1]
- else
- @error "boundary '$boundary' not understood"
- end
-end
-
-## Update ghost nodes for g from current values
-## BC: Dirichlet (g = 0)
-function set_bc_dirichlet(g_ghost, boundary; value=0.0, idx_offset=0)
- if boundary == "-z"
- g_ghost[1+idx_offset] = value
- elseif boundary == "+z"
- g_ghost[end-idx_offset] = value
- else
- @error "boundary '$boundary' not understood"
- end
-end
-
-## Compute shear stress from velocity gradient using finite differences
-# function shear_stress(v, Δz)
-# τ = zero(v)
-#
-# # compute inner values with central differences
-# for i=2:length(v)-1
-# τ[i] = (v[i+1] - v[i-1])/(2.0*Δz)
-# end
-#
-# # use forward/backward finite differences at edges
-# τ[1] = (v[2] - v[1])/Δz
-# τ[end] = (v[end] - v[end-1])/Δz
-#
-# return τ
-# end
-
-#friction(τ, p) = τ./(p .+ 1e-16)
-
-## A single iteration for solving a Poisson equation Laplace(phi) = f on a
-## Cartesian grid. The function returns the normalized residual value
-function poisson_solver_1d_iteration(g_in, g_out, r_norm,
- μ, p, i, Δz,
- verbose=false)
-
- coorp_term = Δz^2.0/(2.0*ξ(μ[i])^2.0)
- g_out[i+1] = 1.0/(1.0 + coorp_term) * (coorp_term*g_local(p[i], μ[i])
- + g_in[i+1+1]/2.0 + g_in[i-1+1]/2.0)
- r_norm[i] = (g_out[i+1] - g_in[i+1])^2.0 / (g_out[i+1]^2.0 + 1e-16)
-
- if verbose
- println("-- $i -----------------")
- println("coorp_term: $coorp_term")
- println(" g_local: $(g_local(p[i], μ[i]))")
- println(" g_in: $(g_in[i+1])")
- println(" g_out: $(g_out[i+1])")
- println(" r_norm: $(r_norm[i])")
- end
-end
-
-## Iteratively solve the system laplace(phi) = f
-function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
- rel_tol=1e-5,
- max_iter=10_000,
- verbose=false)
-
- # allocate second array of g for Jacobian solution procedure
- g_out = zero(g)
-
- if verbose
- println("g: ")
- println(g)
- println()
- println("g_out: ")
- println(g_out)
- end
-
- # array of normalized residuals
- r_norm = zero(p)
- r_norm_max = 0.0
-
- for iter=1:max_iter
- #println("\n@@@ ITERATION $iter @@@")
-
- set_bc_dirichlet(g, "-z")
- set_bc_dirichlet(g, "+z")
- #set_bc_neumann(g, "+z")
-
- if verbose
- println("g after BC: ")
- println(g)
- end
-
- # perform a single jacobi iteration in each cell
- for iz=1:length(p)
- poisson_solver_1d_iteration(g, g_out, r_norm,
- μ, p, iz, Δz)
- end
- r_norm_max = maximum(r_norm)
-
- # Flip-flop arrays
- tmp = g
- g = g_out
- g_out = tmp
-
- if verbose
- @info ".. Relative normalized error: $r_norm_max"
- end
-
- # stop iterating if the relative tolerance is satisfied
- if r_norm_max <= rel_tol
- @info ".. Solution converged after $iter iterations"
- return
- end
- end
- @error "Solution didn't converge after $max_iter iterations ($r_norm_max)"
-end
-
-function shear_velocity(γ_dot, Δz, v_x_bot)
-
- v_x = zero(γ_dot)
-
- # BC
- v_x[1] = v_x_bot
-
- for i=2:length(v_x)
- v_x[i] = v_x[i-1] + γ_dot[i]*Δz
- end
-
- return v_x
-end
-
-function plot_profile(z, v, label, filename)
- PyPlot.figure(figsize=[4,4])
- PyPlot.plot(v, z, "+k")
- PyPlot.xlabel(label)
- PyPlot.ylabel("\$z\$ [m]")
- PyPlot.tight_layout()
- PyPlot.savefig(filename)
- PyPlot.close()
-end
-
-init_μ(μ_wall, ϕ, ρ_s, G, z, P_wall) =
- μ_wall./(1.0 .+ (1.0-ϕ)*ρ_s*G.*(L_z .- z)./P_wall)
-
-
-# Main
-
-for P_wall in P_wall_
-
- ## calculate stresses
- p = p_lithostatic(P_wall, z)
- μ = init_μ(μ_wall, ϕ, ρ_s, G, z, P_wall)
-
- ## solve for fluidity
- implicit_1d_jacobian_poisson_solver(g_ghost, p, μ, Δz)
-
- ## calculate shear velocitiesj
- γ_dot = γ_dot_p(g_ghost[2:end-1], μ)
- v_x = shear_velocity(γ_dot, Δz, v_x_bot)
-
- ## plot results
- P = Int(round(P_wall/1e3))
- PyPlot.plot(v_x/maximum(v_x), z, "+-", label="\$P_{wall}\$ = $P kPa")
- # plot_profile(z, v_x, "Shear velocity, \$v_x\$ [m/s]",
- # "1d_fd_simple_shear_v_x_P$(P)kPa.png")
- # plot_profile(z, μ, "Stress ratio, μ [-]",
- # "1d_fd_simple_shear_mu_P$(P)kPa.png")
- # plot_profile(z, p, "Normal stress, \$p\$ [Pa]",
- # "1d_fd_simple_shear_p_P$(P)kPa.png")
- # plot_profile(z, g_ghost[2:end-1], "Fluidity, \$g\$",
- # "1d_fd_simple_shear_g_P$(P)kPa.png")
-end
-
-PyPlot.xlabel("Normalized shear displacement, [m]")
-PyPlot.ylabel("Vertical position, \$z\$ [m]")
-PyPlot.legend()
-PyPlot.tight_layout()
-PyPlot.savefig("1d_fd_simple_shear.png")
-PyPlot.close()
-
-end # end let
diff --git a/julia/1d_fd_simple_shear.png b/julia/1d_fd_simple_shear.png
Binary files differ.
diff --git a/julia/Makefile b/julia/Makefile
@@ -1,7 +0,0 @@
-JULIA=julia --banner=no --color=yes
-.PHONY: run-julia
-run-julia: 1d_fd_simple_shear.jl
- echo "$<" | entr -s '$(JULIA) "$<"'
-
-1d_fd_simple_shear.png: 1d_fd_simple_shear.jl
- $(JULIA) $<
diff --git a/main.c b/main.c
@@ -4,33 +4,47 @@
#include <getopt.h>
#include "simulation.h"
+#include "fluid.h"
-#define VERSION "0.1"
+#define VERSION "0.2"
#define PROGNAME "1d_fd_simple_shear"
-/* set default simulation parameter values */
-#include "1d_fd_simple_shear_damsgaard2013.h"
+#include "parameter_defaults.h"
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 [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"
- " -b, --rate-dependence VAL rate dependence beyond yield [-]\n"
- " -f, --friction-coefficient VAL grain friction coefficient [-]\n"
- " -p, --porosity VAL porosity fraction [-]\n"
- " -d, --grain-size VAL representative grain size [m]\n"
- " -r, --density VAL grain material density [kg/m^3]\n"
- " -n, --resolution VAL number of cells in domain [-]\n"
- " -o, --origo VAL coordinate system origo [m]\n"
- " -L, --length VAL domain length [m]\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]\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"
+ " -b, --rate-dependence VAL rate dependence beyond yield [-]\n"
+ " -f, --friction-coefficient VAL grain friction coefficient [-]\n"
+ " -p, --porosity VAL porosity fraction [-]\n"
+ " -d, --grain-size VAL representative grain size [m]\n"
+ " -r, --density VAL grain material density [kg/m^3]\n"
+ " -n, --resolution VAL number of cells in domain [-]\n"
+ " -o, --origo VAL coordinate system origo [m]\n"
+ " -L, --length VAL domain length [m]\n"
+ " -F, --fluid enable pore fluid computations\n"
+ " -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 intrinsic 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"
+ " -I, --file-interval VAL interval between output files [s]\n"
+ " -v, --version show version information\n"
+ " -h, --help show this message\n"
, __func__, PROGNAME);
}
@@ -52,12 +66,12 @@ 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: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'},
@@ -69,9 +83,24 @@ int main(int argc, char* argv[])
{"resolution", required_argument, NULL, 'n'},
{"origo", required_argument, NULL, 'o'},
{"length", required_argument, NULL, 'L'},
+ {"fluid", no_argument, NULL, 'F'},
+ {"fluid-compressiblity", required_argument, NULL, 'c'},
+ {"fluid-viscosity", required_argument, NULL, 'i'},
+ {"fluid-density", 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'},
+ {"file-interval", required_argument, NULL, 'I'},
{NULL, 0, NULL, 0}
};
+ double new_phi = sim.phi[0];
+ double new_k = sim.k[0];
while ((opt = getopt_long(argc, argv, optstring, longopts, NULL)) != -1) {
switch (opt) {
case -1: /* no more arguments */
@@ -84,6 +113,9 @@ int main(int argc, char* argv[])
case 'v':
version();
return 0;
+ case 'n':
+ sim.nz = atoi(optarg);
+ break;
case 'N':
normalize = 1;
break;
@@ -109,7 +141,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);
@@ -117,15 +149,51 @@ 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;
case 'L':
sim.L_z = atof(optarg);
break;
+ case 'F':
+ sim.fluid = 1;
+ break;
+ case 'c':
+ sim.beta_f = atof(optarg);
+ break;
+ case 'i':
+ sim.mu_f = atof(optarg);
+ break;
+ case 'R':
+ sim.rho_f = atof(optarg);
+ break;
+ 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;
+ 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", argv[0], opt);
@@ -134,29 +202,90 @@ int main(int argc, char* argv[])
return -2;
}
}
+ for (int i=optind; i<argc; ++i) {
+ if (i>optind) {
+ fprintf(stderr, "error: more than one simulation name specified\n");
+ return 1;
+ }
+ sprintf(sim.name, "%s", argv[i]);
+ }
prepare_arrays(&sim);
- init_pressure(&sim);
- init_friction(&sim);
- compute_cooperativity_length(&sim);
+ if (!isnan(new_phi))
+ for (int i=0; i<sim.nz; ++i)
+ sim.phi[i] = new_phi;
+ if (!isnan(new_k))
+ 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);
+
+#ifdef DEBUG
+ 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("");
+#endif
+
+ double filetimeclock = 0.0;
+ unsigned long iter = 0;
+ while (sim.t <= sim.t_end) {
+
+ if (sim.fluid) {
+ if (darcy_solver_1d(&sim, 10000, 1e-3))
+ exit(1);
+#ifdef DEBUG
+ puts(".. p_f_ghost:"); print_array(sim.p_f_ghost, sim.nz+2);
+ puts("");
+#endif
+ }
+
+ compute_effective_stress(&sim);
+ compute_friction(&sim);
+ compute_cooperativity_length(&sim);
+
+ if (iter == 0)
+ check_simulation_parameters(&sim);
#ifdef DEBUG
- puts("\n## Before solver");
- puts(".. p:"); print_array(sim.p, sim.nz);
- puts(".. mu:"); print_array(sim.mu, sim.nz);
+ puts("\n## Before solver");
+ puts(".. sigma_n_eff:"); print_array(sim.sigma_n_eff, sim.nz);
+ puts(".. mu:"); print_array(sim.mu, sim.nz);
#endif
- if (implicit_1d_jacobian_poisson_solver(&sim, 10000, 1e-5))
- exit(1);
+ if (implicit_1d_jacobian_poisson_solver(&sim, 10000, 1e-3))
+ exit(1);
- compute_shear_strain_rate_plastic(&sim);
- compute_shear_velocity(&sim);
+ compute_shear_strain_rate_plastic(&sim);
+ compute_shear_velocity(&sim);
+
+ sim.t += sim.dt;
+ filetimeclock += sim.dt;
+ iter++;
+
+ if (filetimeclock >= sim.file_dt) {
+ write_output_file(&sim);
+ filetimeclock = 0.0;
+ }
+ }
if (normalize)
print_arrays_2nd_normalized(sim.z, sim.v_x, sim.nz);
else
- print_arrays(sim.z, sim.v_x, sim.nz);
+ if (sim.fluid)
+ for (int i=0; i<sim.nz; ++i)
+ printf("%.17g\t%.17g\t%.17g\t%.17g\n",
+ sim.z[i],
+ sim.v_x[i],
+ sim.sigma_n_eff[i],
+ sim.mu[i]);
+ /* sim.p_f_ghost[idx1g(i)]); */
+ else
+ print_three_arrays(sim.z, sim.v_x, sim.sigma_n_eff, sim.nz);
free_arrays(&sim);
return 0;
diff --git a/parameter_defaults.h b/parameter_defaults.h
@@ -0,0 +1,71 @@
+#ifndef ONED_FD_SIMPLE_SHEAR_
+#define ONED_FD_SIMPLE_SHEAR_
+
+#include <math.h>
+#include <stdio.h>
+#include "arrays.h"
+#include "simulation.h"
+
+/* Simulation settings */
+struct simulation init_sim(void)
+{
+ struct simulation sim;
+
+ sprintf(sim.name, "unnamed");
+
+ sim.G = 9.81;
+
+ sim.P_wall = 120e3; /* larger normal stress deepens the shear depth */
+ sim.mu_wall = 0.40;
+ sim.v_x_bot = 0.0;
+
+ sim.nz = 100;
+
+ /* lower values of A mean that the velocity curve can have sharper curves,
+ * e.g. at the transition from μ ≈ μ_s */
+ sim.A = 0.40; /* Loose fit to Damsgaard et al 2013 */
+
+ /* lower values of b mean larger shear velocity for a given stress ratio
+ * above yield */
+ sim.b = 0.9377; /* Henann and Kamrin 2016 */
+
+ sim.mu_s = atan(DEG2RAD(22.0)); /* Damsgaard et al 2013 */
+
+ sim.phi = initval(0.25, 1); /* 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 */
+
+ /* grain material density [kg/m^3] */
+ 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.dt = 2.0;
+ sim.t_end = 1.0;
+ sim.file_dt = 0.1;
+ sim.n_file = 0;
+
+ /* fluid settings */
+ sim.fluid = 0;
+
+ sim.beta_f = 4.5e-10; /* Water, Goren et al 2011 */
+ sim.mu_f = 1e-3; /* Water, Goren et al 2011 */
+ sim.rho_f = 1e3; /* Water */
+ sim.k = initval(1.9e-15, 1); /* Damsgaard et al 2015 */
+
+ /* 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;
+}
+
+#endif
diff --git a/simulation.c b/simulation.c
@@ -1,5 +1,6 @@
#include <stdio.h>
#include <stdlib.h>
+#include <math.h>
#include "arrays.h"
#include "simulation.h"
@@ -9,25 +10,213 @@ 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_f_ghost = zeros(sim->nz+2); /* fluid pressure with ghost nodes */
+ free(sim->phi);
+ sim->phi = zeros(sim->nz); /* porosity */
+ free(sim->k);
+ sim->k = zeros(sim->nz); /* permeability */
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_f_ghost);
+ free(sim->k);
+ free(sim->phi);
free(sim->xi);
free(sim->gamma_dot_p);
free(sim->v_x);
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.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.0)
+ warn_parameter_value("sim.dz is not a positive number", sim->dz,
+ &return_status);
+
+ check_float("sim.t", sim->t, &return_status);
+ if (sim->t < 0.0)
+ warn_parameter_value("sim.t is a negative number",
+ 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 (sim->file_dt < 0.0)
+ warn_parameter_value("sim.file_dt is a negative number",
+ sim->file_dt, &return_status);
+
+ check_float("sim.phi[0]", sim->phi[0], &return_status);
+ if (sim->phi[0] < 0.0 || sim->phi[0] > 1.0)
+ warn_parameter_value("sim.phi[0] is not within [0;1]",
+ sim->phi[0], &return_status);
+
+ if (sim->fluid != 0 && sim->fluid != 1)
+ warn_parameter_value("sim.fluid has an invalid value",
+ (double)sim->fluid, &return_status);
+
+ if (sim->fluid) {
+
+ check_float("sim.p_f_mod_ampl", sim->p_f_mod_ampl, &return_status);
+ if (sim->p_f_mod_ampl < 0.0)
+ warn_parameter_value("sim.p_f_mod_ampl is not a zero or positive",
+ sim->p_f_mod_ampl, &return_status);
+
+ check_float("sim.p_f_mod_freq", sim->p_f_mod_freq, &return_status);
+ if (sim->p_f_mod_freq < 0.0)
+ warn_parameter_value("sim.p_f_mod_freq is not a zero or positive",
+ sim->p_f_mod_freq, &return_status);
+
+ check_float("sim.beta_f", sim->beta_f, &return_status);
+ if (sim->beta_f <= 0.0)
+ warn_parameter_value("sim.beta_f is not positive",
+ sim->beta_f, &return_status);
+
+ check_float("sim.mu_f", sim->mu_f, &return_status);
+ if (sim->mu_f <= 0.0)
+ warn_parameter_value("sim.mu_f is not positive",
+ sim->mu_f, &return_status);
+
+ check_float("sim.rho_f", sim->rho_f, &return_status);
+ if (sim->rho_f <= 0.0)
+ warn_parameter_value("sim.rho_f is not positive",
+ sim->rho_f, &return_status);
+
+ check_float("sim.k[0]", sim->k[0], &return_status);
+ if (sim->k[0] <= 0.0)
+ warn_parameter_value("sim.k[0] is not positive",
+ sim->k[0], &return_status);
+
+ }
+
+ if (return_status != 0) {
+ fprintf(stderr, "error: aborting due to invalid parameter choices\n");
+ exit(return_status);
+ }
+}
+
+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]);
+}
+
+void compute_friction(struct simulation* sim)
+{
+ if (sim->fluid)
+ for (int i=0; i<sim->nz; ++i)
+ sim->mu[i] = sim->mu_wall/
+ (sim->sigma_n_eff[i]/(sim->P_wall - sim->p_f_top));
+
+ else
+ for (int i=0; i<sim->nz; ++i)
+ sim->mu[i] = sim->mu_wall /
+ (1.0 + (1.0 - sim->phi[i])*sim->rho_s*sim->G*(sim->L_z -
+ sim->z[i])
+ /sim->P_wall);
+
+}
+
double shear_strain_rate_plastic(
const double fluidity,
const double friction)
@@ -52,6 +241,16 @@ 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)
+{
+ if (sim->fluid)
+ for (int i=0; i<sim->nz; ++i)
+ sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_ghost[idx1g(i)];
+ else
+ for (int i=0; i<sim->nz; ++i)
+ sim->sigma_n_eff[i] = sim->sigma_n[i];
+}
+
double cooperativity_length(
const double A,
const double d,
@@ -86,7 +285,12 @@ 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)
@@ -180,7 +384,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,
@@ -192,7 +396,7 @@ int implicit_1d_jacobian_poisson_solver(
if (r_norm_max <= rel_tol) {
set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0);
- set_bc_neumann(sim->g_ghost, sim->nz, +1);
+ set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0);
free(g_ghost_out);
free(r_norm);
/* printf(".. Solution converged after %d iterations\n", iter); */
@@ -207,3 +411,25 @@ int implicit_1d_jacobian_poisson_solver(
fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max);
return 1;
}
+
+void write_output_file(struct simulation* sim)
+{
+
+ char outfile[200];
+ FILE *fp;
+ sprintf(outfile, "%s.output%05d.txt", sim->name, sim->n_file++);
+
+ fp = fopen(outfile, "w");
+ if (sim->fluid)
+ for (int i=0; i<sim->nz; ++i)
+ fprintf(fp, "%.17g\t%.17g\t%.17g\t%.17g\t%.17g\n",
+ sim->z[i],
+ sim->v_x[i],
+ sim->sigma_n_eff[i],
+ sim->p_f_ghost[idx1g(i)],
+ sim->mu[i]);
+ else
+ fprint_three_arrays(fp, sim->z, sim->v_x, sim->sigma_n_eff, sim->nz);
+
+ fclose(fp);
+}
diff --git a/simulation.h b/simulation.h
@@ -1,12 +1,17 @@
#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 {
+ /* simulation name to use for output files */
+ char name[100];
+
/* gravitational acceleration magnitude [m/s^2] */
double G;
@@ -28,9 +33,6 @@ struct simulation {
/* bulk and critical state static yield friction coefficient [-] */
double mu_s;
- /* porosity [-] */
- double phi;
-
/* representative grain size [m] */
double d;
@@ -52,9 +54,38 @@ struct simulation {
/* cell spacing [m] */
double dz;
- /* other arrays */
+ /* 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;
+
+ /* Fluid parameters */
+ int fluid; /* flag to switch fluid on (1) or off (0) */
+ 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] */
+
+ /* 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_f_ghost; /* fluid pressure [Pa] */
+ double* k; /* hydraulic permeability [m^2] */
+ double* phi; /* porosity [-] */
double* xi; /* cooperativity length */
double* gamma_dot_p; /* plastic shear strain rate [1/s] */
double* v_x; /* shear velocity [m/s] */
@@ -65,6 +96,10 @@ struct simulation {
void prepare_arrays(struct simulation* sim);
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,
@@ -79,10 +114,14 @@ 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);
+void compute_friction(struct simulation* sim);
int implicit_1d_jacobian_poisson_solver(
struct simulation* sim,
const int max_iter,
const double rel_tol);
+void write_output_file(struct simulation* sim);
+
#endif
