tFix Jacobi method in iterative 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 6c48f55a2b5bc7078851062650f5e4e31f90d471
DIR parent 599c08baed18df47fd5a6fac8b587e5d51553519
HTML Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 1 Jul 2019 21:42:16 +0200
Fix Jacobi method in iterative solver
Diffstat:
M fluid.c | 19 +++++++------------
1 file changed, 7 insertions(+), 12 deletions(-)
---
DIR diff --git a/fluid.c b/fluid.c
t@@ -124,9 +124,7 @@ darcy_solver_1d(struct simulation* sim,
* epsilon = 0.0: explicit
* epsilon = 0.5: Crank-Nicolson
* epsilon = 1.0: implicit */
- /* epsilon = 0.5; */
- /* epsilon = 0.0; */
- epsilon = 1.0;
+ epsilon = 0.5;
/* choose relaxation factor, parameter in ]0.0; 1.0]
* theta in ]0.0; 1.0]: underrelaxation
t@@ -192,23 +190,19 @@ darcy_solver_1d(struct simulation* sim,
i, dp_f_dt_expl[i], i, dp_f_dt_impl[i]);
#endif
- p_f_old = sim->p_f_ghost[idx1g(i)];
-
- /* I keep adding dp/dt with every iteration
- * p_f_old is not the value from the previous time step,
- * but the value from the previous iteration. */
- p_f_ghost_new[idx1g(i)] = p_f_old
+ p_f_ghost_new[idx1g(i)] = p_f_ghost_old[idx1g(i)]
+ epsilon*dp_f_dt_impl[i]*sim->dt;
if (epsilon < 1.0)
p_f_ghost_new[idx1g(i)] += (1.0 - epsilon)
*dp_f_dt_expl[i]*sim->dt;
- p_f_ghost_new[idx1g(i)] = p_f_old*(1.0 - theta)
+ p_f_ghost_new[idx1g(i)] = p_f_ghost_old[idx1g(i)]*(1.0 - theta)
+ p_f_ghost_new[idx1g(i)]*theta;
- r_norm[i] = (p_f_ghost_new[idx1g(i)] - p_f_old)
- /(p_f_old + 1e-16);
+ r_norm[i] = (p_f_ghost_new[idx1g(i)]
+ - sim->p_f_ghost[idx1g(i)])
+ /(sim->p_f_ghost[idx1g(i)] + 1e-16);
}
r_norm_max = max(r_norm, sim->nz);
t@@ -252,6 +246,7 @@ darcy_solver_1d(struct simulation* sim,
print_array(sim->p_f_ghost, sim->nz+2);
#endif
}
+ set_fluid_bcs(sim, p_f_top);
if (epsilon < 1.0)
free(dp_f_dt_expl);