From db1507388791333f303275eda5029c80c62eb988 Mon Sep 17 00:00:00 2001 From: Logan Alexandria White Date: Wed, 11 Mar 2026 09:43:59 -0700 Subject: [PATCH 01/12] todos plus initial context --- src/fluid/con2prim_robust.hpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/fluid/con2prim_robust.hpp b/src/fluid/con2prim_robust.hpp index 2bbe8de1..f2759cad 100644 --- a/src/fluid/con2prim_robust.hpp +++ b/src/fluid/con2prim_robust.hpp @@ -371,6 +371,17 @@ class ConToPrim { // Doesn't seem to be at a quick glance. // const Real mu_r = res.compute_upper_bound(h0sq_); // solve + + /** + * TODO: implement method to find lower enthalpy bound h0 (upper root find bound) + * + * Katsaun et al. 2021 only assumes h > 0, not h ≥ 1; we shouldn't have a hardcoded upper bound here. + * stellar collapse eos can have h < 1 due to included nuclear binding energies (e.g. Steiner et al. + * 2012), this breaks our root find when mapping onto the gr grid at initialization. + * + * can probably be resolved with an eos-conscious method using primitives (P, rho, eps) passed into here, solution tbd. + */ + root_find::RootFind root(max_iter); const Real mu = root.regula_falsi(res, 0.0, 1.0, rel_tolerance, v(c2p_mu)); v(c2p_mu) = mu; From 44c59d0d49d2f30347cac6448d0e324f325cc24d Mon Sep 17 00:00:00 2001 From: Logan Alexandria White Date: Fri, 13 Mar 2026 10:10:00 -0700 Subject: [PATCH 02/12] first draft of h0 calculation --- src/fluid/con2prim_robust.hpp | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/src/fluid/con2prim_robust.hpp b/src/fluid/con2prim_robust.hpp index f2759cad..05ba21ee 100644 --- a/src/fluid/con2prim_robust.hpp +++ b/src/fluid/con2prim_robust.hpp @@ -296,6 +296,20 @@ class ConToPrim { const bool fail_on_floors_; const bool fail_on_ceilings_; + KOKKOS_INLINE_FUNCTION + Real calc_h0sq( const Microphsyics::EOS::EOS &eos, Real Ye ) { + Real tmin, rhomin, epsmin, pmin, h0; + Real eos_lambda[2] = {Ye, 0.0}; // probably shouldn't hardcode this in the future; following current method below + + tmin = eos.MinimumTemperature(); + rhomin = eos.MinimumDensity(); + epsmin = eos.InternalEnergyFromDensityTemperature(rhomin, tmin, eos_lambda); + pmin = eos.PressureFromDensityTemperature(rhomin, tmin); + + h0 = 1 + epsmin + robust::ratio(pmin, rhomin); // lowest bound for enthalpy in eos at given ye + return h0 * h0; + } + KOKKOS_INLINE_FUNCTION ConToPrimStatus solve(const VarAccessor &v, const CellGeom &g, const Microphysics::EOS::EOS &eos, const Real x1, const Real x2, @@ -360,6 +374,10 @@ class ConToPrim { rbsq = bdotr * bdotr; bsq_rpsq = bsq * rsq - rbsq; } + + // finding a correct lower bound for enthalpy + h0sq_ = calc_h0sq( eos, ye_local); //TODO: check to see if this works since h0sq_ is defined as const. + const Real zsq = rsq / h0sq_; Real v0sq = std::min(zsq / (1.0 + zsq), 1.0 - 1.0 / (gam_max * gam_max)); @@ -369,7 +387,7 @@ class ConToPrim { // find the upper bound // TODO(JCD): revisit this. is it worth it to find the upper bound? // Doesn't seem to be at a quick glance. - // const Real mu_r = res.compute_upper_bound(h0sq_); + const Real mu_r = res.compute_upper_bound(h0sq_); // solve /** @@ -383,7 +401,7 @@ class ConToPrim { */ root_find::RootFind root(max_iter); - const Real mu = root.regula_falsi(res, 0.0, 1.0, rel_tolerance, v(c2p_mu)); + const Real mu = root.regula_falsi(res, 0.0, mu_r, rel_tolerance, v(c2p_mu)); v(c2p_mu) = mu; #if CON2PRIM_STATISTICS con2prim_statistics::Stats::add(root.iteration_count); From f70c93208e1cb78969e09a273e6f6c6e3e8b901d Mon Sep 17 00:00:00 2001 From: Logan Alexandria White Date: Fri, 13 Mar 2026 12:05:03 -0700 Subject: [PATCH 03/12] fixed scope; added method to Residual --- src/fluid/con2prim_robust.hpp | 94 +++++++++++++++++++++++++---------- 1 file changed, 67 insertions(+), 27 deletions(-) diff --git a/src/fluid/con2prim_robust.hpp b/src/fluid/con2prim_robust.hpp index 05ba21ee..bd61040d 100644 --- a/src/fluid/con2prim_robust.hpp +++ b/src/fluid/con2prim_robust.hpp @@ -45,6 +45,8 @@ struct FailFlags { class Residual { public: + + // TODO: check if this causes compile error (leaving const variables in constructor def.) KOKKOS_FUNCTION Residual(const Real D, const Real q, const Real bsq, const Real bsq_rpsq, const Real rsq, const Real rbsq, const Real v0sq, const Real Ye, @@ -61,6 +63,28 @@ class Residual { rho_floor_ *= floor_scale_fac_; } + // NEW: overloaded constructor that handles calculation of h0sq, v0sq internally + KOKKOS_FUNCTION + Residual(const Real D, const Real q, const Real bsq, const Real bsq_rpsq, + const Real rsq, const Real rbsq, const Real Ye, + const Microphysics::EOS::EOS &eos, const fixup::Bounds &bnds, const Real x1, + const Real x2, const Real x3, const Real floor_scale_fac) + : D_(D), q_(q), bsq_(bsq), bsq_rpsq_(bsq_rpsq), rsq_(rsq), rbsq_(rbsq), + eos_(eos), bounds_(bnds), x1_(x1), x2_(x2), x3_(x3), + floor_scale_fac_(floor_scale_fac) { + + lambda_[0] = Ye; + Real garbage = 0.0; + bounds_.GetFloors(x1_, x2_, x3_, rho_floor_, garbage); + bounds_.GetCeilings(x1_, x2_, x3_, gam_max_, e_max_); + + h0sq_ = calc_h0sq(); + Real zsq_ = rsq_ / h0sq_; // TODO: check that nothing breaks in this normalization. + v0sq_ = std::min(zsq_ / (1.0 + zsq_), 1.0 - 1.0 / (gam_max_ * gam_max_)); + + rho_floor_ *= floor_scale_fac_; + } + KOKKOS_FORCEINLINE_FUNCTION Real x_mu(const Real mu) { return robust::ratio(1.0, 1.0 + mu * bsq_); } KOKKOS_FORCEINLINE_FUNCTION @@ -140,11 +164,18 @@ class Residual { return mu - muhat; } + // KOKKOS_INLINE_FUNCTION + // Real compute_upper_bound(const Real h0sq) { + // auto func = [=](const Real x) { return aux_func(x, h0sq); }; + // root_find::RootFind root; + // return root.itp(func, 1.e-16, 1.0 / sqrt(h0sq), 1.e-3, -1.0); + // } + KOKKOS_INLINE_FUNCTION - Real compute_upper_bound(const Real h0sq) { - auto func = [=](const Real x) { return aux_func(x, h0sq); }; + Real compute_upper_bound() { + auto func = [=](const Real x) { return aux_func(x); }; root_find::RootFind root; - return root.itp(func, 1.e-16, 1.0 / sqrt(h0sq), 1.e-3, -1.0); + return root.itp(func, 1.e-16, 1.0 / sqrt(h0sq_), 1.e-3, -1.0); } KOKKOS_INLINE_FUNCTION @@ -162,7 +193,9 @@ class Residual { } private: - const Real D_, q_, bsq_, bsq_rpsq_, rsq_, rbsq_, v0sq_; + // const Real D_, q_, bsq_, bsq_rpsq_, rsq_, rbsq_, v0sq_; + const Real D_, q_, bsq_, bsq_rpsq_, rsq_, rbsq_; + Real h0sq_, v0sq_; // these cannot be const, we need to set them after initialization. const Microphysics::EOS::EOS &eos_; const fixup::Bounds bounds_; const Real x1_, x2_, x3_; @@ -171,12 +204,33 @@ class Residual { Real rho_floor_, e_floor_, gam_max_, e_max_; bool used_density_floor_, used_energy_floor_, used_energy_max_, used_gamma_max_; + // KOKKOS_FORCEINLINE_FUNCTION + // Real aux_func(const Real mu, const Real h0sq) { + // const Real x = 1.0 / (1.0 + mu * bsq_); + // const Real rbarsq = x * (rsq_ * x + mu * (1.0 + x) * rbsq_); + // return mu * std::sqrt(h0sq + rbarsq) - 1.0; + // } + KOKKOS_FORCEINLINE_FUNCTION - Real aux_func(const Real mu, const Real h0sq) { + Real aux_func(const Real mu) { const Real x = 1.0 / (1.0 + mu * bsq_); const Real rbarsq = x * (rsq_ * x + mu * (1.0 + x) * rbsq_); - return mu * std::sqrt(h0sq + rbarsq) - 1.0; + return mu * std::sqrt(h0sq_ + rbarsq) - 1.0; } + + KOKKOS_INLINE_FUNCTION + Real calc_h0sq() { + Real tmin, rhomin, epsmin, pmin, h0; + + tmin = eos_.MinimumTemperature(); + rhomin = eos_.MinimumDensity(); + epsmin = eos_.InternalEnergyFromDensityTemperature(rhomin, tmin, lambda_); + pmin = eos_.PressureFromDensityTemperature(rhomin, tmin); + + h0 = 1 + epsmin + robust::ratio(pmin, rhomin); // lowest bound for enthalpy in eos at given ye + return h0 * h0; + } + }; template @@ -296,20 +350,6 @@ class ConToPrim { const bool fail_on_floors_; const bool fail_on_ceilings_; - KOKKOS_INLINE_FUNCTION - Real calc_h0sq( const Microphsyics::EOS::EOS &eos, Real Ye ) { - Real tmin, rhomin, epsmin, pmin, h0; - Real eos_lambda[2] = {Ye, 0.0}; // probably shouldn't hardcode this in the future; following current method below - - tmin = eos.MinimumTemperature(); - rhomin = eos.MinimumDensity(); - epsmin = eos.InternalEnergyFromDensityTemperature(rhomin, tmin, eos_lambda); - pmin = eos.PressureFromDensityTemperature(rhomin, tmin); - - h0 = 1 + epsmin + robust::ratio(pmin, rhomin); // lowest bound for enthalpy in eos at given ye - return h0 * h0; - } - KOKKOS_INLINE_FUNCTION ConToPrimStatus solve(const VarAccessor &v, const CellGeom &g, const Microphysics::EOS::EOS &eos, const Real x1, const Real x2, @@ -374,20 +414,20 @@ class ConToPrim { rbsq = bdotr * bdotr; bsq_rpsq = bsq * rsq - rbsq; } - - // finding a correct lower bound for enthalpy - h0sq_ = calc_h0sq( eos, ye_local); //TODO: check to see if this works since h0sq_ is defined as const. - const Real zsq = rsq / h0sq_; - Real v0sq = std::min(zsq / (1.0 + zsq), 1.0 - 1.0 / (gam_max * gam_max)); + // const Real zsq = rsq / h0sq_; // TODO: check that nothing breaks in this normalization. + // Real v0sq = std::min(zsq / (1.0 + zsq), 1.0 - 1.0 / (gam_max * gam_max)); + + // Residual res(D, q, bsq, bsq_rpsq, rsq, rbsq, v0sq, ye_local, eos, bounds, x1, x2, x3, + // floor_scale_fac_); - Residual res(D, q, bsq, bsq_rpsq, rsq, rbsq, v0sq, ye_local, eos, bounds, x1, x2, x3, + Residual res(D, q, bsq, bsq_rpsq, rsq, rbsq, ye_local, eos, bounds, x1, x2, x3, floor_scale_fac_); // find the upper bound // TODO(JCD): revisit this. is it worth it to find the upper bound? // Doesn't seem to be at a quick glance. - const Real mu_r = res.compute_upper_bound(h0sq_); + const Real mu_r = res.compute_upper_bound(); // solve /** From 37ee67ef065046c5ca631c534ae04ab43b60f70b Mon Sep 17 00:00:00 2001 From: Logan Alexandria White Date: Fri, 13 Mar 2026 12:26:04 -0700 Subject: [PATCH 04/12] verbose basics for testing --- src/fluid/con2prim_robust.hpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/fluid/con2prim_robust.hpp b/src/fluid/con2prim_robust.hpp index bd61040d..02dcece9 100644 --- a/src/fluid/con2prim_robust.hpp +++ b/src/fluid/con2prim_robust.hpp @@ -192,6 +192,12 @@ class Residual { used_gamma_max_); } + // FOR VERBOSE TESTING ONLY + KOKKOS_INLINE_FUNCTION + Real get_h0sq() { + return h0sq_; + } + private: // const Real D_, q_, bsq_, bsq_rpsq_, rsq_, rbsq_, v0sq_; const Real D_, q_, bsq_, bsq_rpsq_, rsq_, rbsq_; @@ -430,6 +436,9 @@ class ConToPrim { const Real mu_r = res.compute_upper_bound(); // solve + // TESTING, VERBOSE OUTPUT + printf("%-8.5e %-8.5e\n", res.get_h0sq(), mu_r); + /** * TODO: implement method to find lower enthalpy bound h0 (upper root find bound) * From 8a535d06be8c4ea265c2788d66166342718c0871 Mon Sep 17 00:00:00 2001 From: Logan Alexandria White Date: Mon, 16 Mar 2026 16:21:06 -0700 Subject: [PATCH 05/12] works to initialize ccsne (yay) --- src/fluid/con2prim_robust.hpp | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) diff --git a/src/fluid/con2prim_robust.hpp b/src/fluid/con2prim_robust.hpp index 02dcece9..bf08f556 100644 --- a/src/fluid/con2prim_robust.hpp +++ b/src/fluid/con2prim_robust.hpp @@ -46,7 +46,6 @@ struct FailFlags { class Residual { public: - // TODO: check if this causes compile error (leaving const variables in constructor def.) KOKKOS_FUNCTION Residual(const Real D, const Real q, const Real bsq, const Real bsq_rpsq, const Real rsq, const Real rbsq, const Real v0sq, const Real Ye, @@ -164,13 +163,6 @@ class Residual { return mu - muhat; } - // KOKKOS_INLINE_FUNCTION - // Real compute_upper_bound(const Real h0sq) { - // auto func = [=](const Real x) { return aux_func(x, h0sq); }; - // root_find::RootFind root; - // return root.itp(func, 1.e-16, 1.0 / sqrt(h0sq), 1.e-3, -1.0); - // } - KOKKOS_INLINE_FUNCTION Real compute_upper_bound() { auto func = [=](const Real x) { return aux_func(x); }; @@ -231,8 +223,8 @@ class Residual { tmin = eos_.MinimumTemperature(); rhomin = eos_.MinimumDensity(); epsmin = eos_.InternalEnergyFromDensityTemperature(rhomin, tmin, lambda_); - pmin = eos_.PressureFromDensityTemperature(rhomin, tmin); - + pmin = eos_.PressureFromDensityTemperature(rhomin, tmin, lambda_); + h0 = 1 + epsmin + robust::ratio(pmin, rhomin); // lowest bound for enthalpy in eos at given ye return h0 * h0; } @@ -421,12 +413,6 @@ class ConToPrim { bsq_rpsq = bsq * rsq - rbsq; } - // const Real zsq = rsq / h0sq_; // TODO: check that nothing breaks in this normalization. - // Real v0sq = std::min(zsq / (1.0 + zsq), 1.0 - 1.0 / (gam_max * gam_max)); - - // Residual res(D, q, bsq, bsq_rpsq, rsq, rbsq, v0sq, ye_local, eos, bounds, x1, x2, x3, - // floor_scale_fac_); - Residual res(D, q, bsq, bsq_rpsq, rsq, rbsq, ye_local, eos, bounds, x1, x2, x3, floor_scale_fac_); From 0768965e7a19e6c0ffd17ce644782f45b2132842 Mon Sep 17 00:00:00 2001 From: Logan Alexandria White Date: Thu, 19 Mar 2026 13:23:58 -0700 Subject: [PATCH 06/12] new bounds condition from kastaun+21, highly rel. case --- src/fluid/con2prim_robust.hpp | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/src/fluid/con2prim_robust.hpp b/src/fluid/con2prim_robust.hpp index bf08f556..87eef830 100644 --- a/src/fluid/con2prim_robust.hpp +++ b/src/fluid/con2prim_robust.hpp @@ -202,13 +202,6 @@ class Residual { Real rho_floor_, e_floor_, gam_max_, e_max_; bool used_density_floor_, used_energy_floor_, used_energy_max_, used_gamma_max_; - // KOKKOS_FORCEINLINE_FUNCTION - // Real aux_func(const Real mu, const Real h0sq) { - // const Real x = 1.0 / (1.0 + mu * bsq_); - // const Real rbarsq = x * (rsq_ * x + mu * (1.0 + x) * rbsq_); - // return mu * std::sqrt(h0sq + rbarsq) - 1.0; - // } - KOKKOS_FORCEINLINE_FUNCTION Real aux_func(const Real mu) { const Real x = 1.0 / (1.0 + mu * bsq_); @@ -419,11 +412,14 @@ class ConToPrim { // find the upper bound // TODO(JCD): revisit this. is it worth it to find the upper bound? // Doesn't seem to be at a quick glance. - const Real mu_r = res.compute_upper_bound(); - // solve - // TESTING, VERBOSE OUTPUT - printf("%-8.5e %-8.5e\n", res.get_h0sq(), mu_r); + // conditional from Kastaun et al. 2021, we need a tighter upper bound if r > h0 due to sharp kink in lorentz factor. + Real mu_r; + if ( rsq >= res.get_h0sq() ) + mu_r = res.compute_upper_bound(); + else + mu_r = 1 / sqrt(res.get_h0sq()); + // solve /** * TODO: implement method to find lower enthalpy bound h0 (upper root find bound) @@ -454,8 +450,10 @@ class ConToPrim { if (pye > 0) eos_lambda[0] = v(pye); eos_lambda[1] = std::log10(v(tmp)); // initial guess v(peng) = res.ehat_mu(mu, qbar, rbarsq, vsq, W); + // TESTING, VERBOSE OUTPUT + printf("con2prim: %-8.5e %-8.5e %-8.5e %-8.5e %-8.5e\n", res.get_h0sq(), rsq, mu_r, v(peng), ye_local); v(tmp) = eos.TemperatureFromDensityInternalEnergy(v(prho), v(peng), eos_lambda); - v(peng) *= v(prho); + v(peng) *= v(prho); // conversion to u v(prs) = eos.PressureFromDensityTemperature(v(prho), v(tmp), eos_lambda); v(gm1) = eos.BulkModulusFromDensityTemperature(v(prho), v(tmp), eos_lambda) / v(prs); PARTHENON_DEBUG_REQUIRE(v(prs) > robust::SMALL(), "Pressure must be positive"); From b3bcb5a7b785bcd07d3de34a603135306f9ce5ac Mon Sep 17 00:00:00 2001 From: Logan Alexandria White Date: Tue, 24 Mar 2026 23:37:05 -0700 Subject: [PATCH 07/12] refactor, global lower bound for enthalpy --- src/fixup/fixup.cpp | 5 +- src/fluid/con2prim_robust.hpp | 130 ++++++++++++++----- src/fluid/fluid.cpp | 3 +- src/microphysics/eos_phoebus/eos_phoebus.cpp | 47 +++++++ 4 files changed, 150 insertions(+), 35 deletions(-) diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index cec80322..dc9565bb 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -349,6 +349,7 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { auto eos = eos_pkg->Param("d.EOS"); auto geom = Geometry::GetCoordinateSystem(rc); Bounds *pbounds = fix_pkg->MutableParam("bounds"); + const Real h_min = eos_pkg->Param("h_min"); // BLB: Setting EOS bnds for Ceilings/Floors here. pbounds->SetEOSBnds(eos_pkg); @@ -439,6 +440,8 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { du = std::max(du, sie_floor * drho); } + // printf("\tfixup.cpp: %-8.5e %-8.5e %-8.5e %-8.5e\n", sie_floor, u_floor_max, du, drho); + Real dcrho, dS[3], dBcons[3], dtau, dyecons; Real bp[3] = {0}; if (pb_hi > 0) { @@ -469,7 +472,7 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { } // fluid c2p - auto status = invert(geom, eos, coords, k, j, i); + auto status = invert(geom, eos, coords, h_min, k, j, i); if (status == con2prim_robust::ConToPrimStatus::failure) { // If fluid c2p fails, set to floors v(b, prho, k, j, i) = drho; diff --git a/src/fluid/con2prim_robust.hpp b/src/fluid/con2prim_robust.hpp index 87eef830..291668c2 100644 --- a/src/fluid/con2prim_robust.hpp +++ b/src/fluid/con2prim_robust.hpp @@ -46,29 +46,50 @@ struct FailFlags { class Residual { public: - KOKKOS_FUNCTION - Residual(const Real D, const Real q, const Real bsq, const Real bsq_rpsq, - const Real rsq, const Real rbsq, const Real v0sq, const Real Ye, - const Microphysics::EOS::EOS &eos, const fixup::Bounds &bnds, const Real x1, - const Real x2, const Real x3, const Real floor_scale_fac) - : D_(D), q_(q), bsq_(bsq), bsq_rpsq_(bsq_rpsq), rsq_(rsq), rbsq_(rbsq), v0sq_(v0sq), - eos_(eos), bounds_(bnds), x1_(x1), x2_(x2), x3_(x3), - floor_scale_fac_(floor_scale_fac) { - lambda_[0] = Ye; - Real garbage = 0.0; - bounds_.GetFloors(x1_, x2_, x3_, rho_floor_, garbage); - bounds_.GetCeilings(x1_, x2_, x3_, gam_max_, e_max_); - - rho_floor_ *= floor_scale_fac_; - } + // KOKKOS_FUNCTION + // Residual(const Real D, const Real q, const Real bsq, const Real bsq_rpsq, + // const Real rsq, const Real rbsq, const Real v0sq, const Real Ye, + // const Microphysics::EOS::EOS &eos, const fixup::Bounds &bnds, const Real x1, + // const Real x2, const Real x3, const Real floor_scale_fac) + // : D_(D), q_(q), bsq_(bsq), bsq_rpsq_(bsq_rpsq), rsq_(rsq), rbsq_(rbsq), v0sq_(v0sq), + // eos_(eos), bounds_(bnds), x1_(x1), x2_(x2), x3_(x3), + // floor_scale_fac_(floor_scale_fac) { + // lambda_[0] = Ye; + // Real garbage = 0.0; + // bounds_.GetFloors(x1_, x2_, x3_, rho_floor_, garbage); + // bounds_.GetCeilings(x1_, x2_, x3_, gam_max_, e_max_); + + // rho_floor_ *= floor_scale_fac_; + // } // NEW: overloaded constructor that handles calculation of h0sq, v0sq internally + // KOKKOS_FUNCTION + // Residual(const Real D, const Real q, const Real bsq, const Real bsq_rpsq, + // const Real rsq, const Real rbsq, const Real Ye, + // const Microphysics::EOS::EOS &eos, const fixup::Bounds &bnds, const Real x1, + // const Real x2, const Real x3, const Real floor_scale_fac) + // : D_(D), q_(q), bsq_(bsq), bsq_rpsq_(bsq_rpsq), rsq_(rsq), rbsq_(rbsq), + // eos_(eos), bounds_(bnds), x1_(x1), x2_(x2), x3_(x3), + // floor_scale_fac_(floor_scale_fac) { + + // lambda_[0] = Ye; + // Real garbage = 0.0; + // bounds_.GetFloors(x1_, x2_, x3_, rho_floor_, garbage); + // bounds_.GetCeilings(x1_, x2_, x3_, gam_max_, e_max_); + + // h0sq_ = calc_h0sq(); + // Real zsq_ = rsq_ / h0sq_; // TODO: check that nothing breaks in this normalization. + // v0sq_ = std::min(zsq_ / (1.0 + zsq_), 1.0 - 1.0 / (gam_max_ * gam_max_)); + + // rho_floor_ *= floor_scale_fac_; + // } + KOKKOS_FUNCTION Residual(const Real D, const Real q, const Real bsq, const Real bsq_rpsq, - const Real rsq, const Real rbsq, const Real Ye, + const Real rsq, const Real rbsq, const Real h0, const Real Ye, const Microphysics::EOS::EOS &eos, const fixup::Bounds &bnds, const Real x1, const Real x2, const Real x3, const Real floor_scale_fac) - : D_(D), q_(q), bsq_(bsq), bsq_rpsq_(bsq_rpsq), rsq_(rsq), rbsq_(rbsq), + : D_(D), q_(q), bsq_(bsq), bsq_rpsq_(bsq_rpsq), rsq_(rsq), rbsq_(rbsq), h0sq_(h0 * h0), eos_(eos), bounds_(bnds), x1_(x1), x2_(x2), x3_(x3), floor_scale_fac_(floor_scale_fac) { @@ -77,7 +98,6 @@ class Residual { bounds_.GetFloors(x1_, x2_, x3_, rho_floor_, garbage); bounds_.GetCeilings(x1_, x2_, x3_, gam_max_, e_max_); - h0sq_ = calc_h0sq(); Real zsq_ = rsq_ / h0sq_; // TODO: check that nothing breaks in this normalization. v0sq_ = std::min(zsq_ / (1.0 + zsq_), 1.0 - 1.0 / (gam_max_ * gam_max_)); @@ -209,18 +229,48 @@ class Residual { return mu * std::sqrt(h0sq_ + rbarsq) - 1.0; } - KOKKOS_INLINE_FUNCTION - Real calc_h0sq() { - Real tmin, rhomin, epsmin, pmin, h0; - - tmin = eos_.MinimumTemperature(); - rhomin = eos_.MinimumDensity(); - epsmin = eos_.InternalEnergyFromDensityTemperature(rhomin, tmin, lambda_); - pmin = eos_.PressureFromDensityTemperature(rhomin, tmin, lambda_); - - h0 = 1 + epsmin + robust::ratio(pmin, rhomin); // lowest bound for enthalpy in eos at given ye - return h0 * h0; - } + // KOKKOS_INLINE_FUNCTION + // Real calc_h0sq() { // TODO: add option here for global search or edge case... + // int n; // should be set to resolve the highest dimension of the eos table, usually density. + // Real T, rho, ye, dT, drho, dye; + // Real eps, P, h0; + // Real lambda[2]; + + // n = 250; // maybe this shouldn't be hardcoded in the future? + // h0 = 1.0; // initial guess for minimum enthalpy, sufficient in ideal cases. + + // // spacing for each dimension (i.e. np.linspace) + // dT = (eos_.TMax() - eos_.TMin()) / n; + // drho = (eos_.rhoMax() - eos_.rhoMin()) / n; + // dye = (eos_.YeMax() - eos_.YeMin()) / (n / 2); // we don't need to resolve ye as much + + // T = eos_.TMin(); + // rho = eos_.rhoMin(); + // ye = eos_.YeMin(); + + // // WIP: update this to find a global lower bound (still assuming it lies along the minimum edge of the SC-EOS table) + // // realistically this should be something that happens once in something like singularity-EOS and is then callable from there... + // // is there a way to refactor this to be cleaner? this is bad readability. + + // LOOP(y, n/2) { + // lambda[0] = ye; + + // LOOP(r, n) { + // LOOP(t, n) { + + // eps = eos_.InternalEnergyFromDensityTemperature(r, T, lambda); + // P = eos_.PressureFromDensityTemperature(r, T, lambda); + // h0 = std::min(h0, 1 + eps + robust::ratio(P, r)); + + // T += dT; + // } + // rho += drho; + // } + // ye += dye; + // } + + // return h0 * h0; + // } }; @@ -320,6 +370,19 @@ class ConToPrim { return solve(v, g, eos, x1, x2, x3); } + template + KOKKOS_INLINE_FUNCTION ConToPrimStatus operator()(const CoordinateSystem &geom, + const Microphysics::EOS::EOS &eos, + const Coordinates_t &coords, const Real h0, + Args &&...args) const { + VarAccessor v(var, std::forward(args)...); + CellGeom g(geom, std::forward(args)...); + Real x1 = coords.Xc<1>(std::forward(args)...); + Real x2 = coords.Xc<2>(std::forward(args)...); + Real x3 = coords.Xc<3>(std::forward(args)...); + return solve(v, g, eos, x1, x2, x3, h0); + } + int NumBlocks() { return var.GetDim(5); } private: @@ -343,8 +406,9 @@ class ConToPrim { KOKKOS_INLINE_FUNCTION ConToPrimStatus solve(const VarAccessor &v, const CellGeom &g, - const Microphysics::EOS::EOS &eos, const Real x1, const Real x2, - const Real x3) const { + const Microphysics::EOS::EOS &eos, + const Real x1, const Real x2, + const Real x3, const Real h0 = 1.0) const { int num_nans = std::isnan(v(crho)) + std::isnan(v(cmom_lo)) + std::isnan(v(ceng)); if (num_nans > 0) return ConToPrimStatus::failure; const Real igdet = 1.0 / g.gdet; @@ -406,7 +470,7 @@ class ConToPrim { bsq_rpsq = bsq * rsq - rbsq; } - Residual res(D, q, bsq, bsq_rpsq, rsq, rbsq, ye_local, eos, bounds, x1, x2, x3, + Residual res(D, q, bsq, bsq_rpsq, rsq, rbsq, h0, ye_local, eos, bounds, x1, x2, x3, floor_scale_fac_); // find the upper bound diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp index ab3117fa..0161758d 100644 --- a/src/fluid/fluid.cpp +++ b/src/fluid/fluid.cpp @@ -499,6 +499,7 @@ TaskStatus ConservedToPrimitiveRobust(T *rc, const IndexRange &ib, const IndexRa auto eos = eos_pkg->Param("d.EOS"); auto geom = Geometry::GetCoordinateSystem(rc); auto coords = pmb->coords; + const Real h_min = eos_pkg->Param("h_min"); // TODO(JCD): move the setting of this into the solver so we can call this on MeshData auto fail = rc->Get(internal_variables::fail::name()).data; @@ -507,7 +508,7 @@ TaskStatus ConservedToPrimitiveRobust(T *rc, const IndexRange &ib, const IndexRa DEFAULT_LOOP_PATTERN, "ConToPrim::Solve", DevExecSpace(), 0, invert.NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { - auto status = invert(geom, eos, coords, k, j, i); + auto status = invert(geom, eos, coords, h_min, k, j, i); fail(k, j, i) = (status == con2prim_robust::ConToPrimStatus::success ? con2prim_robust::FailFlags::success : con2prim_robust::FailFlags::fail); diff --git a/src/microphysics/eos_phoebus/eos_phoebus.cpp b/src/microphysics/eos_phoebus/eos_phoebus.cpp index 76372b58..1b2d45bb 100644 --- a/src/microphysics/eos_phoebus/eos_phoebus.cpp +++ b/src/microphysics/eos_phoebus/eos_phoebus.cpp @@ -25,6 +25,9 @@ #include "phoebus_utils/unit_conversions.hpp" #include "phoebus_utils/variables.hpp" +// inspired by SPACELOOP in geometry utils, util for global lower bound for enthalpy +#define LOOP(i, n) for (int i = 0; i < n; i++) + using namespace singularity; namespace Microphysics { @@ -65,6 +68,9 @@ std::shared_ptr Initialize(ParameterInput *pin) { Real T_max; Real ye_min; Real ye_max; + Real h_min; // lower enthalpy bound + Real T, dT, rho, drho, ye, dye, eps, P; // for our bound search, temps + int n; std::string eos_type = pin->GetString(block_name, std::string("type")); params.Add("type", eos_type); @@ -93,6 +99,8 @@ std::shared_ptr Initialize(ParameterInput *pin) { T_max = eos_host.TemperatureFromDensityInternalEnergy(rho_max, sie_max, lambda); ye_min = 0.01; ye_max = 1.0; + h_min = 1.0; // is this the right guess for an ideal gas case? + #ifdef SPINER_USE_HDF } else if (eos_type == StellarCollapse::EosType()) { // We request that Ye and temperature exist, but do not provide them. @@ -154,6 +162,42 @@ std::shared_ptr Initialize(ParameterInput *pin) { rho_max = eos_sc.rhoMax() / rho_unit; ye_min = eos_sc.YeMin(); ye_max = eos_sc.YeMax(); + + + // new calculation to find lower bound for enthalpy! + n = 250; // maybe this shouldn't be hardcoded in the future? + h_min = 1.0; // initial guess for minimum enthalpy, sufficient in ideal cases. + eps = 0.0; + + // spacing for each dimension (i.e. np.linspace) + dT = (T_max - T_min) / n; + drho = (rho_max - rho_min) / n; + dye = (ye_max - ye_min) / (n / 2); // we don't need to resolve ye as much + + T = T_min; + rho = rho_min; + ye = ye_min; + + // WIP: update this to find a global lower bound (still assuming it lies along the minimum edge of the SC-EOS table) + // realistically this should be something that happens once in something like singularity-EOS and is then callable from there... + // is there a way to refactor this to be cleaner? this is bad readability. + + LOOP(y, n/2) { + lambda[0] = ye; + LOOP(r, n) { + LOOP(t, n) { + + eps = eos_sc.InternalEnergyFromDensityTemperature(rho, T, lambda) / sie_unit; + P = eos_sc.PressureFromDensityTemperature(rho, T, lambda) / press_unit; + h_min = std::min(h_min, 1 + eps + robust::ratio(P, rho)); + + T += dT; + } + rho += drho; + } + ye += dye; + } + #endif } else { std::stringstream error_mesg; @@ -166,6 +210,8 @@ std::shared_ptr Initialize(ParameterInput *pin) { PARTHENON_THROW(error_mesg); } + printf("h0: %5.8e\n\n", h_min); // VERBOSE + params.Add("needs_ye", needs_ye); params.Add("provides_entropy", provides_entropy); @@ -177,6 +223,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { params.Add("rho_max", rho_max); params.Add("ye_min", ye_min); params.Add("ye_max", ye_max); + params.Add("h_min", h_min); return pkg; } From 64eca66123989214a646191ee776aa935f1f08d9 Mon Sep 17 00:00:00 2001 From: Logan Alexandria White Date: Wed, 25 Mar 2026 08:49:53 -0700 Subject: [PATCH 08/12] general clean up, etc --- src/fluid/con2prim_robust.hpp | 94 +------------------- src/microphysics/eos_phoebus/eos_phoebus.cpp | 12 +-- 2 files changed, 7 insertions(+), 99 deletions(-) diff --git a/src/fluid/con2prim_robust.hpp b/src/fluid/con2prim_robust.hpp index 291668c2..e11670b1 100644 --- a/src/fluid/con2prim_robust.hpp +++ b/src/fluid/con2prim_robust.hpp @@ -46,44 +46,7 @@ struct FailFlags { class Residual { public: - // KOKKOS_FUNCTION - // Residual(const Real D, const Real q, const Real bsq, const Real bsq_rpsq, - // const Real rsq, const Real rbsq, const Real v0sq, const Real Ye, - // const Microphysics::EOS::EOS &eos, const fixup::Bounds &bnds, const Real x1, - // const Real x2, const Real x3, const Real floor_scale_fac) - // : D_(D), q_(q), bsq_(bsq), bsq_rpsq_(bsq_rpsq), rsq_(rsq), rbsq_(rbsq), v0sq_(v0sq), - // eos_(eos), bounds_(bnds), x1_(x1), x2_(x2), x3_(x3), - // floor_scale_fac_(floor_scale_fac) { - // lambda_[0] = Ye; - // Real garbage = 0.0; - // bounds_.GetFloors(x1_, x2_, x3_, rho_floor_, garbage); - // bounds_.GetCeilings(x1_, x2_, x3_, gam_max_, e_max_); - - // rho_floor_ *= floor_scale_fac_; - // } - - // NEW: overloaded constructor that handles calculation of h0sq, v0sq internally - // KOKKOS_FUNCTION - // Residual(const Real D, const Real q, const Real bsq, const Real bsq_rpsq, - // const Real rsq, const Real rbsq, const Real Ye, - // const Microphysics::EOS::EOS &eos, const fixup::Bounds &bnds, const Real x1, - // const Real x2, const Real x3, const Real floor_scale_fac) - // : D_(D), q_(q), bsq_(bsq), bsq_rpsq_(bsq_rpsq), rsq_(rsq), rbsq_(rbsq), - // eos_(eos), bounds_(bnds), x1_(x1), x2_(x2), x3_(x3), - // floor_scale_fac_(floor_scale_fac) { - - // lambda_[0] = Ye; - // Real garbage = 0.0; - // bounds_.GetFloors(x1_, x2_, x3_, rho_floor_, garbage); - // bounds_.GetCeilings(x1_, x2_, x3_, gam_max_, e_max_); - - // h0sq_ = calc_h0sq(); - // Real zsq_ = rsq_ / h0sq_; // TODO: check that nothing breaks in this normalization. - // v0sq_ = std::min(zsq_ / (1.0 + zsq_), 1.0 - 1.0 / (gam_max_ * gam_max_)); - - // rho_floor_ *= floor_scale_fac_; - // } - + // new constructor that allows for a value of h0 to be passed in KOKKOS_FUNCTION Residual(const Real D, const Real q, const Real bsq, const Real bsq_rpsq, const Real rsq, const Real rbsq, const Real h0, const Real Ye, @@ -98,7 +61,7 @@ class Residual { bounds_.GetFloors(x1_, x2_, x3_, rho_floor_, garbage); bounds_.GetCeilings(x1_, x2_, x3_, gam_max_, e_max_); - Real zsq_ = rsq_ / h0sq_; // TODO: check that nothing breaks in this normalization. + Real zsq_ = rsq_ / h0sq_; v0sq_ = std::min(zsq_ / (1.0 + zsq_), 1.0 - 1.0 / (gam_max_ * gam_max_)); rho_floor_ *= floor_scale_fac_; @@ -204,14 +167,13 @@ class Residual { used_gamma_max_); } - // FOR VERBOSE TESTING ONLY + // new accessor for enthalpy lower bound KOKKOS_INLINE_FUNCTION Real get_h0sq() { return h0sq_; } private: - // const Real D_, q_, bsq_, bsq_rpsq_, rsq_, rbsq_, v0sq_; const Real D_, q_, bsq_, bsq_rpsq_, rsq_, rbsq_; Real h0sq_, v0sq_; // these cannot be const, we need to set them after initialization. const Microphysics::EOS::EOS &eos_; @@ -229,49 +191,6 @@ class Residual { return mu * std::sqrt(h0sq_ + rbarsq) - 1.0; } - // KOKKOS_INLINE_FUNCTION - // Real calc_h0sq() { // TODO: add option here for global search or edge case... - // int n; // should be set to resolve the highest dimension of the eos table, usually density. - // Real T, rho, ye, dT, drho, dye; - // Real eps, P, h0; - // Real lambda[2]; - - // n = 250; // maybe this shouldn't be hardcoded in the future? - // h0 = 1.0; // initial guess for minimum enthalpy, sufficient in ideal cases. - - // // spacing for each dimension (i.e. np.linspace) - // dT = (eos_.TMax() - eos_.TMin()) / n; - // drho = (eos_.rhoMax() - eos_.rhoMin()) / n; - // dye = (eos_.YeMax() - eos_.YeMin()) / (n / 2); // we don't need to resolve ye as much - - // T = eos_.TMin(); - // rho = eos_.rhoMin(); - // ye = eos_.YeMin(); - - // // WIP: update this to find a global lower bound (still assuming it lies along the minimum edge of the SC-EOS table) - // // realistically this should be something that happens once in something like singularity-EOS and is then callable from there... - // // is there a way to refactor this to be cleaner? this is bad readability. - - // LOOP(y, n/2) { - // lambda[0] = ye; - - // LOOP(r, n) { - // LOOP(t, n) { - - // eps = eos_.InternalEnergyFromDensityTemperature(r, T, lambda); - // P = eos_.PressureFromDensityTemperature(r, T, lambda); - // h0 = std::min(h0, 1 + eps + robust::ratio(P, r)); - - // T += dT; - // } - // rho += drho; - // } - // ye += dye; - // } - - // return h0 * h0; - // } - }; template @@ -473,17 +392,12 @@ class ConToPrim { Residual res(D, q, bsq, bsq_rpsq, rsq, rbsq, h0, ye_local, eos, bounds, x1, x2, x3, floor_scale_fac_); - // find the upper bound - // TODO(JCD): revisit this. is it worth it to find the upper bound? - // Doesn't seem to be at a quick glance. - // conditional from Kastaun et al. 2021, we need a tighter upper bound if r > h0 due to sharp kink in lorentz factor. Real mu_r; if ( rsq >= res.get_h0sq() ) - mu_r = res.compute_upper_bound(); + mu_r = res.compute_upper_bound(); // we don't always need a second root find else mu_r = 1 / sqrt(res.get_h0sq()); - // solve /** * TODO: implement method to find lower enthalpy bound h0 (upper root find bound) diff --git a/src/microphysics/eos_phoebus/eos_phoebus.cpp b/src/microphysics/eos_phoebus/eos_phoebus.cpp index 1b2d45bb..6b176eb2 100644 --- a/src/microphysics/eos_phoebus/eos_phoebus.cpp +++ b/src/microphysics/eos_phoebus/eos_phoebus.cpp @@ -165,7 +165,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { // new calculation to find lower bound for enthalpy! - n = 250; // maybe this shouldn't be hardcoded in the future? + n = 500; // maybe this shouldn't be hardcoded in the future? h_min = 1.0; // initial guess for minimum enthalpy, sufficient in ideal cases. eps = 0.0; @@ -173,24 +173,18 @@ std::shared_ptr Initialize(ParameterInput *pin) { dT = (T_max - T_min) / n; drho = (rho_max - rho_min) / n; dye = (ye_max - ye_min) / (n / 2); // we don't need to resolve ye as much - + // initial conditions T = T_min; rho = rho_min; ye = ye_min; - // WIP: update this to find a global lower bound (still assuming it lies along the minimum edge of the SC-EOS table) - // realistically this should be something that happens once in something like singularity-EOS and is then callable from there... - // is there a way to refactor this to be cleaner? this is bad readability. - LOOP(y, n/2) { lambda[0] = ye; LOOP(r, n) { LOOP(t, n) { - eps = eos_sc.InternalEnergyFromDensityTemperature(rho, T, lambda) / sie_unit; P = eos_sc.PressureFromDensityTemperature(rho, T, lambda) / press_unit; h_min = std::min(h_min, 1 + eps + robust::ratio(P, rho)); - T += dT; } rho += drho; @@ -210,7 +204,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { PARTHENON_THROW(error_mesg); } - printf("h0: %5.8e\n\n", h_min); // VERBOSE + printf("h0: %5.8e\n\n", h_min); // TESTING, VERBOSE params.Add("needs_ye", needs_ye); params.Add("provides_entropy", provides_entropy); From 9d223314860503547dab96409775ea18af46edf1 Mon Sep 17 00:00:00 2001 From: Logan Alexandria White Date: Wed, 25 Mar 2026 15:28:43 -0700 Subject: [PATCH 09/12] removing verbose testing, comments --- src/fixup/fixup.cpp | 5 ++--- src/fluid/con2prim_robust.hpp | 12 ------------ src/microphysics/eos_phoebus/eos_phoebus.cpp | 2 -- 3 files changed, 2 insertions(+), 17 deletions(-) diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index dc9565bb..921faac7 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -434,14 +434,13 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { Real drho = rho_floor_max - v(b, prho, k, j, i); Real du = u_floor_max - v(b, peng, k, j, i); - if (drho > 0. || du > 0.) { + + if (drho > 0. || du > 0.) { floor_applied = true; drho = std::max(drho, du / sie_floor); du = std::max(du, sie_floor * drho); } - // printf("\tfixup.cpp: %-8.5e %-8.5e %-8.5e %-8.5e\n", sie_floor, u_floor_max, du, drho); - Real dcrho, dS[3], dBcons[3], dtau, dyecons; Real bp[3] = {0}; if (pb_hi > 0) { diff --git a/src/fluid/con2prim_robust.hpp b/src/fluid/con2prim_robust.hpp index e11670b1..134be4a1 100644 --- a/src/fluid/con2prim_robust.hpp +++ b/src/fluid/con2prim_robust.hpp @@ -399,16 +399,6 @@ class ConToPrim { else mu_r = 1 / sqrt(res.get_h0sq()); - /** - * TODO: implement method to find lower enthalpy bound h0 (upper root find bound) - * - * Katsaun et al. 2021 only assumes h > 0, not h ≥ 1; we shouldn't have a hardcoded upper bound here. - * stellar collapse eos can have h < 1 due to included nuclear binding energies (e.g. Steiner et al. - * 2012), this breaks our root find when mapping onto the gr grid at initialization. - * - * can probably be resolved with an eos-conscious method using primitives (P, rho, eps) passed into here, solution tbd. - */ - root_find::RootFind root(max_iter); const Real mu = root.regula_falsi(res, 0.0, mu_r, rel_tolerance, v(c2p_mu)); v(c2p_mu) = mu; @@ -428,8 +418,6 @@ class ConToPrim { if (pye > 0) eos_lambda[0] = v(pye); eos_lambda[1] = std::log10(v(tmp)); // initial guess v(peng) = res.ehat_mu(mu, qbar, rbarsq, vsq, W); - // TESTING, VERBOSE OUTPUT - printf("con2prim: %-8.5e %-8.5e %-8.5e %-8.5e %-8.5e\n", res.get_h0sq(), rsq, mu_r, v(peng), ye_local); v(tmp) = eos.TemperatureFromDensityInternalEnergy(v(prho), v(peng), eos_lambda); v(peng) *= v(prho); // conversion to u v(prs) = eos.PressureFromDensityTemperature(v(prho), v(tmp), eos_lambda); diff --git a/src/microphysics/eos_phoebus/eos_phoebus.cpp b/src/microphysics/eos_phoebus/eos_phoebus.cpp index 6b176eb2..8bac54fe 100644 --- a/src/microphysics/eos_phoebus/eos_phoebus.cpp +++ b/src/microphysics/eos_phoebus/eos_phoebus.cpp @@ -204,8 +204,6 @@ std::shared_ptr Initialize(ParameterInput *pin) { PARTHENON_THROW(error_mesg); } - printf("h0: %5.8e\n\n", h_min); // TESTING, VERBOSE - params.Add("needs_ye", needs_ye); params.Add("provides_entropy", provides_entropy); From 8ef8a954a14cbbd32f2280e0f8a2708ae0ba91a3 Mon Sep 17 00:00:00 2001 From: Brandon Barker Date: Thu, 26 Mar 2026 09:26:11 -0600 Subject: [PATCH 10/12] fix: pin sphinx version for docs --- .github/workflows/docs.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 5f40ed8d..b8e7da6c 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -23,7 +23,7 @@ jobs: run: export DEBIAN_FRONTEND=noninteractive - name: install dependencies run: | - pip install sphinx + pip install "sphinx<8.2" # pinned to make multiversion work pip install sphinx-sitemap pip install sphinx-rtd-theme pip install sphinx-multiversion From bba4e304ebfbb573c536b7f7514ada447fcfe6e5 Mon Sep 17 00:00:00 2001 From: Logan Alexandria White Date: Thu, 26 Mar 2026 08:49:13 -0700 Subject: [PATCH 11/12] clang formatting --- src/fixup/fixup.cpp | 2 +- src/fluid/con2prim_robust.hpp | 24 ++++++++------------ src/microphysics/eos_phoebus/eos_phoebus.cpp | 17 +++++++------- 3 files changed, 19 insertions(+), 24 deletions(-) diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index 921faac7..ff3a5727 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -435,7 +435,7 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { Real drho = rho_floor_max - v(b, prho, k, j, i); Real du = u_floor_max - v(b, peng, k, j, i); - if (drho > 0. || du > 0.) { + if (drho > 0. || du > 0.) { floor_applied = true; drho = std::max(drho, du / sie_floor); du = std::max(du, sie_floor * drho); diff --git a/src/fluid/con2prim_robust.hpp b/src/fluid/con2prim_robust.hpp index 134be4a1..62df9ad3 100644 --- a/src/fluid/con2prim_robust.hpp +++ b/src/fluid/con2prim_robust.hpp @@ -45,15 +45,14 @@ struct FailFlags { class Residual { public: - // new constructor that allows for a value of h0 to be passed in KOKKOS_FUNCTION Residual(const Real D, const Real q, const Real bsq, const Real bsq_rpsq, const Real rsq, const Real rbsq, const Real h0, const Real Ye, const Microphysics::EOS::EOS &eos, const fixup::Bounds &bnds, const Real x1, const Real x2, const Real x3, const Real floor_scale_fac) - : D_(D), q_(q), bsq_(bsq), bsq_rpsq_(bsq_rpsq), rsq_(rsq), rbsq_(rbsq), h0sq_(h0 * h0), - eos_(eos), bounds_(bnds), x1_(x1), x2_(x2), x3_(x3), + : D_(D), q_(q), bsq_(bsq), bsq_rpsq_(bsq_rpsq), rsq_(rsq), rbsq_(rbsq), + h0sq_(h0 * h0), eos_(eos), bounds_(bnds), x1_(x1), x2_(x2), x3_(x3), floor_scale_fac_(floor_scale_fac) { lambda_[0] = Ye; @@ -61,7 +60,7 @@ class Residual { bounds_.GetFloors(x1_, x2_, x3_, rho_floor_, garbage); bounds_.GetCeilings(x1_, x2_, x3_, gam_max_, e_max_); - Real zsq_ = rsq_ / h0sq_; + Real zsq_ = rsq_ / h0sq_; v0sq_ = std::min(zsq_ / (1.0 + zsq_), 1.0 - 1.0 / (gam_max_ * gam_max_)); rho_floor_ *= floor_scale_fac_; @@ -169,9 +168,7 @@ class Residual { // new accessor for enthalpy lower bound KOKKOS_INLINE_FUNCTION - Real get_h0sq() { - return h0sq_; - } + Real get_h0sq() { return h0sq_; } private: const Real D_, q_, bsq_, bsq_rpsq_, rsq_, rbsq_; @@ -190,7 +187,6 @@ class Residual { const Real rbarsq = x * (rsq_ * x + mu * (1.0 + x) * rbsq_); return mu * std::sqrt(h0sq_ + rbarsq) - 1.0; } - }; template @@ -292,8 +288,8 @@ class ConToPrim { template KOKKOS_INLINE_FUNCTION ConToPrimStatus operator()(const CoordinateSystem &geom, const Microphysics::EOS::EOS &eos, - const Coordinates_t &coords, const Real h0, - Args &&...args) const { + const Coordinates_t &coords, + const Real h0, Args &&...args) const { VarAccessor v(var, std::forward(args)...); CellGeom g(geom, std::forward(args)...); Real x1 = coords.Xc<1>(std::forward(args)...); @@ -325,8 +321,7 @@ class ConToPrim { KOKKOS_INLINE_FUNCTION ConToPrimStatus solve(const VarAccessor &v, const CellGeom &g, - const Microphysics::EOS::EOS &eos, - const Real x1, const Real x2, + const Microphysics::EOS::EOS &eos, const Real x1, const Real x2, const Real x3, const Real h0 = 1.0) const { int num_nans = std::isnan(v(crho)) + std::isnan(v(cmom_lo)) + std::isnan(v(ceng)); if (num_nans > 0) return ConToPrimStatus::failure; @@ -392,9 +387,10 @@ class ConToPrim { Residual res(D, q, bsq, bsq_rpsq, rsq, rbsq, h0, ye_local, eos, bounds, x1, x2, x3, floor_scale_fac_); - // conditional from Kastaun et al. 2021, we need a tighter upper bound if r > h0 due to sharp kink in lorentz factor. + // conditional from Kastaun et al. 2021, we need a tighter upper bound if r > h0 due + // to sharp kink in lorentz factor. Real mu_r; - if ( rsq >= res.get_h0sq() ) + if (rsq >= res.get_h0sq()) mu_r = res.compute_upper_bound(); // we don't always need a second root find else mu_r = 1 / sqrt(res.get_h0sq()); diff --git a/src/microphysics/eos_phoebus/eos_phoebus.cpp b/src/microphysics/eos_phoebus/eos_phoebus.cpp index 8bac54fe..e31b3ddc 100644 --- a/src/microphysics/eos_phoebus/eos_phoebus.cpp +++ b/src/microphysics/eos_phoebus/eos_phoebus.cpp @@ -68,7 +68,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { Real T_max; Real ye_min; Real ye_max; - Real h_min; // lower enthalpy bound + Real h_min; // lower enthalpy bound Real T, dT, rho, drho, ye, dye, eps, P; // for our bound search, temps int n; @@ -163,9 +163,8 @@ std::shared_ptr Initialize(ParameterInput *pin) { ye_min = eos_sc.YeMin(); ye_max = eos_sc.YeMax(); - // new calculation to find lower bound for enthalpy! - n = 500; // maybe this shouldn't be hardcoded in the future? + n = 500; // maybe this shouldn't be hardcoded in the future? h_min = 1.0; // initial guess for minimum enthalpy, sufficient in ideal cases. eps = 0.0; @@ -177,15 +176,15 @@ std::shared_ptr Initialize(ParameterInput *pin) { T = T_min; rho = rho_min; ye = ye_min; - - LOOP(y, n/2) { + + LOOP(y, n / 2) { lambda[0] = ye; LOOP(r, n) { LOOP(t, n) { - eps = eos_sc.InternalEnergyFromDensityTemperature(rho, T, lambda) / sie_unit; - P = eos_sc.PressureFromDensityTemperature(rho, T, lambda) / press_unit; - h_min = std::min(h_min, 1 + eps + robust::ratio(P, rho)); - T += dT; + eps = eos_sc.InternalEnergyFromDensityTemperature(rho, T, lambda) / sie_unit; + P = eos_sc.PressureFromDensityTemperature(rho, T, lambda) / press_unit; + h_min = std::min(h_min, 1 + eps + robust::ratio(P, rho)); + T += dT; } rho += drho; } From f3cb5845cdbbea1e07dd64f68159780ae3b99c3c Mon Sep 17 00:00:00 2001 From: Logan Alexandria White Date: Wed, 1 Apr 2026 12:35:51 -0700 Subject: [PATCH 12/12] minor bug fix, additional format cleanup --- src/fluid/con2prim_robust.hpp | 5 +-- src/microphysics/eos_phoebus/eos_phoebus.cpp | 33 ++++++++++---------- 2 files changed, 20 insertions(+), 18 deletions(-) diff --git a/src/fluid/con2prim_robust.hpp b/src/fluid/con2prim_robust.hpp index 62df9ad3..8a6b2c99 100644 --- a/src/fluid/con2prim_robust.hpp +++ b/src/fluid/con2prim_robust.hpp @@ -390,10 +390,11 @@ class ConToPrim { // conditional from Kastaun et al. 2021, we need a tighter upper bound if r > h0 due // to sharp kink in lorentz factor. Real mu_r; - if (rsq >= res.get_h0sq()) + if (rsq >= res.get_h0sq()) { mu_r = res.compute_upper_bound(); // we don't always need a second root find - else + } else { mu_r = 1 / sqrt(res.get_h0sq()); + } root_find::RootFind root(max_iter); const Real mu = root.regula_falsi(res, 0.0, mu_r, rel_tolerance, v(c2p_mu)); diff --git a/src/microphysics/eos_phoebus/eos_phoebus.cpp b/src/microphysics/eos_phoebus/eos_phoebus.cpp index e31b3ddc..984b2948 100644 --- a/src/microphysics/eos_phoebus/eos_phoebus.cpp +++ b/src/microphysics/eos_phoebus/eos_phoebus.cpp @@ -25,9 +25,6 @@ #include "phoebus_utils/unit_conversions.hpp" #include "phoebus_utils/variables.hpp" -// inspired by SPACELOOP in geometry utils, util for global lower bound for enthalpy -#define LOOP(i, n) for (int i = 0; i < n; i++) - using namespace singularity; namespace Microphysics { @@ -164,31 +161,35 @@ std::shared_ptr Initialize(ParameterInput *pin) { ye_max = eos_sc.YeMax(); // new calculation to find lower bound for enthalpy! - n = 500; // maybe this shouldn't be hardcoded in the future? + n = 200; // maybe this shouldn't be hardcoded in the future? h_min = 1.0; // initial guess for minimum enthalpy, sufficient in ideal cases. eps = 0.0; - // spacing for each dimension (i.e. np.linspace) - dT = (T_max - T_min) / n; - drho = (rho_max - rho_min) / n; - dye = (ye_max - ye_min) / (n / 2); // we don't need to resolve ye as much + // spacing for each dimension (logspace for rho, T; linspace for ye) + // we need to use physical (i.e. not scaled) units here for eos_sc calcs + dT = (log10(eos_sc.TMax()) - log10(eos_sc.TMin())) / (n - 1); + drho = (log10(eos_sc.rhoMax()) - log10(eos_sc.rhoMin())) / (n - 1); + dye = (ye_max - ye_min) / ((n / 2) - 1); // we don't need to resolve ye as much + // initial conditions - T = T_min; - rho = rho_min; + T = eos_sc.TMin(); + rho = eos_sc.rhoMin(); ye = ye_min; - LOOP(y, n / 2) { + for (int y = 0; y < n / 2; y++) { lambda[0] = ye; - LOOP(r, n) { - LOOP(t, n) { + for (int r = 0; r < n; r++) { + for (int t = 0; t < n; t++) { eps = eos_sc.InternalEnergyFromDensityTemperature(rho, T, lambda) / sie_unit; P = eos_sc.PressureFromDensityTemperature(rho, T, lambda) / press_unit; h_min = std::min(h_min, 1 + eps + robust::ratio(P, rho)); - T += dT; + T *= pow(10.0, dT); // log spacing } - rho += drho; + T = eos_sc.TMin(); // "reset" + rho *= pow(10.0, drho); // log spacing } - ye += dye; + rho = eos_sc.rhoMin(); // "reset" + ye += dye; // linear spacing } #endif