Merge branch 'main' into forrestglines/parthenon-pr930

parthenon-hpc-lab · Sep 26, 2023 · cee8eff · cee8eff
2 parents 60551ed + 6dc42ec
commit cee8eff
Show file tree

Hide file tree

Showing 9 changed files with 28 additions and 35 deletions.
diff --git a/docs/cluster.md b/docs/cluster.md
@@ -449,7 +449,7 @@ The parameters $\alpha$ and $\ell$ can be changed with
 ```
 <problem/cluster/magnetic_tower>
 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.
 ```
 <problem/cluster/magnetic_tower>
 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
@@ -182,7 +182,7 @@ kinetic_jet_offset = 0.05
 
 <problem/cluster/magnetic_tower>
 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
@@ -110,7 +110,7 @@ thermal_fraction = 0
 
 <problem/cluster/magnetic_tower>
 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
@@ -161,32 +161,13 @@ TaskStatus AddUnsplitSources(MeshData<Real> *md, const SimTime &tm, const Real b
 TaskStatus AddSplitSourcesFirstOrder(MeshData<Real> *md, const SimTime &tm) {
   auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
 
-  // const auto &enable_cooling = hydro_pkg->Param<Cooling>("enable_cooling");
-
-  // if (enable_cooling == Cooling::tabular) {
-  //   const TabularCooling &tabular_cooling =
-  //       hydro_pkg->Param<TabularCooling>("tabular_cooling");
-
-  //   tabular_cooling.SrcTerm(md, tm.dt);
-  // }
   if (ProblemSourceFirstOrder != nullptr) {
     ProblemSourceFirstOrder(md, tm, tm.dt);
   }
   return TaskStatus::complete;
 }
 
 TaskStatus AddSplitSourcesStrang(MeshData<Real> *md, const SimTime &tm) {
-  // auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
-
-  // const auto &enable_cooling = hydro_pkg->Param<Cooling>("enable_cooling");
-
-  // if (enable_cooling == Cooling::tabular) {
-  //   const TabularCooling &tabular_cooling =
-  //       hydro_pkg->Param<TabularCooling>("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
@@ -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<Real>("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
@@ -89,7 +89,7 @@ parthenon::Real LocalReduceAGNExtent(parthenon::MeshData<parthenon::Real> *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
@@ -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<StellarFeedback>("stellar_feedback", *this);
@@ -98,6 +109,7 @@ void StellarFeedback::FeedbackSrcTerm(parthenon::MeshData<parthenon::Real> *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<parthenon::Real> *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
@@ -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<parthenon::Real> *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 <typename EOS>
   void FeedbackSrcTerm(parthenon::MeshData<parthenon::Real> *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
@@ -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}",