Added exclusion radius to stellar feedback

parthenon-hpc-lab · Sep 10, 2023 · 2bdf1f8 · 2bdf1f8
1 parent 160d10a
commit 2bdf1f8
Show file tree

Hide file tree

Showing 2 changed files with 15 additions and 4 deletions.
diff --git a/src/pgen/cluster/stellar_feedback.cpp b/src/pgen/cluster/stellar_feedback.cpp
@@ -33,19 +33,28 @@ 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 &&
+  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 &&
+
+  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.");
@@ -98,6 +107,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 +135,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,