tCreated implicit scheme for Darcy Solver, waiting for review, validation, and optimization from Anders. - cngf-pf - continuum model for granular flows with pore-pressure dynamics (renamed from 1d_fd_simple_shear)
HTML git clone git://src.adamsgaard.dk/cngf-pf
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
DIR commit e416cbbbc45454688a91db2ec0ebedd4ed7afd71
DIR parent c0ecbde7b32404ef869fa3ca9959cbee18777309
HTML Author: Ian Madden <iamadden@stanford.edu>
Date: Mon, 10 May 2021 15:28:03 -0700
Created implicit scheme for Darcy Solver, waiting for review, validation, and optimization from Anders.
Diffstat:
M fluid.c | 104 ++++++++++++++++++++++++++++---
M simulation.c | 24 ++++++++++++++++--------
M simulation.h | 1 +
3 files changed, 111 insertions(+), 18 deletions(-)
---
DIR diff --git a/fluid.c b/fluid.c
t@@ -164,6 +164,68 @@ darcy_pressure_change_1d(const int i,
}
}
+static double
+darcy_pressure_change_1d_impl(const int i,
+ const int nz,
+ const double dt,
+ const double *p_f_old_val,
+ const double *p_f_ghost_in,
+ double *p_f_ghost_out,
+ const double *phi,
+ const double *phi_dot,
+ const double *k,
+ const double dz,
+ const double beta_f,
+ const double alpha,
+ const double mu_f,
+ const double D)
+{
+ double k_, div_k_grad_p, k_zn, k_zp,rhs_term;
+ double omega = 1.0;
+
+ if (D > 0.0)
+ return D * (p_f_ghost_in[i + 2]
+ - 2.0 * p_f_ghost_in[i + 1]
+ + p_f_ghost_in[i]) / (dz * dz);
+ else {
+ k_ = k[i];
+ if (i == 0)
+ k_zn = k_;
+ else
+ k_zn = k[i - 1];
+ if (i == nz - 1)
+ k_zp = k_;
+ else
+ k_zp = k[i + 1];
+ rhs_term = dt / ((alpha + beta_f * phi[i]) * mu_f) * ((2.0 * k_zp * k_ / (k_zp + k_) /(dz*dz)) + (2.0 * k_zn * k_ / (k_zn + k_) /(dz*dz)));
+ p_f_ghost_out[i+1] = (1 / (1.0 + rhs_term)) * (p_f_old_val[i+1] + dt* (1.0 / ((alpha + beta_f * phi[i]) * mu_f) *
+ (2.0 * k_zp * k_ / (k_zp + k_)
+ * (p_f_ghost_in[i + 2]) / dz
+ - 2.0 * k_zn * k_ / (k_zn + k_)
+ * ( - p_f_ghost_in[i]) / dz)/dz
+ - 1.0 / ((alpha + beta_f * phi[i]) * (1.0 - phi[i])) * phi_dot[i]));
+ p_f_ghost_out[i+1] = omega * p_f_ghost_out[i+1] + (1.0-omega)*p_f_ghost_in[i+1];
+
+ div_k_grad_p = (2.0 * k_zp * k_ / (k_zp + k_)
+ * (p_f_ghost_out[i + 2] - p_f_ghost_out[i + 1]) / dz
+ - 2.0 * k_zn * k_ / (k_zn + k_)
+ * (p_f_ghost_out[i + 1] - p_f_ghost_out[i]) / dz
+) / dz;
+#ifdef DEBUG
+ printf("%s [%d]: phi=%g\tdiv_k_grad_p=%g\tphi_dot=%g\n",
+ __func__, i, phi[i], div_k_grad_p, phi_dot[i]);
+
+ printf(" p[%d, %d, %d] = [%g, %g, %g]\tk=[%g, %g, %g]\n",
+ i, i + 1, i + 2,
+ p_f_ghost_in[i], p_f_ghost_in[i + 1], p_f_ghost_in[i + 2],
+ k_zn, k_, k_zp);
+#endif
+ /* use the values from the next time step as the time derivative for this iteration */
+ return 1.0 / ((alpha + beta_f * phi[i]) * mu_f) * div_k_grad_p
+ - 1.0 / ((alpha + beta_f * phi[i]) * (1.0 - phi[i])) * phi_dot[i];
+ }
+}
+
int
darcy_solver_1d(struct simulation *sim,
const int max_iter,
t@@ -208,6 +270,7 @@ darcy_solver_1d(struct simulation *sim,
/* set fluid BCs (1 of 2) */
set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
+ set_fluid_bcs(sim->p_f_next, sim, p_f_top);
/* explicit solution to pressure change */
if (epsilon < 1.0) {
t@@ -231,11 +294,21 @@ darcy_solver_1d(struct simulation *sim,
/* implicit solution with Jacobian iterations */
if (epsilon > 0.0) {
-
+ /* grabbing the n + 1 iteration values for k and phi */
+ double * k_n;
+ double * phi_n;
+ k_n = zeros(sim->nz);
+ phi_n = zeros(sim->nz);
+ for (i = 0; i < sim->nz; ++i) {
+ phi_n[i] = sim->phi[i] + sim->dt*sim->phi_dot[i];
+ k_n[i] = sim->d * sim->d / 180.0
+ * phi_n[i] * phi_n[i] * phi_n[i]
+ / pow(1.0 - phi_n[i], 2.0);
+ }
#ifdef DEBUG
printf("\nEXPLICIT SOLVER IN %s\n", __func__);
#endif
- copy_values(sim->p_f_ghost, sim->tmp_ghost, sim->nz + 2);
+ copy_values(sim->p_f_next, sim->tmp_ghost, sim->nz + 2);
for (iter = 0; iter < max_iter; ++iter) {
copy_values(sim->p_f_dot_impl, sim->fluid_old_val, sim->nz);
t@@ -254,12 +327,15 @@ darcy_solver_1d(struct simulation *sim,
#endif
for (i = 0; i < sim->nz - 1; ++i)
- sim->p_f_dot_impl[i] = darcy_pressure_change_1d(i,
+ sim->p_f_dot_impl[i] = darcy_pressure_change_1d_impl(i,
sim->nz,
+ sim->dt,
+ sim->p_f_ghost,
sim->tmp_ghost,
- sim->phi,
+ sim->p_f_next,
+ phi_n,
sim->phi_dot,
- sim->k,
+ k_n,
sim->dz,
sim->beta_f,
sim->alpha,
t@@ -271,10 +347,12 @@ darcy_solver_1d(struct simulation *sim,
errx(1, "NaN at sim->p_f_dot_impl[%d] (t = %g s, iter = %d)",
i, sim->t, iter);
- for (i = 0; i < sim->nz - 1; ++i)
- sim->p_f_dot_impl_r_norm[i] = fabs(residual(sim->p_f_dot_impl[i],
- sim->fluid_old_val[i]));
+ set_fluid_bcs(sim->p_f_next, sim, p_f_top);
+ for (i = 0; i < sim->nz-1; ++i)
+ sim->p_f_dot_impl_r_norm[i] = fabs(residual(sim->p_f_next[i],
+ sim->tmp_ghost[i]));
r_norm_max = max(sim->p_f_dot_impl_r_norm, sim->nz - 1);
+ copy_values(sim->p_f_next, sim->tmp_ghost, sim->nz + 2);
#ifdef DEBUG
puts(".. p_f_ghost_new:");
t@@ -289,6 +367,8 @@ darcy_solver_1d(struct simulation *sim,
break;
}
}
+ free(k_n);
+ free(phi_n);
if (!solved) {
fprintf(stderr, "darcy_solver_1d: ");
fprintf(stderr, "Solution did not converge after %d iterations\n",
t@@ -303,11 +383,15 @@ darcy_solver_1d(struct simulation *sim,
sim->p_f_dot[i] = epsilon * sim->p_f_dot_impl[i]
+ (1.0 - epsilon) * sim->p_f_dot_expl[i];
- set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
-
for (i = 0; i < sim->nz; ++i)
sim->p_f_dot[i] = omega * sim->p_f_dot[i]
+ (1.0 - omega) * sim->p_f_dot_old[i];
+
+ for (i = 0; i < sim->nz-1; ++i)
+ sim->p_f_next[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, sim, p_f_top);
#ifdef DEBUG
printf(".. epsilon = %g\n", epsilon);
puts(".. p_f_dot_expl:");
DIR diff --git a/simulation.c b/simulation.c
t@@ -9,7 +9,7 @@
/* iteration limits for solvers */
#define MAX_ITER_GRANULAR 100000
-#define MAX_ITER_DARCY 10000
+#define MAX_ITER_DARCY 1000000
#define MAX_ITER_STRESS 20000
/* tolerance criteria when in velocity driven or velocity limited mode */
t@@ -47,9 +47,10 @@ init_sim(struct simulation *sim)
/* Henann and Kamrin 2016 */
/* sim->mu_s = 0.3819; */
/* sim->C = 0.0; */
-
+ /* Testing */
+ sim->mu_s = 0.30;
/* Damsgaard et al 2013 */
- sim->mu_s = tan(DEG2RAD(22.0));
+ /*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 */
t@@ -149,6 +150,7 @@ prepare_arrays(struct simulation *sim)
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 = 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);
t@@ -181,6 +183,7 @@ free_arrays(struct simulation *sim)
free(sim->sigma_n_eff);
free(sim->sigma_n);
free(sim->p_f_ghost);
+ free(sim->p_f_next);
free(sim->p_f_dot);
free(sim->p_f_dot_expl);
free(sim->p_f_dot_impl);
t@@ -541,7 +544,9 @@ compute_effective_stress(struct simulation *sim)
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];
+ /* average of current and next step pressure values for effective stress - may not be optimal */
+ sim->sigma_n_eff[i] = sim->sigma_n[i] - ((sim->p_f_ghost[i + 1] + sim->p_f_next[i+1])/2.0);
+ //sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_ghost[i + 1];
if (sim->sigma_n_eff[i] < 0)
errx(1, "%s: sigma_n_eff[%d] is negative with value %g\n",
__func__, i, sim->sigma_n_eff[i]);
t@@ -803,6 +808,9 @@ coupled_shear_solver(struct simulation *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,sim->nz+2);
+ compute_effective_stress(sim); /* Eq. 9 */
+
do { /* stress iterations */
coupled_iter = 0;
do { /* coupled iterations */
t@@ -827,8 +835,8 @@ coupled_shear_solver(struct simulation *sim,
compute_friction(sim); /* Eq. 4 */
/* step 5, Eq. 13 */
- if (sim->fluid)
- if (darcy_solver_1d(sim, MAX_ITER_DARCY, rel_tol))
+ if (sim->fluid && (sim->t > 0) )
+ if (darcy_solver_1d(sim, MAX_ITER_DARCY, 1e-3 * rel_tol))
exit(11);
/* step 6 */
t@@ -869,7 +877,7 @@ coupled_shear_solver(struct simulation *sim,
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(sim->g_r_norm, sim->nz);
- if (r_norm_max <= rel_tol)
+ if (r_norm_max <= rel_tol && coupled_iter > 0)
break;
if (coupled_iter++ >= max_iter) {
t@@ -893,7 +901,7 @@ coupled_shear_solver(struct simulation *sim,
vel_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 + (vel_res_norm * 1e-2);
+ sim->mu_wall *= 1.0 + (vel_res_norm * 1e-3);
}
if (++stress_iter > MAX_ITER_STRESS) {
fprintf(stderr, "error: stress solution did not converge:\n");
DIR diff --git a/simulation.h b/simulation.h
t@@ -106,6 +106,7 @@ struct simulation {
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; /* 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] */