Merge pull request #210 from lanl/eos_floors

Add eos information to floors, misc UserWorkBeforeOutput fix

src/fixup/fixup.cpp

@@ -49,6 +49,7 @@ namespace fixup {
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   auto fix = std::make_shared<StateDescriptor>("fixup");
   Params &params = fix->AllParams();
+  auto mutable_param = parthenon::Params::Mutability::Mutable;
 
   const bool enable_mhd = pin->GetOrAddBoolean("fluid", "mhd", false);
   const bool enable_rad = pin->GetOrAddBoolean("physics", "rad", false);
@@ -273,7 +274,8 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
              Bounds(params.Get<Floors>("floor"), params.Get<Ceilings>("ceiling"),
                     params.Get<MHDCeilings>("mhd_ceiling"),
                     params.Get<RadiationFloors>("rad_floor"),
-                    params.Get<RadiationCeilings>("rad_ceiling")));
+                    params.Get<RadiationCeilings>("rad_ceiling")),
+             mutable_param);
 
   return fix;
 }
@@ -307,6 +309,9 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) {
 
   if (!enable_floors) return TaskStatus::complete;
 
+  const Real ye_min = eos_pkg->Param<Real>("ye_min");
+  const Real ye_max = eos_pkg->Param<Real>("ye_max");
+
   const std::vector<std::string> vars(
       {p::density::name(), c::density::name(), p::velocity::name(), c::momentum::name(),
        p::energy::name(), c::energy::name(), p::bfield::name(), p::ye::name(),
@@ -344,7 +349,13 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) {
 
   auto eos = eos_pkg->Param<EOS>("d.EOS");
   auto geom = Geometry::GetCoordinateSystem(rc);
-  auto bounds = fix_pkg->Param<Bounds>("bounds");
+  Bounds *pbounds = fix_pkg->MutableParam<Bounds>("bounds");
+
+  // BLB: Setting EOS bnds for Ceilings/Floors here.
+  pbounds->SetEOSBnds(eos_pkg);
+
+  // copy bounds by value for kernel
+  Bounds bounds = *pbounds;
 
   const Real c2p_tol = fluid_pkg->Param<Real>("c2p_tol");
   const int c2p_max_iter = fluid_pkg->Param<int>("c2p_max_iter");
@@ -437,16 +448,18 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) {
 
         if (floor_applied) {
           Real vp_normalobs[3] = {0}; // Inject floors at rest in normal observer frame
-          Real ye_prim_default = 0.5;
-          eos_lambda[0] = ye_prim_default;
+          Real ye = 0.5;
+          if (pye > 0) {
+            ye = v(b, cye, k, j, i);
+          }
+          eos_lambda[0] = ye;
           Real dprs =
               eos.PressureFromDensityInternalEnergy(drho, ratio(du, drho), eos_lambda);
           Real dgm1 = ratio(
               eos.BulkModulusFromDensityInternalEnergy(drho, ratio(du, drho), eos_lambda),
               dprs);
-          prim2con::p2c(drho, vp_normalobs, bp, du, ye_prim_default, dprs, dgm1, gcov,
-                        gammacon, betacon, alpha, sdetgam, dcrho, dS, dBcons, dtau,
-                        dyecons);
+          prim2con::p2c(drho, vp_normalobs, bp, du, ye, dprs, dgm1, gcov, gammacon,
+                        betacon, alpha, sdetgam, dcrho, dS, dBcons, dtau, dyecons);
 
           // Update cons vars (not B field)
           v(b, crho, k, j, i) += dcrho;
@@ -464,7 +477,7 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) {
             SPACELOOP(ii) { v(b, pvel_lo + ii, k, j, i) = vp_normalobs[ii]; }
             v(b, peng, k, j, i) = du;
             if (pye > 0) {
-              v(b, pye, k, j, i) = ye_prim_default;
+              v(b, pye, k, j, i) = ye_min;
             }
 
             // Update auxiliary primitives
@@ -476,9 +489,9 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) {
                 eos_lambda);
 
             // Update cons vars (not B field)
-            prim2con::p2c(drho, vp_normalobs, bp, du, ye_prim_default, dprs, dgm1, gcov,
-                          gammacon, betacon, alpha, sdetgam, v(b, crho, k, j, i), dS,
-                          dBcons, v(b, ceng, k, j, i), dyecons);
+            prim2con::p2c(drho, vp_normalobs, bp, du, ye, dprs, dgm1, gcov, gammacon,
+                          betacon, alpha, sdetgam, v(b, crho, k, j, i), dS, dBcons,
+                          v(b, ceng, k, j, i), dyecons);
             SPACELOOP(ii) { v(b, cmom_lo + ii, k, j, i) = dS[ii]; }
             if (pye > 0) {
               v(b, cye, k, j, i) = dyecons;

src/fixup/fixup.hpp

@@ -91,37 +91,60 @@ class Floors {
   KOKKOS_INLINE_FUNCTION
   void GetFloors(const Real x1, const Real x2, const Real x3, Real &rflr,
                  Real &sflr) const {
+    if (!eos_bnds_set_) {
+      PARTHENON_FAIL("EOS bounds not set in floors.");
+    }
+
     switch (floor_flag_) {
     case FloorFlag::ConstantRhoSie:
       rflr = r0_;
       sflr = s0_;
+      rflr = std::max(rflr, rho_min_eos_);
+      sflr = std::max(sflr, sie_min_eos_);
       break;
     case FloorFlag::ExpX1RhoSie:
       rflr = r0_ * std::exp(ralpha_ * x1);
       sflr = s0_ * std::exp(salpha_ * x1);
+      rflr = std::max(rflr, rho_min_eos_);
+      sflr = std::max(sflr, sie_min_eos_);
       break;
     case FloorFlag::ExpX1RhoU: {
       Real scratch = r0_ * std::exp(ralpha_ * x1);
       sflr = s0_ * std::exp(salpha_ * x1) / std::max(rflr, scratch);
       rflr = scratch;
+      rflr = std::max(rflr, rho_min_eos_);
+      sflr = std::max(sflr, sie_min_eos_);
     } break;
     case FloorFlag::X1RhoSie:
       rflr = r0_ * std::min(1.0, std::pow(x1, ralpha_));
       sflr = s0_ * std::min(1.0, std::pow(x1, salpha_));
+      rflr = std::max(rflr, rho_min_eos_);
+      sflr = std::max(sflr, sie_min_eos_);
       break;
     case FloorFlag::RRhoSie: {
       Real r = std::sqrt(x1 * x1 + x2 * x2 + x3 * x3);
       rflr = r0_ * std::min(1.0, std::pow(r, ralpha_));
       sflr = s0_ * std::min(1.0, std::pow(r, salpha_));
+      rflr = std::max(rflr, rho_min_eos_);
+      sflr = std::max(sflr, sie_min_eos_);
     } break;
     default:
       PARTHENON_FAIL("No valid floor set.");
     }
   }
 
+  void SetEOSBnds(StateDescriptor *eos_pkg) {
+    if (!eos_bnds_set_) {
+      rho_min_eos_ = eos_pkg->Param<Real>("rho_min");
+      sie_min_eos_ = eos_pkg->Param<Real>("sie_min");
+      eos_bnds_set_ = true;
+    }
+  }
+
  private:
-  Real r0_, s0_, ralpha_, salpha_;
+  Real r0_, s0_, ralpha_, salpha_, rho_min_eos_, sie_min_eos_;
   const FloorFlag floor_flag_;
+  bool eos_bnds_set_;
 };
 
 static struct ConstantJFloor {
@@ -173,6 +196,7 @@ class Ceilings {
   KOKKOS_INLINE_FUNCTION
   void GetCeilings(const Real x1, const Real x2, const Real x3, Real &gmax,
                    Real &smax) const {
+
     switch (ceiling_flag_) {
     case 1:
       gmax = g0_;
@@ -257,13 +281,18 @@ class Bounds {
          const RadiationFloors &rfl, const RadiationCeilings &rcl)
       : floors_(fl), ceilings_(cl), mhd_ceilings_(mcl), radiation_floors_(rfl),
         radiation_ceilings_(rcl) {}
-  explicit Bounds(const Floors &fl)
+  explicit Bounds(Floors &fl)
       : floors_(fl), ceilings_(Ceilings()), mhd_ceilings_(MHDCeilings()),
         radiation_floors_(RadiationFloors()), radiation_ceilings_(RadiationCeilings()) {}
-  explicit Bounds(const Ceilings &cl)
+  explicit Bounds(Ceilings &cl)
       : floors_(Floors()), ceilings_(cl), mhd_ceilings_(MHDCeilings()),
         radiation_floors_(RadiationFloors()), radiation_ceilings_(RadiationCeilings()) {}
 
+  template <class... Args>
+  void SetEOSBnds(Args &&...args) {
+    floors_.SetEOSBnds(std::forward<Args>(args)...);
+  }
+
   template <class... Args>
   KOKKOS_INLINE_FUNCTION void GetFloors(Args &&...args) const {
     floors_.GetFloors(std::forward<Args>(args)...);
@@ -290,8 +319,8 @@ class Bounds {
   }
 
  private:
-  const Floors floors_;
-  const Ceilings ceilings_;
+  Floors floors_;
+  Ceilings ceilings_;
   const MHDCeilings mhd_ceilings_;
   const RadiationFloors radiation_floors_;
   const RadiationCeilings radiation_ceilings_;

src/fixup/fixup_c2p.cpp

@@ -155,7 +155,8 @@ TaskStatus ConservedToPrimitiveFixupImpl(T *rc) {
 
   auto eos = eos_pkg->Param<Microphysics::EOS::EOS>("d.EOS");
   auto geom = Geometry::GetCoordinateSystem(rc);
-  auto bounds = fix_pkg->Param<Bounds>("bounds");
+  Bounds *pbounds = fix_pkg->MutableParam<Bounds>("bounds");
+  Bounds bounds = *pbounds;
 
   Coordinates_t coords = rc->GetParentPointer()->coords;
 

src/fixup/fixup_radc2p.cpp

@@ -113,7 +113,7 @@ TaskStatus RadConservedToPrimitiveFixupImpl(T *rc) {
   }
 
   auto geom = Geometry::GetCoordinateSystem(rc);
-  auto bounds = fix_pkg->Param<Bounds>("bounds");
+  Bounds *bounds = fix_pkg->MutableParam<Bounds>("bounds");
 
   Coordinates_t coords = rc->GetParentPointer()->coords;
 
@@ -129,8 +129,8 @@ TaskStatus RadConservedToPrimitiveFixupImpl(T *rc) {
       KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
         Real xi_max;
         Real garbage;
-        bounds.GetRadiationCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
-                                    coords.Xc<3>(k, j, i), xi_max, garbage);
+        bounds->GetRadiationCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
+                                     coords.Xc<3>(k, j, i), xi_max, garbage);
 
         // It is assumed that the fluid is already fixed up
         auto fail = [&](const int k, const int j, const int i) {

src/fixup/fixup_src.cpp

@@ -73,7 +73,7 @@ TaskStatus SourceFixupImpl(T *rc) {
   }
 
   auto eos = eos_pkg->Param<EOS>("d.EOS");
-  auto bounds = fix_pkg->Param<Bounds>("bounds");
+  Bounds *bounds = fix_pkg->MutableParam<Bounds>("bounds");
 
   const std::vector<std::string> vars(
       {p::density::name(), c::density::name(), p::velocity::name(), c::momentum::name(),
@@ -151,12 +151,12 @@ TaskStatus SourceFixupImpl(T *rc) {
       kb.e, jb.s, jb.e, ib.s, ib.e,
       KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
         double gamma_max, e_max;
-        bounds.GetCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
-                           coords.Xc<3>(k, j, i), gamma_max, e_max);
+        bounds->GetCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
+                            coords.Xc<3>(k, j, i), gamma_max, e_max);
         Real xi_max;
         Real garbage;
-        bounds.GetRadiationCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
-                                    coords.Xc<3>(k, j, i), xi_max, garbage);
+        bounds->GetRadiationCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
+                                     coords.Xc<3>(k, j, i), xi_max, garbage);
 
         double eos_lambda[2]; // used for stellarcollapse eos and other EOS's that require
                               // root finding.

src/fluid/con2prim_robust.hpp

@@ -164,7 +164,7 @@ class Residual {
  private:
   const Real D_, q_, bsq_, bsq_rpsq_, rsq_, rbsq_, v0sq_;
   const Microphysics::EOS::EOS &eos_;
-  const fixup::Bounds &bounds_;
+  const fixup::Bounds bounds_;
   const Real x1_, x2_, x3_;
   const Real floor_scale_fac_;
   Real lambda_[2];

src/fluid/fluid.cpp

@@ -476,7 +476,10 @@ TaskStatus ConservedToPrimitiveRobust(T *rc, const IndexRange &ib, const IndexRa
   auto *pmb = rc->GetParentPointer();
 
   StateDescriptor *fix_pkg = pmb->packages.Get("fixup").get();
-  auto bounds = fix_pkg->Param<fixup::Bounds>("bounds");
+  StateDescriptor *eos_pkg = pmb->packages.Get("eos").get();
+  fixup::Bounds *pbounds = fix_pkg->MutableParam<fixup::Bounds>("bounds");
+  pbounds->SetEOSBnds(eos_pkg);
+  fixup::Bounds bounds = *pbounds;
 
   StateDescriptor *pkg = pmb->packages.Get("fluid").get();
   const Real c2p_tol = pkg->Param<Real>("c2p_tol");
@@ -488,7 +491,6 @@ TaskStatus ConservedToPrimitiveRobust(T *rc, const IndexRange &ib, const IndexRa
                                                 c2p_floor_scale_fac, c2p_fail_on_floors,
                                                 c2p_fail_on_ceilings);
 
-  StateDescriptor *eos_pkg = pmb->packages.Get("eos").get();
   auto eos = eos_pkg->Param<Microphysics::EOS::EOS>("d.EOS");
   auto geom = Geometry::GetCoordinateSystem(rc);
   auto coords = pmb->coords;
@@ -516,7 +518,7 @@ TaskStatus ConservedToPrimitiveClassic(T *rc, const IndexRange &ib, const IndexR
   auto *pmesh = rc->GetMeshPointer();
 
   StateDescriptor *fix_pkg = pmesh->packages.Get("fixup").get();
-  auto bounds = fix_pkg->Param<fixup::Bounds>("bounds");
+  fixup::Bounds *pbounds = fix_pkg->MutableParam<fixup::Bounds>("bounds");
 
   StateDescriptor *pkg = pmesh->packages.Get("fluid").get();
   const Real c2p_tol = pkg->Param<Real>("c2p_tol");

src/fluid/riemann.hpp

@@ -205,6 +205,7 @@ class FluxState {
   const ParArrayND<Real> qr;
   const Geometry::CoordSysMeshBlock geom;
   const Coordinates_t coords;
+  fixup::Bounds *pbounds;
   fixup::Bounds bounds;
 
  private:
@@ -216,9 +217,11 @@ class FluxState {
         ql(rc->Get("ql").data), qr(rc->Get("qr").data),
         geom(Geometry::GetCoordinateSystem(rc)),
         coords(rc->GetParentPointer()->coords), // problem for packs
-        bounds(rc->GetParentPointer()->packages.Get("fixup").get()->Param<fixup::Bounds>(
-            "bounds")),
-        prho(imap[fluid_prim::density::name()].first),
+        pbounds(rc->GetParentPointer()
+                    ->packages.Get("fixup")
+                    .get()
+                    ->MutableParam<fixup::Bounds>("bounds")),
+        bounds(*pbounds), prho(imap[fluid_prim::density::name()].first),
         pvel_lo(imap[fluid_prim::velocity::name()].first),
         peng(imap[fluid_prim::energy::name()].first),
         pb_lo(imap[fluid_prim::bfield::name()].first),

src/microphysics/eos_phoebus/eos_phoebus.cpp

@@ -63,6 +63,8 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   Real rho_max;
   Real sie_max;
   Real T_max;
+  Real ye_min;
+  Real ye_max;
 
   std::string eos_type = pin->GetString(block_name, std::string("type"));
   params.Add("type", eos_type);
@@ -82,13 +84,16 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     params.Add("d.EOS", eos_device);
     params.Add("h.EOS", eos_host);
 
-    rho_min = pin->GetOrAddReal("fixup", "rho0_floor", 0.0);
-    sie_min = pin->GetOrAddReal("fixup", "sie0_floor", 0.0);
+    // Can specify rho_min, etc, in <eos>
+    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("fixup", "rho0_ceiling", 1e18);
-    sie_max = pin->GetOrAddReal("fixup", "sie0_ceiling", 1e35);
+    rho_max = pin->GetOrAddReal("eos", "rho_max", 1e18);
+    sie_max = pin->GetOrAddReal("eos", "sie_max", 1e35);
     T_max = eos_host.TemperatureFromDensityInternalEnergy(rho_max, sie_max, lambda);
+    ye_min = 0.01;
+    ye_max = 1.0;
 #ifdef SPINER_USE_HDF
   } else if (eos_type == StellarCollapse::EosType()) {
     // We request that Ye and temperature exist, but do not provide them.
@@ -148,6 +153,8 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     T_max = eos_sc.TMax() / T_unit;
     rho_min = eos_sc.rhoMin() / rho_unit;
     rho_max = eos_sc.rhoMax() / rho_unit;
+    ye_min = eos_sc.YeMin();
+    ye_max = eos_sc.YeMax();
 #endif
   } else {
     std::stringstream error_mesg;
@@ -169,6 +176,8 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   params.Add("T_max", T_max);
   params.Add("rho_min", rho_min);
   params.Add("rho_max", rho_max);
+  params.Add("ye_min", ye_min);
+  params.Add("ye_max", ye_max);
 
   return pkg;
 }

src/pgen/torus.cpp

@@ -198,9 +198,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   const Real &minx_k = pmb->coords.Xf<3>(kb.s);
 
   auto coords = pmb->coords;
-  auto eos = pmb->packages.Get("eos")->Param<EOS>("d.EOS");
+  StateDescriptor *eos_pkg = pmb->packages.Get("eos").get();
+  auto eos = eos_pkg->Param<EOS>("d.EOS");
   auto eos_h = pmb->packages.Get("eos")->Param<EOS>("h.EOS");
   auto floor = pmb->packages.Get("fixup")->Param<fixup::Floors>("floor");
+  floor.SetEOSBnds(eos_pkg);
   auto &unit_conv =
       pmb->packages.Get("phoebus")->Param<phoebus::UnitConversions>("unit_conv");
   S *= unit_conv.GetEntropyCGSToCode();

src/phoebus_driver.cpp

@@ -1049,6 +1049,8 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) {
   IndexRange kb = rc->GetBoundsK(IndexDomain::interior);
 
   auto eos = pmb->packages.Get("eos")->Param<Microphysics::EOS::EOS>("d.EOS");
+  auto analysis = pmb->packages.Get("analysis").get();
+  const Real sigma = analysis->Param<Real>("sigma");
   parthenon::par_for(
       DEFAULT_LOOP_PATTERN, "UserWorkBeforeOutput::H5", DevExecSpace(), kb.s, kb.e, jb.s,
       jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) {
@@ -1071,9 +1073,7 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) {
         Real gam1[3][3];
         Real gam2[3][3];
         Real gam3[3][3];
-        auto analysis = pmb->packages.Get("analysis").get();
         const Real z = coords.Xc<3>(k, j, i);
-        const Real sigma = analysis->Param<Real>("sigma");
         const Real pi = 3.14;
         const Real s0 =
             s * std::exp(-z * z / sigma / sigma) / std::sqrt(pi) / sigma; // sigma > 0

src/radiation/moments_source.cpp

@@ -269,7 +269,9 @@ TaskStatus MomentFluidSourceImpl(T *rc, Real dt, bool update_fluid) {
     species_d[s] = species[s];
   }
 
-  auto bounds = fix_pkg->Param<fixup::Bounds>("bounds");
+  Bounds *pbounds = fix_pkg->MutableParam<fixup::Bounds>("bounds");
+  pbounds->SetEOSBnds(eos_pkg);
+  Bounds bounds = *pbounds;
 
   const parthenon::Coordinates_t &coords = pmb->coords;