From be2f5b9c7c01cbf162b1e2ea5d5bf369a03bb687 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 12:08:54 -0700 Subject: [PATCH] ok... lets try this --- singularity-eos/eos/eos_base.hpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 74fc286085..600cace9cf 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -856,10 +856,10 @@ class EosBase { using RootFinding1D::findRoot; // more robust but slower. Better default. using RootFinding1D::Status; - // Pressure is not monotone in density at low densities, due to a - // phase transition, which can prevent convergence. We want to - // approach tension from above, not below. Choose close to, but - // above, normal density for a metal like copper. + // Pressure is not monotone in density at low densities, which can + // prevent convergence. We want to approach tension from above, + // not below. Choose close to, but above, normal density for a + // metal like copper. constexpr Real DEFAULT_RHO_GUESS = 12; CRTP copy = *(static_cast(this)); @@ -870,14 +870,13 @@ class EosBase { return copy.PressureFromDensityTemperature(r, temp, lambda); }; Real rhoguess = rho; // use input density - if ((rhoguess < rhomin) || (rhoguess > rhomax)) { - if ((rhomin < DEFAULT_RHO_GUESS) && (DEFAULT_RHO_GUESS < rhomax)) { + if ((rhoguess <= rhomin) || (rhoguess >= rhomax)) { + if ((rhomin <= DEFAULT_RHO_GUESS) && (DEFAULT_RHO_GUESS <= rhomax)) { rhoguess = DEFAULT_RHO_GUESS; } else { rhoguess = 0.5 * (rhomin + rhomax); } } - printf("%.14e %.14e %.14e\n", rhomin, rhomax, rhoguess); auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, robust::EPS(), robust::EPS(), rho); // JMM: This needs to not fail and instead return something sane