From 7ed2f9a394a82537ec8959e071f70e9ef6a7f8e0 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Thu, 1 Feb 2024 01:44:31 +0100 Subject: [PATCH] src/equilibrium.f90: fix floating point exception MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit If frhopol is called with a very small ρ_p, less than then step size of ε_machine^⅓, the two sided finite step differentiation evaluates frhotor at ρ_p < 0, producing a floating point exception. In this case we use the single- sided definition. --- src/equilibrium.f90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index ad8cb61..43dedbb 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -1155,7 +1155,12 @@ contains f(1) = frhotor(x(1)) - rho_t else ! return f'(x), computed numerically - df(1,1) = (frhotor(x(1) + e) - frhotor(x(1) - e)) / (2*e) + if (x(1) - e > 0) then + df(1,1) = (frhotor(x(1) + e) - frhotor(x(1) - e)) / (2*e) + else + ! single-sided when close to ρ=0 + df(1,1) = (frhotor(x(1) + e) - frhotor(x(1))) / e + end if end if end subroutine