commit 2ac6bc191520c22444978580b39a69749d6bb804
parent 9bfaf834812f4c720008acdaba8c75db2f254915
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Thu, 16 Apr 2020 12:56:42 +0200
Move simulation struct out of global scope for main program
Diffstat:
| M | 1d_fd_simple_shear.c | | | 28 | ++++++++++++++-------------- |
| M | fluid.c | | | 161 | ++++++++++++++++++++++++++++++++++++++++--------------------------------------- |
| M | fluid.h | | | 10 | +++++----- |
| M | parameter_defaults.h | | | 117 | ++++++++++++++++++++++++++++++++++++++++--------------------------------------- |
| M | simulation.c | | | 666 | ++++++++++++++++++++++++++++++++++++++++--------------------------------------- |
| M | simulation.h | | | 43 | +++++++++++++++++++++++-------------------- |
6 files changed, 517 insertions(+), 508 deletions(-)
diff --git a/1d_fd_simple_shear.c b/1d_fd_simple_shear.c
@@ -2,10 +2,11 @@
#include <math.h>
#include <string.h>
#include <time.h>
+#include <unistd.h>
#include "simulation.h"
-#include "fluid.h"
#include "arg.h"
+#include "fluid.h"
#include "parameter_defaults.h"
@@ -17,7 +18,6 @@
/*#define BENCHMARK_PERFORMANCE*/
char *argv0;
-struct simulation sim;
static void
usage(void)
@@ -69,6 +69,7 @@ main(int argc, char* argv[])
int i, normalize;
unsigned long iter;
double new_phi, new_k, filetimeclock;
+ struct simulation sim;
#ifdef BENCHMARK_PERFORMANCE
clock_t t_begin, t_end;
double t_elapsed;
@@ -81,9 +82,7 @@ main(int argc, char* argv[])
}
#endif
- atexit(free_arrays);
-
- init_sim();
+ init_sim(&sim);
normalize = 0;
new_phi = sim.phi[0];
new_k = sim.k[0];
@@ -223,7 +222,7 @@ main(int argc, char* argv[])
if (sim.nz < 1)
sim.nz = (int)ceil(sim.L_z/sim.d);
- prepare_arrays();
+ prepare_arrays(&sim);
if (!isnan(new_phi))
for (i=0; i<sim.nz; ++i)
@@ -232,15 +231,15 @@ main(int argc, char* argv[])
for (i=0; i<sim.nz; ++i)
sim.k[i] = new_k;
- lithostatic_pressure_distribution();
+ lithostatic_pressure_distribution(&sim);
if (sim.fluid) {
- hydrostatic_fluid_pressure_distribution();
- set_largest_fluid_timestep(0.5);
+ hydrostatic_fluid_pressure_distribution(&sim);
+ set_largest_fluid_timestep(&sim, 0.5);
}
- compute_effective_stress(); /* Eq. 9 */
+ compute_effective_stress(&sim);
- check_simulation_parameters();
+ check_simulation_parameters(&sim);
filetimeclock = 0.0;
iter = 0;
@@ -250,7 +249,7 @@ main(int argc, char* argv[])
t_begin = clock();
#endif
- if (coupled_shear_solver(MAX_ITER_1D_FD_SIMPLE_SHEAR, RTOL))
+ if (coupled_shear_solver(&sim, MAX_ITER_1D_FD_SIMPLE_SHEAR, RTOL))
exit(10);
#ifdef BENCHMARK_PERFORMANCE
@@ -266,14 +265,15 @@ main(int argc, char* argv[])
if ((filetimeclock - sim.dt >= sim.file_dt || iter == 1) &&
strncmp(sim.name, DEFAULT_SIMULATION_NAME,
sizeof(DEFAULT_SIMULATION_NAME)) != 0) {
- write_output_file(normalize);
+ write_output_file(&sim, normalize);
filetimeclock = 0.0;
if (iter == 1)
sim.t = 0.0;
}
}
- print_output(stdout, normalize);
+ print_output(&sim, stdout, normalize);
+ free_arrays(&sim);
return 0;
}
diff --git a/fluid.c b/fluid.c
@@ -6,12 +6,12 @@
/* #define DEBUG true */
void
-hydrostatic_fluid_pressure_distribution()
+hydrostatic_fluid_pressure_distribution(struct simulation *sim)
{
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->L_z - sim->z[i]);
}
/* Determines the largest time step for the current simulation state. Note
@@ -19,15 +19,15 @@ hydrostatic_fluid_pressure_distribution()
* (i.e., permeabilities, porosities, viscosities, or compressibilities)
* change. The safety factor should be in ]0;1] */
void
-set_largest_fluid_timestep(const double safety)
+set_largest_fluid_timestep(struct simulation *sim, const double safety)
{
int i;
double dx_min, diff, diff_max;
double *dx;
- dx = spacing(sim.z, sim.nz);
+ dx = spacing(sim->z, sim->nz);
dx_min = INFINITY;
- for (i=0; i<sim.nz-1; ++i) {
+ 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);
@@ -39,14 +39,14 @@ set_largest_fluid_timestep(const double safety)
free(dx);
diff_max = -INFINITY;
- for (i=0; i<sim.nz; ++i) {
- diff = sim.k[i]/(sim.phi[i]*sim.beta_f*sim.mu_f);
+ for (i=0; i<sim->nz; ++i) {
+ diff = sim->k[i]/(sim->phi[i]*sim->beta_f*sim->mu_f);
if (diff > diff_max) diff_max = diff;
}
- 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 = safety*0.5*dx_min*dx_min/diff_max;
+ if (sim->file_dt*0.5 < sim->dt)
+ sim->dt = sim->file_dt;
}
static double
@@ -86,15 +86,15 @@ square_pulse(const double time,
}
static void
-set_fluid_bcs(const double p_f_top)
+set_fluid_bcs(struct simulation *sim, const double p_f_top)
{
- set_bc_dirichlet(sim.p_f_ghost, sim.nz, +1, p_f_top);
- sim.p_f_ghost[sim.nz] = p_f_top; /* Include top node in BC */
- set_bc_neumann(sim.p_f_ghost,
- sim.nz,
+ set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top);
+ sim->p_f_ghost[sim->nz] = p_f_top; /* Include top node in BC */
+ set_bc_neumann(sim->p_f_ghost,
+ sim->nz,
-1,
- sim.phi[0]*sim.rho_f*sim.G,
- sim.dz);
+ sim->phi[0]*sim->rho_f*sim->G,
+ sim->dz);
}
static double
@@ -135,11 +135,12 @@ darcy_pressure_change_1d(const int i,
/* TODO: add advective term */
return 1.0/(beta_f*phi[i]*mu_f)*div_k_grad_p
- - 1.0/(beta_f*(1.0 - phi[i]))*phi_dot[i];
+ - 1.0/(beta_f*(1.0 - phi[i]))*phi_dot[i];
}
int
-darcy_solver_1d(const int max_iter,
+darcy_solver_1d(struct simulation *sim,
+ const int max_iter,
const double rel_tol)
{
int i, iter, solved;
@@ -162,104 +163,104 @@ darcy_solver_1d(const int max_iter,
/* TODO: values other than 1.0 do not work! */
theta = 1.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);
+ 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);
+ 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);
+ 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(p_f_top);
+ set_fluid_bcs(sim, p_f_top);
solved = 0;
if (epsilon < 1.0) {
/* compute explicit solution to pressure change */
- dp_f_dt_expl = zeros(sim.nz);
- for (i=0; i<sim.nz; ++i)
+ dp_f_dt_expl = zeros(sim->nz);
+ for (i=0; i<sim->nz; ++i)
dp_f_dt_expl[i] = darcy_pressure_change_1d(i,
- sim.nz,
- sim.p_f_ghost,
- sim.phi,
- sim.phi_dot,
- sim.k,
- sim.dz,
- sim.beta_f,
- sim.mu_f);
+ sim->nz,
+ sim->p_f_ghost,
+ sim->phi,
+ sim->phi_dot,
+ sim->k,
+ sim->dz,
+ sim->beta_f,
+ sim->mu_f);
}
if (epsilon > 0.0) {
/* compute implicit solution with Jacobian iterations */
- dp_f_dt_impl = zeros(sim.nz);
- p_f_ghost_new = zeros(sim.nz+2);
- r_norm = zeros(sim.nz);
- p_f_ghost_old = empty(sim.nz+2);
- copy_values(sim.p_f_ghost, p_f_ghost_old, sim.nz+2);
+ dp_f_dt_impl = zeros(sim->nz);
+ p_f_ghost_new = zeros(sim->nz+2);
+ r_norm = zeros(sim->nz);
+ p_f_ghost_old = empty(sim->nz+2);
+ copy_values(sim->p_f_ghost, p_f_ghost_old, sim->nz+2);
for (iter=0; iter<max_iter; ++iter) {
/* set fluid BCs (2 of 2) */
- set_fluid_bcs( p_f_top);
+ set_fluid_bcs(sim, p_f_top);
#ifdef DEBUG
puts(".. p_f_ghost after BC:");
- print_array(sim.p_f_ghost, sim.nz+2);
+ print_array(sim->p_f_ghost, sim->nz+2);
#endif
- for (i=0; i<sim.nz-1; ++i)
+ for (i=0; i<sim->nz-1; ++i)
dp_f_dt_impl[i] = darcy_pressure_change_1d(i,
- sim.nz,
- sim.p_f_ghost,
- sim.phi,
- sim.phi_dot,
- sim.k,
- sim.dz,
- sim.beta_f,
- sim.mu_f);
-
- for (i=0; i<sim.nz-1; ++i) {
+ sim->nz,
+ sim->p_f_ghost,
+ sim->phi,
+ sim->phi_dot,
+ sim->k,
+ sim->dz,
+ sim->beta_f,
+ sim->mu_f);
+
+ for (i=0; i<sim->nz-1; ++i) {
#ifdef DEBUG
printf("dp_f_expl[%d] = %g\ndp_f_impl[%d] = %g\n",
i, dp_f_dt_expl[i], i, dp_f_dt_impl[i]);
#endif
p_f_ghost_new[i+1] = p_f_ghost_old[i+1]
- + epsilon*dp_f_dt_impl[i]*sim.dt;
+ + epsilon*dp_f_dt_impl[i]*sim->dt;
if (epsilon < 1.0)
p_f_ghost_new[i+1] += (1.0 - epsilon)
- *dp_f_dt_expl[i]*sim.dt;
+ *dp_f_dt_expl[i]*sim->dt;
p_f_ghost_new[i+1] = p_f_ghost_old[i+1]*(1.0 - theta)
+ p_f_ghost_new[i+1]*theta;
- r_norm[i] = fabs((p_f_ghost_new[i+1] - sim.p_f_ghost[i+1])
- /(sim.p_f_ghost[i+1] + 1e-16));
+ r_norm[i] = fabs((p_f_ghost_new[i+1] - sim->p_f_ghost[i+1])
+ /(sim->p_f_ghost[i+1] + 1e-16));
}
- r_norm_max = max(r_norm, sim.nz);
+ r_norm_max = max(r_norm, sim->nz);
#ifdef DEBUG
puts(".. p_f_ghost_new:");
- print_array(p_f_ghost_new, sim.nz+2);
+ print_array(p_f_ghost_new, sim->nz+2);
#endif
- copy_values(p_f_ghost_new, sim.p_f_ghost, sim.nz+2);
+ copy_values(p_f_ghost_new, sim->p_f_ghost, sim->nz+2);
#ifdef DEBUG
puts(".. p_f_ghost after update:");
- print_array(sim.p_f_ghost, sim.nz+2);
+ print_array(sim->p_f_ghost, sim->nz+2);
#endif
if (r_norm_max <= rel_tol) {
@@ -281,17 +282,17 @@ darcy_solver_1d(const int max_iter,
fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max);
}
} else {
- for (i=0; i<sim.nz; ++i)
- sim.p_f_ghost[i+1] += dp_f_dt_expl[i]*sim.dt;
+ for (i=0; i<sim->nz; ++i)
+ sim->p_f_ghost[i+1] += dp_f_dt_expl[i]*sim->dt;
solved = 1;
#ifdef DEBUG
puts(".. dp_f_dt_expl:");
- print_array(dp_f_dt_expl, sim.nz);
+ print_array(dp_f_dt_expl, sim->nz);
puts(".. p_f_ghost after explicit solution:");
- print_array(sim.p_f_ghost, sim.nz+2);
+ print_array(sim->p_f_ghost, sim->nz+2);
#endif
}
- set_fluid_bcs(p_f_top);
+ set_fluid_bcs(sim, p_f_top);
if (epsilon < 1.0)
free(dp_f_dt_expl);
diff --git a/fluid.h b/fluid.h
@@ -5,12 +5,12 @@
extern struct simulation sim;
-void hydrostatic_fluid_pressure_distribution();
+void hydrostatic_fluid_pressure_distribution(struct simulation *sim);
-void set_largest_fluid_timestep(const double safety);
+void set_largest_fluid_timestep(struct simulation *sim, const double safety);
-int darcy_solver_1d(
- const int max_iter,
- const double rel_tol);
+int darcy_solver_1d(struct simulation *sim,
+ const int max_iter,
+ const double rel_tol);
#endif
diff --git a/parameter_defaults.h b/parameter_defaults.h
@@ -1,5 +1,5 @@
-#ifndef ONED_FD_SIMPLE_SHEAR_
-#define ONED_FD_SIMPLE_SHEAR_
+#ifndef PARAMETER_DEFAULTS_H_
+#define PARAMETER_DEFAULTS_H_
#include <math.h>
#include <stdio.h>
@@ -10,106 +10,107 @@
/* Simulation settings */
void
-init_sim()
+init_sim(struct simulation *sim)
{
- snprintf(sim.name, sizeof(sim.name), DEFAULT_SIMULATION_NAME);
+ snprintf(sim->name, sizeof(sim->name), DEFAULT_SIMULATION_NAME);
- sim.G = 9.81;
+ 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->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->nz = -1; /* cell size equals grain size if negative */
/* lower values of A mean that the velocity curve can have sharper curves,
* e.g. at the transition from μ ≈ μ_s */
- sim.A = 0.40; /* Loose fit to Damsgaard et al 2013 */
+ sim->A = 0.40; /* Loose fit to Damsgaard et al 2013 */
/* lower values of b mean larger shear velocity for a given stress ratio
* above yield */
- sim.b = 0.9377; /* Henann and Kamrin 2016 */
+ sim->b = 0.9377; /* Henann and Kamrin 2016 */
/* Henann and Kamrin 2016 */
- /* sim.mu_s = 0.3819; */
- /* sim.C = 0.0; */
+ /* sim->mu_s = 0.3819; */
+ /* sim->C = 0.0; */
/* Damsgaard et al 2013 */
- sim.mu_s = tan(DEG2RAD(22.0));
- sim.C = 0.0;
- sim.phi = initval(0.25, 1);
- sim.d = 0.04; /* Damsgaard et al 2013 */
+ sim->mu_s = tan(DEG2RAD(22.0));
+ sim->C = 0.0;
+ sim->phi = initval(0.25, 1);
+ sim->d = 0.04; /* Damsgaard et al 2013 */
- sim.phi_min = 0.25;
- sim.phi_max = 0.55;
- sim.dilatancy_angle = 1.0;
+ sim->transient = 0;
+ sim->phi_min = 0.25;
+ sim->phi_max = 0.55;
+ sim->dilatancy_angle = 1.0;
/* 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 = ??; */
+ /* 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 = ??; */
+ /* 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->mu_s = tan(DEG2RAD(23.9)); */
+ /* sim->C = 3.0e3; */
+ /* sim->phi = initval(0.35, 1); */
+ /* sim->d = ??; */
/* grain material density [kg/m^3] */
- sim.rho_s = 2.6e3; /* Damsgaard et al 2013 */
+ sim->rho_s = 2.6e3; /* Damsgaard et al 2013 */
/* spatial settings */
- sim.origo_z = 0.0;
- sim.L_z = 1.0;
+ sim->origo_z = 0.0;
+ sim->L_z = 1.0;
/* temporal settings */
- sim.t = 0.0;
- sim.dt = 2.0;
- sim.t_end = 1.0;
- sim.file_dt = 0.1;
- sim.n_file = 0;
+ sim->t = 0.0;
+ sim->dt = 2.0;
+ sim->t_end = 1.0;
+ sim->file_dt = 0.1;
+ sim->n_file = 0;
/* fluid settings */
- sim.fluid = 0;
+ sim->fluid = 0;
- sim.rho_f = 1e3;
+ sim->rho_f = 1e3;
/* Water at 20 deg C, Goren et al 2011 */
- /* sim.beta_f = 4.5e-10; */
- /* sim.mu_f = 1.0-3; */
+ /* sim->beta_f = 4.5e-10; */
+ /* sim->mu_f = 1.0-3; */
/* 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->beta_f = 3.9e-10; /* doi:10.1063/1.1679903 */
+ sim->mu_f = 1.787e-3; /* Cuffey and Paterson 2010 */
/* Damsgaard et al 2015 */
- sim.k = initval(1.9e-15, 1);
+ sim->k = initval(1.9e-15, 1);
/* Iverson et al 1994: Storglaciaren */
- /* sim.k = initval(1.3e-14, 1); */
+ /* sim->k = initval(1.3e-14, 1); */
/* Engelhardt et al 1990: Upstream B */
- /* sim.k = initval(2.0e-16, 1); */
+ /* sim->k = initval(2.0e-16, 1); */
/* Leeman et al 2016: Upstream B till */
- /* sim.k = initval(4.9e-17, 1); */
+ /* 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;
+ 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;
}
#endif
diff --git a/simulation.c b/simulation.c
@@ -20,62 +20,59 @@
#define MAX_ITER_STRESS 20000
void
-prepare_arrays()
+prepare_arrays(struct simulation *sim)
{
- if (sim.nz < 2) {
+ 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.phi_c);
- free(sim.phi_dot);
- free(sim.k);
- free(sim.g_local);
-
- sim.z = linspace(sim.origo_z, /* spatial coordinates */
- sim.origo_z + sim.L_z,
- sim.nz);
- sim.dz = sim.z[1] - sim.z[0]; /* cell spacing */
- sim.mu = zeros(sim.nz); /* stress ratio */
- sim.mu_c = zeros(sim.nz); /* critical-state stress ratio */
- sim.sigma_n_eff = zeros(sim.nz); /* effective normal stress */
- sim.sigma_n = zeros(sim.nz); /* normal stess */
- sim.p_f_ghost = zeros(sim.nz+2); /* fluid pressure with ghost nodes */
- sim.phi = zeros(sim.nz); /* porosity */
- sim.phi_c = zeros(sim.nz); /* critical-state porosity */
- sim.phi_dot = zeros(sim.nz); /* rate of porosity change */
- sim.k = zeros(sim.nz); /* permeability */
- sim.xi = zeros(sim.nz); /* cooperativity length */
- sim.gamma_dot_p = zeros(sim.nz); /* shear velocity */
- sim.v_x = zeros(sim.nz); /* shear velocity */
- sim.g_local = zeros(sim.nz); /* local fluidity */
- sim.g_ghost = zeros(sim.nz+2); /* fluidity with ghost nodes */
- sim.I = zeros(sim.nz); /* inertia number */
- sim.tan_psi = zeros(sim.nz); /* tan(dilatancy_angle) */
+ free(sim->phi);
+ free(sim->k);
+
+ sim->z = linspace(sim->origo_z, /* spatial coordinates */
+ sim->origo_z + sim->L_z,
+ sim->nz);
+ sim->dz = sim->z[1] - sim->z[0]; /* cell spacing */
+ sim->mu = zeros(sim->nz); /* stress ratio */
+ sim->mu_c = zeros(sim->nz); /* critical-state stress ratio */
+ sim->sigma_n_eff = zeros(sim->nz); /* effective normal stress */
+ sim->sigma_n = zeros(sim->nz); /* normal stess */
+ sim->p_f_ghost = zeros(sim->nz+2); /* fluid pressure with ghost nodes */
+ sim->phi = zeros(sim->nz); /* porosity */
+ sim->phi_c = zeros(sim->nz); /* critical-state porosity */
+ sim->phi_dot = zeros(sim->nz); /* rate of porosity change */
+ sim->k = zeros(sim->nz); /* permeability */
+ sim->xi = zeros(sim->nz); /* cooperativity length */
+ sim->gamma_dot_p = zeros(sim->nz); /* shear velocity */
+ sim->v_x = zeros(sim->nz); /* shear velocity */
+ sim->g_local = zeros(sim->nz); /* local fluidity */
+ sim->g_ghost = zeros(sim->nz+2); /* fluidity with ghost nodes */
+ sim->I = zeros(sim->nz); /* inertia number */
+ sim->tan_psi = zeros(sim->nz); /* tan(dilatancy_angle) */
}
void
-free_arrays()
+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.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.g_local);
- free(sim.g_ghost);
- free(sim.I);
- free(sim.tan_psi);
+ 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->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->g_local);
+ free(sim->g_ghost);
+ free(sim->I);
+ free(sim->tan_psi);
}
static void
@@ -110,156 +107,156 @@ check_float(const char name[], const double value, int* return_status)
}
void
-check_simulation_parameters()
+check_simulation_parameters(struct simulation *sim)
{
int return_status;
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->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,
+ 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->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,
+ 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->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->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,
+ 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->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,
+ 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,
+ 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,
+ 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);
+ 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,
+ 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,
+ 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 not a positive number",
- sim.dt, &return_status);
-
- check_float("sim.file_dt", sim.file_dt, &return_status);
- if (sim.file_dt < 0.0)
- warn_parameter_value("sim.file_dt is a negative number",
- sim.file_dt, &return_status);
-
- check_float("sim.phi[0]", sim.phi[0], &return_status);
- if (sim.phi[0] < 0.0 || sim.phi[0] > 1.0)
- warn_parameter_value("sim.phi[0] is not within [0;1]",
- sim.phi[0], &return_status);
-
- 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_angle", sim.dilatancy_angle, &return_status);
- if (sim.dilatancy_angle < 0.0 || sim.dilatancy_angle > 100.0)
- warn_parameter_value("sim.dilatancy_angle is not within [0;100]",
- sim.dilatancy_angle, &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.fluid has an invalid value",
- (double)sim.fluid, &return_status);
-
- if (sim.fluid) {
-
- check_float("sim.p_f_mod_ampl", sim.p_f_mod_ampl, &return_status);
- if (sim.p_f_mod_ampl < 0.0)
- warn_parameter_value("sim.p_f_mod_ampl is not a zero or positive",
- sim.p_f_mod_ampl, &return_status);
-
- 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,
+ 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 not a positive number",
+ sim->dt, &return_status);
+
+ check_float("sim->file_dt", sim->file_dt, &return_status);
+ if (sim->file_dt < 0.0)
+ warn_parameter_value("sim->file_dt is a negative number",
+ sim->file_dt, &return_status);
+
+ check_float("sim->phi[0]", sim->phi[0], &return_status);
+ if (sim->phi[0] < 0.0 || sim->phi[0] > 1.0)
+ warn_parameter_value("sim->phi[0] is not within [0;1]",
+ sim->phi[0], &return_status);
+
+ 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_angle", sim->dilatancy_angle, &return_status);
+ if (sim->dilatancy_angle < 0.0 || sim->dilatancy_angle > 100.0)
+ warn_parameter_value("sim->dilatancy_angle is not within [0;100]",
+ sim->dilatancy_angle, &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);
+
+ 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);
- check_float("sim.p_f_mod_freq", sim.p_f_mod_freq, &return_status);
- if (sim.p_f_mod_freq < 0.0)
- warn_parameter_value("sim.p_f_mod_freq is not a zero or positive",
- sim.p_f_mod_freq, &return_status);
-
- check_float("sim.beta_f", sim.beta_f, &return_status);
- if (sim.beta_f <= 0.0)
- warn_parameter_value("sim.beta_f is not positive",
- sim.beta_f, &return_status);
-
- check_float("sim.mu_f", sim.mu_f, &return_status);
- if (sim.mu_f <= 0.0)
- warn_parameter_value("sim.mu_f is not positive",
- sim.mu_f, &return_status);
-
- check_float("sim.rho_f", sim.rho_f, &return_status);
- if (sim.rho_f <= 0.0)
- warn_parameter_value("sim.rho_f is not positive",
- sim.rho_f, &return_status);
-
- check_float("sim.k[0]", sim.k[0], &return_status);
- if (sim.k[0] <= 0.0)
- warn_parameter_value("sim.k[0] is not positive",
- sim.k[0], &return_status);
+ check_float("sim->p_f_mod_freq", sim->p_f_mod_freq, &return_status);
+ if (sim->p_f_mod_freq < 0.0)
+ warn_parameter_value("sim->p_f_mod_freq is not a zero or positive",
+ sim->p_f_mod_freq, &return_status);
+
+ check_float("sim->beta_f", sim->beta_f, &return_status);
+ if (sim->beta_f <= 0.0)
+ warn_parameter_value("sim->beta_f is not positive",
+ sim->beta_f, &return_status);
+
+ check_float("sim->mu_f", sim->mu_f, &return_status);
+ if (sim->mu_f <= 0.0)
+ warn_parameter_value("sim->mu_f is not positive",
+ sim->mu_f, &return_status);
+
+ check_float("sim->rho_f", sim->rho_f, &return_status);
+ if (sim->rho_f <= 0.0)
+ warn_parameter_value("sim->rho_f is not positive",
+ sim->rho_f, &return_status);
+
+ check_float("sim->k[0]", sim->k[0], &return_status);
+ if (sim->k[0] <= 0.0)
+ warn_parameter_value("sim->k[0] is not positive",
+ sim->k[0], &return_status);
}
if (return_status != 0) {
@@ -269,102 +266,102 @@ check_simulation_parameters()
}
void
-lithostatic_pressure_distribution()
+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]);
+ 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]);
}
/* NEW FUNCTIONS START */
void
-compute_inertia_number()
+compute_inertia_number(struct simulation *sim)
{
int i;
- 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);
+ 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);
}
void
-compute_critical_state_porosity()
+compute_critical_state_porosity(struct simulation *sim)
{
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()
+compute_critical_state_friction(struct simulation *sim)
{
int i;
- 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));
+ 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));
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);
+ 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);
}
void
-compute_friction()
+compute_friction(struct simulation *sim)
{
int i;
- if (sim.transient)
- for (i=0; i<sim.nz; ++i)
- sim.mu[i] = sim.tan_psi[i] + sim.mu_c[i];
+ if (sim->transient)
+ for (i=0; i<sim->nz; ++i)
+ sim->mu[i] = sim->tan_psi[i] + sim->mu_c[i];
else
- for (i=0; i<sim.nz; ++i)
- sim.mu[i] = sim.mu_c[i];
+ for (i=0; i<sim->nz; ++i)
+ sim->mu[i] = sim->mu_c[i];
}
void
-compute_tan_dilatancy_angle()
+compute_tan_dilatancy_angle(struct simulation *sim)
{
int i;
- for (i=0; i<sim.nz; ++i)
- sim.tan_psi[i] = sim.dilatancy_angle*(sim.phi_c[i] - sim.phi[i]);
+ for (i=0; i<sim->nz; ++i)
+ sim->tan_psi[i] = sim->dilatancy_angle*(sim->phi_c[i] - sim->phi[i]);
}
void
-compute_porosity_change()
+compute_porosity_change(struct simulation *sim)
{
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];
- sim.phi[i] += sim.phi_dot[i]*sim.dt;
+ for (i=0; i<sim->nz; ++i) {
+ sim->phi_dot[i] = sim->tan_psi[i]*sim->gamma_dot_p[i]*sim->phi[i];
+ sim->phi[i] += sim->phi_dot[i]*sim->dt;
}
}
void
-compute_permeability()
+compute_permeability(struct simulation *sim)
{
int i;
- 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);
+ 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);
}
/* NEW FUNCTIONS END */
/*void
-compute_friction()
+compute_friction(struct simulation *sim)
{
int i;
- if (sim.fluid)
- for (i=0; i<sim.nz; ++i)
- sim.mu[i] = sim.mu_wall/
- (sim.sigma_n_eff[i]/(sim.P_wall - sim.p_f_top));
+ if (sim->fluid)
+ for (i=0; i<sim->nz; ++i)
+ sim->mu[i] = sim->mu_wall/
+ (sim->sigma_n_eff[i]/(sim->P_wall - sim->p_f_top));
else
- for (i=0; i<sim.nz; ++i)
- sim.mu[i] = sim.mu_wall/
- (1.0 + (1.0 - sim.phi[i])*sim.rho_s*sim.G*
- (sim.L_z - sim.z[i])/sim.P_wall);
+ for (i=0; i<sim->nz; ++i)
+ sim->mu[i] = sim->mu_wall/
+ (1.0 + (1.0 - sim->phi[i])*sim->rho_s*sim->G*
+ (sim->L_z - sim->z[i])/sim->P_wall);
}*/
double
@@ -374,37 +371,37 @@ shear_strain_rate_plastic(const double fluidity, const double friction)
}
void
-compute_shear_strain_rate_plastic()
+compute_shear_strain_rate_plastic(struct simulation *sim)
{
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]);
}
void
-compute_shear_velocity()
+compute_shear_velocity(struct simulation *sim)
{
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;
+ 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()
+compute_effective_stress(struct simulation *sim)
{
int i;
- if (sim.fluid)
- for (i=0; i<sim.nz; ++i)
- sim.sigma_n_eff[i] = sim.sigma_n[i] - sim.p_f_ghost[i+1];
+ if (sim->fluid)
+ for (i=0; i<sim->nz; ++i)
+ sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_ghost[i+1];
else
- for (i=0; i<sim.nz; ++i)
- sim.sigma_n_eff[i] = sim.sigma_n[i];
+ for (i=0; i<sim->nz; ++i)
+ sim->sigma_n_eff[i] = sim->sigma_n[i];
}
static double
@@ -419,16 +416,16 @@ cooperativity_length(const double A,
}
void
-compute_cooperativity_length()
+compute_cooperativity_length(struct simulation *sim)
{
int i;
- for (i=0; i<sim.nz; ++i)
- sim.xi[i] = cooperativity_length(sim.A,
- sim.d,
- sim.mu[i],
- sim.sigma_n_eff[i],
- sim.mu_s,
- sim.C);
+ for (i=0; i<sim->nz; ++i)
+ sim->xi[i] = cooperativity_length(sim->A,
+ sim->d,
+ sim->mu[i],
+ sim->sigma_n_eff[i],
+ sim->mu_s,
+ sim->C);
}
static double
@@ -447,17 +444,17 @@ local_fluidity(const double p,
}
void
-compute_local_fluidity()
+compute_local_fluidity(struct simulation *sim)
{
int i;
- for (i=0; i<sim.nz; ++i)
- sim.g_local[i] = local_fluidity(sim.sigma_n_eff[i],
- 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(sim->sigma_n_eff[i],
+ sim->mu[i],
+ sim->mu_s,
+ sim->C,
+ sim->b,
+ sim->rho_s,
+ sim->d);
}
void
@@ -529,7 +526,8 @@ poisson_solver_1d_cell_update(int i,
}
int
-implicit_1d_jacobian_poisson_solver(const int max_iter,
+implicit_1d_jacobian_poisson_solver(struct simulation *sim,
+ const int max_iter,
const double rel_tol)
{
double r_norm_max;
@@ -537,42 +535,42 @@ implicit_1d_jacobian_poisson_solver(const int max_iter,
int iter, i;
r_norm_max = NAN;
- g_ghost_out = empty(sim.nz+2);
- r_norm = empty(sim.nz);
+ g_ghost_out = empty(sim->nz+2);
+ r_norm = empty(sim->nz);
for (iter=0; iter<max_iter; ++iter) {
#ifdef DEBUG
printf("\n@@@ ITERATION %d @@@\n", iter);
#endif
/* Dirichlet BCs resemble fixed particle velocities */
- set_bc_dirichlet(sim.g_ghost, sim.nz, -1, 0.0);
- set_bc_dirichlet(sim.g_ghost, sim.nz, +1, 0.0);
+ set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0);
+ set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0);
/* Neumann BCs resemble free surfaces */
- /* set_bc_neumann(sim.g_ghost, sim.nz, +1, 0.0, sim.dz); */
+ /* set_bc_neumann(sim->g_ghost, sim->nz, +1, 0.0, sim->dz); */
- for (i=0; i<sim.nz; ++i)
+ for (i=0; i<sim->nz; ++i)
poisson_solver_1d_cell_update(i,
- sim.g_ghost,
- sim.g_local,
+ sim->g_ghost,
+ sim->g_local,
g_ghost_out,
r_norm,
- sim.dz,
- sim.mu,
- sim.sigma_n_eff,
- sim.xi,
- sim.mu_s,
- sim.C,
- sim.b,
- sim.rho_s,
- sim.d);
- r_norm_max = max(r_norm, sim.nz);
-
- copy_values(g_ghost_out, sim.g_ghost, sim.nz+2);
+ sim->dz,
+ sim->mu,
+ sim->sigma_n_eff,
+ sim->xi,
+ sim->mu_s,
+ sim->C,
+ sim->b,
+ sim->rho_s,
+ sim->d);
+ r_norm_max = max(r_norm, sim->nz);
+
+ copy_values(g_ghost_out, sim->g_ghost, sim->nz+2);
if (r_norm_max <= rel_tol) {
- set_bc_dirichlet(sim.g_ghost, sim.nz, -1, 0.0);
- set_bc_dirichlet(sim.g_ghost, sim.nz, +1, 0.0);
+ set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0);
+ set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0);
free(g_ghost_out);
free(r_norm);
/* printf(".. Solution converged after %d iterations\n", iter); */
@@ -589,16 +587,16 @@ implicit_1d_jacobian_poisson_solver(const int max_iter,
}
void
-write_output_file(const int normalize)
+write_output_file(struct simulation *sim, const int normalize)
{
char outfile[200];
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(fp, normalize);
+ print_output(sim, fp, normalize);
fclose(fp);
} else {
fprintf(stderr, "could not open output file: %s", outfile);
@@ -607,43 +605,46 @@ write_output_file(const int normalize)
}
void
-print_output(FILE* fp, const int norm)
+print_output(struct simulation *sim, FILE* fp, const int norm)
{
int i;
double *v_x_out;
if (norm)
- v_x_out = normalize(sim.v_x, sim.nz);
+ v_x_out = normalize(sim->v_x, sim->nz);
else
- v_x_out = copy(sim.v_x, sim.nz);
+ v_x_out = copy(sim->v_x, sim->nz);
- for (i=0; i<sim.nz; ++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\n",
- sim.z[i],
+ "%.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]*sim.sigma_n_eff[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]*sim->sigma_n_eff[i]);
free(v_x_out);
}
int
-coupled_shear_solver(const int max_iter,
+coupled_shear_solver(struct simulation *sim,
+ const int max_iter,
const double rel_tol)
{
int coupled_iter, stress_iter;
- double res_norm, mu_wall_orig;
+ double res_norm, r_norm_max, mu_wall_orig;
+ double *r_norm;
- res_norm = NAN;
- mu_wall_orig = sim.mu_wall;
+ r_norm_max = NAN;
+ r_norm = empty(sim->nz);
+ mu_wall_orig = sim->mu_wall;
stress_iter = 0;
do { /* stress iterations */
@@ -651,86 +652,89 @@ coupled_shear_solver(const int max_iter,
coupled_iter = 0;
do { /* coupled iterations */
- if (sim.transient) {
+ if (sim->transient) {
/* step 1 */
- compute_inertia_number(); /* Eq. 1 */
+ compute_inertia_number(sim); /* Eq. 1 */
/* step 2 */
- compute_critical_state_porosity(); /* Eq. 2 */
+ compute_critical_state_porosity(sim); /* Eq. 2 */
/* step 3 */
- compute_tan_dilatancy_angle(); /* Eq. 5 */
+ compute_tan_dilatancy_angle(sim); /* Eq. 5 */
}
- compute_critical_state_friction(); /* Eq. 7 */
+ compute_critical_state_friction(sim); /* Eq. 7 */
/* step 4 */
- if (sim.transient) {
- compute_porosity_change(); /* Eq. 3 */
- compute_permeability(); /* Eq. 6 */
+ if (sim->transient) {
+ compute_porosity_change(sim); /* Eq. 3 */
+ compute_permeability(sim); /* Eq. 6 */
}
- compute_friction(); /* Eq. 4 */
+ compute_friction(sim); /* Eq. 4 */
/* step 5, Eq. 13 */
- if (sim.fluid)
- if (darcy_solver_1d(MAX_ITER_DARCY, RTOL_DARCY))
+ if (sim->fluid)
+ if (darcy_solver_1d(sim, MAX_ITER_DARCY, RTOL_DARCY))
exit(11);
/* step 6 */
- compute_effective_stress(); /* Eq. 9 */
+ compute_effective_stress(sim); /* Eq. 9 */
/* step 7 */
- compute_local_fluidity(); /* Eq. 10 */
- compute_cooperativity_length(); /* Eq. 12 */
+ compute_local_fluidity(sim); /* Eq. 10 */
+ compute_cooperativity_length(sim); /* Eq. 12 */
/* step 8, Eq. 11 */
- if (implicit_1d_jacobian_poisson_solver(MAX_ITER_GRANULAR,
+ if (implicit_1d_jacobian_poisson_solver(sim,
+ MAX_ITER_GRANULAR,
RTOL_GRANULAR))
exit(12);
/* step 9 */
- compute_shear_strain_rate_plastic(); /* Eq. 8 */
- compute_inertia_number();
- compute_shear_velocity();
+ compute_shear_strain_rate_plastic(sim); /* Eq. 8 */
+ compute_inertia_number(sim);
+ compute_shear_velocity(sim);
#ifdef DEBUG
- for (i=0; i<sim.nz) {
- 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_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) {
+ 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_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
- if (!isnan(sim.v_x_limit) || !isnan(sim.v_x_fix)) {
- if (!isnan(sim.v_x_limit)) {
- res_norm = (sim.v_x_limit - sim.v_x[sim.nz-1])
- /(sim.v_x[sim.nz-1] + 1e-12);
+ if (!isnan(sim->v_x_limit) || !isnan(sim->v_x_fix)) {
+ if (!isnan(sim->v_x_limit)) {
+ res_norm = (sim->v_x_limit - sim->v_x[sim->nz-1])
+ /(sim->v_x[sim->nz-1] + 1e-12);
if (res_norm > 0.0)
res_norm = 0.0;
} else {
- res_norm = (sim.v_x_fix - sim.v_x[sim.nz-1])
- /(sim.v_x[sim.nz-1] + 1e-12);
+ res_norm = (sim->v_x_fix - sim->v_x[sim->nz-1])
+ /(sim->v_x[sim->nz-1] + 1e-12);
}
- sim.mu_wall *= 1.0 + (res_norm*1e-2);
+ sim->mu_wall *= 1.0 + (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, "
- "res_norm=%g, mu_wall=%g\n",
- sim.v_x[sim.nz-1], sim.v_x_fix, sim.v_x_limit,
- res_norm, sim.mu_wall);
+ "res_norm=%g, mu_wall=%g\n",
+ sim->v_x[sim->nz-1], sim->v_x_fix, sim->v_x_limit,
+ res_norm, sim->mu_wall);
return 10;
}
} while (0);
- } while ((!isnan(sim.v_x_fix) || !isnan(sim.v_x_limit))
+ } while ((!isnan(sim->v_x_fix) || !isnan(sim->v_x_limit))
&& fabs(res_norm) > RTOL_STRESS);
- if (!isnan(sim.v_x_limit))
- sim.mu_wall = mu_wall_orig;
+ if (!isnan(sim->v_x_limit))
+ sim->mu_wall = mu_wall_orig;
+
+ free(r_norm);
return 0;
}
diff --git a/simulation.h b/simulation.h
@@ -115,19 +115,19 @@ struct simulation {
double* tan_psi; /* tan(dilatancy_angle) [-] */
};
-void prepare_arrays();
-void free_arrays();
+void prepare_arrays(struct simulation *sim);
+void free_arrays(struct simulation *sim);
-void check_simulation_parameters();
+void check_simulation_parameters(struct simulation *sim);
-void lithostatic_pressure_distribution();
+void lithostatic_pressure_distribution(struct simulation *sim);
-void compute_inertia_number();
-void compute_critical_state_porosity();
-void compute_critical_state_friction();
-void compute_porosity_change();
-void compute_permeability();
-void compute_tan_dilatancy_angle();
+void compute_inertia_number(struct simulation *sim);
+void compute_critical_state_porosity(struct simulation *sim);
+void compute_critical_state_friction(struct simulation *sim);
+void compute_porosity_change(struct simulation *sim);
+void compute_permeability(struct simulation *sim);
+void compute_tan_dilatancy_angle(struct simulation *sim);
void set_bc_neumann(double* g_ghost,
const int nz,
@@ -140,19 +140,22 @@ void set_bc_dirichlet(double* g_ghost,
const int boundary,
const double value);
-void compute_cooperativity_length();
-void compute_local_fluidity();
-void compute_shear_strain_rate_plastic();
-void compute_shear_velocity();
-void compute_effective_stress();
-void compute_friction();
+void compute_cooperativity_length(struct simulation *sim);
+void compute_local_fluidity(struct simulation *sim);
+void compute_shear_strain_rate_plastic(struct simulation *sim);
+void compute_shear_velocity(struct simulation *sim);
+void compute_effective_stress(struct simulation *sim);
+void compute_friction(struct simulation *sim);
-int implicit_1d_jacobian_poisson_solver(const int max_iter,
+int implicit_1d_jacobian_poisson_solver(struct simulation *sim,
+ const int max_iter,
const double rel_tol);
-void write_output_file(const int normalize);
-void print_output(FILE* fp, const int normalize);
+void write_output_file(struct simulation *sim, const int normalize);
+void print_output(struct simulation *sim, FILE* fp, const int normalize);
-int coupled_shear_solver(const int max_iter, const double rel_tol);
+int coupled_shear_solver(struct simulation *sim,
+ const int max_iter,
+ const double rel_tol);
#endif
