diff --git a/external/singularity-eos b/external/singularity-eos index 1ff87551..d443cb54 160000 --- a/external/singularity-eos +++ b/external/singularity-eos @@ -1 +1 @@ -Subproject commit 1ff8755149a6ea6e1af44983f9af2b149391a718 +Subproject commit d443cb54a09234ac2fa0507210d6f978680124f4 diff --git a/src/geometry/mckinney_gammie_ryan.hpp b/src/geometry/mckinney_gammie_ryan.hpp index 10c695e3..a175ee77 100644 --- a/src/geometry/mckinney_gammie_ryan.hpp +++ b/src/geometry/mckinney_gammie_ryan.hpp @@ -54,7 +54,7 @@ class McKinneyGammieRyan { Real smooth, Real hexp_br, Real hexp_nsq, Real hexp_csq) : derefine_poles_(derefine_poles), h_(h), xt_(xt), alpha_(alpha), x0_(x0), smooth_(smooth), norm_(GetNorm_(alpha_, xt_)), hexp_br_(hexp_br), - hexp_nsq_(hexp_nsq), hexp_csq_(hexp_csq_) {} + hexp_nsq_(hexp_nsq), hexp_csq_(hexp_csq) {} McKinneyGammieRyan(bool derefine_poles, Real h, Real xt, Real alpha, Real x0, Real smooth) : derefine_poles_(derefine_poles), h_(h), xt_(xt), alpha_(alpha), x0_(x0), diff --git a/src/microphysics/eos_phoebus/eos_phoebus.cpp b/src/microphysics/eos_phoebus/eos_phoebus.cpp index 8fa9def8..76372b58 100644 --- a/src/microphysics/eos_phoebus/eos_phoebus.cpp +++ b/src/microphysics/eos_phoebus/eos_phoebus.cpp @@ -87,7 +87,6 @@ std::shared_ptr Initialize(ParameterInput *pin) { // Can specify rho_min, etc, in rho_min = pin->GetOrAddReal("eos", "rho_min", 0.0); sie_min = pin->GetOrAddReal("eos", "sie_min", 0.0); - lambda[2] = {0.}; T_min = eos_host.TemperatureFromDensityInternalEnergy(rho_min, sie_min, lambda); rho_max = pin->GetOrAddReal("eos", "rho_max", 1e18); sie_max = pin->GetOrAddReal("eos", "sie_max", 1e35); diff --git a/src/pgen/pgen.cpp b/src/pgen/pgen.cpp index 67a90336..d53bd4a5 100644 --- a/src/pgen/pgen.cpp +++ b/src/pgen/pgen.cpp @@ -14,7 +14,6 @@ #include "geometry/geometry.hpp" #include "pgen.hpp" -#include "phoebus_utils/root_find.hpp" namespace phoebus { @@ -56,32 +55,4 @@ void PostInitializationModifier(ParameterInput *pin, Mesh *pmesh) { f(pin, pmesh); } -class PressResidual { - public: - KOKKOS_INLINE_FUNCTION - PressResidual(const EOS &eos, const Real rho, const Real P, const Real Ye) - : eos_(eos), rho_(rho), P_(P) { - lambda_[0] = Ye; - } - KOKKOS_INLINE_FUNCTION - Real operator()(const Real e) { - return eos_.PressureFromDensityInternalEnergy(rho_, e, lambda_) - P_; - } - - private: - const EOS &eos_; - Real rho_, P_; - Real lambda_[2]; -}; - -KOKKOS_FUNCTION -Real energy_from_rho_P(const EOS &eos, const Real rho, const Real P, const Real emin, - const Real emax, const Real Ye) { - PARTHENON_REQUIRE(P >= 0, "Pressure is negative!"); - PressResidual res(eos, rho, P, Ye); - root_find::RootFind root; - Real eroot = root.regula_falsi(res, emin, emax, 1.e-10 * P, emin - 1.e10); - return rho * eroot; -} - } // namespace phoebus diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 1491f3fd..5ae89d6d 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -22,6 +22,7 @@ using namespace parthenon::package::prelude; // internal includes #include "fluid/fluid.hpp" #include "geometry/geometry.hpp" +#include "phoebus_utils/root_find.hpp" #include "phoebus_utils/variables.hpp" #include "radiation/radiation.hpp" #include "tracers/tracers.hpp" @@ -117,9 +118,34 @@ static std::map +KOKKOS_INLINE_FUNCTION Real energy_from_rho_P(T &eos, const Real rho, const Real P, + const Real emin, const Real emax, + const Real Ye = 0.0) { + PARTHENON_REQUIRE(P >= 0, "Pressure is negative!"); + PressResidual res(eos, rho, P, Ye); + root_find::RootFind root; + Real eroot = root.regula_falsi(res, emin, emax, 1.e-10 * P, emin - 1.e10); + return rho * eroot; +} } // namespace phoebus diff --git a/src/pgen/sedov.cpp b/src/pgen/sedov.cpp index c317e199..a16bd0e7 100644 --- a/src/pgen/sedov.cpp +++ b/src/pgen/sedov.cpp @@ -51,7 +51,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real rhoa = pin->GetOrAddReal("sedov", "rho_ambient", 1.0); const Real rinner = pin->GetOrAddReal("sedov", "rinner", 0.01); - const bool spherical = pin->GetOrAddReal("sedov", "spherical_coords", true); + const bool spherical = pin->GetOrAddBoolean("sedov", "spherical_coords", true); auto &coords = pmb->coords; auto pmesh = pmb->pmy_mesh;