commit 7a0546b0c631a1241defb79e95c29cc3c7c91dce
parent 57ec0c00034da7759ef7a84e9e9ba33efd98b1ad
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 6 Apr 2020 10:53:09 +0200
Fix function definitions and rearrange calls according to analysis
Diffstat:
3 files changed, 61 insertions(+), 48 deletions(-)
diff --git a/1d_fd_simple_shear.c b/1d_fd_simple_shear.c
@@ -248,6 +248,7 @@ main(int argc, char* argv[])
hydrostatic_fluid_pressure_distribution();
set_largest_fluid_timestep(0.5);
}
+ compute_effective_stress(); /* Eq. 9 */
check_simulation_parameters();
@@ -263,29 +264,44 @@ main(int argc, char* argv[])
stressiter = 0;
do {
+
+ if (sim.transient) {
+ /* step 1 */
+ compute_inertia_number(); /* Eq. 1 */
+ /* step 2 */
+ compute_critical_state_porosity(); /* Eq. 2 */
+ /* step 3 */
+ compute_tan_dilatancy_angle(); /* Eq. 5 */
+ compute_critical_state_friction(); /* Eq. 7 */
+ }
+
+ /* step 4 */
+ if (sim.transient) {
+ compute_porosity_change(); /* Eq. 3 */
+ compute_permeability(); /* Eq. 6 */
+ }
+ compute_friction(); /* Eq. 4 */
+
+ /* step 5, Eq. 13 */
if (sim.fluid)
if (darcy_solver_1d(MAX_ITER_DARCY, RTOL_DARCY))
exit(1);
- compute_effective_stress();
- if (sim.transient) {
- compute_inertia_number(); /* new */
- compute_critical_state_porosity(); /* new */
- compute_tan_dilatancy_angle(); /* new */
- compute_critical_state_shear_stress(); /* new */
- compute_shear_stress(); /* new */
- }
- compute_friction();
- compute_cooperativity_length();
+ /* step 6 */
+ compute_effective_stress(); /* Eq. 9 */
+
+ /* step 7 */
+ compute_local_fluidity(); /* Eq. 10 */
+ compute_cooperativity_length(); /* Eq. 12 */
+ /* step 8, Eq. 11 */
if (implicit_1d_jacobian_poisson_solver(MAX_ITER_GRANULAR,
RTOL_GRANULAR))
exit(1);
- compute_shear_strain_rate_plastic();
+ /* step 9 */
+ compute_shear_strain_rate_plastic(); /* Eq. 8 */
compute_shear_velocity();
- if (sim.transient)
- compute_porosity_change(); /* new */
if (!isnan(sim.v_x_limit) || !isnan(sim.v_x_fix)) {
if (!isnan(sim.v_x_limit)) {
diff --git a/simulation.c b/simulation.c
@@ -19,6 +19,7 @@ prepare_arrays()
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,
@@ -36,6 +37,7 @@ prepare_arrays()
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.tau = zeros(sim.nz); /* shear stress */
@@ -58,6 +60,7 @@ free_arrays()
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.tau);
@@ -267,7 +270,7 @@ lithostatic_pressure_distribution()
/* NEW FUNCTIONS START */
void
-compute_inetria_number()
+compute_inertia_number()
{
int i;
for (i=0; i<sim.nz; ++i)
@@ -306,27 +309,11 @@ compute_friction()
}
void
-compute_critical_state_shear_stress()
-{
- int i;
- for (i=0; i<sim.nz; ++i)
- sim.tau_c[i] = sim.mu_c[i]*sim.sigma_n_eff[i];
-}
-
-void
-compute_shear_stress()
-{
- int i;
- for (i=0; i<sim.nz; ++i)
- sim.tau[i] = sim.tan_dilatancy_angle[i]*sim.sigma_n_eff[i] + sim.tau_c[i];
-}
-
-void
compute_tan_dilatancy_angle()
{
int i;
for (i=0; i<sim.nz; ++i)
- sim.tan_dilatancy_angle[i] = sim.K*(sim.phi_c[i] - sim.phi[i]);
+ sim.tan_psi[i] = sim.dilatancy_angle*(sim.phi_c[i] - sim.phi[i]);
}
void
@@ -334,7 +321,7 @@ compute_porosity_change()
{
int i;
for (i=0; i<sim.nz; ++i)
- sim.phi[i] += sim.tan_dilatancy_angle[i]*sim.gamma_dot_p[i]*sim.phi[i]*sim.dt;
+ sim.phi[i] += sim.tan_psi[i]*sim.gamma_dot_p[i]*sim.phi[i]*sim.dt;
}
void
@@ -342,7 +329,8 @@ compute_permeability()
{
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);
+ 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 */
@@ -419,11 +407,11 @@ compute_cooperativity_length()
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);
+ sim.d,
+ sim.mu[i],
+ sim.sigma_n_eff[i],
+ sim.mu_s,
+ sim.C);
}
static double
@@ -446,13 +434,13 @@ compute_local_fluidity()
{
int i;
for (i=0; i<sim.nz; ++i)
- sim.g_ghost[i+1] = local_fluidity(sim.sigma_n_eff[i],
- sim.mu[i],
- sim.mu_s,
- sim.C,
- sim.b,
- sim.rho_s,
- sim.d);
+ 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
@@ -491,6 +479,7 @@ set_bc_dirichlet(double* g_ghost,
static void
poisson_solver_1d_cell_update(int i,
const double* g_in,
+ const double* g_local,
double* g_out,
double* r_norm,
const double dz,
@@ -507,9 +496,7 @@ poisson_solver_1d_cell_update(int i,
coorp_term = dz*dz/(2.0*pow(xi[i], 2.0));
g_out[i+1] = 1.0/(1.0 + coorp_term)*(coorp_term*
- local_fluidity(p[i], mu[i], mu_s, C, b, rho_s, d)
- + g_in[i+2]/2.0
- + g_in[i]/2.0);
+ g_local[i] + g_in[i+2]/2.0 + g_in[i]/2.0);
r_norm[i] = pow(g_out[i+1] - g_in[i+1], 2.0) / (pow(g_out[i+1], 2.0) + 1e-16);
@@ -549,6 +536,7 @@ implicit_1d_jacobian_poisson_solver(const int max_iter,
for (i=0; i<sim.nz; ++i)
poisson_solver_1d_cell_update(i,
sim.g_ghost,
+ sim.g_local,
g_ghost_out,
r_norm,
sim.dz,
diff --git a/simulation.h b/simulation.h
@@ -109,6 +109,7 @@ struct simulation {
double* xi; /* cooperativity length */
double* gamma_dot_p; /* plastic shear strain rate [s^-1] */
double* v_x; /* shear velocity [m/s] */
+ double* g_local; /* local fluidity */
double* g_ghost; /* fluidity with ghost nodes */
double* I; /* inertia number [-] */
double* tau; /* shear stress [Pa] */
@@ -122,6 +123,13 @@ void check_simulation_parameters();
void lithostatic_pressure_distribution();
+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 set_bc_neumann(double* g_ghost,
const int nz,
const int boundary,
@@ -134,6 +142,7 @@ void set_bc_dirichlet(double* g_ghost,
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();
