Add magic heating from stellar feedback

docs/cluster.md

@@ -390,7 +390,7 @@ with kinetic jet feedback.
 Material launched by the kinetic jet component can be traced by setting
 ```
 <problem/cluster/agn_feedback>
-enable_tracer = true
+enable_tracer = true   # disabled by default
 ```
 
 At the moment, this also requires to enable a single passive scalar by setting
@@ -483,3 +483,27 @@ where `power_per_bcg_mass` and `mass_rate_per_bcg_mass` is the power and mass
 per time respectively injected per BCG mass at a given radius. This SNIA
 feedback is otherwise fixed in time, spherically symmetric, and dependant on
 the BCG specified in `<problem/cluster/gravity>`.
+
+## Stellar feedback
+
+Cold, dense, potentially magnetically supported gas may accumulate in a small
+disk-like structure around the AGN depending on the specific setup.
+In reality, this gas would form stars.
+
+In absence of star particles (or similar) in the current setup, we use a simple
+formulation to convert some fraction of this gas to thermal energy.
+More specifically, gas within a given radius, above a critical number density
+and below a temperature threshold, will be converted to thermal energy with a
+given efficiency.
+The parameters can be set as follows (numerical values here just for illustration purpose -- by default they are all 0, i.e., stellar feedback is disabled).
+
+```
+<problem/cluster/stellar_feedback>
+stellar_radius = 0.025 # in code length
+efficiency = 5e-6
+number_density_threshold = 2.93799894e+74 # in code_length**-3
+temperature_threshold = 2e4 # in K
+```
+
+Note that all parameters need to be specified explicitly for the feedback to work
+(i.e., no hidden default values).

src/pgen/CMakeLists.txt

