diff --git a/docs/cluster.md b/docs/cluster.md index 7e650cfa..b4d83983 100644 --- a/docs/cluster.md +++ b/docs/cluster.md @@ -449,7 +449,7 @@ The parameters $\alpha$ and $\ell$ can be changed with ``` potential_type = li -alpha = 20 +li_alpha = 20 l_scale = 0.001 ``` When injected as a fraction of @@ -498,9 +498,9 @@ with all other components being zero. ``` potential_type = donut -l_scale = 0.0005 # in code length -offset = 0.001 # in code length -thickness = 0.0001 # in code length +l_scale = 0.0005 # in code length +donut_offset = 0.001 # in code length +donut_thickness = 0.0001 # in code length ``` It is recommended to match the donut thickness to the thickness of the kinetic jet launching region diff --git a/inputs/cluster/cluster.in b/inputs/cluster/cluster.in index fe1ff02b..e3192524 100644 --- a/inputs/cluster/cluster.in +++ b/inputs/cluster/cluster.in @@ -182,7 +182,7 @@ kinetic_jet_offset = 0.05 potential_type = li -alpha = 20 +li_alpha = 20 l_scale = 0.001 initial_field = 0.12431560000204142 #NOTE: Change to this line to disable initial magnetic tower diff --git a/inputs/cluster/magnetic_tower.in b/inputs/cluster/magnetic_tower.in index 15b5647d..44ba58a0 100644 --- a/inputs/cluster/magnetic_tower.in +++ b/inputs/cluster/magnetic_tower.in @@ -110,7 +110,7 @@ thermal_fraction = 0 potential_type = li -alpha = 20 +li_alpha = 20 l_scale = 0.01 initial_field = 0.12431560000204142 fixed_field_rate = 12.431560000204144 diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index b46db6b2..e5402454 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -161,14 +161,6 @@ TaskStatus AddUnsplitSources(MeshData *md, const SimTime &tm, const Real b TaskStatus AddSplitSourcesFirstOrder(MeshData *md, const SimTime &tm) { auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); - // const auto &enable_cooling = hydro_pkg->Param("enable_cooling"); - - // if (enable_cooling == Cooling::tabular) { - // const TabularCooling &tabular_cooling = - // hydro_pkg->Param("tabular_cooling"); - - // tabular_cooling.SrcTerm(md, tm.dt); - // } if (ProblemSourceFirstOrder != nullptr) { ProblemSourceFirstOrder(md, tm, tm.dt); } @@ -176,17 +168,6 @@ TaskStatus AddSplitSourcesFirstOrder(MeshData *md, const SimTime &tm) { } TaskStatus AddSplitSourcesStrang(MeshData *md, const SimTime &tm) { - // auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); - - // const auto &enable_cooling = hydro_pkg->Param("enable_cooling"); - - // if (enable_cooling == Cooling::tabular) { - // const TabularCooling &tabular_cooling = - // hydro_pkg->Param("tabular_cooling"); - - // tabular_cooling.SrcTerm(md, 0.5 * tm.dt); - // } - if (ProblemSourceStrangSplit != nullptr) { ProblemSourceStrangSplit(md, tm, tm.dt); } diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index d0f981d2..a26b1da6 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -327,7 +327,6 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd const Real agn_tracer_thresh = pin->GetOrAddReal("problem/cluster/reductions", "agn_tracer_thresh", -1.0); if (agn_tracer_thresh >= 0) { - auto mbar_over_kb = hydro_pkg->Param("mbar_over_kb"); PARTHENON_REQUIRE( pin->GetOrAddBoolean("problem/cluster/agn_feedback", "enable_tracer", false), "AGN Tracer must be enabled to reduce AGN tracer extent"); diff --git a/src/pgen/cluster/cluster_reductions.cpp b/src/pgen/cluster/cluster_reductions.cpp index 0e8c88ed..9e9edb4b 100644 --- a/src/pgen/cluster/cluster_reductions.cpp +++ b/src/pgen/cluster/cluster_reductions.cpp @@ -89,7 +89,7 @@ parthenon::Real LocalReduceAGNExtent(parthenon::MeshData *md) { const auto r2 = SQR(coords.Xc<1>(k, j, i)) + SQR(coords.Xc<2>(k, j, i)) + SQR(coords.Xc<3>(k, j, i)); - if (cons(nhydro, k, j, i) > tracer_thresh && r2 > max_r2) { + if (cons(nhydro, k, j, i) / cons(IDN, k, j, i) > tracer_thresh && r2 > max_r2) { max_r2_team = r2; } }, diff --git a/src/pgen/cluster/stellar_feedback.cpp b/src/pgen/cluster/stellar_feedback.cpp index f4eb774e..aa130d65 100644 --- a/src/pgen/cluster/stellar_feedback.cpp +++ b/src/pgen/cluster/stellar_feedback.cpp @@ -33,21 +33,32 @@ StellarFeedback::StellarFeedback(parthenon::ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg) : stellar_radius_( pin->GetOrAddReal("problem/cluster/stellar_feedback", "stellar_radius", 0.0)), + exclusion_radius_( + pin->GetOrAddReal("problem/cluster/stellar_feedback", "exclusion_radius", 0.0)), efficiency_( pin->GetOrAddReal("problem/cluster/stellar_feedback", "efficiency", 0.0)), number_density_threshold_(pin->GetOrAddReal("problem/cluster/stellar_feedback", "number_density_threshold", 0.0)), temperatue_threshold_(pin->GetOrAddReal("problem/cluster/stellar_feedback", "temperature_threshold", 0.0)) { - if (stellar_radius_ == 0.0 && efficiency_ == 0.0 && number_density_threshold_ == 0.0 && - temperatue_threshold_ == 0.0) { + if (stellar_radius_ == 0.0 && exclusion_radius_ == 0.0 && efficiency_ == 0.0 && + number_density_threshold_ == 0.0 && temperatue_threshold_ == 0.0) { disabled_ = true; } else { disabled_ = false; } - PARTHENON_REQUIRE(disabled_ || (stellar_radius_ != 0.0 && efficiency_ != 0.00 && - number_density_threshold_ != 0.0 && - temperatue_threshold_ != 0.0), + + if (!disabled_ && exclusion_radius_ == 0.0) { + // If exclusion_radius_ is not specified, use AGN triggering accretion radius. + // If both are zero, the PARTHENON_REQUIRE will fail + exclusion_radius_ = + pin->GetOrAddReal("problem/cluster/agn_triggering", "accretion_radius", 0); + } + + PARTHENON_REQUIRE(disabled_ || + (stellar_radius_ != 0.0 && exclusion_radius_ != 0.0 && + efficiency_ != 0.00 && number_density_threshold_ != 0.0 && + temperatue_threshold_ != 0.0), "Enabling stellar feedback requires setting all parameters."); hydro_pkg->AddParam("stellar_feedback", *this); @@ -98,6 +109,7 @@ void StellarFeedback::FeedbackSrcTerm(parthenon::MeshData *md, const auto mass_to_energy = efficiency_ * SQR(units.speed_of_light()); const auto stellar_radius = stellar_radius_; + const auto exclusion_radius = exclusion_radius_; const auto temperature_threshold = temperatue_threshold_; const auto number_density_threshold = number_density_threshold_; @@ -125,7 +137,7 @@ void StellarFeedback::FeedbackSrcTerm(parthenon::MeshData *md, const auto z = coords.Xc<3>(k); const auto r = sqrt(x * x + y * y + z * z); - if (r > stellar_radius) { + if (r > stellar_radius || r <= exclusion_radius) { return; } diff --git a/src/pgen/cluster/stellar_feedback.hpp b/src/pgen/cluster/stellar_feedback.hpp index d91dff6a..e3d4b20c 100644 --- a/src/pgen/cluster/stellar_feedback.hpp +++ b/src/pgen/cluster/stellar_feedback.hpp @@ -26,6 +26,7 @@ class StellarFeedback { private: // feedback parameters in code units const parthenon::Real stellar_radius_; // length + parthenon::Real exclusion_radius_; // length const parthenon::Real efficiency_; // dimless const parthenon::Real number_density_threshold_; // 1/(length^3) const parthenon::Real temperatue_threshold_; // K @@ -38,7 +39,7 @@ class StellarFeedback { void FeedbackSrcTerm(parthenon::MeshData *md, const parthenon::Real beta_dt, const parthenon::SimTime &tm) const; - // Apply the feedback from stellare tied to the BCG density + // Apply stellar feedback following cold gas density above a density threshold template void FeedbackSrcTerm(parthenon::MeshData *md, const parthenon::Real beta_dt, const parthenon::SimTime &tm, diff --git a/tst/regression/test_suites/cluster_magnetic_tower/cluster_magnetic_tower.py b/tst/regression/test_suites/cluster_magnetic_tower/cluster_magnetic_tower.py index 4da1eeb1..c80ef354 100644 --- a/tst/regression/test_suites/cluster_magnetic_tower/cluster_magnetic_tower.py +++ b/tst/regression/test_suites/cluster_magnetic_tower/cluster_magnetic_tower.py @@ -275,7 +275,7 @@ def Prepare(self, parameters, step): f"problem/cluster/agn_feedback/magnetic_fraction=1", f"problem/cluster/agn_feedback/kinetic_fraction=0", f"problem/cluster/agn_feedback/thermal_fraction=0", - f"problem/cluster/magnetic_tower/alpha={self.magnetic_tower_alpha}", + f"problem/cluster/magnetic_tower/li_alpha={self.magnetic_tower_alpha}", f"problem/cluster/magnetic_tower/l_scale={self.magnetic_tower_l_scale.in_units('code_length').v}", f"problem/cluster/magnetic_tower/initial_field={self.initial_magnetic_tower_field.in_units('sqrt(code_mass)/sqrt(code_length)/code_time').v}", f"problem/cluster/magnetic_tower/fixed_field_rate={fixed_field_rate}",