tFix solver by disabling under/overrelaxation - cngf-pf - continuum model for granular flows with pore-pressure dynamics (renamed from 1d_fd_simple_shear)
(HTM) git clone git://src.adamsgaard.dk/cngf-pf
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
(DIR) commit acfcb445ba8d101ce61446f4500f6269d7e74af0
(DIR) parent 6c1091d29c2c4b5519cd270a35c3dc427696bd50
(HTM) Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Fri, 28 Jun 2019 16:03:20 +0200
Fix solver by disabling under/overrelaxation
Diffstat:
M fluid.c | 22 +++++++++++-----------
1 file changed, 11 insertions(+), 11 deletions(-)
---
(DIR) diff --git a/fluid.c b/fluid.c
t@@ -47,7 +47,7 @@ set_largest_fluid_timestep(struct simulation* sim, const double safety)
sim->dt = safety*0.5*dx_min*dx_min/diff_max;
if (sim->file_dt*0.5 < sim->dt)
- sim->dt = sim->file_dt*0.5;
+ sim->dt = sim->file_dt;
}
static double
t@@ -87,16 +87,16 @@ darcy_pressure_change_1d(const int i,
k_zn, k_, k_zp);
#endif
- div_k_grad_p = (2.0*k_zp*k_/(k_zp + k_) * (p_zp - p)/dz -
- 2.0*k_zn*k_/(k_zn + k_) * (p - p_zn)/dz
- )/dz;
+ div_k_grad_p = (2.0*k_zp*k_/(k_zp + k_) * (p_zp - p)/dz -
+ 2.0*k_zn*k_/(k_zn + k_) * (p - p_zn)/dz
+ )/dz;
+
#ifdef DEBUG
printf("phi[%d]=%g\tdiv_k_grad_p[%d]=%g\n",
i, phi[i], i, div_k_grad_p);
#endif
- /* return delta p */
- return dt/(beta_f*phi[i]*mu_f)*div_k_grad_p;
+ return dt/(beta_f*phi[i]*mu_f)*div_k_grad_p;
}
int
t@@ -112,14 +112,14 @@ darcy_solver_1d(struct simulation* sim,
* epsilon = 0.0: explicit
* epsilon = 0.5: Crank-Nicolson
* epsilon = 1.0: implicit */
- epsilon = 0.5;
+ 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 */
- theta = 0.05;
- /* theta = 1.7; */
+ /* TODO: values other than 1.0 do not work! */
+ theta = 1.0;
p_f_top = sine_wave(sim->t,
sim->p_f_mod_ampl,
t@@ -200,7 +200,7 @@ darcy_solver_1d(struct simulation* sim,
print_array(sim->p_f_ghost, sim->nz+2);
#endif
- if (r_norm_max <= rel_tol) {
+ if (r_norm_max <= rel_tol) {
set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top);
sim->p_f_ghost[idx1g(sim->nz-1)] = p_f_top; /* top node in BC */
set_bc_neumann(sim->p_f_ghost, sim->nz, -1);
t@@ -212,7 +212,7 @@ darcy_solver_1d(struct simulation* sim,
printf(".. Solution converged after %d iterations\n", iter);
#endif
return 0;
- }
+ }
}
free(dp_f_expl);