@@ -11,6 +11,7 @@ target_sources(athenaPK PRIVATE
     cluster/hydrostatic_equilibrium_sphere.cpp
     cluster/magnetic_tower.cpp
     cluster/snia_feedback.cpp
+    cluster/stellar_feedback.cpp
     cpaw.cpp
     diffusion.cpp
     field_loop.cpp

src/pgen/cluster.cpp

@@ -9,7 +9,7 @@
 // Setups up an idealized galaxy cluster with an ACCEPT-like entropy profile in
 // hydrostatic equilbrium with an NFW+BCG+SMBH gravitational profile,
 // optionally with an initial magnetic tower field. Includes AGN feedback, AGN
-// triggering via cold gas, simple SNIA Feedback(TODO)
+// triggering via cold gas, simple SNIA Feedback, and simple stellar feedback
 //========================================================================================
 
 // C headers
@@ -28,6 +28,8 @@
 #include "kokkos_abstraction.hpp"
 #include "mesh/domain.hpp"
 #include "mesh/mesh.hpp"
+#include "parthenon_array_generic.hpp"
+#include "utils/error_checking.hpp"
 #include <parthenon/driver.hpp>
 #include <parthenon/package.hpp>
 
@@ -48,8 +50,7 @@
 #include "cluster/hydrostatic_equilibrium_sphere.hpp"
 #include "cluster/magnetic_tower.hpp"
 #include "cluster/snia_feedback.hpp"
-#include "parthenon_array_generic.hpp"
-#include "utils/error_checking.hpp"
+#include "cluster/stellar_feedback.hpp"
 
 namespace cluster {
 using namespace parthenon::driver::prelude;
@@ -194,6 +195,9 @@ void ClusterSrcTerm(MeshData<Real> *md, const parthenon::SimTime &tm,
   const auto &snia_feedback = hydro_pkg->Param<SNIAFeedback>("snia_feedback");
   snia_feedback.FeedbackSrcTerm(md, beta_dt, tm);
 
+  const auto &stellar_feedback = hydro_pkg->Param<StellarFeedback>("stellar_feedback");
+  stellar_feedback.FeedbackSrcTerm(md, beta_dt, tm);
+
   ApplyClusterClips(md, tm, beta_dt);
 };
 
@@ -341,6 +345,12 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd
 
   SNIAFeedback snia_feedback(pin, hydro_pkg);
 
+  /************************************************************
+   * Read Stellar Feedback
+   ************************************************************/
+
+  StellarFeedback stellar_feedback(pin, hydro_pkg);
+
   /************************************************************
    * Read Clips  (ceilings and floors)
    ************************************************************/

src/pgen/cluster/agn_feedback.cpp

@@ -164,7 +164,7 @@ AGNFeedback::AGNFeedback(parthenon::ParameterInput *pin,
   hydro_pkg->UpdateParam(parthenon::hist_param_key, hst_vars);
 
   // Double check that tracers are also enabled in fluid solver
-  PARTHENON_REQUIRE_THROWS(enable_tracer_ && hydro_pkg->Param<int>("nscalars") == 1,
+  PARTHENON_REQUIRE_THROWS(!enable_tracer_ || hydro_pkg->Param<int>("nscalars") == 1,
                            "Enabling tracer for AGN feedback requires hydro/nscalars=1");
 
   hydro_pkg->AddParam<>("agn_feedback", *this);

src/pgen/cluster/snia_feedback.hpp

@@ -47,4 +47,4 @@ class SNIAFeedback {
 
 } // namespace cluster
 
-#endif // CLUSTER_SNIAE_FEEDBACK_HPP_
+#endif // CLUSTER_SNIA_FEEDBACK_HPP_

src/pgen/cluster/stellar_feedback.cpp

@@ -0,0 +1,151 @@
+//========================================================================================
+// AthenaPK - a performance portable block structured AMR astrophysical MHD code.
+// Copyright (c) 2023, Athena-Parthenon Collaboration. All rights reserved.
+// Licensed under the 3-clause BSD License, see LICENSE file for details
+//========================================================================================
+//! \file stellar_feedback.cpp
+//  \brief  Class for magic heating modeling star formation
+
+#include <cmath>
+
+// Parthenon headers
+#include <coordinates/uniform_cartesian.hpp>
+#include <globals.hpp>
+#include <interface/state_descriptor.hpp>
+#include <mesh/domain.hpp>
+#include <parameter_input.hpp>
+#include <parthenon/package.hpp>
+
+// Athena headers
+#include "../../eos/adiabatic_glmmhd.hpp"
+#include "../../eos/adiabatic_hydro.hpp"
+#include "../../main.hpp"
+#include "../../units.hpp"
+#include "cluster_gravity.hpp"
+#include "cluster_utils.hpp"
+#include "stellar_feedback.hpp"
+#include "utils/error_checking.hpp"
+
+namespace cluster {
+using namespace parthenon;
+
+StellarFeedback::StellarFeedback(parthenon::ParameterInput *pin,
+                                 parthenon::StateDescriptor *hydro_pkg)
+    : stellar_radius_(
+          pin->GetOrAddReal("problem/cluster/stellar_feedback", "stellar_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) {
+    disabled_ = true;
+  }
+  PARTHENON_REQUIRE(disabled_ || (stellar_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);
+}
+
+void StellarFeedback::FeedbackSrcTerm(parthenon::MeshData<parthenon::Real> *md,
+                                      const parthenon::Real beta_dt,
+                                      const parthenon::SimTime &tm) const {
+  auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
+  auto fluid = hydro_pkg->Param<Fluid>("fluid");
+  if (fluid == Fluid::euler) {
+    FeedbackSrcTerm(md, beta_dt, tm, hydro_pkg->Param<AdiabaticHydroEOS>("eos"));
+  } else if (fluid == Fluid::glmmhd) {
+    FeedbackSrcTerm(md, beta_dt, tm, hydro_pkg->Param<AdiabaticGLMMHDEOS>("eos"));
+  } else {
+    PARTHENON_FAIL("StellarFeedback::FeedbackSrcTerm: Unknown EOS");
+  }
+}
+template <typename EOS>
+void StellarFeedback::FeedbackSrcTerm(parthenon::MeshData<parthenon::Real> *md,
+                                      const parthenon::Real beta_dt,
+                                      const parthenon::SimTime &tm,
+                                      const EOS &eos) const {
+  using parthenon::IndexDomain;
+  using parthenon::IndexRange;
+  using parthenon::Real;
+
+  auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
+
+  if (disabled_) {
+    // No stellar feedback, return
+    return;
+  }
+
+  // Grab some necessary variables
+  const auto &prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
+  const auto &cons_pack = md->PackVariables(std::vector<std::string>{"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<int>("nhydro");
+  const auto nscalars = hydro_pkg->Param<int>("nscalars");
+
+  // const auto gm1 = (hydro_pkg->Param<Real>("AdiabaticIndex") - 1.0);
+  const auto units = hydro_pkg->Param<Units>("units");
+  const auto mbar = hydro_pkg->Param<Real>("mu") * units.mh();
+  const auto mbar_over_kb = hydro_pkg->Param<Real>("mbar_over_kb");
+
+  const auto mass_to_energy = efficiency_ * SQR(units.speed_of_light());
+  const auto stellar_radius = stellar_radius_;
+  const auto temperature_threshold = temperatue_threshold_;
+  const auto number_density_threshold = number_density_threshold_;
+
+  ////////////////////////////////////////////////////////////////////////////////
+
+  // Constant volumetric heating
+  parthenon::par_for(
+      DEFAULT_LOOP_PATTERN, "StellarFeedback::FeedbackSrcTerm", 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 auto x = coords.Xc<1>(i);
+        const auto y = coords.Xc<2>(j);
+        const auto z = coords.Xc<3>(k);
+
+        const auto r = sqrt(x * x + y * y + z * z);
+        if (r > stellar_radius) {
+          return;
+        }
+
+        auto number_density = prim(IDN, k, j, i) / mbar;
+        if (number_density < number_density_threshold) {
+          return;
+        }
+
+        auto temp = mbar_over_kb * prim(IPR, k, j, i) / prim(IDN, k, j, i);
+        if (temp > temperature_threshold) {
+          return;
+        }
+
+        // All conditions to convert mass to energy are met
+        const auto cell_delta_rho = number_density_threshold * mbar - prim(IDN, k, j, i);
+
+        // First remove density at fixed temperature
+        AddDensityToConsAtFixedVelTemp(cell_delta_rho, cons, prim, eos.GetGamma(), k, j,
+                                       i);
+        //  Then add thermal energy
+        const auto cell_delta_mass = -cell_delta_rho * coords.CellVolume(k, j, i);
+        const auto cell_delta_energy_density =
+            mass_to_energy * cell_delta_mass / coords.CellVolume(k, j, i);
+        PARTHENON_REQUIRE(cell_delta_energy_density > 0.0,
+                          "Sanity check failed. Added thermal energy should be positive.")
+        cons(IEN, k, j, i) += cell_delta_energy_density;
+
+        // Update prims
+        eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
+      });
+}
+
+} // namespace cluster

src/pgen/cluster/stellar_feedback.hpp

@@ -0,0 +1,50 @@
+#ifndef CLUSTER_STELLAR_FEEDBACK_HPP_
+#define CLUSTER_STELLAR_FEEDBACK_HPP_
+//========================================================================================
+// AthenaPK - a performance portable block structured AMR astrophysical MHD code.
+// Copyright (c) 2023, Athena-Parthenon Collaboration. All rights reserved.
+// Licensed under the 3-clause BSD License, see LICENSE file for details
+//========================================================================================
+//! \file stellar_feedback.hpp
+//  \brief  Class for injecting Stellar feedback following BCG density
+
+// parthenon headers
+#include <basic_types.hpp>
+#include <mesh/domain.hpp>
+#include <mesh/mesh.hpp>
+#include <parameter_input.hpp>
+#include <parthenon/package.hpp>
+
+#include "jet_coords.hpp"
+
+namespace cluster {
+
+/************************************************************
+ *  StellarFeedback
+ ************************************************************/
+class StellarFeedback {
+ private:
+  // feedback parameters in code units
+  const parthenon::Real stellar_radius_;           // length
+  const parthenon::Real efficiency_;               // dimless
+  const parthenon::Real number_density_threshold_; // 1/(length^3)
+  const parthenon::Real temperatue_threshold_;     // K
+
+  bool disabled_;
+
+ public:
+  StellarFeedback(parthenon::ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg);
+
+  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
+  template <typename EOS>
+  void FeedbackSrcTerm(parthenon::MeshData<parthenon::Real> *md,
+                       const parthenon::Real beta_dt, const parthenon::SimTime &tm,
+                       const EOS &eos) const;
+};
+
+} // namespace cluster
+
+#endif // CLUSTER_STELLAR_FEEDBACK_HPP_