tCreated implicit scheme for Darcy Solver, waiting for review, validation, and … | |
git clone git://src.adamsgaard.dk/cngf-pf | |
Log | |
Files | |
Refs | |
README | |
LICENSE | |
--- | |
commit e416cbbbc45454688a91db2ec0ebedd4ed7afd71 | |
parent c0ecbde7b32404ef869fa3ca9959cbee18777309 | |
Author: Ian Madden <[email protected]> | |
Date: Mon, 10 May 2021 15:28:03 -0700 | |
Created implicit scheme for Darcy Solver, waiting for review, validation, and o… | |
Diffstat: | |
M fluid.c | 104 ++++++++++++++++++++++++++++-… | |
M simulation.c | 24 ++++++++++++++++-------- | |
M simulation.h | 1 + | |
3 files changed, 111 insertions(+), 18 deletions(-) | |
--- | |
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_ / (… | |
+ p_f_ghost_out[i+1] = (1 / (1.0 + rhs_term)) * (p_f_old_val[i+1] + dt* (1.… | |
+ (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])) * … | |
+ p_f_ghost_out[i+1] = omega * p_f_ghost_out[i+1] + (1.0-omega)*… | |
+ | |
+ div_k_grad_p = (2.0 * k_zp * k_ / (k_zp + k_) | |
+ * (p_f_ghost_out[i + 2] - p_f_ghost_out[i + 1]… | |
+ - 2.0 * k_zn * k_ / (k_zn + k_) | |
+ * (p_f_ghost_out[i + 1] - p_f_ghost_out[i]) / … | |
+) / 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 + … | |
+ k_zn, k_, k_zp); | |
+#endif | |
+ /* use the values from the next time step as the time derivative for this … | |
+ 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])) * … | |
+ } | |
+} | |
+ | |
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… | |
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_1… | |
+ sim->p_f_dot_impl[i] = darcy_pressure_change_1… | |
… | |
+ … | |
+ sim-… | |
… | |
- … | |
+ … | |
+ … | |
… | |
- … | |
+ … | |
… | |
… | |
… | |
t@@ -271,10 +347,12 @@ darcy_solver_1d(struct simulation *sim, | |
errx(1, "NaN at sim->p_f_dot_impl[%d] … | |
i, sim->t, iter); | |
- for (i = 0; i < sim->nz - 1; ++i) | |
- sim->p_f_dot_impl_r_norm[i] = fabs(residual(si… | |
- si… | |
+ 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(si… | |
+ si… | |
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 it… | |
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… | |
+ | |
+ 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:"); | |
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… | |
+ /* average of current and next step pressure values for effe… | |
+ sim->sigma_n_eff[i] = sim->sigma_n[i] - ((sim->p_f_gho… | |
+ //sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_gho… | |
if (sim->sigma_n_eff[i] < 0) | |
errx(1, "%s: sigma_n_eff[%d] is negative with … | |
__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_t… | |
+ if (sim->fluid && (sim->t > 0) ) | |
+ if (darcy_solver_1d(sim, MAX_ITER_DARCY, 1e-3 … | |
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->… | |
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->n… | |
/ (sim->v_x[sim->nz - 1] + 1e-1… | |
} | |
- 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 conver… | |
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… | |
double *p_f_dot_impl; /* fluid pressure change (implicit solution) [Pa… |