Velocity and temperature ceilings

parthenon-hpc-lab · Jul 17, 2023 · e695577 · e695577
1 parent 08ab35c
commit e695577
Show file tree

Hide file tree

Showing 4 changed files with 95 additions and 11 deletions.
diff --git a/src/eos/adiabatic_glmmhd.hpp b/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;
+      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;

diff --git a/src/eos/adiabatic_hydro.hpp b/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;

diff --git a/src/eos/eos.hpp b/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_
diff --git a/src/hydro/hydro.cpp b/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, gamma, vceil, eceil);
       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, gamma, vceil, eceil);
       pkg->AddParam<>("eos", eos);
       pkg->FillDerivedMesh = ConsToPrim<AdiabaticGLMMHDEOS>;
       pkg->EstimateTimestepMesh = EstimateTimestep<Fluid::glmmhd>;