calculate con2prim enthalpy bound

.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

src/fixup/fixup.cpp

@@ -349,6 +349,7 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) {
   auto eos = eos_pkg->Param<EOS>("d.EOS");
   auto geom = Geometry::GetCoordinateSystem(rc);
   Bounds *pbounds = fix_pkg->MutableParam<Bounds>("bounds");
+  const Real h_min = eos_pkg->Param<Real>("h_min");
 
   // BLB: Setting EOS bnds for Ceilings/Floors here.
   pbounds->SetEOSBnds(eos_pkg);
@@ -433,6 +434,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.) {
           floor_applied = true;
           drho = std::max<Real>(drho, du / sie_floor);
@@ -469,7 +471,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;

src/fluid/con2prim_robust.hpp

@@ -45,19 +45,24 @@ 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 v0sq, 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), v0sq_(v0sq),
-        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;
     Real garbage = 0.0;
     bounds_.GetFloors(x1_, x2_, x3_, rho_floor_, garbage);
     bounds_.GetCeilings(x1_, x2_, x3_, gam_max_, e_max_);
 
+    Real zsq_ = rsq_ / h0sq_;
+    v0sq_ = std::min(zsq_ / (1.0 + zsq_), 1.0 - 1.0 / (gam_max_ * gam_max_));
+
     rho_floor_ *= floor_scale_fac_;
   }
 
@@ -141,10 +146,10 @@ class Residual {
   }
 
   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
@@ -161,8 +166,13 @@ class Residual {
             used_gamma_max_);
   }
 
+  // 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_;
   const fixup::Bounds bounds_;
   const Real x1_, x2_, x3_;
@@ -172,10 +182,10 @@ class Residual {
   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) {
+  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;
   }
 };
 
@@ -275,6 +285,19 @@ class ConToPrim {
     return solve(v, g, eos, x1, x2, x3);
   }
 
+  template <typename CoordinateSystem, class... Args>
+  KOKKOS_INLINE_FUNCTION ConToPrimStatus operator()(const CoordinateSystem &geom,
+                                                    const Microphysics::EOS::EOS &eos,
+                                                    const Coordinates_t &coords,
+                                                    const Real h0, Args &&...args) const {
+    VarAccessor<T> v(var, std::forward<Args>(args)...);
+    CellGeom g(geom, std::forward<Args>(args)...);
+    Real x1 = coords.Xc<1>(std::forward<Args>(args)...);
+    Real x2 = coords.Xc<2>(std::forward<Args>(args)...);
+    Real x3 = coords.Xc<3>(std::forward<Args>(args)...);
+    return solve(v, g, eos, x1, x2, x3, h0);
+  }
+
   int NumBlocks() { return var.GetDim(5); }
 
  private:
@@ -299,7 +322,7 @@ class ConToPrim {
   KOKKOS_INLINE_FUNCTION
   ConToPrimStatus solve(const VarAccessor<T> &v, const CellGeom &g,
                         const Microphysics::EOS::EOS &eos, const Real x1, const Real x2,
-                        const Real x3) const {
+                        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;
@@ -360,19 +383,21 @@ class ConToPrim {
       rbsq = bdotr * bdotr;
       bsq_rpsq = bsq * rsq - rbsq;
     }
-    const Real zsq = rsq / h0sq_;
-    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,
+    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.
-    // const Real mu_r = res.compute_upper_bound(h0sq_);
-    // solve
+    // 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(); // we don't always need a second root find
+    } else {
+      mu_r = 1 / sqrt(res.get_h0sq());
+    }
+
     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);
@@ -391,7 +416,7 @@ class ConToPrim {
     eos_lambda[1] = std::log10(v(tmp)); // initial guess
     v(peng) = res.ehat_mu(mu, qbar, rbarsq, vsq, W);
     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");

src/fluid/fluid.cpp

@@ -499,6 +499,7 @@ TaskStatus ConservedToPrimitiveRobust(T *rc, const IndexRange &ib, const IndexRa
   auto eos = eos_pkg->Param<Microphysics::EOS::EOS>("d.EOS");
   auto geom = Geometry::GetCoordinateSystem(rc);
   auto coords = pmb->coords;
+  const Real h_min = eos_pkg->Param<Real>("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);

src/microphysics/eos_phoebus/eos_phoebus.cpp

@@ -65,6 +65,9 @@ std::shared_ptr<StateDescriptor> 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 +96,8 @@ std::shared_ptr<StateDescriptor> 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 +159,39 @@ std::shared_ptr<StateDescriptor> 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 = 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 (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 = eos_sc.TMin();
+    rho = eos_sc.rhoMin();
+    ye = ye_min;
+
+    for (int y = 0; y < n / 2; y++) {
+      lambda[0] = ye;
+      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 *= pow(10.0, dT); // log spacing
+        }
+        T = eos_sc.TMin();      // "reset"
+        rho *= pow(10.0, drho); // log spacing
+      }
+      rho = eos_sc.rhoMin(); // "reset"
+      ye += dye;             // linear spacing
+    }
+
 #endif
   } else {
     std::stringstream error_mesg;
@@ -177,6 +215,7 @@ std::shared_ptr<StateDescriptor> 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;
 }