tfix code style and correct top fluid pressure with dirichlet BC - cngf-pf - continuum model for granular flows with pore-pressure dynamics (renamed from 1d_fd_simple_shear)
(HTM) git clone git://src.adamsgaard.dk/cngf-pf
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
(DIR) commit 58aabfabbd2e3c01e27e30d908514cebef245b2a
(DIR) parent 74fadc86b83bfaa0064ecf36dc68689c7db2036f
(HTM) Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 23 Nov 2020 13:20:43 +0100
fix code style and correct top fluid pressure with dirichlet BC
Diffstat:
M cngf-pf.1 | 4 ++--
M cngf-pf.c | 87 +++++++++++++++----------------
M fluid.c | 161 +++++++++++++++++--------------
M max_depth_simple_shear.c | 54 +++++++++++++++----------------
M shear_flux.c | 15 +++++++--------
M simulation.c | 161 ++++++++++++++++---------------
M simulation.h | 2 +-
7 files changed, 255 insertions(+), 229 deletions(-)
---
(DIR) diff --git a/cngf-pf.1 b/cngf-pf.1
t@@ -75,9 +75,9 @@ Only relevant with fluid dynamics enabled
.It Fl c Ar grain-cohesion
Granular material cohesion [Pa] (default 0).
.It Fl D Ar fluid-diffusivity
-Fluid diffusion coefficient [m^2/s] (default -1). Overrides fluid
+Fluid diffusion coefficient [m^2/s] (default -1). Overrides fluid
permeability (-k), grain compressibility (-P), fluid compressibility
-(-C), and fluid viscosity (-i). Do not use for transient simulations
+(-C), and fluid viscosity (-i). Do not use for transient simulations
(-T). Disabled when set to a negative value.
.It Fl d Ar grain-size
Granular material representative grain size [m] (default 0.04).
(DIR) diff --git a/cngf-pf.c b/cngf-pf.c
t@@ -3,6 +3,7 @@
#include <string.h>
#include <time.h>
#include <unistd.h>
+#include <err.h>
#include "simulation.h"
#include "fluid.h"
t@@ -22,45 +23,45 @@ static void
usage(void)
{
fprintf(stderr, "usage: %s "
- "[-A grain-nonlocal-ampl] "
- "[-a fluid-pressure-ampl] "
- "[-b grain-rate-dependence] "
- "[-C fluid-compressibility] "
- "[-c grain-cohesion] "
- "[-d grain-size] "
- "[-D fluid-diffusivity] "
- "[-e end-time] "
- "[-F] "
- "[-f applied-shear-friction] "
- "[-g gravity-accel] "
- "[-H fluid-pressure-phase] "
- "[-h] "
- "[-I file-interval] "
- "[-i fluid-viscosity] "
- "[-K dilatancy-constant] "
- "[-k fluid-permeability] "
- "[-L length] "
- "[-l applied-shear-vel-limit] "
- "[-m grain-friction] "
- "[-N] "
- "[-n normal-stress] "
- "[-O fluid-pressure-top] "
- "[-o origo] "
- "[-P grain-compressibility] "
- "[-p grain-porosity] "
- "[-q fluid-pressure-freq] "
- "[-R fluid-density] "
- "[-r grain-density] "
- "[-S fluid-pressure-pulse-shape] "
- "[-s applied-shear-vel] "
- "[-T] "
- "[-t curr-time] "
- "[-U resolution] "
- "[-u fluid-pulse-time] "
- "[-v] "
- "[-Y max-porosity] "
- "[-y min-porosity] "
- "[name]\n", argv0);
+ "[-A grain-nonlocal-ampl] "
+ "[-a fluid-pressure-ampl] "
+ "[-b grain-rate-dependence] "
+ "[-C fluid-compressibility] "
+ "[-c grain-cohesion] "
+ "[-d grain-size] "
+ "[-D fluid-diffusivity] "
+ "[-e end-time] "
+ "[-F] "
+ "[-f applied-shear-friction] "
+ "[-g gravity-accel] "
+ "[-H fluid-pressure-phase] "
+ "[-h] "
+ "[-I file-interval] "
+ "[-i fluid-viscosity] "
+ "[-K dilatancy-constant] "
+ "[-k fluid-permeability] "
+ "[-L length] "
+ "[-l applied-shear-vel-limit] "
+ "[-m grain-friction] "
+ "[-N] "
+ "[-n normal-stress] "
+ "[-O fluid-pressure-top] "
+ "[-o origo] "
+ "[-P grain-compressibility] "
+ "[-p grain-porosity] "
+ "[-q fluid-pressure-freq] "
+ "[-R fluid-density] "
+ "[-r grain-density] "
+ "[-S fluid-pressure-pulse-shape] "
+ "[-s applied-shear-vel] "
+ "[-T] "
+ "[-t curr-time] "
+ "[-U resolution] "
+ "[-u fluid-pulse-time] "
+ "[-v] "
+ "[-Y max-porosity] "
+ "[-y min-porosity] "
+ "[name]\n", argv0);
exit(1);
}
t@@ -79,10 +80,8 @@ main(int argc, char *argv[])
#endif
#ifdef __OpenBSD__
- if (pledge("stdio wpath cpath", NULL) == -1) {
- fprintf(stderr, "error: pledge failed");
- exit(2);
- }
+ if (pledge("stdio wpath cpath", NULL) == -1)
+ err(2, "pledge failed");
#endif
init_sim(&sim);
t@@ -276,7 +275,7 @@ main(int argc, char *argv[])
if ((filetimeclock >= sim.file_dt || iter == 1) &&
strncmp(sim.name, DEFAULT_SIMULATION_NAME,
- sizeof(DEFAULT_SIMULATION_NAME)) != 0) {
+ sizeof(DEFAULT_SIMULATION_NAME)) != 0) {
write_output_file(&sim, normalize);
filetimeclock = 0.0;
}
(DIR) diff --git a/fluid.c b/fluid.c
t@@ -1,5 +1,6 @@
#include <stdlib.h>
#include <math.h>
+#include <err.h>
#include "simulation.h"
#include "arrays.h"
t@@ -10,8 +11,8 @@ hydrostatic_fluid_pressure_distribution(struct simulation *sim)
for (i = 0; i < sim->nz; ++i)
sim->p_f_ghost[i + 1] = sim->p_f_top
- + sim->phi[i] * sim->rho_f * sim->G
- * (sim->L_z - sim->z[i]);
+ + sim->phi[i] * sim->rho_f * sim->G
+ * (sim->L_z - sim->z[i]);
}
static double
t@@ -19,7 +20,8 @@ diffusivity(struct simulation *sim, int i) {
if (sim->D > 0.0)
return sim->D;
else
- return sim->k[i]/((sim->alpha + sim->phi[i] * sim->beta_f) * sim->mu_f);
+ return sim->k[i]
+ / ((sim->alpha + sim->phi[i] * sim->beta_f) * sim->mu_f);
}
/* Determines the largest time step for the current simulation state. Note
t@@ -56,7 +58,7 @@ set_largest_fluid_timestep(struct simulation *sim, const double safety)
sim->dt = safety * 0.5 * dx_min * dx_min / diff_max;
if (sim->file_dt * 0.5 < sim->dt)
- sim->dt = sim->file_dt;
+ sim->dt = sim->file_dt;
return 0;
}
t@@ -100,8 +102,10 @@ square_pulse(const double time,
static void
set_fluid_bcs(double *p_f_ghost, struct simulation *sim, const double p_f_top)
{
- set_bc_dirichlet(p_f_ghost, sim->nz, +1, p_f_top);
- p_f_ghost[sim->nz] = p_f_top; /* Include top node in BC */
+ /* correct ghost node at top BC for hydrostatic pressure gradient */
+ set_bc_dirichlet(p_f_ghost, sim->nz, +1,
+ p_f_top - sim->phi[0] * sim->rho_f * sim->G * sim->dz);
+ p_f_ghost[sim->nz] = p_f_top; /* include top node in BC */
set_bc_neumann(p_f_ghost,
sim->nz,
-1,
t@@ -127,11 +131,8 @@ darcy_pressure_change_1d(const int i,
if (D > 0.0)
return D * (p_f_ghost_in[i + 2]
- 2.0 * p_f_ghost_in[i + 1]
- + p_f_ghost_in[i])
- / (dz * dz);
-
+ + p_f_ghost_in[i]) / (dz * dz);
else {
-
k_ = k[i];
if (i == 0)
k_zn = k_;
t@@ -141,29 +142,25 @@ darcy_pressure_change_1d(const int i,
k_zp = k_;
else
k_zp = k[i + 1];
-#ifdef DEBUG
- printf("%d->%d: p=[%g, %g, %g]\tk=[%g, %g, %g]\n",
- i, i + 1,
- p_f_ghost_in[i], p_f_ghost_in[i + 1], p_f_ghost_in[i + 2],
- k_zn, k_, k_zp);
-#endif
div_k_grad_p = (2.0 * k_zp * k_ / (k_zp + k_)
- * (p_f_ghost_in[i + 2]
- - p_f_ghost_in[i + 1]) / dz
- - 2.0 * k_zn * k_ / (k_zn + k_)
- * (p_f_ghost_in[i + 1]
- - p_f_ghost_in[i]) / dz
- ) / dz;
+ * (p_f_ghost_in[i + 2] - p_f_ghost_in[i + 1]) / dz
+ - 2.0 * k_zn * k_ / (k_zn + k_)
+ * (p_f_ghost_in[i + 1] - p_f_ghost_in[i]) / dz
+) / dz;
#ifdef DEBUG
- printf("darcy_pressure_change [%d]: phi=%g\tdiv_k_grad_p=%g\tphi_dot=%g\n",
- i, phi[i], div_k_grad_p, phi_dot[i]);
+ printf("%s [%d]: phi=%g\tdiv_k_grad_p=%g\tphi_dot=%g\n",
+ __func__, i, phi[i], div_k_grad_p, phi_dot[i]);
+
+ printf(" p[%d, %d, %d] = [%g, %g, %g]\tk=[%g, %g, %g]\n",
+ i, i + 1, i + 2,
+ p_f_ghost_in[i], p_f_ghost_in[i + 1], p_f_ghost_in[i + 2],
+ k_zn, k_, k_zp);
#endif
- /* TODO: add advective term */
return 1.0 / ((alpha + beta_f * phi[i]) * mu_f) * div_k_grad_p
- - 1.0 / ((alpha + beta_f * phi[i]) * (1.0 - phi[i])) * phi_dot[i];
+ - 1.0 / ((alpha + beta_f * phi[i]) * (1.0 - phi[i])) * phi_dot[i];
}
}
t@@ -174,41 +171,44 @@ darcy_solver_1d(struct simulation *sim,
{
int i, iter, solved = 0;
double epsilon, p_f_top, r_norm_max = NAN;
- double *tmp;
/* choose integration method, parameter in [0.0; 1.0]
* epsilon = 0.0: explicit
* epsilon = 0.5: Crank-Nicolson
* epsilon = 1.0: implicit */
- epsilon = 0.0;
+ epsilon = 0.5;
+
+ for (i = 0; i < sim->nz; ++i)
+ sim->p_f_dot_impl[i] = sim->p_f_dot_expl[i] = 0.0;
if (isnan(sim->p_f_mod_pulse_time))
p_f_top = sim->p_f_top
- + sine_wave(sim->t,
- sim->p_f_mod_ampl,
- sim->p_f_mod_freq,
- sim->p_f_mod_phase);
+ + sine_wave(sim->t,
+ sim->p_f_mod_ampl,
+ sim->p_f_mod_freq,
+ sim->p_f_mod_phase);
else if (sim->p_f_mod_pulse_shape == 1)
p_f_top = sim->p_f_top
- + square_pulse(sim->t,
- sim->p_f_mod_ampl,
- sim->p_f_mod_freq,
- sim->p_f_mod_pulse_time);
+ + square_pulse(sim->t,
+ sim->p_f_mod_ampl,
+ sim->p_f_mod_freq,
+ sim->p_f_mod_pulse_time);
else
p_f_top = sim->p_f_top
- + triangular_pulse(sim->t,
- sim->p_f_mod_ampl,
- sim->p_f_mod_freq,
- sim->p_f_mod_pulse_time);
-
+ + triangular_pulse(sim->t,
+ sim->p_f_mod_ampl,
+ sim->p_f_mod_freq,
+ sim->p_f_mod_pulse_time);
/* set fluid BCs (1 of 2) */
set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
- copy_values(sim->p_f_ghost, sim->new_ghost, sim->nz + 2);
/* explicit solution to pressure change */
if (epsilon < 1.0) {
- for (i = 0; i < sim->nz; ++i)
+#ifdef DEBUG
+ printf("\nEXPLICIT SOLVER IN %s\n", __func__);
+#endif
+ for (i = 0; i < sim->nz - 1; ++i)
sim->p_f_dot_expl[i] = darcy_pressure_change_1d(i,
sim->nz,
sim->p_f_ghost,
t@@ -219,26 +219,38 @@ darcy_solver_1d(struct simulation *sim,
sim->beta_f,
sim->alpha,
sim->mu_f,
- sim->D);
+ sim->D);
+
}
/* implicit solution with Jacobian iterations */
if (epsilon > 0.0) {
+#ifdef DEBUG
+ printf("\nEXPLICIT SOLVER IN %s\n", __func__);
+#endif
+ copy_values(sim->p_f_ghost, sim->tmp_ghost, sim->nz + 2);
+
for (iter = 0; iter < max_iter; ++iter) {
copy_values(sim->p_f_dot_impl, sim->old_val, sim->nz);
+#ifdef DEBUG
+ puts(".. p_f_ghost bfore BC:");
+ print_array(sim->tmp_ghost, sim->nz + 2);
+#endif
+
/* set fluid BCs (2 of 2) */
- set_fluid_bcs(sim->new_ghost, sim, p_f_top);
+ set_fluid_bcs(sim->tmp_ghost, sim, p_f_top);
+
#ifdef DEBUG
puts(".. p_f_ghost after BC:");
- print_array(sim->new_ghost, sim->nz + 2);
+ print_array(sim->tmp_ghost, sim->nz + 2);
#endif
for (i = 0; i < sim->nz - 1; ++i)
sim->p_f_dot_impl[i] = darcy_pressure_change_1d(i,
sim->nz,
- sim->new_ghost,
+ sim->tmp_ghost,
sim->phi,
sim->phi_dot,
sim->k,
t@@ -246,27 +258,26 @@ darcy_solver_1d(struct simulation *sim,
sim->beta_f,
sim->alpha,
sim->mu_f,
- sim->D);
+ sim->D);
- for (i = 0; i < sim->nz - 1; ++i) {
- sim->new_ghost[i + 1] = sim->p_f_ghost[i + 1]
- + sim->p_f_dot_impl[i] * sim->dt;
- sim->p_f_dot_impl_r_norm[i] = fabs(residual(
- sim->p_f_dot_impl[i], sim->old_val[i]));
- }
+ /*for (i = 0; i < sim->nz; ++i)
+ printf("sim->p_f_dot_impl[%d] = %g (t = %g s, iter = %d)\n",
+ i, sim->p_f_dot_impl[i], sim->t, iter);*/
- r_norm_max = max(sim->p_f_dot_impl_r_norm, sim->nz);
-#ifdef DEBUG
- puts(".. p_f_ghost_new:");
- print_array(sim->new_ghost, sim->nz + 2);
-#endif
+ /*for (i = 0; i < sim->nz; ++i)
+ if (fabs(sim->p_f_dot_impl[i]) > 1e-12)
+ errx(1, "sim->p_f_dot_impl[%d] = %g (t = %g s, iter = %d)",
+ i, sim->p_f_dot_impl[i], sim->t, iter);*/
+
+ for (i = 0; i < sim->nz; ++i)
+ if (isnan(sim->p_f_dot_impl[i]))
+ errx(1, "NaN at sim->p_f_dot_impl[%d] (t = %g s, iter = %d)",
+ i, sim->t, iter);
- tmp = sim->new_ghost;
- sim->new_ghost = sim->p_f_ghost;
- sim->p_f_ghost = tmp;
+ r_norm_max = max(sim->p_f_dot_impl_r_norm, sim->nz - 1);
#ifdef DEBUG
- puts(".. p_f_ghost after update:");
- print_array(sim->new_ghost, sim->nz + 2);
+ puts(".. p_f_ghost_new:");
+ print_array(sim->tmp_ghost, sim->nz + 2);
#endif
if (r_norm_max <= rel_tol) {
t@@ -291,15 +302,25 @@ darcy_solver_1d(struct simulation *sim,
sim->p_f_dot[i] = epsilon * sim->p_f_dot_impl[i]
+ (1.0 - epsilon) * sim->p_f_dot_expl[i];
+ set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
+
#ifdef DEBUG
- printf(".. epsilon = %s\n", epsilon);
- puts(".. p_f_dot_impl:");
- print_array(sim->p_f_dot_impl, sim->nz);
- puts(".. p_f_dot_expl:");
- print_array(sim->p_f_dot_expl, sim->nz);
+ printf(".. epsilon = %g\n", epsilon);
+ puts(".. p_f_dot_expl:");
+ print_array(sim->p_f_dot_expl, sim->nz);
+ puts(".. p_f_dot_impl:");
+ print_array(sim->p_f_dot_impl, sim->nz);
#endif
- set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
+ for (i = 0; i < sim->nz; ++i)
+ if (isnan(sim->p_f_dot_expl[i]) || isinf(sim->p_f_dot_expl[i]))
+ errx(1, "invalid: sim->p_f_dot_expl[%d] = %g (t = %g s)",
+ i, sim->p_f_dot_expl[i], sim->t);
+
+ for (i = 0; i < sim->nz; ++i)
+ if (isnan(sim->p_f_dot_impl[i]) || isinf(sim->p_f_dot_impl[i]))
+ errx(1, "invalid: sim->p_f_dot_impl[%d] = %g (t = %g s)",
+ i, sim->p_f_dot_impl[i], sim->t);
return solved - 1;
}
(DIR) diff --git a/max_depth_simple_shear.c b/max_depth_simple_shear.c
t@@ -4,6 +4,7 @@
#include <math.h>
#include <getopt.h>
#include <time.h>
+#include <err.h>
#include "simulation.h"
#include "arg.h"
t@@ -27,21 +28,21 @@ static void
usage(void)
{
fprintf(stderr, "usage: %s "
- "[-a fluid-pressure-ampl] "
- "[-C fluid-compressibility] "
- "[-D fluid-diffusivity] "
- "[-g gravity-accel] "
- "[-h] "
- "[-i fluid-viscosity] "
- "[-k fluid-permeability] "
- "[-O fluid-pressure-top] "
- "[-P grain-compressibility] "
- "[-p grain-porosity] "
- "[-q fluid-pressure-freq] "
- "[-R fluid-density] "
- "[-r grain-density] "
- "[-v] "
- "\n", argv0);
+ "[-a fluid-pressure-ampl] "
+ "[-C fluid-compressibility] "
+ "[-D fluid-diffusivity] "
+ "[-g gravity-accel] "
+ "[-h] "
+ "[-i fluid-viscosity] "
+ "[-k fluid-permeability] "
+ "[-O fluid-pressure-top] "
+ "[-P grain-compressibility] "
+ "[-p grain-porosity] "
+ "[-q fluid-pressure-freq] "
+ "[-R fluid-density] "
+ "[-r grain-density] "
+ "[-v] "
+ "\n", argv0);
exit(1);
}
t@@ -51,7 +52,9 @@ skin_depth(const struct simulation *sim)
if (sim->D > 0.0)
return sqrt(sim->D / (M_PI * sim->p_f_mod_freq));
else
- return sqrt(sim->k[0] / (sim->mu_f * (sim->alpha + sim->phi[0] * sim->beta_f) * M_PI * sim->p_f_mod_freq));
+ return sqrt(sim->k[0] / (sim->mu_f
+ * (sim->alpha + sim->phi[0] * sim->beta_f)
+ * M_PI * sim->p_f_mod_freq));
}
/* using alternate form: sin(x) + cos(x) = sqrt(2)*sin(x + pi/4) */
t@@ -65,8 +68,8 @@ eff_normal_stress_gradient(struct simulation * sim, double d_s, double z_)
exit(10);
}
return sqrt(2.0) * sin((3.0 * M_PI / 2.0 - z_ / d_s) + M_PI / 4.0)
- + (sim->rho_s - sim->rho_f) * sim->G * d_s
- / sim->p_f_mod_ampl * exp(z_ / d_s);
+ + (sim->rho_s - sim->rho_f) * sim->G * d_s
+ / sim->p_f_mod_ampl * exp(z_ / d_s);
}
/* use Brent's method for finding roots (Press et al., 2007) */
t@@ -88,9 +91,8 @@ zbrent(struct simulation *sim,
fb = (*f) (sim, d_s, b);
if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) {
- fprintf(stderr,
- "error: %s: no root in range %g,%g m\n",
- __func__, x1, x2);
+ warnx("%s: no root in range %g,%g m\n",
+ __func__, x1, x2);
free_arrays(sim);
exit(11);
}
t@@ -149,9 +151,7 @@ zbrent(struct simulation *sim,
b += ((xm) >= 0.0 ? fabs(tol1) : -fabs(tol1));
fb = (*f) (sim, d_s, b);
}
- fprintf(stderr,
- "error: %s: exceeded maximum number of iterations",
- __func__);
+ warnx("error: %s: exceeded maximum number of iterations", __func__);
free_arrays(sim);
exit(12);
return NAN;
t@@ -172,10 +172,8 @@ main(int argc, char *argv[])
#endif
#ifdef __OpenBSD__
- if (pledge("stdio", NULL) == -1) {
- fprintf(stderr, "error: pledge failed");
- exit(2);
- }
+ if (pledge("stdio", NULL) == -1)
+ err(2, "pledge failed");
#endif
init_sim(&sim);
(DIR) diff --git a/shear_flux.c b/shear_flux.c
t@@ -2,6 +2,7 @@
#include <unistd.h>
#include <stdio.h>
#include <stdlib.h>
+#include <err.h>
#include "arg.h"
t@@ -11,10 +12,10 @@ static void
usage(void)
{
fprintf(stderr, "usage: %s "
- "[-v] "
- "[-h] "
- "[file ...] "
- "\n", argv0);
+ "[-v] "
+ "[-h] "
+ "[file ...] "
+ "\n", argv0);
exit(1);
}
t@@ -46,10 +47,8 @@ main(int argc, char *argv[])
FILE *fp;
#ifdef __OpenBSD__
- if (pledge("stdio rpath", NULL) == -1) {
- fprintf(stderr, "error: pledge failed");
- exit(2);
- }
+ if (pledge("stdio rpath", NULL) == -1)
+ err(2, "pledge failed");
#endif
ARGBEGIN {
(DIR) diff --git a/simulation.c b/simulation.c
t@@ -1,6 +1,7 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
+#include <err.h>
#include "arrays.h"
#include "simulation.h"
#include "fluid.h"
t@@ -135,15 +136,15 @@ prepare_arrays(struct simulation *sim)
{
if (sim->nz < 2) {
fprintf(stderr, "error: grid size (nz) must be at least 2 but is %d\n",
- sim->nz);
+ sim->nz);
exit(1);
}
free(sim->phi);
free(sim->k);
sim->z = linspace(sim->origo_z,
- sim->origo_z + sim->L_z,
- sim->nz);
+ sim->origo_z + sim->L_z,
+ sim->nz);
sim->dz = sim->z[1] - sim->z[0];
sim->mu = zeros(sim->nz);
sim->mu_c = zeros(sim->nz);
t@@ -167,7 +168,7 @@ prepare_arrays(struct simulation *sim)
sim->I = zeros(sim->nz);
sim->tan_psi = zeros(sim->nz);
sim->old_val = empty(sim->nz);
- sim->new_ghost = empty(sim->nz + 2);
+ sim->tmp_ghost = empty(sim->nz + 2);
}
void
t@@ -196,13 +197,13 @@ free_arrays(struct simulation *sim)
free(sim->I);
free(sim->tan_psi);
free(sim->old_val);
- free(sim->new_ghost);
+ free(sim->tmp_ghost);
}
static void
warn_parameter_value(const char message[],
- const double value,
- int *return_status)
+ const double value,
+ int *return_status)
{
fprintf(stderr,
"check_simulation_parameters: %s (%.17g)\n",
t@@ -243,14 +244,14 @@ check_simulation_parameters(struct simulation *sim)
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);
+ &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);
+ &return_status);
check_float("sim->A", sim->A, &return_status);
if (sim->A < 0.0)
t@@ -263,128 +264,134 @@ check_simulation_parameters(struct simulation *sim)
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);
+ &return_status);
check_float("sim->C", sim->C, &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);
+ &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);
+ &return_status);
if (sim->nz <= 0)
warn_parameter_value("sim->nz is not a positive number", sim->nz,
- &return_status);
+ &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);
+ sim->L_z, &return_status);
if (sim->nz <= 0)
warn_parameter_value("sim->nz is not a positive number", sim->nz,
- &return_status);
+ &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);
+ &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);
+ 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);
+ sim->t, &return_status);
check_float("sim->dt", sim->dt, &return_status);
if (sim->dt <= 0.0)
warn_parameter_value("sim->dt is not a positive number",
- sim->dt, &return_status);
+ 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);
+ 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);
+ sim->phi[0], &return_status);
check_float("sim->phi_min", sim->phi_min, &return_status);
if (sim->phi_min < 0.0 || sim->phi_min > 1.0)
warn_parameter_value("sim->phi_min is not within [0;1]",
- sim->phi_min, &return_status);
+ sim->phi_min, &return_status);
check_float("sim->phi_max", sim->phi_max, &return_status);
if (sim->phi_max < 0.0 || sim->phi_max > 1.0)
warn_parameter_value("sim->phi_max is not within [0;1]",
- sim->phi_max, &return_status);
+ sim->phi_max, &return_status);
check_float("sim->dilatancy_constant", sim->dilatancy_constant, &return_status);
if (sim->dilatancy_constant < 0.0 || sim->dilatancy_constant > 100.0)
warn_parameter_value("sim->dilatancy_constant is not within [0;100]",
- sim->dilatancy_constant, &return_status);
+ sim->dilatancy_constant, &return_status);
if (sim->fluid != 0 && sim->fluid != 1)
warn_parameter_value("sim->fluid has an invalid value",
- (double) sim->fluid, &return_status);
+ (double) sim->fluid, &return_status);
if (sim->transient != 0 && sim->transient != 1)
warn_parameter_value("sim->transient has an invalid value",
- (double) sim->transient, &return_status);
+ (double) sim->transient, &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);
+ sim->p_f_mod_ampl, &return_status);
if (sim->P_wall - sim->p_f_mod_ampl < 0.0)
warn_parameter_value("sim->P_wall - sim->p_f_mod_ampl is negative",
- sim->P_wall - sim->p_f_mod_ampl,
- &return_status);
+ sim->P_wall - 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);
+ 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);
+ sim->beta_f, &return_status);
check_float("sim->alpha", sim->alpha, &return_status);
if (sim->alpha <= 0.0)
warn_parameter_value("sim->alpha is not positive",
- sim->alpha, &return_status);
+ sim->alpha, &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);
+ 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);
+ 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);
+ sim->k[0], &return_status);
+
+ check_float("sim->D", sim->D, &return_status);
+ if (sim->transient && sim->D > 0.0)
+ warn_parameter_value("constant diffusivity does not work in "
+ "transient simulations",
+ sim->D, &return_status);
+
}
if (return_status != 0) {
fprintf(stderr, "error: aborting due to invalid parameter choices\n");
t@@ -398,9 +405,9 @@ lithostatic_pressure_distribution(struct simulation *sim)
int i;
for (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]);
+ sim->sigma_n[i] = sim->P_wall
+ + (1.0 - sim->phi[i]) * sim->rho_s * sim->G
+ * (sim->L_z - sim->z[i]);
}
/* NEW FUNCTIONS START */
t@@ -412,7 +419,7 @@ compute_inertia_number(struct simulation *sim)
for (i = 0; i < sim->nz; ++i)
sim->I[i] = sim->gamma_dot_p[i] * sim->d
- / sqrt(sim->sigma_n_eff[i] / sim->rho_s);
+ / sqrt(sim->sigma_n_eff[i] / sim->rho_s);
}
void
t@@ -432,13 +439,13 @@ compute_critical_state_friction(struct simulation *sim)
if (sim->fluid)
for (i = 0; i < sim->nz; ++i)
- sim->mu_c[i] = sim->mu_wall /
- (sim->sigma_n_eff[i] / (sim->P_wall - sim->p_f_top));
+ sim->mu_c[i] = sim->mu_wall / (sim->sigma_n_eff[i]
+ / (sim->P_wall - sim->p_f_top));
else
for (i = 0; i < sim->nz; ++i)
- sim->mu_c[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);
+ sim->mu_c[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);
}
static void
t@@ -479,8 +486,8 @@ compute_permeability(struct simulation *sim)
for (i = 0; i < sim->nz; ++i)
sim->k[i] = pow(sim->d, 2.0) / 180.0
- * pow(sim->phi[i], 3.0)
- / pow(1.0 - sim->phi[i], 2.0);
+ * pow(sim->phi[i], 3.0)
+ / pow(1.0 - sim->phi[i], 2.0);
}
/* NEW FUNCTIONS END */
t@@ -498,7 +505,7 @@ compute_shear_strain_rate_plastic(struct simulation *sim)
for (i = 0; i < sim->nz; ++i)
sim->gamma_dot_p[i] = shear_strain_rate_plastic(sim->g_ghost[i + 1],
- sim->mu[i]);
+ sim->mu[i]);
}
static void
t@@ -593,10 +600,8 @@ set_bc_neumann(double *a,
a[0] = a[1] + df * dx;
else if (boundary == +1)
a[nz + 1] = a[nz] - df * dx;
- else {
- fprintf(stderr, "set_bc_neumann: Unknown boundary %d\n", boundary);
- exit(1);
- }
+ else
+ errx(1, "%s: Unknown boundary %d\n", __func__, boundary);
}
void
t@@ -609,10 +614,8 @@ set_bc_dirichlet(double *a,
a[0] = value;
else if (boundary == +1)
a[nz + 1] = value;
- else {
- fprintf(stderr, "set_bc_dirichlet: Unknown boundary %d\n", boundary);
- exit(1);
- }
+ else
+ errx(1, "%s: Unknown boundary %d\n", __func__, boundary);
}
double
t@@ -633,8 +636,8 @@ poisson_solver_1d_cell_update(int i,
double coorp_term = dz * dz / (2.0 * pow(xi[i], 2.0));
g_out[i + 1] = 1.0 / (1.0 + coorp_term)
- * (coorp_term * g_local[i] + g_in[i + 2] / 2.0 + g_in[i] / 2.0);
-
+ * (coorp_term * g_local[i] + g_in[i + 2] / 2.0
+ + g_in[i] / 2.0);
r_norm[i] = residual(g_out[i + 1], g_in[i + 1]);
#ifdef DEBUG
t@@ -671,14 +674,14 @@ implicit_1d_jacobian_poisson_solver(struct simulation *sim,
poisson_solver_1d_cell_update(i,
sim->g_ghost,
sim->g_local,
- sim->new_ghost,
+ sim->tmp_ghost,
sim->g_r_norm,
sim->dz,
sim->xi);
r_norm_max = max(sim->g_r_norm, sim->nz);
- tmp = sim->new_ghost;
- sim->new_ghost = sim->g_ghost;
+ tmp = sim->tmp_ghost;
+ sim->tmp_ghost = sim->g_ghost;
sim->g_ghost = tmp;
if (r_norm_max <= rel_tol) {
t@@ -704,7 +707,7 @@ write_output_file(struct simulation *sim, const int normalize)
FILE *fp;
snprintf(outfile, sizeof(outfile), "%s.output%05d.txt",
- sim->name, sim->n_file++);
+ sim->name, sim->n_file++);
if ((fp = fopen(outfile, "w")) != NULL) {
print_output(sim, fp, normalize);
t@@ -753,9 +756,16 @@ temporal_increment(struct simulation *sim)
for (i = 0; i < sim->nz; ++i)
sim->phi[i] += sim->phi_dot[i] * sim->dt;
- if (sim->fluid)
- for (i = 0; i < sim->nz; ++i)
- sim->p_f_ghost[i + 1] += sim->p_f_dot[i] * sim->dt;
+ if (sim->fluid) {
+ for (i = 0; i < sim->nz; ++i) {
+ if (isnan(sim->p_f_dot[i])) {
+ errx(1, "encountered NaN at sim->p_f_dot[%d] (t = %g s)",
+ i, sim->t);
+ } else {
+ sim->p_f_ghost[i + 1] += sim->p_f_dot[i] * sim->dt;
+ }
+ }
+ }
sim->t += sim->dt;
}
t@@ -769,7 +779,6 @@ coupled_shear_solver(struct simulation *sim,
double r_norm_max, vel_res_norm = NAN, mu_wall_orig = sim->mu_wall;
do { /* stress iterations */
-
coupled_iter = 0;
do { /* coupled iterations */
t@@ -806,8 +815,8 @@ coupled_shear_solver(struct simulation *sim,
/* step 8, Eq. 11 */
if (implicit_1d_jacobian_poisson_solver(sim,
- MAX_ITER_GRANULAR,
- RTOL_GRANULAR))
+ MAX_ITER_GRANULAR,
+ RTOL_GRANULAR))
exit(12);
/* step 9 */
t@@ -840,9 +849,9 @@ coupled_shear_solver(struct simulation *sim,
if (coupled_iter++ >= max_iter) {
fprintf(stderr, "coupled_shear_solver: ");
fprintf(stderr, "Transient solution did not converge "
- "after %d iterations\n", coupled_iter);
+ "after %d iterations\n", coupled_iter);
fprintf(stderr, ".. Residual normalized error: %g\n",
- r_norm_max);
+ r_norm_max);
return 1;
}
}
t@@ -851,25 +860,25 @@ coupled_shear_solver(struct simulation *sim,
if (!isnan(sim->v_x_limit) || !isnan(sim->v_x_fix)) {
if (!isnan(sim->v_x_limit)) {
vel_res_norm = (sim->v_x_limit - sim->v_x[sim->nz - 1])
- / (sim->v_x[sim->nz - 1] + 1e-12);
+ / (sim->v_x[sim->nz - 1] + 1e-12);
if (vel_res_norm > 0.0)
vel_res_norm = 0.0;
} else {
vel_res_norm = (sim->v_x_fix - sim->v_x[sim->nz - 1])
- / (sim->v_x[sim->nz - 1] + 1e-12);
+ / (sim->v_x[sim->nz - 1] + 1e-12);
}
sim->mu_wall *= 1.0 + (vel_res_norm * 1e-2);
}
if (++stress_iter > MAX_ITER_STRESS) {
fprintf(stderr, "error: stress solution did not converge:\n");
fprintf(stderr, "v_x=%g, v_x_fix=%g, v_x_limit=%g, "
- "vel_res_norm=%g, mu_wall=%g\n",
- sim->v_x[sim->nz - 1], sim->v_x_fix, sim->v_x_limit,
- vel_res_norm, sim->mu_wall);
+ "vel_res_norm=%g, mu_wall=%g\n",
+ sim->v_x[sim->nz - 1], sim->v_x_fix, sim->v_x_limit,
+ vel_res_norm, sim->mu_wall);
return 10;
}
} while ((!isnan(sim->v_x_fix) || !isnan(sim->v_x_limit))
- && fabs(vel_res_norm) > RTOL_STRESS);
+ && fabs(vel_res_norm) > RTOL_STRESS);
if (!isnan(sim->v_x_limit))
sim->mu_wall = mu_wall_orig;
(DIR) diff --git a/simulation.h b/simulation.h
t@@ -123,7 +123,7 @@ struct simulation {
double *I; /* inertia number [-] */
double *tan_psi; /* tan(dilatancy_angle) [-] */
double *old_val; /* temporary storage for iterative solvers */
- double *new_ghost; /* temporary storage for iterative solvers */
+ double *tmp_ghost; /* temporary storage for iterative solvers */
};
void init_sim(struct simulation *sim);