commit c95b1939d5e265279eea39af7c47dbb863c21c22 parent 169a3e9cb29a2010a6d25dba71998e869bf7a446 Author: Anders Damsgaard <anders@adamsgaard.dk> Date: Wed, 14 Jan 2026 09:52:14 +0100 fix(fluid): add safety floor to permeability harmonic mean Prevent division by zero when k_zp + k_ or k_zn + k_ approaches zero by adding 1e-30 floor value to denominator. Diffstat:
| M | fluid.c | | | 8 | ++++---- |
1 file changed, 4 insertions(+), 4 deletions(-)
diff --git a/fluid.c b/fluid.c @@ -143,10 +143,10 @@ darcy_pressure_change_1d(const int i, else k_zp = k[i + 1]; - div_k_grad_p = (2.0 * k_zp * k_ / (k_zp + k_) - * (p_f_ghost_in[i + 2] - p_f_ghost_in[i + 1]) / dz - - 2.0 * k_zn * k_ / (k_zn + k_) - * (p_f_ghost_in[i + 1] - p_f_ghost_in[i]) / dz + double k_harm_p = 2.0 * k_zp * k_ / fmax(k_zp + k_, 1e-30); + double k_harm_n = 2.0 * k_zn * k_ / fmax(k_zn + k_, 1e-30); + div_k_grad_p = (k_harm_p * (p_f_ghost_in[i + 2] - p_f_ghost_in[i + 1]) / dz + - k_harm_n * (p_f_ghost_in[i + 1] - p_f_ghost_in[i]) / dz ) / dz; #ifdef DEBUG