commit 92fb592d0bd18bef9855ce1b5a596107c0fa40bb
parent 3dbb7bbe713574ed78923bcc2bddc1a3046139bc
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Fri, 3 Jul 2026 11:24:54 +0200
fix(simulation): harden fluid BC, origo depth, and solver error handling
Prepare cngf-pf for use as an embedded per-column till solver:
- Correct the origo_z offset in all four depth terms (hydrostatic,
lithostatic, critical-state friction, fused transient friction) so
physical depth below the top surface is independent of origo_z; the
ice-till interface node now equals p_f_top.
- Enforce the top Dirichlet pressure as an identity row in the Darcy
TDMA system and use top porosity phi[nz-1] in the top ghost correction.
- Return error codes instead of exiting from library solver paths:
darcy_solver_1d and temporal_increment on NaN/inf, and
coupled_shear_solver (11 Darcy, 12 Poisson, 13 temporal increment);
prepare_arrays and check_simulation_parameters now return int. CLI
drivers translate codes back to exit/errx.
- Add reset_column() for in-place per-column reuse without reallocation.
- Reindent simulation.c and fluid.c to tabs (UseTab: ForIndentation).
Diffstat:
| M | cngf-pf.c | | | 12 | ++++++++---- |
| M | deborah.c | | | 8 | ++++++-- |
| M | fluid.c | | | 425 | ++++++++++++++++++++++++++++++++++++++++--------------------------------------- |
| M | max_depth_simple_shear.c | | | 8 | ++++++-- |
| M | simulation.c | | | 1554 | +++++++++++++++++++++++++++++++++++++++++-------------------------------------- |
| M | simulation.h | | | 263 | +++++++++++++++++++++++++++++++++++++++++-------------------------------------- |
6 files changed, 1176 insertions(+), 1094 deletions(-)
diff --git a/cngf-pf.c b/cngf-pf.c
@@ -247,7 +247,8 @@ main(int argc, char *argv[])
if (sim.nz < 1)
sim.nz = (int) ceil(sim.L_z / sim.d);
- prepare_arrays(&sim);
+ if (prepare_arrays(&sim))
+ errx(1, "prepare_arrays failed");
if (!isnan(new_phi))
for (i = 0; i < sim.nz; ++i)
@@ -279,15 +280,18 @@ main(int argc, char *argv[])
sim.dt = sim.t_end;
compute_effective_stress(&sim);
- check_simulation_parameters(&sim);
+ if (check_simulation_parameters(&sim)) {
+ free_arrays(&sim);
+ exit(1);
+ }
filetimeclock = 0.0;
iter = 0;
t_begin = clock();
do {
- if (coupled_shear_solver(&sim, max_iter, rtol)) {
+ if ((ret = coupled_shear_solver(&sim, max_iter, rtol))) {
free_arrays(&sim);
- exit(10);
+ exit(ret);
}
if ((filetimeclock >= sim.file_dt || iter == 1) &&
diff --git a/deborah.c b/deborah.c
@@ -103,7 +103,8 @@ main(int argc, char *argv[])
usage();
sim.nz = 2;
- prepare_arrays(&sim);
+ if (prepare_arrays(&sim))
+ errx(1, "prepare_arrays failed");
if (!isnan(new_phi))
for (i = 0; i < sim.nz; ++i)
@@ -112,7 +113,10 @@ main(int argc, char *argv[])
for (i = 0; i < sim.nz; ++i)
sim.k[i] = new_k;
- check_simulation_parameters(&sim);
+ if (check_simulation_parameters(&sim)) {
+ free_arrays(&sim);
+ exit(1);
+ }
#ifdef BENCHMARK_PERFORMANCE
t_begin = clock();
diff --git a/fluid.c b/fluid.c
@@ -1,23 +1,23 @@
#include "fluid.h"
#include "arrays.h"
#include "simulation.h"
-#include <err.h>
#include <math.h>
#include <stdlib.h>
void hydrostatic_fluid_pressure_distribution(struct simulation *sim) {
- int i;
+ int i;
- 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]);
+ 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->origo_z + sim->L_z - sim->z[i]);
}
static double 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);
+ if (sim->D > 0.0)
+ return sim->D;
+ else
+ 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
@@ -25,229 +25,234 @@ static double diffusivity(struct simulation *sim, int i) {
* diffusivities (i.e., permeabilities, porosities, viscosities, or
* compressibilities) change. The safety factor should be in ]0;1] */
int set_largest_fluid_timestep(struct simulation *sim, const double safety) {
- int i;
- double dx_min, diff, diff_max, *dx;
-
- dx = spacing(sim->z, sim->nz);
- dx_min = INFINITY;
- for (i = 0; i < sim->nz - 1; ++i) {
- if (dx[i] < 0.0) {
- fprintf(stderr, "error: cell spacing negative (%g) in cell %d\n", dx[i],
- i);
- free(dx);
- return 1;
- }
- if (dx[i] < dx_min)
- dx_min = dx[i];
- }
- free(dx);
-
- diff_max = -INFINITY;
- for (i = 0; i < sim->nz; ++i) {
- diff = diffusivity(sim, i);
- if (diff > diff_max)
- diff_max = diff;
- }
-
- sim->dt = safety * 0.5 * dx_min * dx_min / diff_max;
- if (sim->file_dt < sim->dt)
- sim->dt = sim->file_dt;
-
- return 0;
+ int i;
+ double dx_min, diff, diff_max, *dx;
+
+ dx = spacing(sim->z, sim->nz);
+ dx_min = INFINITY;
+ for (i = 0; i < sim->nz - 1; ++i) {
+ if (dx[i] < 0.0) {
+ fprintf(stderr, "error: cell spacing negative (%g) in cell %d\n", dx[i],
+ i);
+ free(dx);
+ return 1;
+ }
+ if (dx[i] < dx_min)
+ dx_min = dx[i];
+ }
+ free(dx);
+
+ diff_max = -INFINITY;
+ for (i = 0; i < sim->nz; ++i) {
+ diff = diffusivity(sim, i);
+ if (diff > diff_max)
+ diff_max = diff;
+ }
+
+ sim->dt = safety * 0.5 * dx_min * dx_min / diff_max;
+ if (sim->file_dt < sim->dt)
+ sim->dt = sim->file_dt;
+
+ return 0;
}
static double sine_wave(const double time, const double ampl, const double freq,
const double phase) {
- return ampl * sin(2.0 * PI * freq * time + phase);
+ return ampl * sin(2.0 * PI * freq * time + phase);
}
static double triangular_pulse(const double time, const double peak_ampl,
const double freq, const double peak_time) {
- if (peak_time - 1.0 / (2.0 * freq) < time && time <= peak_time)
- return peak_ampl * 2.0 * freq * (time - peak_time) + peak_ampl;
- else if (peak_time < time && time < peak_time + 1.0 / (2.0 * freq))
- return peak_ampl * 2.0 * freq * (peak_time - time) + peak_ampl;
- else
- return 0.0;
+ if (peak_time - 1.0 / (2.0 * freq) < time && time <= peak_time)
+ return peak_ampl * 2.0 * freq * (time - peak_time) + peak_ampl;
+ else if (peak_time < time && time < peak_time + 1.0 / (2.0 * freq))
+ return peak_ampl * 2.0 * freq * (peak_time - time) + peak_ampl;
+ else
+ return 0.0;
}
static double square_pulse(const double time, const double peak_ampl,
const double freq, const double peak_time) {
- if (peak_time - 1.0 / (2.0 * freq) < time &&
- time < peak_time + 1.0 / (2.0 * freq))
- return peak_ampl;
- else
- return 0.0;
+ if (peak_time - 1.0 / (2.0 * freq) < time &&
+ time < peak_time + 1.0 / (2.0 * freq))
+ return peak_ampl;
+ else
+ return 0.0;
}
static void set_fluid_bcs(double *p_f_ghost, struct simulation *sim,
const double p_f_top) {
- /* 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; */ /* REMOVED: Do not pin top physical node,
- let solver handle it */
- set_bc_neumann(p_f_ghost, sim->nz, -1, sim->phi[0] * sim->rho_f * sim->G,
- sim->dz);
+ /* correct ghost node at top BC for hydrostatic pressure gradient (uses top
+ porosity phi[nz-1] at the ice-till interface) */
+ set_bc_dirichlet(p_f_ghost, sim->nz, +1,
+ p_f_top -
+ sim->phi[sim->nz - 1] * sim->rho_f * sim->G * sim->dz);
+ p_f_ghost[sim->nz] = p_f_top; /* pin top physical node */
+ set_bc_neumann(p_f_ghost, sim->nz, -1, sim->phi[0] * sim->rho_f * sim->G,
+ sim->dz);
}
int darcy_solver_1d(struct simulation *sim) {
- int i, solved = 0;
- double p_f_top;
- double *a = sim->tdma_a;
- double *b = sim->tdma_b;
- double *c = sim->tdma_c;
- double *d = sim->tdma_d;
- double *x = sim->tdma_x;
- double *k_n = sim->darcy_k_n;
- double *phi_n = sim->darcy_phi_n;
-
- 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);
- 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);
- 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);
-
- /* set fluid BCs (1 of 2) */
- set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
- set_fluid_bcs(sim->p_f_next_ghost, sim, p_f_top);
-
- /* implicit solution with TDMA */
- {
+ int i, solved = 0;
+ double p_f_top;
+ double *a = sim->tdma_a;
+ double *b = sim->tdma_b;
+ double *c = sim->tdma_c;
+ double *d = sim->tdma_d;
+ double *x = sim->tdma_x;
+ double *k_n = sim->darcy_k_n;
+ double *phi_n = sim->darcy_phi_n;
+
+ 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);
+ 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);
+ 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);
+
+ /* set fluid BCs (1 of 2) */
+ set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
+ set_fluid_bcs(sim->p_f_next_ghost, sim, p_f_top);
+
+ /* implicit solution with TDMA */
+ {
#ifdef DEBUG
- printf("\nIMPLICIT SOLVER IN %s\n", __func__);
+ printf("\nIMPLICIT SOLVER IN %s\n", __func__);
#endif
- /* Predictor step for nonlinear coefficients */
- if (sim->transient)
- for (i = 0; i < sim->nz; ++i) {
- phi_n[i] = sim->phi[i] + sim->dt * sim->phi_dot[i];
- k_n[i] = kozeny_carman(sim->d, phi_n[i]);
- }
- else
- for (i = 0; i < sim->nz; ++i) {
- phi_n[i] = sim->phi[i];
- k_n[i] = sim->k[i];
- }
-
- /* Build Tridiagonal System */
- for (i = 0; i < sim->nz; ++i) {
- if (sim->D > 0.0) {
- /* Constant diffusivity mode */
- double coeff = sim->D * sim->dt / (sim->dz * sim->dz);
- a[i] = -coeff;
- b[i] = 1.0 + 2.0 * coeff;
- c[i] = -coeff;
- /* RHS is just p_old (plus potential source terms if applicable,
- but explicit function for D > 0 uses only diffusion term) */
- d[i] = sim->p_f_ghost[i + 1];
- } else {
- /* Coefficients calculation matches the discretized equation */
- double k_i = k_n[i];
- double k_zn, k_zp;
-
- if (i == 0)
- k_zn = k_i;
- else
- k_zn = k_n[i - 1];
-
- if (i == sim->nz - 1)
- k_zp = k_i;
- else
- k_zp = k_n[i + 1];
-
- /* Harmonic means for conductivity at faces */
- double k_harm_p = 2.0 * k_zp * k_i / fmax(k_zp + k_i, 1e-30);
- double k_harm_n = 2.0 * k_zn * k_i / fmax(k_zn + k_i, 1e-30);
-
- /* Diffusion terms */
- double coupling =
- 1.0 / ((sim->alpha + sim->beta_f * phi_n[i]) * sim->mu_f);
- double porosity_term =
- -1.0 / ((sim->alpha + sim->beta_f * phi_n[i]) * (1.0 - phi_n[i]));
-
- /* Matrix coefficients (LHS) */
- /* term: - dt * coupling * (k_n * (p_i - p_{i-1}) / dz^2) */
- double alpha_i = -sim->dt * coupling * k_harm_n /
- (sim->dz * sim->dz); // coeff for p_{i-1}
- double gamma_i = -sim->dt * coupling * k_harm_p /
- (sim->dz * sim->dz); // coeff for p_{i+1}
- double beta_i =
- 1.0 -
- (alpha_i + gamma_i); // coeff for p_i (sum of abs(off-diags) + 1)
-
- a[i] = alpha_i;
- b[i] = beta_i;
- c[i] = gamma_i;
-
- /* RHS: p_old + source term */
- d[i] =
- sim->p_f_ghost[i + 1] + sim->dt * porosity_term * sim->phi_dot[i];
- }
- }
-
- /* Apply Boundary Conditions to Linear System */
- /* Bottom (i=0): Neumann. p_{-1} = p_0 + C. */
- /* Bottom (i=0): Neumann. p_{-1} = p_0 + C. */
- double bc_neumann_val = sim->phi[0] * sim->rho_f * sim->G * sim->dz;
-
- b[0] += a[0];
- d[0] -= a[0] * bc_neumann_val;
- a[0] = 0.0;
-
- /* Top (i=nz-1): Dirichlet. p_{n} = p_ghost_top.
- Term c[n-1] * p_{n} becomes c[n-1] * p_ghost_top.
- Subtract from d[n-1].
- FIX: Use sim->p_f_ghost[sim->nz + 1] which has correct BC correction.
- */
- d[sim->nz - 1] -= c[sim->nz - 1] * sim->p_f_ghost[sim->nz + 1];
- c[sim->nz - 1] = 0.0;
-
- /* Solve */
- tridiagonal_solver(x, a, b, c, d, sim->tdma_c_prime, sim->tdma_d_prime,
- sim->nz);
-
- /* Store result in p_f_dot (rate) */
- for (i = 0; i < sim->nz; ++i) {
- sim->p_f_dot[i] = (x[i] - sim->p_f_ghost[i + 1]) / sim->dt;
- /* Store in impl array too for consistency/debug */
- sim->p_f_dot_impl[i] = sim->p_f_dot[i];
- }
-
- add_darcy_iters(1);
- solved = 1;
- }
-
- for (i = 0; i < sim->nz; ++i)
- sim->p_f_next_ghost[i + 1] =
- sim->p_f_dot[i] * sim->dt + sim->p_f_ghost[i + 1];
-
- set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
- set_fluid_bcs(sim->p_f_next_ghost, sim, p_f_top);
+ /* Predictor step for nonlinear coefficients */
+ if (sim->transient)
+ for (i = 0; i < sim->nz; ++i) {
+ phi_n[i] = sim->phi[i] + sim->dt * sim->phi_dot[i];
+ k_n[i] = kozeny_carman(sim->d, phi_n[i]);
+ }
+ else
+ for (i = 0; i < sim->nz; ++i) {
+ phi_n[i] = sim->phi[i];
+ k_n[i] = sim->k[i];
+ }
+
+ /* Build Tridiagonal System */
+ for (i = 0; i < sim->nz; ++i) {
+ if (sim->D > 0.0) {
+ /* Constant diffusivity mode */
+ double coeff = sim->D * sim->dt / (sim->dz * sim->dz);
+ a[i] = -coeff;
+ b[i] = 1.0 + 2.0 * coeff;
+ c[i] = -coeff;
+ /* RHS is just p_old (plus potential source terms if applicable,
+ but explicit function for D > 0 uses only diffusion term) */
+ d[i] = sim->p_f_ghost[i + 1];
+ } else {
+ /* Coefficients calculation matches the discretized equation */
+ double k_i = k_n[i];
+ double k_zn, k_zp;
+
+ if (i == 0)
+ k_zn = k_i;
+ else
+ k_zn = k_n[i - 1];
+
+ if (i == sim->nz - 1)
+ k_zp = k_i;
+ else
+ k_zp = k_n[i + 1];
+
+ /* Harmonic means for conductivity at faces */
+ double k_harm_p = 2.0 * k_zp * k_i / fmax(k_zp + k_i, 1e-30);
+ double k_harm_n = 2.0 * k_zn * k_i / fmax(k_zn + k_i, 1e-30);
+
+ /* Diffusion terms */
+ double coupling =
+ 1.0 / ((sim->alpha + sim->beta_f * phi_n[i]) * sim->mu_f);
+ double porosity_term =
+ -1.0 / ((sim->alpha + sim->beta_f * phi_n[i]) * (1.0 - phi_n[i]));
+
+ /* Matrix coefficients (LHS) */
+ /* term: - dt * coupling * (k_n * (p_i - p_{i-1}) / dz^2) */
+ double alpha_i = -sim->dt * coupling * k_harm_n /
+ (sim->dz * sim->dz); // coeff for p_{i-1}
+ double gamma_i = -sim->dt * coupling * k_harm_p /
+ (sim->dz * sim->dz); // coeff for p_{i+1}
+ double beta_i =
+ 1.0 -
+ (alpha_i + gamma_i); // coeff for p_i (sum of abs(off-diags) + 1)
+
+ a[i] = alpha_i;
+ b[i] = beta_i;
+ c[i] = gamma_i;
+
+ /* RHS: p_old + source term */
+ d[i] =
+ sim->p_f_ghost[i + 1] + sim->dt * porosity_term * sim->phi_dot[i];
+ }
+ }
+
+ /* Apply Boundary Conditions to Linear System */
+ /* Bottom (i=0): Neumann. p_{-1} = p_0 + C. */
+ /* Bottom (i=0): Neumann. p_{-1} = p_0 + C. */
+ double bc_neumann_val = sim->phi[0] * sim->rho_f * sim->G * sim->dz;
+
+ b[0] += a[0];
+ d[0] -= a[0] * bc_neumann_val;
+ a[0] = 0.0;
+
+ /* Top (i=nz-1): Dirichlet at the physical top node (ice-till interface):
+ p[nz-1] = p_f_top(t). Identity row so the solved top-node pressure
+ equals the imposed value exactly. */
+ a[sim->nz - 1] = 0.0;
+ b[sim->nz - 1] = 1.0;
+ c[sim->nz - 1] = 0.0;
+ d[sim->nz - 1] = p_f_top;
+
+ /* Solve */
+ tridiagonal_solver(x, a, b, c, d, sim->tdma_c_prime, sim->tdma_d_prime,
+ sim->nz);
+
+ /* Store result in p_f_dot (rate) */
+ for (i = 0; i < sim->nz; ++i) {
+ sim->p_f_dot[i] = (x[i] - sim->p_f_ghost[i + 1]) / sim->dt;
+ /* Store in impl array too for consistency/debug */
+ sim->p_f_dot_impl[i] = sim->p_f_dot[i];
+ }
+
+ add_darcy_iters(1);
+ solved = 1;
+ }
+
+ for (i = 0; i < sim->nz; ++i)
+ sim->p_f_next_ghost[i + 1] =
+ sim->p_f_dot[i] * sim->dt + sim->p_f_ghost[i + 1];
+
+ set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
+ set_fluid_bcs(sim->p_f_next_ghost, sim, p_f_top);
#ifdef DEBUG
- 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);
+ 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
- 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;
+ for (i = 0; i < sim->nz; ++i)
+ if (isnan(sim->p_f_dot_expl[i]) || isinf(sim->p_f_dot_expl[i])) {
+ fprintf(stderr, "invalid: sim->p_f_dot_expl[%d] = %g (t = %g s)\n", i,
+ sim->p_f_dot_expl[i], sim->t);
+ return 1;
+ }
+
+ for (i = 0; i < sim->nz; ++i)
+ if (isnan(sim->p_f_dot_impl[i]) || isinf(sim->p_f_dot_impl[i])) {
+ fprintf(stderr, "invalid: sim->p_f_dot_impl[%d] = %g (t = %g s)\n", i,
+ sim->p_f_dot_impl[i], sim->t);
+ return 1;
+ }
+
+ return solved - 1;
}
diff --git a/max_depth_simple_shear.c b/max_depth_simple_shear.c
@@ -229,7 +229,8 @@ main(int argc, char *argv[])
usage();
sim.nz = 2;
- prepare_arrays(&sim);
+ if (prepare_arrays(&sim))
+ errx(1, "prepare_arrays failed");
if (!isnan(new_phi))
for (i = 0; i < sim.nz; ++i)
@@ -238,7 +239,10 @@ main(int argc, char *argv[])
for (i = 0; i < sim.nz; ++i)
sim.k[i] = new_k;
- check_simulation_parameters(&sim);
+ if (check_simulation_parameters(&sim)) {
+ free_arrays(&sim);
+ exit(1);
+ }
depth = 0.0;
d_s = skin_depth(&sim);
diff --git a/simulation.c b/simulation.c
@@ -5,6 +5,7 @@
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
+#include <string.h>
/* iteration limits for solvers */
#define MAX_ITER_GRANULAR 100000
@@ -15,11 +16,11 @@
/* solver statistics for benchmarking */
struct solver_stats {
- long poisson_iters;
- long darcy_iters;
- long coupled_iters;
- long stress_iters;
- long timesteps;
+ long poisson_iters;
+ long darcy_iters;
+ long coupled_iters;
+ long stress_iters;
+ long timesteps;
};
static struct solver_stats g_stats = {0, 0, 0, 0, 0};
@@ -29,910 +30,961 @@ static struct solver_stats g_stats = {0, 0, 0, 0, 0};
/* Simulation settings */
void init_sim(struct simulation *sim) {
- int ret;
-
- ret = snprintf(sim->name, sizeof(sim->name), DEFAULT_SIMULATION_NAME);
- if (ret < 0 || (size_t)ret == sizeof(sim->name))
- err(1, "%s: could not write simulation name", __func__);
-
- sim->G = 9.81;
- sim->P_wall = 120e3;
- sim->mu_wall = 0.45;
- sim->v_x_bot = 0.0;
- sim->v_x_fix = NAN;
- sim->v_x_limit = NAN;
- sim->nz = -1; /* cell size equals grain size if negative */
-
- sim->A = 0.40; /* Loose fit to Damsgaard et al 2013 */
- sim->b = 0.9377; /* Henann and Kamrin 2016 */
- /* sim->mu_s = 0.3819; */ /* Henann and Kamrin 2016 */
- /* sim->C = 0.0; */ /* Henann and Kamrin 2016 */
- sim->mu_s = tan(DEG2RAD(22.0)); /* Damsgaard et al 2013 */
- sim->C = 0.0; /* Damsgaard et al 2013 */
- sim->phi = initval(0.25, 1);
- sim->d = 0.04; /* Damsgaard et al 2013 */
- sim->transient = 0;
- sim->phi_min = 0.20;
- sim->phi_max = 0.55;
- sim->dilatancy_constant = 4.09; /* Pailha & Pouliquen 2009 */
-
- /* Iverson et al 1997, 1998: Storglaciaren till */
- /* sim->mu_s = tan(DEG2RAD(26.3)); */
- /* sim->C = 5.0e3; */
- /* sim->phi = initval(0.22, 1); */
- /* sim->d = ??; */
-
- /* Iverson et al 1997, 1998: Two Rivers till */
- /* sim->mu_s = tan(DEG2RAD(17.8)); */
- /* sim->C = 14.0e3; */
- /* sim->phi = initval(0.37, 1); */
- /* sim->d = ??; */
-
- /* Tulaczyk et al 2000a: Upstream B till */
- /* sim->mu_s = tan(DEG2RAD(23.9)); */
- /* sim->C = 3.0e3; */
- /* sim->phi = initval(0.35, 1); */
- /* sim->d = ??; */
-
- sim->rho_s = 2.6e3; /* Damsgaard et al 2013 */
- sim->origo_z = 0.0;
- sim->L_z = 1.0;
- sim->t = 0.0;
- sim->dt = 1.0;
- sim->t_end = 1.0;
- sim->file_dt = 1.0;
- sim->n_file = 0;
- sim->fluid = 0;
- sim->rho_f = 1e3;
-
- /* Water at 20 deg C */
- /* sim->beta_f = 4.5e-10; */ /* Goren et al 2011 */
- /* sim->mu_f = 1.0-3; */ /* Goren et al 2011 */
-
- /* Water at 0 deg C */
- sim->beta_f = 3.9e-10; /* doi:10.1063/1.1679903 */
- sim->mu_f = 1.787e-3; /* Cuffey and Paterson 2010 */
-
- sim->alpha = 1e-8;
- sim->D = -1.0; /* disabled when negative */
-
- sim->k = initval(1.9e-15, 1); /* Damsgaard et al 2015 */
-
- /* Iverson et al 1994: Storglaciaren */
- /* sim->k = initval(1.3e-14, 1); */
-
- /* Engelhardt et al 1990: Upstream B */
- /* sim->k = initval(2.0e-16, 1); */
-
- /* Leeman et al 2016: Upstream B till */
- /* sim->k = initval(4.9e-17, 1); */
-
- /* 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;
- sim->p_f_mod_pulse_time = NAN;
- sim->p_f_mod_pulse_shape = 0;
-}
-
-void 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);
- exit(1);
- }
- free(sim->phi);
- free(sim->k);
-
- sim->z = linspace(sim->origo_z, 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);
- sim->sigma_n_eff = zeros(sim->nz);
- sim->sigma_n = zeros(sim->nz);
- sim->p_f_ghost = zeros(sim->nz + 2);
- sim->p_f_next_ghost = zeros(sim->nz + 2);
- sim->p_f_dot = zeros(sim->nz);
- sim->p_f_dot_expl = zeros(sim->nz);
- sim->p_f_dot_impl = zeros(sim->nz);
- sim->phi = zeros(sim->nz);
- sim->phi_c = zeros(sim->nz);
- sim->phi_dot = zeros(sim->nz);
- sim->k = zeros(sim->nz);
- sim->xi = zeros(sim->nz);
- sim->gamma_dot_p = zeros(sim->nz);
- sim->v_x = zeros(sim->nz);
- sim->d_x = zeros(sim->nz);
- sim->g_local = zeros(sim->nz);
- sim->g_ghost = zeros(sim->nz + 2);
- sim->g_r_norm = zeros(sim->nz);
- sim->I = zeros(sim->nz);
- sim->tan_psi = zeros(sim->nz);
- sim->old_val = empty(sim->nz);
- sim->tdma_a = empty(sim->nz);
- sim->tdma_b = empty(sim->nz);
- sim->tdma_c = empty(sim->nz);
- sim->tdma_d = empty(sim->nz);
- sim->tdma_x = empty(sim->nz);
- sim->tdma_c_prime = empty(sim->nz);
- sim->tdma_d_prime = empty(sim->nz);
- sim->darcy_k_n = empty(sim->nz);
- sim->darcy_phi_n = empty(sim->nz);
+ int ret;
+
+ ret = snprintf(sim->name, sizeof(sim->name), DEFAULT_SIMULATION_NAME);
+ if (ret < 0 || (size_t)ret == sizeof(sim->name))
+ err(1, "%s: could not write simulation name", __func__);
+
+ sim->G = 9.81;
+ sim->P_wall = 120e3;
+ sim->mu_wall = 0.45;
+ sim->v_x_bot = 0.0;
+ sim->v_x_fix = NAN;
+ sim->v_x_limit = NAN;
+ sim->nz = -1; /* cell size equals grain size if negative */
+
+ sim->A = 0.40; /* Loose fit to Damsgaard et al 2013 */
+ sim->b = 0.9377; /* Henann and Kamrin 2016 */
+ /* sim->mu_s = 0.3819; */ /* Henann and Kamrin 2016 */
+ /* sim->C = 0.0; */ /* Henann and Kamrin 2016 */
+ sim->mu_s = tan(DEG2RAD(22.0)); /* Damsgaard et al 2013 */
+ sim->C = 0.0; /* Damsgaard et al 2013 */
+ sim->phi = initval(0.25, 1);
+ sim->d = 0.04; /* Damsgaard et al 2013 */
+ sim->transient = 0;
+ sim->phi_min = 0.20;
+ sim->phi_max = 0.55;
+ sim->dilatancy_constant = 4.09; /* Pailha & Pouliquen 2009 */
+
+ /* Iverson et al 1997, 1998: Storglaciaren till */
+ /* sim->mu_s = tan(DEG2RAD(26.3)); */
+ /* sim->C = 5.0e3; */
+ /* sim->phi = initval(0.22, 1); */
+ /* sim->d = ??; */
+
+ /* Iverson et al 1997, 1998: Two Rivers till */
+ /* sim->mu_s = tan(DEG2RAD(17.8)); */
+ /* sim->C = 14.0e3; */
+ /* sim->phi = initval(0.37, 1); */
+ /* sim->d = ??; */
+
+ /* Tulaczyk et al 2000a: Upstream B till */
+ /* sim->mu_s = tan(DEG2RAD(23.9)); */
+ /* sim->C = 3.0e3; */
+ /* sim->phi = initval(0.35, 1); */
+ /* sim->d = ??; */
+
+ sim->rho_s = 2.6e3; /* Damsgaard et al 2013 */
+ sim->origo_z = 0.0;
+ sim->L_z = 1.0;
+ sim->t = 0.0;
+ sim->dt = 1.0;
+ sim->t_end = 1.0;
+ sim->file_dt = 1.0;
+ sim->n_file = 0;
+ sim->fluid = 0;
+ sim->rho_f = 1e3;
+
+ /* Water at 20 deg C */
+ /* sim->beta_f = 4.5e-10; */ /* Goren et al 2011 */
+ /* sim->mu_f = 1.0-3; */ /* Goren et al 2011 */
+
+ /* Water at 0 deg C */
+ sim->beta_f = 3.9e-10; /* doi:10.1063/1.1679903 */
+ sim->mu_f = 1.787e-3; /* Cuffey and Paterson 2010 */
+
+ sim->alpha = 1e-8;
+ sim->D = -1.0; /* disabled when negative */
+
+ sim->k = initval(1.9e-15, 1); /* Damsgaard et al 2015 */
+
+ /* Iverson et al 1994: Storglaciaren */
+ /* sim->k = initval(1.3e-14, 1); */
+
+ /* Engelhardt et al 1990: Upstream B */
+ /* sim->k = initval(2.0e-16, 1); */
+
+ /* Leeman et al 2016: Upstream B till */
+ /* sim->k = initval(4.9e-17, 1); */
+
+ /* 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;
+ sim->p_f_mod_pulse_time = NAN;
+ sim->p_f_mod_pulse_shape = 0;
+}
+
+int 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);
+ return 1;
+ }
+ free(sim->phi);
+ free(sim->k);
+
+ sim->z = linspace(sim->origo_z, 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);
+ sim->sigma_n_eff = zeros(sim->nz);
+ sim->sigma_n = zeros(sim->nz);
+ sim->p_f_ghost = zeros(sim->nz + 2);
+ sim->p_f_next_ghost = zeros(sim->nz + 2);
+ sim->p_f_dot = zeros(sim->nz);
+ sim->p_f_dot_expl = zeros(sim->nz);
+ sim->p_f_dot_impl = zeros(sim->nz);
+ sim->phi = zeros(sim->nz);
+ sim->phi_c = zeros(sim->nz);
+ sim->phi_dot = zeros(sim->nz);
+ sim->k = zeros(sim->nz);
+ sim->xi = zeros(sim->nz);
+ sim->gamma_dot_p = zeros(sim->nz);
+ sim->v_x = zeros(sim->nz);
+ sim->d_x = zeros(sim->nz);
+ sim->g_local = zeros(sim->nz);
+ sim->g_ghost = zeros(sim->nz + 2);
+ sim->g_r_norm = zeros(sim->nz);
+ sim->I = zeros(sim->nz);
+ sim->tan_psi = zeros(sim->nz);
+ sim->old_val = empty(sim->nz);
+ sim->tdma_a = empty(sim->nz);
+ sim->tdma_b = empty(sim->nz);
+ sim->tdma_c = empty(sim->nz);
+ sim->tdma_d = empty(sim->nz);
+ sim->tdma_x = empty(sim->nz);
+ sim->tdma_c_prime = empty(sim->nz);
+ sim->tdma_d_prime = empty(sim->nz);
+ sim->darcy_k_n = empty(sim->nz);
+ sim->darcy_phi_n = empty(sim->nz);
+
+ return 0;
}
void free_arrays(struct simulation *sim) {
- free(sim->z);
- free(sim->mu);
- free(sim->mu_c);
- free(sim->sigma_n_eff);
- free(sim->sigma_n);
- free(sim->p_f_ghost);
- free(sim->p_f_next_ghost);
- free(sim->p_f_dot);
- free(sim->p_f_dot_expl);
- free(sim->p_f_dot_impl);
- free(sim->k);
- free(sim->phi);
- free(sim->phi_c);
- free(sim->phi_dot);
- free(sim->xi);
- free(sim->gamma_dot_p);
- free(sim->v_x);
- free(sim->d_x);
- free(sim->g_local);
- free(sim->g_ghost);
- free(sim->g_r_norm);
- free(sim->I);
- free(sim->tan_psi);
- free(sim->old_val);
- free(sim->tdma_a);
- free(sim->tdma_b);
- free(sim->tdma_c);
- free(sim->tdma_d);
- free(sim->tdma_x);
- free(sim->tdma_c_prime);
- free(sim->tdma_d_prime);
- free(sim->darcy_k_n);
- free(sim->darcy_phi_n);
+ free(sim->z);
+ free(sim->mu);
+ free(sim->mu_c);
+ free(sim->sigma_n_eff);
+ free(sim->sigma_n);
+ free(sim->p_f_ghost);
+ free(sim->p_f_next_ghost);
+ free(sim->p_f_dot);
+ free(sim->p_f_dot_expl);
+ free(sim->p_f_dot_impl);
+ free(sim->k);
+ free(sim->phi);
+ free(sim->phi_c);
+ free(sim->phi_dot);
+ free(sim->xi);
+ free(sim->gamma_dot_p);
+ free(sim->v_x);
+ free(sim->d_x);
+ free(sim->g_local);
+ free(sim->g_ghost);
+ free(sim->g_r_norm);
+ free(sim->I);
+ free(sim->tan_psi);
+ free(sim->old_val);
+ free(sim->tdma_a);
+ free(sim->tdma_b);
+ free(sim->tdma_c);
+ free(sim->tdma_d);
+ free(sim->tdma_x);
+ free(sim->tdma_c_prime);
+ free(sim->tdma_d_prime);
+ free(sim->darcy_k_n);
+ free(sim->darcy_phi_n);
+}
+
+void reset_column(struct simulation *sim) {
+ int i;
+ const int nz = sim->nz;
+ const size_t nz_bytes = (size_t)nz * sizeof(double);
+ const size_t ghost_bytes = (size_t)(nz + 2) * sizeof(double);
+
+ /* zero accumulated state fields (size nz) */
+ memset(sim->mu, 0, nz_bytes);
+ memset(sim->mu_c, 0, nz_bytes);
+ memset(sim->sigma_n_eff, 0, nz_bytes);
+ memset(sim->sigma_n, 0, nz_bytes);
+ memset(sim->p_f_dot, 0, nz_bytes);
+ memset(sim->p_f_dot_expl, 0, nz_bytes);
+ memset(sim->p_f_dot_impl, 0, nz_bytes);
+ memset(sim->phi_c, 0, nz_bytes);
+ memset(sim->phi_dot, 0, nz_bytes);
+ memset(sim->xi, 0, nz_bytes);
+ memset(sim->gamma_dot_p, 0, nz_bytes);
+ memset(sim->v_x, 0, nz_bytes);
+ memset(sim->d_x, 0, nz_bytes);
+ memset(sim->g_local, 0, nz_bytes);
+ memset(sim->g_r_norm, 0, nz_bytes);
+ memset(sim->I, 0, nz_bytes);
+ memset(sim->tan_psi, 0, nz_bytes);
+
+ /* zero ghost fields (size nz + 2) */
+ memset(sim->p_f_ghost, 0, ghost_bytes);
+ memset(sim->p_f_next_ghost, 0, ghost_bytes);
+ memset(sim->g_ghost, 0, ghost_bytes);
+
+ /* restore init_sim() defaults for porosity and permeability */
+ for (i = 0; i < nz; ++i) {
+ sim->phi[i] = 0.25; /* init_sim default */
+ sim->k[i] = 1.9e-15; /* init_sim default */
+ }
+
+ /* reset per-column scalars to init_sim() defaults */
+ sim->t = 0.0;
+ sim->mu_wall = 0.45;
+ sim->v_x_bot = 0.0;
}
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;
+ 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) {
- int ret;
- char message[100];
+ int ret;
+ char message[100];
#ifdef SHOW_PARAMETERS
- printf("%30s: %.17g\n", name, value);
+ printf("%30s: %.17g\n", name, value);
#endif
- if (isnan(value)) {
- ret = snprintf(message, sizeof(message), "%s is NaN", name);
- if (ret < 0 || (size_t)ret >= sizeof(message))
- err(1, "%s: message parsing", __func__);
- warn_parameter_value(message, value, return_status);
- } else if (isinf(value)) {
- ret = snprintf(message, sizeof(message), "%s is infinite", name);
- if (ret < 0 || (size_t)ret >= sizeof(message))
- err(1, "%s: message parsing", __func__);
- warn_parameter_value(message, value, return_status);
- }
-}
-
-void check_simulation_parameters(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->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);
-
- 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->dt, &return_status);
- if (sim->dt < 0.0)
- warn_parameter_value("sim->dt is less than zero", 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);
-
- 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);
-
- 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);
-
- 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);
-
- if (sim->fluid != 0 && sim->fluid != 1)
- warn_parameter_value("sim->fluid has an invalid value", (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);
-
- 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->alpha", sim->alpha, &return_status);
- if (sim->alpha <= 0.0)
- warn_parameter_value("sim->alpha is not positive", 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);
-
- 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);
-
- 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");
- exit(return_status);
- }
+ if (isnan(value)) {
+ ret = snprintf(message, sizeof(message), "%s is NaN", name);
+ if (ret < 0 || (size_t)ret >= sizeof(message))
+ err(1, "%s: message parsing", __func__);
+ warn_parameter_value(message, value, return_status);
+ } else if (isinf(value)) {
+ ret = snprintf(message, sizeof(message), "%s is infinite", name);
+ if (ret < 0 || (size_t)ret >= sizeof(message))
+ err(1, "%s: message parsing", __func__);
+ warn_parameter_value(message, value, return_status);
+ }
+}
+
+int check_simulation_parameters(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->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);
+
+ 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->dt, &return_status);
+ if (sim->dt < 0.0)
+ warn_parameter_value("sim->dt is less than zero", 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);
+
+ 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);
+
+ 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);
+
+ 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);
+
+ if (sim->fluid != 0 && sim->fluid != 1)
+ warn_parameter_value("sim->fluid has an invalid value", (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);
+
+ 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->alpha", sim->alpha, &return_status);
+ if (sim->alpha <= 0.0)
+ warn_parameter_value("sim->alpha is not positive", 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);
+
+ 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);
+
+ 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");
+
+ return return_status;
}
void lithostatic_pressure_distribution(struct simulation *sim) {
- int i;
+ 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]);
+ 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->origo_z + sim->L_z - sim->z[i]);
}
inline static double inertia_number(double gamma_dot_p, double d,
double sigma_n_eff, double rho_s) {
- return fabs(gamma_dot_p) * d / sqrt(sigma_n_eff / rho_s);
+ return fabs(gamma_dot_p) * d / sqrt(sigma_n_eff / rho_s);
}
void compute_inertia_number(struct simulation *sim) {
- int i;
+ int i;
- for (i = 0; i < sim->nz; ++i)
- sim->I[i] =
- inertia_number(sim->gamma_dot_p[i], sim->d,
- fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN), sim->rho_s);
+ for (i = 0; i < sim->nz; ++i)
+ sim->I[i] =
+ inertia_number(sim->gamma_dot_p[i], sim->d,
+ fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN), sim->rho_s);
}
void compute_critical_state_porosity(struct simulation *sim) {
- int i;
+ int i;
- for (i = 0; i < sim->nz; ++i)
- sim->phi_c[i] = sim->phi_min + (sim->phi_max - sim->phi_min) * sim->I[i];
+ for (i = 0; i < sim->nz; ++i)
+ sim->phi_c[i] = sim->phi_min + (sim->phi_max - sim->phi_min) * sim->I[i];
}
void compute_critical_state_friction(struct simulation *sim) {
- int i;
+ int i;
- if (sim->fluid)
- for (i = 0; i < sim->nz; ++i)
- sim->mu_c[i] =
- sim->mu_wall / (fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN) /
- (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);
+ if (sim->fluid)
+ for (i = 0; i < sim->nz; ++i)
+ sim->mu_c[i] =
+ sim->mu_wall / (fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN) /
+ (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->origo_z + sim->L_z - sim->z[i]) / sim->P_wall);
}
static void compute_friction(struct simulation *sim) {
- int i;
+ int i;
- if (sim->transient)
- for (i = 0; i < sim->nz; ++i)
- sim->mu[i] = sim->mu_c[i] + sim->tan_psi[i];
- else
- for (i = 0; i < sim->nz; ++i)
- sim->mu[i] = sim->mu_c[i];
+ if (sim->transient)
+ for (i = 0; i < sim->nz; ++i)
+ sim->mu[i] = sim->mu_c[i] + sim->tan_psi[i];
+ else
+ for (i = 0; i < sim->nz; ++i)
+ sim->mu[i] = sim->mu_c[i];
}
static void compute_tan_dilatancy_angle(struct simulation *sim) {
- int i;
+ int i;
- for (i = 0; i < sim->nz; ++i)
- sim->tan_psi[i] = sim->dilatancy_constant * (sim->phi_c[i] - sim->phi[i]);
+ for (i = 0; i < sim->nz; ++i)
+ sim->tan_psi[i] = sim->dilatancy_constant * (sim->phi_c[i] - sim->phi[i]);
}
static void compute_transient_fields(struct simulation *sim) {
- int i;
+ int i;
- /* Fused loop: compute I, phi_c, and tan_psi in single pass */
- for (i = 0; i < sim->nz; ++i) {
- /* Eq. 1: Inertia number */
- sim->I[i] =
- inertia_number(sim->gamma_dot_p[i], sim->d,
- fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN), sim->rho_s);
+ /* Fused loop: compute I, phi_c, and tan_psi in single pass */
+ for (i = 0; i < sim->nz; ++i) {
+ /* Eq. 1: Inertia number */
+ sim->I[i] =
+ inertia_number(sim->gamma_dot_p[i], sim->d,
+ fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN), sim->rho_s);
- /* Eq. 2: Critical state porosity */
- sim->phi_c[i] = sim->phi_min + (sim->phi_max - sim->phi_min) * sim->I[i];
+ /* Eq. 2: Critical state porosity */
+ sim->phi_c[i] = sim->phi_min + (sim->phi_max - sim->phi_min) * sim->I[i];
- /* Eq. 5: Dilatancy angle */
- sim->tan_psi[i] = sim->dilatancy_constant * (sim->phi_c[i] - sim->phi[i]);
- }
+ /* Eq. 5: Dilatancy angle */
+ sim->tan_psi[i] = sim->dilatancy_constant * (sim->phi_c[i] - sim->phi[i]);
+ }
}
static void compute_porosity_change(struct simulation *sim) {
- int i;
+ int i;
- for (i = 0; i < sim->nz; ++i)
- sim->phi_dot[i] = sim->tan_psi[i] * sim->gamma_dot_p[i] * sim->phi[i];
+ for (i = 0; i < sim->nz; ++i)
+ sim->phi_dot[i] = sim->tan_psi[i] * sim->gamma_dot_p[i] * sim->phi[i];
}
double kozeny_carman(const double diameter, const double porosity) {
- return (diameter * diameter) / 180.0 * (porosity * porosity * porosity) /
- ((1.0 - porosity) * (1.0 - porosity));
+ return (diameter * diameter) / 180.0 * (porosity * porosity * porosity) /
+ ((1.0 - porosity) * (1.0 - porosity));
}
static void compute_permeability(struct simulation *sim) {
- int i;
+ int i;
- for (i = 0; i < sim->nz; ++i)
- sim->k[i] = kozeny_carman(sim->d, sim->phi[i]);
+ for (i = 0; i < sim->nz; ++i)
+ sim->k[i] = kozeny_carman(sim->d, sim->phi[i]);
}
static double shear_strain_rate_plastic(const double fluidity,
const double friction) {
- return fluidity * friction;
+ return fluidity * friction;
}
static void compute_shear_strain_rate_plastic(struct simulation *sim) {
- int i;
+ int i;
- for (i = 0; i < sim->nz; ++i)
- sim->gamma_dot_p[i] =
- shear_strain_rate_plastic(sim->g_ghost[i + 1], sim->mu[i]);
+ for (i = 0; i < sim->nz; ++i)
+ sim->gamma_dot_p[i] =
+ shear_strain_rate_plastic(sim->g_ghost[i + 1], sim->mu[i]);
}
static void compute_shear_velocity(struct simulation *sim) {
- int i;
+ int i;
- /* TODO: implement iterative solver for v_x from gamma_dot_p field */
- /* Dirichlet BC at bottom */
- sim->v_x[0] = sim->v_x_bot;
+ /* TODO: implement iterative solver for v_x from gamma_dot_p field */
+ /* Dirichlet BC at bottom */
+ sim->v_x[0] = sim->v_x_bot;
- for (i = 1; i < sim->nz; ++i)
- sim->v_x[i] = sim->v_x[i - 1] + sim->gamma_dot_p[i] * sim->dz;
+ for (i = 1; i < sim->nz; ++i)
+ sim->v_x[i] = sim->v_x[i - 1] + sim->gamma_dot_p[i] * sim->dz;
}
void compute_effective_stress(struct simulation *sim) {
- int i;
+ int i;
- if (sim->fluid)
- for (i = 0; i < sim->nz; ++i) {
- /* use implicit (next-step) pressure for tighter coupling */
- sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_next_ghost[i + 1];
- if (sim->sigma_n_eff[i] < 0)
- warnx("%s: sigma_n_eff[%d] is negative with value %g", __func__, i,
- sim->sigma_n_eff[i]);
- }
- else
- for (i = 0; i < sim->nz; ++i)
- sim->sigma_n_eff[i] = sim->sigma_n[i];
+ if (sim->fluid)
+ for (i = 0; i < sim->nz; ++i) {
+ /* use implicit (next-step) pressure for tighter coupling */
+ sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_next_ghost[i + 1];
+ if (sim->sigma_n_eff[i] < 0)
+ warnx("%s: sigma_n_eff[%d] is negative with value %g", __func__, i,
+ sim->sigma_n_eff[i]);
+ }
+ else
+ for (i = 0; i < sim->nz; ++i)
+ sim->sigma_n_eff[i] = sim->sigma_n[i];
}
static double cooperativity_length(const double A, const double d,
const double mu, const double p,
const double mu_s, const double C) {
- double denom = fmax(fabs((mu - C / p) - mu_s), 1e-10);
- return A * d / sqrt(denom);
+ double denom = fmax(fabs((mu - C / p) - mu_s), 1e-10);
+ return A * d / sqrt(denom);
}
static void compute_cooperativity_length(struct simulation *sim) {
- int i;
+ int i;
- for (i = 0; i < sim->nz; ++i)
- sim->xi[i] = cooperativity_length(
- sim->A, sim->d, sim->mu[i], fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN),
- sim->mu_s, sim->C);
+ for (i = 0; i < sim->nz; ++i)
+ sim->xi[i] = cooperativity_length(
+ sim->A, sim->d, sim->mu[i], fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN),
+ sim->mu_s, sim->C);
}
static double local_fluidity(const double p, const double mu, const double mu_s,
const double C, const double b, const double rho_s,
const double d) {
- if (mu - C / p <= mu_s)
- return 0.0;
- else
- return sqrt(p / (rho_s * d * d)) * ((mu - C / p) - mu_s) / (b * mu);
+ if (mu - C / p <= mu_s)
+ return 0.0;
+ else
+ return sqrt(p / (rho_s * d * d)) * ((mu - C / p) - mu_s) / (b * mu);
}
static void compute_local_fluidity(struct simulation *sim) {
- int i;
+ int i;
- for (i = 0; i < sim->nz; ++i)
- sim->g_local[i] =
- local_fluidity(fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN), sim->mu[i],
- sim->mu_s, sim->C, sim->b, sim->rho_s, sim->d);
+ for (i = 0; i < sim->nz; ++i)
+ sim->g_local[i] =
+ local_fluidity(fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN), sim->mu[i],
+ sim->mu_s, sim->C, sim->b, sim->rho_s, sim->d);
}
void set_bc_neumann(double *a, const int nz, const int boundary,
const double df, const double dx) {
- if (boundary == -1)
- a[0] = a[1] + df * dx;
- else if (boundary == +1)
- a[nz + 1] = a[nz] - df * dx;
- else
- errx(1, "%s: Unknown boundary %d\n", __func__, boundary);
+ if (boundary == -1)
+ a[0] = a[1] + df * dx;
+ else if (boundary == +1)
+ a[nz + 1] = a[nz] - df * dx;
+ else
+ errx(1, "%s: Unknown boundary %d\n", __func__, boundary);
}
void set_bc_dirichlet(double *a, const int nz, const int boundary,
const double value) {
- if (boundary == -1)
- a[0] = value;
- else if (boundary == +1)
- a[nz + 1] = value;
- else
- errx(1, "%s: Unknown boundary %d\n", __func__, boundary);
+ if (boundary == -1)
+ a[0] = value;
+ else if (boundary == +1)
+ a[nz + 1] = value;
+ else
+ errx(1, "%s: Unknown boundary %d\n", __func__, boundary);
}
double residual(double new_val, double old_val) {
- return (new_val - old_val) / (fmax(fabs(old_val), fabs(new_val)) + 1e-16);
+ return (new_val - old_val) / (fmax(fabs(old_val), fabs(new_val)) + 1e-16);
}
void tridiagonal_solver(double *x, const double *a, const double *b,
const double *c, const double *d, double *c_prime,
double *d_prime, int n) {
- /*
- * TDMA (Thomas Algorithm) solver for tridiagonal system
- * a: sub-diagonal (means a[i] is coeff for x[i-1])
- * b: main diagonal (means b[i] is coeff for x[i])
- * c: super-diagonal (means c[i] is coeff for x[i+1])
- * d: right hand side
- * x: solution vector
- * n: system size
- *
- * Note: This implementation assumes 0-indexed arrays, where:
- * Eq i (0 <= i < n): a[i]*x[i-1] + b[i]*x[i] + c[i]*x[i+1] = d[i]
- * Boundary conditions: a[0] = 0 and c[n-1] = 0 (assumed to be handled by
- * caller or zeroed)
- */
- int i;
-
- /* Forward sweep */
- c_prime[0] = c[0] / b[0];
- d_prime[0] = d[0] / b[0];
-
- for (i = 1; i < n; i++) {
- double temp = 1.0 / (b[i] - a[i] * c_prime[i - 1]);
- c_prime[i] = c[i] * temp;
- d_prime[i] = (d[i] - a[i] * d_prime[i - 1]) * temp;
- }
-
- /* Back substitution */
- x[n - 1] = d_prime[n - 1];
- for (i = n - 2; i >= 0; i--) {
- x[i] = d_prime[i] - c_prime[i] * x[i + 1];
- }
+ /*
+ * TDMA (Thomas Algorithm) solver for tridiagonal system
+ * a: sub-diagonal (means a[i] is coeff for x[i-1])
+ * b: main diagonal (means b[i] is coeff for x[i])
+ * c: super-diagonal (means c[i] is coeff for x[i+1])
+ * d: right hand side
+ * x: solution vector
+ * n: system size
+ *
+ * Note: This implementation assumes 0-indexed arrays, where:
+ * Eq i (0 <= i < n): a[i]*x[i-1] + b[i]*x[i] + c[i]*x[i+1] = d[i]
+ * Boundary conditions: a[0] = 0 and c[n-1] = 0 (assumed to be handled by
+ * caller or zeroed)
+ */
+ int i;
+
+ /* Forward sweep */
+ c_prime[0] = c[0] / b[0];
+ d_prime[0] = d[0] / b[0];
+
+ for (i = 1; i < n; i++) {
+ double temp = 1.0 / (b[i] - a[i] * c_prime[i - 1]);
+ c_prime[i] = c[i] * temp;
+ d_prime[i] = (d[i] - a[i] * d_prime[i - 1]) * temp;
+ }
+
+ /* Back substitution */
+ x[n - 1] = d_prime[n - 1];
+ for (i = n - 2; i >= 0; i--) {
+ x[i] = d_prime[i] - c_prime[i] * x[i + 1];
+ }
}
static int implicit_1d_sor_poisson_solver(struct simulation *sim) {
- /*
- * Replaced SOR solver with direct TDMA solver.
- * System: -0.5*g[i-1] + (1 + C)*g[i] - 0.5*g[i+1] = C*g_local[i]
- * where C = dz^2 / (2 * xi^2)
- */
- int i;
- int n = sim->nz;
- double *a = sim->tdma_a;
- double *b = sim->tdma_b;
- double *c = sim->tdma_c;
- double *d = sim->tdma_d;
- double *x = sim->tdma_x;
-
- /* Set up TDMA arrays */
- for (i = 0; i < n; i++) {
- double coorp_term = sim->dz * sim->dz / (2.0 * sim->xi[i] * sim->xi[i]);
-
- a[i] = -0.5;
- b[i] = 1.0 + coorp_term;
- c[i] = -0.5;
- d[i] = coorp_term * sim->g_local[i];
- }
-
- /* Boundary conditions adjustment */
- /* g[-1] = 0 (sim->g_ghost[0]) -> a[0]*0 term vanishes */
- /* g[n] = 0 (sim->g_ghost[n+1]) -> c[n-1]*0 term vanishes */
- /* But TDMA solver assumes a[0] and c[n-1] are coefficients in the matrix,
- which should be 0 for the first and last row if we strictly follow standard
- TDMA for isolated systems. However, our equation for i=0 is: -0.5*g{-1} +
- b[0]*g{0} + c[0]*g{1} = d[0] Since g{-1}=0, the term -0.5*g{-1} is 0. So
- effectively a[0]=0 in the matrix sense. Same for i=n-1: c[n-1]*g{n} is 0.
- */
- a[0] = 0.0;
- c[n - 1] = 0.0;
-
- tridiagonal_solver(x, a, b, c, d, sim->tdma_c_prime, sim->tdma_d_prime, n);
-
- /* Copy result back to ghost array (indices 1 to n) */
- set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0);
- set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0);
- for (i = 0; i < n; i++) {
- sim->g_ghost[i + 1] = x[i];
- }
-
- g_stats.poisson_iters += 1; /* Count as 1 iteration for stats */
- return 0;
+ /*
+ * Replaced SOR solver with direct TDMA solver.
+ * System: -0.5*g[i-1] + (1 + C)*g[i] - 0.5*g[i+1] = C*g_local[i]
+ * where C = dz^2 / (2 * xi^2)
+ */
+ int i;
+ int n = sim->nz;
+ double *a = sim->tdma_a;
+ double *b = sim->tdma_b;
+ double *c = sim->tdma_c;
+ double *d = sim->tdma_d;
+ double *x = sim->tdma_x;
+
+ /* Set up TDMA arrays */
+ for (i = 0; i < n; i++) {
+ double coorp_term = sim->dz * sim->dz / (2.0 * sim->xi[i] * sim->xi[i]);
+
+ a[i] = -0.5;
+ b[i] = 1.0 + coorp_term;
+ c[i] = -0.5;
+ d[i] = coorp_term * sim->g_local[i];
+ }
+
+ /* Boundary conditions adjustment */
+ /* g[-1] = 0 (sim->g_ghost[0]) -> a[0]*0 term vanishes */
+ /* g[n] = 0 (sim->g_ghost[n+1]) -> c[n-1]*0 term vanishes */
+ /* But TDMA solver assumes a[0] and c[n-1] are coefficients in the matrix,
+ which should be 0 for the first and last row if we strictly follow standard
+ TDMA for isolated systems. However, our equation for i=0 is: -0.5*g{-1} +
+ b[0]*g{0} + c[0]*g{1} = d[0] Since g{-1}=0, the term -0.5*g{-1} is 0. So
+ effectively a[0]=0 in the matrix sense. Same for i=n-1: c[n-1]*g{n} is 0.
+ */
+ a[0] = 0.0;
+ c[n - 1] = 0.0;
+
+ tridiagonal_solver(x, a, b, c, d, sim->tdma_c_prime, sim->tdma_d_prime, n);
+
+ /* Copy result back to ghost array (indices 1 to n) */
+ set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0);
+ set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0);
+ for (i = 0; i < n; i++) {
+ sim->g_ghost[i + 1] = x[i];
+ }
+
+ g_stats.poisson_iters += 1; /* Count as 1 iteration for stats */
+ return 0;
}
void write_output_file(struct simulation *sim, const int normalize) {
- int ret;
- char outfile[200];
- FILE *fp;
+ int ret;
+ char outfile[200];
+ FILE *fp;
- ret = snprintf(outfile, sizeof(outfile), "%s.output%05d.txt", sim->name,
- sim->n_file++);
- if (ret < 0 || (size_t)ret >= sizeof(outfile))
- err(1, "%s: outfile snprintf", __func__);
+ ret = snprintf(outfile, sizeof(outfile), "%s.output%05d.txt", sim->name,
+ sim->n_file++);
+ if (ret < 0 || (size_t)ret >= sizeof(outfile))
+ err(1, "%s: outfile snprintf", __func__);
- if ((fp = fopen(outfile, "w")) != NULL) {
- print_output(sim, fp, normalize);
- fclose(fp);
- } else {
- fprintf(stderr, "could not open output file: %s", outfile);
- exit(1);
- }
+ if ((fp = fopen(outfile, "w")) != NULL) {
+ print_output(sim, fp, normalize);
+ fclose(fp);
+ } else {
+ fprintf(stderr, "could not open output file: %s", outfile);
+ exit(1);
+ }
}
void print_output(struct simulation *sim, FILE *fp, const int norm) {
- int i;
- double *v_x_out;
+ int i;
+ double *v_x_out;
- if (norm)
- v_x_out = normalize(sim->v_x, sim->nz);
- else
- v_x_out = copy(sim->v_x, sim->nz);
+ if (norm)
+ v_x_out = normalize(sim->v_x, sim->nz);
+ else
+ v_x_out = copy(sim->v_x, sim->nz);
- for (i = 0; i < sim->nz; ++i)
- fprintf(fp,
- "%.17g\t%.17g\t%.17g\t"
- "%.17g\t%.17g\t%.17g\t"
- "%.17g\t%.17g\t%.17g\t%.17g"
- "\n",
- sim->z[i], v_x_out[i], sim->sigma_n_eff[i], sim->p_f_ghost[i + 1],
- sim->mu[i], sim->gamma_dot_p[i], sim->phi[i], sim->I[i],
- sim->mu[i] * fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN),
- sim->d_x[i]);
+ for (i = 0; i < sim->nz; ++i)
+ fprintf(fp,
+ "%.17g\t%.17g\t%.17g\t"
+ "%.17g\t%.17g\t%.17g\t"
+ "%.17g\t%.17g\t%.17g\t%.17g"
+ "\n",
+ sim->z[i], v_x_out[i], sim->sigma_n_eff[i], sim->p_f_ghost[i + 1],
+ sim->mu[i], sim->gamma_dot_p[i], sim->phi[i], sim->I[i],
+ sim->mu[i] * fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN),
+ sim->d_x[i]);
- free(v_x_out);
+ free(v_x_out);
}
-static void temporal_increment(struct simulation *sim) {
- int i;
+static int temporal_increment(struct simulation *sim) {
+ int i;
- if (sim->transient)
- for (i = 0; i < sim->nz; ++i)
- sim->phi[i] += sim->phi_dot[i] * sim->dt;
+ if (sim->transient)
+ 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) {
- 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;
- }
- }
+ if (sim->fluid)
+ for (i = 0; i < sim->nz; ++i) {
+ if (isnan(sim->p_f_dot[i])) {
+ fprintf(stderr, "encountered NaN at sim->p_f_dot[%d] (t = %g s)\n", i,
+ sim->t);
+ return 1;
+ } else {
+ sim->p_f_ghost[i + 1] += sim->p_f_dot[i] * sim->dt;
+ }
+ }
- for (i = 0; i < sim->nz; ++i)
- sim->d_x[i] += sim->v_x[i] * sim->dt;
- sim->t += sim->dt;
+ for (i = 0; i < sim->nz; ++i)
+ sim->d_x[i] += sim->v_x[i] * sim->dt;
+ sim->t += sim->dt;
+
+ return 0;
}
int coupled_shear_solver(struct simulation *sim, const int max_iter,
const double rel_tol) {
- int i, coupled_iter, stress_iter = 0;
- double r_norm_max, vel_res_norm = NAN, mu_wall_orig = sim->mu_wall;
-
- copy_values(sim->p_f_ghost, sim->p_f_next_ghost, sim->nz + 2);
- compute_effective_stress(sim); /* Eq. 9 */
-
- do { /* stress iterations */
- coupled_iter = 0;
- do { /* coupled iterations */
-
- if (sim->transient) {
- copy_values(sim->phi_dot, sim->old_val, sim->nz);
-
- /* Fused loop for Eqs. 1-6 */
- for (i = 0; i < sim->nz; ++i) {
- /* Eq. 1: Inertia number */
- sim->I[i] = inertia_number(sim->gamma_dot_p[i], sim->d,
- fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN),
- sim->rho_s);
-
- /* Eq. 2: Critical state porosity */
- sim->phi_c[i] =
- sim->phi_min + (sim->phi_max - sim->phi_min) * sim->I[i];
-
- /* Eq. 5: Dilatancy angle */
- sim->tan_psi[i] =
- sim->dilatancy_constant * (sim->phi_c[i] - sim->phi[i]);
-
- /* Eq. 7: Critical state friction */
- if (sim->fluid)
- sim->mu_c[i] =
- sim->mu_wall / (fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN) /
- (sim->P_wall - sim->p_f_top));
- else
- 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);
-
- /* Eq. 3: Porosity change */
- sim->phi_dot[i] = sim->tan_psi[i] * sim->gamma_dot_p[i] * sim->phi[i];
-
- /* Eq. 6: Permeability */
- sim->k[i] = kozeny_carman(sim->d, sim->phi[i]);
-
- /* Eq. 4: Friction */
- sim->mu[i] = sim->mu_c[i] + sim->tan_psi[i];
- }
- } else {
- /* Non-transient case */
- compute_critical_state_friction(sim);
- compute_friction(sim);
- }
-
- /* step 5, Eq. 13 */
- if (sim->fluid && (sim->t > 0))
- if (darcy_solver_1d(sim))
- exit(11);
-
- /* step 6 */
- compute_effective_stress(sim); /* Eq. 9 */
-
- /* step 7 */
- compute_local_fluidity(sim); /* Eq. 10 */
- compute_cooperativity_length(sim); /* Eq. 12 */
-
- /* step 8, Eq. 11 */
- if (implicit_1d_sor_poisson_solver(sim))
- exit(12);
-
- /* step 9 */
- compute_shear_strain_rate_plastic(sim); /* Eq. 8 */
- compute_inertia_number(sim); /* Eq. 1 */
- compute_shear_velocity(sim);
+ int i, coupled_iter, stress_iter = 0;
+ double r_norm_max, vel_res_norm = NAN, mu_wall_orig = sim->mu_wall;
+
+ copy_values(sim->p_f_ghost, sim->p_f_next_ghost, sim->nz + 2);
+ compute_effective_stress(sim); /* Eq. 9 */
+
+ do { /* stress iterations */
+ coupled_iter = 0;
+ do { /* coupled iterations */
+
+ if (sim->transient) {
+ copy_values(sim->phi_dot, sim->old_val, sim->nz);
+
+ /* Fused loop for Eqs. 1-6 */
+ for (i = 0; i < sim->nz; ++i) {
+ /* Eq. 1: Inertia number */
+ sim->I[i] = inertia_number(sim->gamma_dot_p[i], sim->d,
+ fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN),
+ sim->rho_s);
+
+ /* Eq. 2: Critical state porosity */
+ sim->phi_c[i] =
+ sim->phi_min + (sim->phi_max - sim->phi_min) * sim->I[i];
+
+ /* Eq. 5: Dilatancy angle */
+ sim->tan_psi[i] =
+ sim->dilatancy_constant * (sim->phi_c[i] - sim->phi[i]);
+
+ /* Eq. 7: Critical state friction */
+ if (sim->fluid)
+ sim->mu_c[i] =
+ sim->mu_wall / (fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN) /
+ (sim->P_wall - sim->p_f_top));
+ else
+ sim->mu_c[i] =
+ sim->mu_wall /
+ (1.0 + (1.0 - sim->phi[i]) * sim->rho_s * sim->G *
+ (sim->origo_z + sim->L_z - sim->z[i]) / sim->P_wall);
+
+ /* Eq. 3: Porosity change */
+ sim->phi_dot[i] = sim->tan_psi[i] * sim->gamma_dot_p[i] * sim->phi[i];
+
+ /* Eq. 6: Permeability */
+ sim->k[i] = kozeny_carman(sim->d, sim->phi[i]);
+
+ /* Eq. 4: Friction */
+ sim->mu[i] = sim->mu_c[i] + sim->tan_psi[i];
+ }
+ } else {
+ /* Non-transient case */
+ compute_critical_state_friction(sim);
+ compute_friction(sim);
+ }
+
+ /* step 5, Eq. 13 */
+ if (sim->fluid && (sim->t > 0))
+ if (darcy_solver_1d(sim))
+ return 11;
+
+ /* step 6 */
+ compute_effective_stress(sim); /* Eq. 9 */
+
+ /* step 7 */
+ compute_local_fluidity(sim); /* Eq. 10 */
+ compute_cooperativity_length(sim); /* Eq. 12 */
+
+ /* step 8, Eq. 11 */
+ if (implicit_1d_sor_poisson_solver(sim))
+ return 12;
+
+ /* step 9 */
+ compute_shear_strain_rate_plastic(sim); /* Eq. 8 */
+ compute_inertia_number(sim); /* Eq. 1 */
+ compute_shear_velocity(sim);
#ifdef DEBUG
- /* for (i = 0; i < sim->nz; ++i) { */
- for (i = sim->nz - 1; i < sim->nz; ++i) {
- printf("\nsim->t = %g s\n", sim->t);
- printf("sim->I[%d] = %g\n", i, sim->I[i]);
- printf("sim->phi_c[%d] = %g\n", i, sim->phi_c[i]);
- printf("sim->tan_psi[%d] = %g\n", i, sim->tan_psi[i]);
- printf("sim->mu_c[%d] = %g\n", i, sim->mu_c[i]);
- printf("sim->phi[%d] = %g\n", i, sim->phi[i]);
- printf("sim->phi_dot[%d] = %g\n", i, sim->phi_dot[i]);
- printf("sim->k[%d] = %g\n", i, sim->k[i]);
- printf("sim->mu[%d] = %g\n", i, sim->mu[i]);
- }
+ /* for (i = 0; i < sim->nz; ++i) { */
+ for (i = sim->nz - 1; i < sim->nz; ++i) {
+ printf("\nsim->t = %g s\n", sim->t);
+ printf("sim->I[%d] = %g\n", i, sim->I[i]);
+ printf("sim->phi_c[%d] = %g\n", i, sim->phi_c[i]);
+ printf("sim->tan_psi[%d] = %g\n", i, sim->tan_psi[i]);
+ printf("sim->mu_c[%d] = %g\n", i, sim->mu_c[i]);
+ printf("sim->phi[%d] = %g\n", i, sim->phi[i]);
+ printf("sim->phi_dot[%d] = %g\n", i, sim->phi_dot[i]);
+ printf("sim->k[%d] = %g\n", i, sim->k[i]);
+ printf("sim->mu[%d] = %g\n", i, sim->mu[i]);
+ }
#endif
- /* stable porosity change field == coupled solution converged */
- if (sim->transient) {
- for (i = 0; i < sim->nz; ++i)
- sim->g_r_norm[i] = fabs(residual(sim->phi_dot[i], sim->old_val[i]));
- r_norm_max = max_with_threshold(sim->g_r_norm, sim->nz, rel_tol);
- if (r_norm_max <= rel_tol && coupled_iter > 0)
- break;
- if (coupled_iter++ >= max_iter) {
- fprintf(stderr, "coupled_shear_solver: ");
- fprintf(stderr,
- "Transient solution did not converge "
- "after %d iterations\n",
- coupled_iter);
- fprintf(stderr, ".. Residual normalized error: %g\n", r_norm_max);
- return 1;
- }
- }
-
- } while (sim->transient);
- if (!isnan(sim->v_x_limit) || !isnan(sim->v_x_fix)) {
- if (!isnan(sim->v_x_limit)) {
- double v_ref =
- fmax(fabs(sim->v_x_limit), fabs(sim->v_x[sim->nz - 1])) + 1e-12;
- vel_res_norm = (sim->v_x_limit - sim->v_x[sim->nz - 1]) / v_ref;
- if (vel_res_norm > 0.0)
- vel_res_norm = 0.0;
- } else {
- double v_ref =
- fmax(fabs(sim->v_x_fix), fabs(sim->v_x[sim->nz - 1])) + 1e-12;
- vel_res_norm = (sim->v_x_fix - sim->v_x[sim->nz - 1]) / v_ref;
- }
- sim->mu_wall *= 1.0 + (vel_res_norm * 1e-3);
- }
- if (++stress_iter > max_iter) {
- 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);
- return 10;
- }
- } while ((!isnan(sim->v_x_fix) || !isnan(sim->v_x_limit)) &&
- fabs(vel_res_norm) > RTOL_VELOCITY);
-
- if (!isnan(sim->v_x_limit) || !isnan(sim->v_x_fix))
- sim->mu_wall = mu_wall_orig;
-
- temporal_increment(sim);
-
- g_stats.coupled_iters += coupled_iter;
- g_stats.stress_iters += stress_iter;
- g_stats.timesteps++;
-
- return 0;
+ /* stable porosity change field == coupled solution converged */
+ if (sim->transient) {
+ for (i = 0; i < sim->nz; ++i)
+ sim->g_r_norm[i] = fabs(residual(sim->phi_dot[i], sim->old_val[i]));
+ r_norm_max = max_with_threshold(sim->g_r_norm, sim->nz, rel_tol);
+ if (r_norm_max <= rel_tol && coupled_iter > 0)
+ break;
+ if (coupled_iter++ >= max_iter) {
+ fprintf(stderr, "coupled_shear_solver: ");
+ fprintf(stderr,
+ "Transient solution did not converge "
+ "after %d iterations\n",
+ coupled_iter);
+ fprintf(stderr, ".. Residual normalized error: %g\n", r_norm_max);
+ return 1;
+ }
+ }
+
+ } while (sim->transient);
+ if (!isnan(sim->v_x_limit) || !isnan(sim->v_x_fix)) {
+ if (!isnan(sim->v_x_limit)) {
+ double v_ref =
+ fmax(fabs(sim->v_x_limit), fabs(sim->v_x[sim->nz - 1])) + 1e-12;
+ vel_res_norm = (sim->v_x_limit - sim->v_x[sim->nz - 1]) / v_ref;
+ if (vel_res_norm > 0.0)
+ vel_res_norm = 0.0;
+ } else {
+ double v_ref =
+ fmax(fabs(sim->v_x_fix), fabs(sim->v_x[sim->nz - 1])) + 1e-12;
+ vel_res_norm = (sim->v_x_fix - sim->v_x[sim->nz - 1]) / v_ref;
+ }
+ sim->mu_wall *= 1.0 + (vel_res_norm * 1e-3);
+ }
+ if (++stress_iter > max_iter) {
+ 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);
+ return 10;
+ }
+ } while ((!isnan(sim->v_x_fix) || !isnan(sim->v_x_limit)) &&
+ fabs(vel_res_norm) > RTOL_VELOCITY);
+
+ if (!isnan(sim->v_x_limit) || !isnan(sim->v_x_fix))
+ sim->mu_wall = mu_wall_orig;
+
+ if (temporal_increment(sim))
+ return 13;
+
+ g_stats.coupled_iters += coupled_iter;
+ g_stats.stress_iters += stress_iter;
+ g_stats.timesteps++;
+
+ return 0;
}
double find_flux(const struct simulation *sim) {
- int i;
- double flux = 0.0;
+ int i;
+ double flux = 0.0;
- for (i = 1; i < sim->nz; ++i)
- flux += (sim->v_x[i] + sim->v_x[i - 1]) / 2.0 * sim->dz;
+ for (i = 1; i < sim->nz; ++i)
+ flux += (sim->v_x[i] + sim->v_x[i - 1]) / 2.0 * sim->dz;
- return flux;
+ return flux;
}
void set_coupled_fluid_transient_timestep(struct simulation *sim,
const double safety) {
- double max_gamma_dot, mu, xi, max_I, dt;
+ double max_gamma_dot, mu, xi, max_I, dt;
- /* max expected strain rate */
- max_gamma_dot = 1.0 / sim->d;
- if (!isnan(sim->v_x_fix))
- max_gamma_dot = sim->v_x_fix / sim->dz;
- if (!isnan(sim->v_x_limit))
- max_gamma_dot = sim->v_x_limit / sim->dz;
+ /* max expected strain rate */
+ max_gamma_dot = 1.0 / sim->d;
+ if (!isnan(sim->v_x_fix))
+ max_gamma_dot = sim->v_x_fix / sim->dz;
+ if (!isnan(sim->v_x_limit))
+ max_gamma_dot = sim->v_x_limit / sim->dz;
- /* estimate for shear friction */
- mu = (sim->mu_wall / ((sim->sigma_n[sim->nz - 1] - sim->p_f_mod_ampl) /
- (sim->P_wall - sim->p_f_top))) +
- sim->dilatancy_constant * sim->phi[sim->nz - 1];
+ /* estimate for shear friction */
+ mu = (sim->mu_wall / ((sim->sigma_n[sim->nz - 1] - sim->p_f_mod_ampl) /
+ (sim->P_wall - sim->p_f_top))) +
+ sim->dilatancy_constant * sim->phi[sim->nz - 1];
- /* estimate for cooperativity length */
- xi = cooperativity_length(sim->A, sim->d, mu,
- (sim->sigma_n[sim->nz - 1] - sim->p_f_mod_ampl),
- sim->mu_s, sim->C);
+ /* estimate for cooperativity length */
+ xi = cooperativity_length(sim->A, sim->d, mu,
+ (sim->sigma_n[sim->nz - 1] - sim->p_f_mod_ampl),
+ sim->mu_s, sim->C);
- /* max expected inertia number */
- max_I =
- inertia_number(max_gamma_dot, sim->d,
- sim->sigma_n[sim->nz - 1] - sim->p_f_mod_ampl, sim->rho_s);
+ /* max expected inertia number */
+ max_I =
+ inertia_number(max_gamma_dot, sim->d,
+ sim->sigma_n[sim->nz - 1] - sim->p_f_mod_ampl, sim->rho_s);
- dt = xi * (sim->alpha + sim->phi[sim->nz - 1] * sim->beta_f) * sim->mu_f /
- (sim->phi[sim->nz - 1] * sim->phi[sim->nz - 1] * sim->phi[sim->nz - 1] *
- sim->L_z * max_I);
+ dt = xi * (sim->alpha + sim->phi[sim->nz - 1] * sim->beta_f) * sim->mu_f /
+ (sim->phi[sim->nz - 1] * sim->phi[sim->nz - 1] * sim->phi[sim->nz - 1] *
+ sim->L_z * max_I);
- if (sim->dt > safety * dt)
- sim->dt = safety * dt;
+ if (sim->dt > safety * dt)
+ sim->dt = safety * dt;
}
void print_solver_stats(FILE *fp) {
- fprintf(fp,
- "solver_stats: timesteps=%ld poisson_iters=%ld "
- "darcy_iters=%ld coupled_iters=%ld stress_iters=%ld\n",
- g_stats.timesteps, g_stats.poisson_iters, g_stats.darcy_iters,
- g_stats.coupled_iters, g_stats.stress_iters);
+ fprintf(fp,
+ "solver_stats: timesteps=%ld poisson_iters=%ld "
+ "darcy_iters=%ld coupled_iters=%ld stress_iters=%ld\n",
+ g_stats.timesteps, g_stats.poisson_iters, g_stats.darcy_iters,
+ g_stats.coupled_iters, g_stats.stress_iters);
}
void reset_solver_stats(void) {
- g_stats.poisson_iters = 0;
- g_stats.darcy_iters = 0;
- g_stats.coupled_iters = 0;
- g_stats.stress_iters = 0;
- g_stats.timesteps = 0;
+ g_stats.poisson_iters = 0;
+ g_stats.darcy_iters = 0;
+ g_stats.coupled_iters = 0;
+ g_stats.stress_iters = 0;
+ g_stats.timesteps = 0;
}
void add_darcy_iters(int iters) { g_stats.darcy_iters += iters; }
diff --git a/simulation.h b/simulation.h
@@ -12,137 +12,147 @@ extern struct simulation sim;
/* Simulation settings */
struct simulation {
- /* simulation name to use for output files */
- char name[100];
+ /* simulation name to use for output files */
+ char name[100];
- /* gravitational acceleration magnitude [m/s^2] */
- double G;
+ /* gravitational acceleration magnitude [m/s^2] */
+ double G;
- /* normal stress from the top wall [Pa] */
- double P_wall;
-
- /* optionally fix top shear velocity to this value [m/s] */
- double v_x_fix;
-
- /* optionally fix top shear velocity to this value [m/s] */
- double v_x_limit;
-
- /* bottom velocity along x [m/s] */
- double v_x_bot;
-
- /* stress ratio at top wall */
- double mu_wall;
-
- /* nonlocal amplitude [-] */
- double A;
-
- /* rate dependence beyond yield [-] */
- double b;
-
- /* bulk and critical state static yield friction coefficient [-] */
- double mu_s;
-
- /* material cohesion [Pa] */
- double C;
-
- /* representative grain size [m] */
- double d; /* ohlala */
-
- /* grain material density [kg/m^3] */
- double rho_s;
-
- /* nodes along z */
- int nz;
-
- /* origo of axis [m] */
- double origo_z;
-
- /* length of domain [m] */
- double L_z;
-
- /* array of cell coordinates */
- double *z;
-
- /* cell spacing [m] */
- double dz;
-
- /* 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;
-
- double transient;
- double phi_min;
- double phi_max;
- double dilatancy_constant;
-
- /* 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 p_f_mod_pulse_time; /* single pressure pulse at this time [s] */
- int p_f_mod_pulse_shape; /* waveform for fluid-pressure pulse */
- double beta_f; /* adiabatic fluid compressibility [Pa^-1] */
- double alpha; /* adiabatic grain compressibility [Pa^-1] */
- double mu_f; /* fluid dynamic viscosity [Pa*s] */
- double rho_f; /* fluid density [kg/m^3] */
- double D; /* diffusivity [m^2/s], overrides k, beta_f, alpha, mu_f */
-
- /* arrays */
- double *mu; /* static yield friction [-] */
- double *mu_c; /* critical-state static yield friction [-] */
- double *sigma_n_eff; /* effective normal pressure [Pa] */
- double *sigma_n; /* normal stress [Pa] */
- double *p_f_ghost; /* fluid pressure [Pa] */
- double *p_f_next_ghost; /* fluid pressure for next iteration [Pa] */
- double *p_f_dot; /* fluid pressure change [Pa/s] */
- double *p_f_dot_expl; /* fluid pressure change (explicit solution) [Pa/s] */
- double *p_f_dot_impl; /* fluid pressure change (implicit solution) [Pa/s] */
-
- double *k; /* hydraulic permeability [m^2] */
- double *phi; /* porosity [-] */
- double *phi_c; /* critical-state porosity [-] */
- double *phi_dot; /* porosity change [s^-1] */
- double *xi; /* cooperativity length */
- double *gamma_dot_p; /* plastic shear strain rate [s^-1] */
- double *v_x; /* shear velocity [m/s] */
- double *d_x; /* cumulative shear displacement [m] */
- double *g_local; /* local fluidity */
- double *g_ghost; /* fluidity with ghost nodes */
- double *g_r_norm; /* normalized residual of fluidity field */
- double *I; /* inertia number [-] */
- double *tan_psi; /* tan(dilatancy_angle) [-] */
- double *old_val; /* temporary storage for iterative solvers */
-
- /* Persistent solver workspace (size nz unless noted otherwise), allocated
- * once in prepare_arrays() and reused by Darcy/Poisson TDMA paths. */
- double *tdma_a; /* shared TDMA sub-diagonal coefficients */
- double *tdma_b; /* shared TDMA diagonal coefficients */
- double *tdma_c; /* shared TDMA super-diagonal coefficients */
- double *tdma_d; /* shared TDMA RHS coefficients */
- double *tdma_x; /* shared TDMA solution vector */
- double *tdma_c_prime; /* shared TDMA forward-sweep scratch */
- double *tdma_d_prime; /* shared TDMA forward-sweep scratch */
- double *darcy_k_n; /* Darcy predictor permeability workspace */
- double *darcy_phi_n; /* Darcy predictor porosity workspace */
+ /* normal stress from the top wall [Pa] */
+ double P_wall;
+
+ /* optionally fix top shear velocity to this value [m/s] */
+ double v_x_fix;
+
+ /* optionally fix top shear velocity to this value [m/s] */
+ double v_x_limit;
+
+ /* bottom velocity along x [m/s] */
+ double v_x_bot;
+
+ /* stress ratio at top wall */
+ double mu_wall;
+
+ /* nonlocal amplitude [-] */
+ double A;
+
+ /* rate dependence beyond yield [-] */
+ double b;
+
+ /* bulk and critical state static yield friction coefficient [-] */
+ double mu_s;
+
+ /* material cohesion [Pa] */
+ double C;
+
+ /* representative grain size [m] */
+ double d; /* ohlala */
+
+ /* grain material density [kg/m^3] */
+ double rho_s;
+
+ /* nodes along z */
+ int nz;
+
+ /* origo of axis [m] */
+ double origo_z;
+
+ /* length of domain [m] */
+ double L_z;
+
+ /* array of cell coordinates */
+ double *z;
+
+ /* cell spacing [m] */
+ double dz;
+
+ /* 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;
+
+ double transient;
+ double phi_min;
+ double phi_max;
+ double dilatancy_constant;
+
+ /* 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 p_f_mod_pulse_time; /* single pressure pulse at this time [s] */
+ int p_f_mod_pulse_shape; /* waveform for fluid-pressure pulse */
+ double beta_f; /* adiabatic fluid compressibility [Pa^-1] */
+ double alpha; /* adiabatic grain compressibility [Pa^-1] */
+ double mu_f; /* fluid dynamic viscosity [Pa*s] */
+ double rho_f; /* fluid density [kg/m^3] */
+ double D; /* diffusivity [m^2/s], overrides k, beta_f, alpha, mu_f */
+
+ /* arrays */
+ double *mu; /* static yield friction [-] */
+ double *mu_c; /* critical-state static yield friction [-] */
+ double *sigma_n_eff; /* effective normal pressure [Pa] */
+ double *sigma_n; /* normal stress [Pa] */
+ double *p_f_ghost; /* fluid pressure [Pa] */
+ double *p_f_next_ghost; /* fluid pressure for next iteration [Pa] */
+ double *p_f_dot; /* fluid pressure change [Pa/s] */
+ double *p_f_dot_expl; /* fluid pressure change (explicit solution) [Pa/s] */
+ double *p_f_dot_impl; /* fluid pressure change (implicit solution) [Pa/s] */
+
+ double *k; /* hydraulic permeability [m^2] */
+ double *phi; /* porosity [-] */
+ double *phi_c; /* critical-state porosity [-] */
+ double *phi_dot; /* porosity change [s^-1] */
+ double *xi; /* cooperativity length */
+ double *gamma_dot_p; /* plastic shear strain rate [s^-1] */
+ double *v_x; /* shear velocity [m/s] */
+ double *d_x; /* cumulative shear displacement [m] */
+ double *g_local; /* local fluidity */
+ double *g_ghost; /* fluidity with ghost nodes */
+ double *g_r_norm; /* normalized residual of fluidity field */
+ double *I; /* inertia number [-] */
+ double *tan_psi; /* tan(dilatancy_angle) [-] */
+ double *old_val; /* temporary storage for iterative solvers */
+
+ /* Persistent solver workspace (size nz unless noted otherwise), allocated
+ * once in prepare_arrays() and reused by Darcy/Poisson TDMA paths. */
+ double *tdma_a; /* shared TDMA sub-diagonal coefficients */
+ double *tdma_b; /* shared TDMA diagonal coefficients */
+ double *tdma_c; /* shared TDMA super-diagonal coefficients */
+ double *tdma_d; /* shared TDMA RHS coefficients */
+ double *tdma_x; /* shared TDMA solution vector */
+ double *tdma_c_prime; /* shared TDMA forward-sweep scratch */
+ double *tdma_d_prime; /* shared TDMA forward-sweep scratch */
+ double *darcy_k_n; /* Darcy predictor permeability workspace */
+ double *darcy_phi_n; /* Darcy predictor porosity workspace */
};
void init_sim(struct simulation *sim);
-void prepare_arrays(struct simulation *sim);
+
+/* returns: 0 ok, 1 grid size (nz) is less than 2 */
+int prepare_arrays(struct simulation *sim);
+
void free_arrays(struct simulation *sim);
-void check_simulation_parameters(struct simulation *sim);
+
+/* Reset per-column state (stresses, velocities, fluidity, porosity, time)
+ * without reallocating; geometry (nz, L_z, z, dz) and workspaces are kept.
+ * For reuse of one prepared simulation across many independent columns. */
+void reset_column(struct simulation *sim);
+
+/* returns: 0 ok, 1 one or more parameters are invalid */
+int check_simulation_parameters(struct simulation *sim);
void lithostatic_pressure_distribution(struct simulation *sim);
void compute_effective_stress(struct simulation *sim);
@@ -158,6 +168,9 @@ double kozeny_carman(const double diameter, const double porosity);
void write_output_file(struct simulation *sim, const int normalize);
void print_output(struct simulation *sim, FILE *fp, const int normalize);
+/* returns: 0 ok, 1 transient solution not converged, 10 stress solution not
+ * converged, 11 fluid (Darcy) solver failed, 12 fluidity (Poisson) solver
+ * failed, 13 NaN in temporal increment */
int coupled_shear_solver(struct simulation *sim, const int max_iter,
const double rel_tol);
