WIP: Limit velocity and temperature in Kinetic jet injection region

src/eos/adiabatic_glmmhd.hpp

@@ -23,9 +23,11 @@ using parthenon::Real;
 
 class AdiabaticGLMMHDEOS : public EquationOfState {
  public:
-  AdiabaticGLMMHDEOS(Real pressure_floor, Real density_floor, Real internal_e_floor,
-                     Real gamma)
-      : EquationOfState(pressure_floor, density_floor, internal_e_floor), gamma_{gamma} {}
+   AdiabaticGLMMHDEOS(Real pressure_floor, Real density_floor,
+       Real internal_e_floor, Real velocity_ceiling, Real internal_e_ceiling,
+       Real gamma)
+     : EquationOfState(pressure_floor, density_floor, internal_e_floor,
+       velocity_ceiling, internal_e_ceiling), gamma_{gamma} {}
 
   void ConservedToPrimitive(MeshData<Real> *md) const override;
 
@@ -66,6 +68,9 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
     auto pressure_floor_ = GetPressureFloor();
     auto e_floor_ = GetInternalEFloor();
 
+    auto velocity_ceiling_ = GetVelocityCeiling();
+    auto e_ceiling_ = GetInternalECeiling();
+
     Real &u_d = cons(IDN, k, j, i);
     Real &u_m1 = cons(IM1, k, j, i);
     Real &u_m2 = cons(IM2, k, j, i);
@@ -109,6 +114,23 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
     Real e_B = 0.5 * (SQR(u_b1) + SQR(u_b2) + SQR(u_b3));
     w_p = gm1 * (u_e - e_k - e_B);
 
+    // apply velocity ceiling. By default ceiling is std::numeric_limits<Real>::infinity()
+    const Real w_v2 = SQR(w_vx) + SQR(w_vy) + SQR(w_vz);
+    if( w_v2 > SQR(velocity_ceiling_) ){
+      const Real w_v = sqrt(w_v2);
+      w_vx *= velocity_ceiling_/w_v;
+      w_vy *= velocity_ceiling_/w_v;
+      w_vz *= velocity_ceiling_/w_v;
+
+      u_m1 *= velocity_ceiling_/w_v;
+      u_m2 *= velocity_ceiling_/w_v;
+      u_m3 *= velocity_ceiling_/w_v;
+
+      Real e_k_new = 0.5 * u_d * SQR(velocity_ceiling_);
+      u_e -= e_k - e_k_new;
+      e_k = e_k_new;
+    }
+
     // Let's apply floors explicitly, i.e., by default floor will be disabled (<=0)
     // and the code will fail if a negative pressure is encountered.
     PARTHENON_REQUIRE(w_p > 0.0 || pressure_floor_ > 0.0 || e_floor_ > 0.0,
@@ -130,6 +152,14 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
       w_p = eff_pressure_floor;
     }
 
+    // temperature (internal energy) based pressure ceiling
+    const Real eff_pressure_ceiling = gm1 * u_d * e_ceiling_;
+    if (w_p > eff_pressure_ceiling) {
+      // apply temperature ceiling, correct total energy
+      u_e = (u_d * e_ceiling_) + e_k + e_B;
+      w_p = eff_pressure_ceiling;
+    }
+
     // Convert passive scalars
     for (auto n = nhydro; n < nhydro + nscalars; ++n) {
       prim(n, k, j, i) = cons(n, k, j, i) * di;

src/eos/adiabatic_hydro.hpp

@@ -24,9 +24,11 @@ using parthenon::Real;
 
 class AdiabaticHydroEOS : public EquationOfState {
  public:
-  AdiabaticHydroEOS(Real pressure_floor, Real density_floor, Real internal_e_floor,
-                    Real gamma)
-      : EquationOfState(pressure_floor, density_floor, internal_e_floor), gamma_{gamma} {}
+   AdiabaticHydroEOS(Real pressure_floor, Real density_floor, 
+       Real internal_e_floor, Real velocity_ceiling, Real internal_e_ceiling,
+       Real gamma)
+     : EquationOfState(pressure_floor, density_floor, internal_e_floor,
+       velocity_ceiling, internal_e_ceiling), gamma_{gamma} {}
 
   void ConservedToPrimitive(MeshData<Real> *md) const override;
 
@@ -55,6 +57,9 @@ class AdiabaticHydroEOS : public EquationOfState {
     auto pressure_floor_ = GetPressureFloor();
     auto e_floor_ = GetInternalEFloor();
 
+    auto velocity_ceiling_ = GetVelocityCeiling();
+    auto e_ceiling_ = GetInternalECeiling();
+
     Real &u_d = cons(IDN, k, j, i);
     Real &u_m1 = cons(IM1, k, j, i);
     Real &u_m2 = cons(IM2, k, j, i);
@@ -84,6 +89,23 @@ class AdiabaticHydroEOS : public EquationOfState {
     Real e_k = 0.5 * di * (SQR(u_m1) + SQR(u_m2) + SQR(u_m3));
     w_p = gm1 * (u_e - e_k);
 
+    // apply velocity ceiling. By default ceiling is std::numeric_limits<Real>::infinity()
+    const Real w_v2 = SQR(w_vx) + SQR(w_vy) + SQR(w_vz);
+    if( w_v2 > SQR(velocity_ceiling_) ){
+      const Real w_v = sqrt(w_v2);
+      w_vx *= velocity_ceiling_/w_v;
+      w_vy *= velocity_ceiling_/w_v;
+      w_vz *= velocity_ceiling_/w_v;
+
+      u_m1 *= velocity_ceiling_/w_v;
+      u_m2 *= velocity_ceiling_/w_v;
+      u_m3 *= velocity_ceiling_/w_v;
+
+      Real e_k_new = 0.5 * u_d * SQR(velocity_ceiling_);
+      u_e -= e_k - e_k_new;
+      e_k = e_k_new;
+    }
+
     // Let's apply floors explicitly, i.e., by default floor will be disabled (<=0)
     // and the code will fail if a negative pressure is encountered.
     PARTHENON_REQUIRE(w_p > 0.0 || pressure_floor_ > 0.0 || e_floor_ > 0.0,
@@ -105,6 +127,14 @@ class AdiabaticHydroEOS : public EquationOfState {
       w_p = eff_pressure_floor;
     }
 
+    // temperature (internal energy) based pressure ceiling
+    const Real eff_pressure_ceiling = gm1 * u_d * e_ceiling_;
+    if (w_p > eff_pressure_ceiling) {
+      // apply temperature ceiling, correct total energy
+      u_e = (u_d * e_ceiling_) + e_k;
+      w_p = eff_pressure_ceiling;
+    }
+
     // Convert passive scalars
     for (auto n = nhydro; n < nhydro + nscalars; ++n) {
       prim(n, k, j, i) = cons(n, k, j, i) * di;

src/eos/eos.hpp

@@ -32,9 +32,11 @@ using parthenon::Real;
 
 class EquationOfState {
  public:
-  EquationOfState(Real pressure_floor, Real density_floor, Real internal_e_floor)
-      : pressure_floor_(pressure_floor), density_floor_(density_floor),
-        internal_e_floor_(internal_e_floor) {}
+   EquationOfState(Real pressure_floor, Real density_floor, Real
+       internal_e_floor, Real velocity_ceiling, Real internal_e_ceiling)
+     : pressure_floor_(pressure_floor), density_floor_(density_floor),
+     internal_e_floor_(internal_e_floor), velocity_ceiling_(velocity_ceiling),
+     internal_e_ceiling_(internal_e_ceiling) {}
   virtual void ConservedToPrimitive(MeshData<Real> *md) const = 0;
 
   KOKKOS_INLINE_FUNCTION
@@ -47,8 +49,15 @@ class EquationOfState {
   KOKKOS_INLINE_FUNCTION
   Real GetInternalEFloor() const { return internal_e_floor_; }
 
+  KOKKOS_INLINE_FUNCTION
+  Real GetVelocityCeiling() const { return velocity_ceiling_; }
+
+  KOKKOS_INLINE_FUNCTION
+  Real GetInternalECeiling() const { return internal_e_ceiling_; }
+
  private:
   Real pressure_floor_, density_floor_, internal_e_floor_;
+  Real velocity_ceiling_, internal_e_ceiling_;
 };
 
 #endif // EOS_EOS_HPP_

src/hydro/hydro.cpp

@@ -439,6 +439,21 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
       efloor = Tfloor / mbar_over_kb / (gamma - 1.0);
     }
 
+    // By default disable ceilings by setting to infinity
+    Real vceil = pin->GetOrAddReal("hydro", "vceil", std::numeric_limits<Real>::infinity());
+    Real Tceil = pin->GetOrAddReal("hydro", "Tceil", std::numeric_limits<Real>::infinity());
+    Real eceil = Tceil;
+    if (eceil < std::numeric_limits<Real>::infinity()) {
+      if (!pkg->AllParams().hasKey("mbar_over_kb")) {
+        PARTHENON_FAIL("Temperature ceiling requires units and gas composition. "
+                       "Either set a 'units' block and the 'hydro/He_mass_fraction' in "
+                       "input file or use a pressure floor "
+                       "(defined code units) instead.");
+      }
+      auto mbar_over_kb = pkg->Param<Real>("mbar_over_kb");
+      eceil = Tceil / mbar_over_kb / (gamma - 1.0);
+    }
+
     auto conduction = Conduction::none;
     auto conduction_str = pin->GetOrAddString("diffusion", "conduction", "none");
     if (conduction_str == "spitzer") {
@@ -472,12 +487,12 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     pkg->AddParam<>("conduction", conduction);
 
     if (fluid == Fluid::euler) {
-      AdiabaticHydroEOS eos(pfloor, dfloor, efloor, gamma);
+      AdiabaticHydroEOS eos(pfloor, dfloor, efloor, vceil, eceil, gamma);
       pkg->AddParam<>("eos", eos);
       pkg->FillDerivedMesh = ConsToPrim<AdiabaticHydroEOS>;
       pkg->EstimateTimestepMesh = EstimateTimestep<Fluid::euler>;
     } else if (fluid == Fluid::glmmhd) {
-      AdiabaticGLMMHDEOS eos(pfloor, dfloor, efloor, gamma);
+      AdiabaticGLMMHDEOS eos(pfloor, dfloor, efloor, vceil, eceil, gamma);
       pkg->AddParam<>("eos", eos);
       pkg->FillDerivedMesh = ConsToPrim<AdiabaticGLMMHDEOS>;
       pkg->EstimateTimestepMesh = EstimateTimestep<Fluid::glmmhd>;

src/pgen/cluster.cpp

@@ -32,6 +32,8 @@
 #include <parthenon/package.hpp>
 
 // AthenaPK headers
+#include "../eos/adiabatic_glmmhd.hpp"
+#include "../eos/adiabatic_hydro.hpp"
 #include "../hydro/hydro.hpp"
 #include "../hydro/srcterms/gravitational_field.hpp"
 #include "../hydro/srcterms/tabular_cooling.hpp"
@@ -54,6 +56,122 @@ using namespace parthenon::driver::prelude;
 using namespace parthenon::package::prelude;
 using utils::few_modes_ft::FewModesFT;
 
+template<class EOS>
+void ApplyClusterClips(MeshData<Real> *md, const parthenon::SimTime &tm,
+                    const Real beta_dt, const EOS eos) {
+
+  auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
+
+  //Apply clips -- ceilings on temperature, velocity, alfven velocity, and
+  //density floor -- within a radius of the AGN
+  const auto &dfloor = hydro_pkg->Param<Real>("cluster_dfloor");
+  const auto &eceil = hydro_pkg->Param<Real>("cluster_eceil");
+  const auto &vceil = hydro_pkg->Param<Real>("cluster_vceil");
+  const auto &vAceil = hydro_pkg->Param<Real>("cluster_vAceil");
+  const auto &clip_r = hydro_pkg->Param<Real>("cluster_clip_r");
+
+  if( clip_r > 0 && (dfloor > 0 || 
+                     eceil  < std::numeric_limits<Real>::infinity() ||
+                     vceil  < std::numeric_limits<Real>::infinity() ||
+                     vAceil < std::numeric_limits<Real>::infinity()) ){
+    // Grab some necessary variables
+    const auto &prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
+    const auto &cons_pack = md->PackVariables(std::vector<std::string>{"cons"});
+    IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior);
+    IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior);
+    IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior);
+    const auto nhydro = hydro_pkg->Param<int>("nhydro");
+    const auto nscalars = hydro_pkg->Param<int>("nscalars");
+
+    const Real clip_r2 = SQR(clip_r);
+    const Real vceil2 = SQR(vceil);
+    const Real vAceil2 = SQR(vAceil);
+    const Real gm1 = (hydro_pkg->Param<Real>("AdiabaticIndex") - 1.0);
+
+    parthenon::par_for(
+        DEFAULT_LOOP_PATTERN, "Cluster::ApplyClusterClips",
+        parthenon::DevExecSpace(), 0, cons_pack.GetDim(5) - 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 &cons = cons_pack(b);
+          auto &prim = prim_pack(b);
+          const auto &coords = cons_pack.GetCoords(b);
+
+          const Real r2 = SQR(coords.Xc<1>(i)) + SQR(coords.Xc<2>(j)) + SQR(coords.Xc<3>(k));
+
+          if ( r2 < clip_r2 ) {
+            //Cell falls within clipping radius
+            eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
+
+            if( dfloor > 0 ){
+              const Real rho = prim(IDN, k, j, i);
+              if( rho < dfloor){
+                cons(IDN, k, j, i) = dfloor;
+                prim(IDN, k, j, i) = dfloor;
+              }
+            }
+
+            if( vceil < std::numeric_limits<Real>::infinity() ){
+              // Apply velocity ceiling
+              const Real v2 =
+                  SQR(prim(IV1, k, j, i)) + SQR(prim(IV2, k, j, i)) + SQR(prim(IV3, k, j, i));
+              if (v2 > vceil2) {
+                // Fix the velocity to the velocity ceiling
+                const Real v = sqrt(v2);
+                cons(IM1, k, j, i) *= vceil / v;
+                cons(IM2, k, j, i) *= vceil / v;
+                cons(IM3, k, j, i) *= vceil / v;
+                prim(IV1, k, j, i) *= vceil / v;
+                prim(IV2, k, j, i) *= vceil / v;
+                prim(IV3, k, j, i) *= vceil / v;
+
+                // Update the internal energy
+                cons(IEN, k, j, i) -= 0.5 * prim(IDN, k, j, i) * (v2 - vceil2);
+              }
+            }
+
+            if( vAceil2 < std::numeric_limits<Real>::infinity() ){
+              // Apply Alfven velocity ceiling by raising density
+              const Real rho = prim(IDN, k, j, i);
+              const Real B2 = (SQR(prim(IB1, k, j, i)) + SQR(prim(IB2, k, j, i)) + SQR(prim(IB3, k, j, i)));
+
+              // compute Alfven mach number
+              const Real va2 = (B2 / rho);
+
+              if (va2 > vAceil2) {
+                //Increase the density to match the alfven velocity ceiling
+                const Real rho_new = std::sqrt(vAceil2/B2);
+                cons(IDN, k, j, i) = rho_new;
+                prim(IDN, k, j, i) = rho_new;
+              }
+            }
+
+            if( eceil < std::numeric_limits<Real>::infinity() ){
+              // Apply  internal energy ceiling as a pressure ceiling
+              const Real internal_e = prim(IPR, k, j, i) / (gm1 * prim(IDN, k, j, i));
+              if (internal_e > eceil) {
+                cons(IEN, k, j, i) -= prim(IDN, k, j, i) * (internal_e - eceil);
+                prim(IPR, k, j, i) = gm1 * prim(IDN, k, j, i) * eceil;
+              }
+            }
+          }
+        });
+  }    
+}
+
+void ApplyClusterClips(MeshData<Real> *md, const parthenon::SimTime &tm,
+                    const Real beta_dt) {
+  auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
+  auto fluid = hydro_pkg->Param<Fluid>("fluid");
+  if (fluid == Fluid::euler) {
+    ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param<AdiabaticHydroEOS>("eos"));
+  } else if (fluid == Fluid::glmmhd) {
+    ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param<AdiabaticGLMMHDEOS>("eos"));
+  } else {
+    PARTHENON_FAIL("Cluster::ApplyClusterClips: Unknown EOS");
+  }
+
+}
+
 void ClusterSrcTerm(MeshData<Real> *md, const parthenon::SimTime &tm,
                     const Real beta_dt) {
   auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
@@ -75,7 +193,9 @@ void ClusterSrcTerm(MeshData<Real> *md, const parthenon::SimTime &tm,
 
   const auto &snia_feedback = hydro_pkg->Param<SNIAFeedback>("snia_feedback");
   snia_feedback.FeedbackSrcTerm(md, beta_dt, tm);
-}
+
+  ApplyClusterClips(md, tm, beta_dt);
+};
 
 Real ClusterEstimateTimestep(MeshData<Real> *md) {
   Real min_dt = std::numeric_limits<Real>::max();
@@ -221,6 +341,38 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd
 
   SNIAFeedback snia_feedback(pin, hydro_pkg);
 
+
+  /************************************************************
+   * Read Clips  (ceilings and floors)
+   ************************************************************/
+
+  //Disable all clips by default with a negative radius clip
+  Real clip_r = pin->GetOrAddReal("cluster/clips", "clip_r", -1.0);
+
+  // By default disable floors by setting a negative value
+  Real dfloor = pin->GetOrAddReal("cluster/clips", "dfloor", -1.0);
+
+  // By default disable ceilings by setting to infinity
+  Real vceil = pin->GetOrAddReal("cluster/clips", "vceil", std::numeric_limits<Real>::infinity());
+  Real vAceil = pin->GetOrAddReal("cluster/clips", "vAceil", std::numeric_limits<Real>::infinity());
+  Real Tceil = pin->GetOrAddReal("cluster/clips", "Tceil", std::numeric_limits<Real>::infinity());
+  Real eceil = Tceil;
+  if (eceil < std::numeric_limits<Real>::infinity()) {
+    if (!hydro_pkg->AllParams().hasKey("mbar_over_kb")) {
+      PARTHENON_FAIL("Temperature ceiling requires units and gas composition. "
+                     "Either set a 'units' block and the 'hydro/He_mass_fraction' in "
+                     "input file or use a pressure floor "
+                     "(defined code units) instead.");
+    }
+    auto mbar_over_kb = hydro_pkg->Param<Real>("mbar_over_kb");
+    eceil = Tceil / mbar_over_kb / (hydro_pkg->Param<Real>("adiabatic_index") - 1.0);
+  }
+  hydro_pkg->AddParam("cluster_dfloor",dfloor);
+  hydro_pkg->AddParam("cluster_eceil",eceil);
+  hydro_pkg->AddParam("cluster_vceil",vceil);
+  hydro_pkg->AddParam("cluster_vAceil",vAceil);
+  hydro_pkg->AddParam("cluster_clip_r",clip_r);
+
   /************************************************************
    * Add derived fields
    * NOTE: these must be filled in UserWorkBeforeOutput
@@ -686,7 +838,7 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) {
   // for computing temperature from primitives
   auto units = pkg->Param<Units>("units");
   auto mbar_over_kb = pkg->Param<Real>("mbar_over_kb");
-  auto mbar = mbar_over_kb*units.k_boltzmann();
+  auto mbar = mbar_over_kb * units.k_boltzmann();
 
   // fill derived vars (*including ghost cells*)
   auto &coords = pmb->coords;
@@ -712,7 +864,7 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) {
         log10_radius(k, j, i) = 0.5 * std::log10(r2);
 
         // compute entropy
-        const Real K = P / std::pow(rho/mbar, gam);
+        const Real K = P / std::pow(rho / mbar, gam);
         entropy(k, j, i) = K;
 
         const Real v_mag = std::sqrt(SQR(v1) + SQR(v2) + SQR(v3));