diff --git a/src/eos/adiabatic_glmmhd.hpp b/src/eos/adiabatic_glmmhd.hpp index 9fa80b76..baebc681 100644 --- a/src/eos/adiabatic_glmmhd.hpp +++ b/src/eos/adiabatic_glmmhd.hpp @@ -24,8 +24,10 @@ 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} {} + 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 *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::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; diff --git a/src/eos/adiabatic_hydro.hpp b/src/eos/adiabatic_hydro.hpp index 4dca9ce2..bf317923 100644 --- a/src/eos/adiabatic_hydro.hpp +++ b/src/eos/adiabatic_hydro.hpp @@ -25,8 +25,10 @@ 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} {} + 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 *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::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 index b6986d7e..b75eedce 100644 --- 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) + 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) {} + internal_e_floor_(internal_e_floor), velocity_ceiling_(velocity_ceiling), + internal_e_ceiling_(internal_e_ceiling) {} virtual void ConservedToPrimitive(MeshData *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 index 40cdd27b..3808bde6 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -439,6 +439,23 @@ std::shared_ptr 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::infinity()); + Real Tceil = + pin->GetOrAddReal("hydro", "Tceil", std::numeric_limits::infinity()); + Real eceil = Tceil; + if (eceil < std::numeric_limits::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("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 +489,12 @@ std::shared_ptr 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; pkg->EstimateTimestepMesh = EstimateTimestep; } 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; pkg->EstimateTimestepMesh = EstimateTimestep; diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index 5452e5d1..ab112be4 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -32,6 +32,8 @@ #include // 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 +void ApplyClusterClips(MeshData *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("cluster_dfloor"); + const auto &eceil = hydro_pkg->Param("cluster_eceil"); + const auto &vceil = hydro_pkg->Param("cluster_vceil"); + const auto &vAceil = hydro_pkg->Param("cluster_vAceil"); + const auto &clip_r = hydro_pkg->Param("cluster_clip_r"); + + if (clip_r > 0 && (dfloor > 0 || eceil < std::numeric_limits::infinity() || + vceil < std::numeric_limits::infinity() || + vAceil < std::numeric_limits::infinity())) { + // Grab some necessary variables + const auto &prim_pack = md->PackVariables(std::vector{"prim"}); + const auto &cons_pack = md->PackVariables(std::vector{"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("nhydro"); + const auto nscalars = hydro_pkg->Param("nscalars"); + + const Real clip_r2 = SQR(clip_r); + const Real vceil2 = SQR(vceil); + const Real vAceil2 = SQR(vAceil); + const Real gm1 = (hydro_pkg->Param("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::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; + + // Remove kinetic energy + cons(IEN, k, j, i) -= 0.5 * prim(IDN, k, j, i) * (v2 - vceil2); + } + } + + if (vAceil2 < std::numeric_limits::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(B2 / vAceil2); + cons(IDN, k, j, i) = rho_new; + prim(IDN, k, j, i) = rho_new; + } + } + + if (eceil < std::numeric_limits::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 *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"); + if (fluid == Fluid::euler) { + ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param("eos")); + } else if (fluid == Fluid::glmmhd) { + ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param("eos")); + } else { + PARTHENON_FAIL("Cluster::ApplyClusterClips: Unknown EOS"); + } +} + void ClusterSrcTerm(MeshData *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 *md, const parthenon::SimTime &tm, const auto &snia_feedback = hydro_pkg->Param("snia_feedback"); snia_feedback.FeedbackSrcTerm(md, beta_dt, tm); -} + + ApplyClusterClips(md, tm, beta_dt); +}; Real ClusterEstimateTimestep(MeshData *md) { Real min_dt = std::numeric_limits::max(); @@ -221,6 +341,40 @@ 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("problem/cluster/clips", "clip_r", -1.0); + + // By default disable floors by setting a negative value + Real dfloor = pin->GetOrAddReal("problem/cluster/clips", "dfloor", -1.0); + + // By default disable ceilings by setting to infinity + Real vceil = pin->GetOrAddReal("problem/cluster/clips", "vceil", + std::numeric_limits::infinity()); + Real vAceil = pin->GetOrAddReal("problem/cluster/clips", "vAceil", + std::numeric_limits::infinity()); + Real Tceil = pin->GetOrAddReal("problem/cluster/clips", "Tceil", + std::numeric_limits::infinity()); + Real eceil = Tceil; + if (eceil < std::numeric_limits::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("mbar_over_kb"); + eceil = Tceil / mbar_over_kb / (hydro_pkg->Param("AdiabaticIndex") - 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 +840,7 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { // for computing temperature from primitives auto units = pkg->Param("units"); auto mbar_over_kb = pkg->Param("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 +866,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)); diff --git a/src/pgen/cluster/agn_feedback.cpp b/src/pgen/cluster/agn_feedback.cpp index 16bb647c..1c76e913 100644 --- a/src/pgen/cluster/agn_feedback.cpp +++ b/src/pgen/cluster/agn_feedback.cpp @@ -17,7 +17,7 @@ #include #include -// Athena headers +// AthenaPK headers #include "../../eos/adiabatic_glmmhd.hpp" #include "../../eos/adiabatic_hydro.hpp" #include "../../main.hpp" @@ -33,6 +33,8 @@ using namespace parthenon; AGNFeedback::AGNFeedback(parthenon::ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg) : fixed_power_(pin->GetOrAddReal("problem/cluster/agn_feedback", "fixed_power", 0.0)), + vceil_(pin->GetOrAddReal("problem/cluster/agn_feedback", "vceil", + std::numeric_limits::infinity())), efficiency_(pin->GetOrAddReal("problem/cluster/agn_feedback", "efficiency", 1e-3)), thermal_fraction_( pin->GetOrAddReal("problem/cluster/agn_feedback", "thermal_fraction", 0.0)), @@ -119,7 +121,8 @@ AGNFeedback::AGNFeedback(parthenon::ParameterInput *pin, (1 - efficiency_) * kinetic_jet_e_))) < 10 * std::numeric_limits::epsilon(), "Specified kinetic jet velocity and temperature are incompatible with mass to " - "energy conversion efficiency. Either the specified velocity, temperature, or efficiency are incompatible"); + "energy conversion efficiency. Either the specified velocity, temperature, or " + "efficiency are incompatible"); PARTHENON_REQUIRE(kinetic_jet_velocity_ <= units.speed_of_light() * sqrt(2 * efficiency_), @@ -134,6 +137,11 @@ AGNFeedback::AGNFeedback(parthenon::ParameterInput *pin, PARTHENON_REQUIRE(kinetic_jet_temperature_ >= 0, "Kinetic jet temperature must be non-negative"); + // Compute the internal energy ceiling from the temperature ceiling + const Real tceil = pin->GetOrAddReal("problem/cluster/agn_feedback", "Tceil", + std::numeric_limits::infinity()); + eceil_ = tceil / mbar_gm1_over_kb; + // Add user history output variable for AGN power auto hst_vars = hydro_pkg->Param(parthenon::hist_param_key); if (!disabled_) { @@ -261,15 +269,18 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData *md, // Velocity of added gas const Real jet_velocity = kinetic_jet_velocity_; -#ifndef NDEBUG const Real jet_specific_internal_e = kinetic_jet_e_; -#endif // Amount of momentum density ( density * velocity) to dump in each cell const Real jet_momentum = jet_density * jet_velocity; // Amount of total energy to dump in each cell const Real jet_feedback = kinetic_fraction_ * power * kinetic_scaling_factor * beta_dt; + + const Real vceil = vceil_; + const Real vceil2 = SQR(vceil); + const Real eceil = eceil_; + const Real gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); //////////////////////////////////////////////////////////////////////////////// const parthenon::Real time = tm.time; @@ -278,7 +289,7 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData *md, hydro_pkg->Param("jet_coords_factory"); const JetCoords jet_coords = jet_coords_factory.CreateJetCoords(time); - // Constant volumetric heating + // Appy kinietic jet and thermal feedback parthenon::par_for( DEFAULT_LOOP_PATTERN, "HydroAGNFeedback::FeedbackSrcTerm", parthenon::DevExecSpace(), 0, cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, @@ -321,17 +332,15 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData *md, const Real sign_jet = (h > 0) ? 1 : -1; // Above or below jet-disk - /////////////////////////////////////////////////////////////////// - // We add the kinetic jet with a fixed jet velocity and specific - // internal energy/temperature of the added gas. The density, - // momentum, and total energy added depend on the triggered power. - /////////////////////////////////////////////////////////////////// + /////////////////////////////////////////////////////////////////// + // We add the kinetic jet with a fixed jet velocity and specific + // internal energy/temperature of the added gas. The density, + // momentum, and total energy added depend on the triggered power. + /////////////////////////////////////////////////////////////////// -#ifndef NDEBUG eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); const Real old_specific_internal_e = prim(IPR, k, j, i) / (prim(IDN, k, j, i) * (eos.GetGamma() - 1.)); -#endif cons(IDN, k, j, i) += jet_density; cons(IM1, k, j, i) += jet_momentum * sign_jet * jet_axis_x; @@ -339,20 +348,43 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData *md, cons(IM3, k, j, i) += jet_momentum * sign_jet * jet_axis_z; cons(IEN, k, j, i) += jet_feedback; -#ifndef NDEBUG eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); const Real new_specific_internal_e = prim(IPR, k, j, i) / (prim(IDN, k, j, i) * (eos.GetGamma() - 1.)); - PARTHENON_DEBUG_REQUIRE( + PARTHENON_REQUIRE( new_specific_internal_e > jet_specific_internal_e || new_specific_internal_e > old_specific_internal_e, "Kinetic injection leads to temperature below jet and existing gas"); -#endif + } + + // 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; + + // Remove kinetic energy + cons(IEN, k, j, i) -= 0.5 * prim(IDN, k, j, i) * (v2 - vceil2); + } + + // 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; } } + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); - PARTHENON_DEBUG_REQUIRE(prim(IPR, k, j, i) > 0, - "Kinetic injection leads to negative pressure"); + PARTHENON_REQUIRE(prim(IPR, k, j, i) > 0, + "Kinetic injection leads to negative pressure"); }); // Apply magnetic tower feedback diff --git a/src/pgen/cluster/agn_feedback.hpp b/src/pgen/cluster/agn_feedback.hpp index 1c88a61e..de516949 100644 --- a/src/pgen/cluster/agn_feedback.hpp +++ b/src/pgen/cluster/agn_feedback.hpp @@ -31,6 +31,9 @@ class AGNFeedback { // Efficiency converting mass to energy const parthenon::Real efficiency_; + // Velocity and temperature ceilings + parthenon::Real vceil_, eceil_; + // Thermal Heating Parameters const parthenon::Real thermal_radius_;