tperform temporal integration of fluid pressure outside of fluid solver - 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 97e88f765d0b872280c99f2cd5753ed1afefe420
DIR parent 3dfa0a14ecc48bdd9021dbe83072c161af72cbf8
HTML Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Wed, 18 Nov 2020 11:31:56 +0100
perform temporal integration of fluid pressure outside of fluid solver
Diffstat:
M fluid.c | 82 +++++++++++++------------------
M simulation.c | 16 ++++++++--------
2 files changed, 41 insertions(+), 57 deletions(-)
---
DIR diff --git a/fluid.c b/fluid.c
t@@ -98,11 +98,11 @@ square_pulse(const double time,
}
static void
-set_fluid_bcs(struct simulation *sim, const double p_f_top)
+set_fluid_bcs(double *p_f_ghost, 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,
+ set_bc_dirichlet(p_f_ghost, sim->nz, +1, p_f_top);
+ p_f_ghost[sim->nz] = p_f_top; /* Include top node in BC */
+ set_bc_neumann(p_f_ghost,
sim->nz,
-1,
sim->phi[0] * sim->rho_f * sim->G,
t@@ -170,9 +170,13 @@ darcy_solver_1d(struct simulation *sim,
const double rel_tol)
{
int i, iter, solved = 0;
- double epsilon, theta, p_f_top, r_norm_max = NAN;
- double *dp_f_dt_expl, *p_f_ghost_old, *dp_f_dt_impl, *p_f_ghost_new;
- double *tmp, *r_norm;
+ double epsilon, p_f_top, r_norm_max = NAN;
+ double *p_f_dot_expl, *p_f_dot_impl, *p_f_ghost_new;
+ double *tmp, *r_norm, *old_val;
+
+ p_f_dot_impl = zeros(sim->nz);
+ p_f_dot_expl = zeros(sim->nz);
+ p_f_ghost_new = zeros(sim->nz + 2);
/* choose integration method, parameter in [0.0; 1.0]
* epsilon = 0.0: explicit
t@@ -180,12 +184,6 @@ darcy_solver_1d(struct simulation *sim,
* epsilon = 1.0: * implicit */
epsilon = 0.5;
- /* choose relaxation factor, parameter in ]0.0; 1.0] theta in
- * ]0.0; 1.0]: underrelaxation theta = 1.0: Gauss-Seidel theta >
- * 1.0: overrelaxation */
- /* 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,
t@@ -205,14 +203,15 @@ darcy_solver_1d(struct simulation *sim,
sim->p_f_mod_freq,
sim->p_f_mod_pulse_time);
+
/* set fluid BCs (1 of 2) */
- set_fluid_bcs(sim, p_f_top);
+ copy_values(sim->p_f_ghost, p_f_ghost_new, sim->nz + 2);
+ set_fluid_bcs(p_f_ghost_new, sim, p_f_top);
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[i] = darcy_pressure_change_1d(i,
+ p_f_dot_expl[i] = darcy_pressure_change_1d(i,
sim->nz,
sim->p_f_ghost,
sim->phi,
t@@ -226,25 +225,23 @@ darcy_solver_1d(struct simulation *sim,
}
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);
+ old_val = empty(sim->nz);
for (iter = 0; iter < max_iter; ++iter) {
+ copy_values(sim->p_f_dot, old_val, sim->nz);
/* set fluid BCs (2 of 2) */
- set_fluid_bcs(sim, p_f_top);
+ set_fluid_bcs(p_f_ghost_new, sim, p_f_top);
#ifdef DEBUG
puts(".. p_f_ghost after BC:");
- print_array(sim->p_f_ghost, sim->nz + 2);
+ print_array(p_f_ghost_new, sim->nz + 2);
#endif
for (i = 0; i < sim->nz - 1; ++i)
- dp_f_dt_impl[i] = darcy_pressure_change_1d(i,
+ p_f_dot_impl[i] = darcy_pressure_change_1d(i,
sim->nz,
- sim->p_f_ghost,
+ p_f_ghost_new,
sim->phi,
sim->phi_dot,
sim->k,
t@@ -255,22 +252,13 @@ darcy_solver_1d(struct simulation *sim,
sim->D);
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;
+ sim->p_f_dot[i] = epsilon * p_f_dot_impl[i];
if (epsilon < 1.0)
- p_f_ghost_new[i + 1] += (1.0 - epsilon)
- * 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;
+ sim->p_f_dot[i + 1] += (1.0 - epsilon) * p_f_dot_expl[i];
- r_norm[i] = fabs(residual(p_f_ghost_new[i + 1],
- sim->p_f_ghost[i + 1]));
+ p_f_ghost_new[i + 1] += sim->p_f_dot[i] * sim->dt;
+ r_norm[i] = fabs(residual(sim->p_f_dot[i], old_val[i]));
}
r_norm_max = max(r_norm, sim->nz);
t@@ -284,7 +272,7 @@ darcy_solver_1d(struct simulation *sim,
sim->p_f_ghost = tmp;
#ifdef DEBUG
puts(".. p_f_ghost after update:");
- print_array(sim->p_f_ghost, sim->nz + 2);
+ print_array(p_f_ghost_new, sim->nz + 2);
#endif
if (r_norm_max <= rel_tol) {
t@@ -295,10 +283,10 @@ darcy_solver_1d(struct simulation *sim,
break;
}
}
- free(dp_f_dt_impl);
+ free(p_f_dot_impl);
+ free(p_f_dot_expl);
free(p_f_ghost_new);
free(r_norm);
- free(p_f_ghost_old);
if (!solved) {
fprintf(stderr, "darcy_solver_1d: ");
fprintf(stderr, "Solution did not converge after %d iterations\n",
t@@ -307,19 +295,15 @@ darcy_solver_1d(struct simulation *sim,
}
} else {
for (i = 0; i < sim->nz; ++i)
- sim->p_f_ghost[i + 1] += dp_f_dt_expl[i] * sim->dt;
+ sim->p_f_dot[i] = p_f_dot_expl[i];
solved = 1;
#ifdef DEBUG
- puts(".. dp_f_dt_expl:");
- print_array(dp_f_dt_expl, sim->nz);
- puts(".. p_f_ghost after explicit solution:");
- print_array(sim->p_f_ghost, sim->nz + 2);
+ puts(".. p_f_dot_expl:");
+ print_array(p_f_dot_expl, sim->nz);
#endif
+ free(p_f_dot_expl);
}
- set_fluid_bcs(sim, p_f_top);
-
- if (epsilon < 1.0)
- free(dp_f_dt_expl);
+ set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
return solved - 1;
}
DIR diff --git a/simulation.c b/simulation.c
t@@ -744,9 +744,9 @@ temporal_increment(struct simulation *sim)
for (i = 0; i < sim->nz; ++i)
sim->phi[i] += sim->phi_dot[i] * sim->dt;
- /*if (sim->fluid)
+ if (sim->fluid)
for (i = 1; i <= sim->nz; ++i)
- sim->p_f_ghost[i] += sim->p_f_dot[i] * sim->dt;*/
+ sim->p_f_ghost[i] += sim->p_f_dot[i] * sim->dt;
sim->t += sim->dt;
}
t@@ -757,12 +757,12 @@ coupled_shear_solver(struct simulation *sim,
const double rel_tol)
{
int i, coupled_iter, stress_iter = 0;
- double *r_norm, *oldval;
+ double *r_norm, *old_val;
double r_norm_max, vel_res_norm = NAN, mu_wall_orig = sim->mu_wall;
if (sim->transient) {
r_norm = empty(sim->nz);
- oldval = empty(sim->nz);
+ old_val = empty(sim->nz);
}
do { /* stress iterations */
t@@ -770,7 +770,7 @@ coupled_shear_solver(struct simulation *sim,
do { /* coupled iterations */
if (sim->transient) {
- copy_values(sim->phi_dot, oldval, sim->nz);
+ copy_values(sim->phi_dot, old_val, sim->nz);
/* step 1 */
compute_inertia_number(sim); /* Eq. 1 */
t@@ -828,14 +828,14 @@ coupled_shear_solver(struct simulation *sim,
/* stable porosity change field == coupled solution converged */
if (sim->transient) {
for (i = 0; i < sim->nz; ++i)
- r_norm[i] = fabs(residual(sim->phi_dot[i], oldval[i]));
+ r_norm[i] = fabs(residual(sim->phi_dot[i], old_val[i]));
r_norm_max = max(r_norm, sim->nz);
if (r_norm_max <= rel_tol)
break;
if (coupled_iter++ >= max_iter) {
free(r_norm);
- free(oldval);
+ free(old_val);
fprintf(stderr, "coupled_shear_solver: ");
fprintf(stderr, "Transient solution did not converge "
"after %d iterations\n", coupled_iter);
t@@ -874,7 +874,7 @@ coupled_shear_solver(struct simulation *sim,
if (sim->transient) {
free(r_norm);
- free(oldval);
+ free(old_val);
}
temporal_increment(sim);