Add jet tracing

docs/cluster.md

@@ -306,6 +306,8 @@ kinetic_fraction = 0.3333
 ```
 These values are automatically normalized to sum to 1.0 at run time.
 
+### Thermal feedback
+
 Thermal feedback is deposited at a flat power density within a sphere of defined radius
 ```
 <problem/cluster/agn_feedback>
@@ -318,6 +320,8 @@ feedback, e.g.
 thermal_injected_mass = mdot * (1 - efficiency) * normalized_thermal_fraction;
 ```
 
+### Kinetic feedback
+
 Kinetic feedback is deposited into two disks along the axis of the jet within a
 defined radius, thickness of each disk, and an offset above and below the plane
 of the AGN disk where each disk begins.
@@ -381,6 +385,31 @@ energy density; changes in kinetic energy density will depend on the velocity
 of the ambient gas. Temperature will also change but should always increase
 with kinetic jet feedback.
 
+#### Tracing jet material
+
+Material launched by the kinetic jet component can be traced by setting
+```
+<problem/cluster/agn_feedback>
+enable_tracer = true
+```
+
+At the moment, this also requires to enable a single passive scalar by setting
+```
+<hydro>
+nscalars = 1
+```
+This may change in the future as the current implementation restricts the
+passive scalar use to tracing jet material.
+This is a design decision motivated by simplicity and not for technical
+reasons (KISS!).
+
+Note that all material launched from within the jet region is traced, i.e.,
+passive scalar concentration does not differentiate between original cell
+material and mass added through the kinetic jet feedback mechanism.
+
+
+### Magnetic feedback
+
 Magnetic feedback is injected following  ([Hui 2006](doi.org/10.1086/501499))
 where the injected magnetic field follows 
 

src/pgen/cluster/agn_feedback.cpp

@@ -26,6 +26,7 @@
 #include "agn_triggering.hpp"
 #include "cluster_utils.hpp"
 #include "magnetic_tower.hpp"
+#include "utils/error_checking.hpp"
 
 namespace cluster {
 using namespace parthenon;
@@ -50,6 +51,8 @@ AGNFeedback::AGNFeedback(parthenon::ParameterInput *pin,
                                                "kinetic_jet_thickness", 0.02)),
       kinetic_jet_offset_(
           pin->GetOrAddReal("problem/cluster/agn_feedback", "kinetic_jet_offset", 0.02)),
+      enable_tracer_(
+          pin->GetOrAddBoolean("problem/cluster/agn_feedback", "enable_tracer", false)),
       disabled_(pin->GetOrAddBoolean("problem/cluster/agn_feedback", "disabled", false)) {
 
   // Normalize the thermal, kinetic, and magnetic fractions to sum to 1.0
@@ -160,6 +163,10 @@ 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,
+                           "Enabling tracer for AGN feedback requires hydro/nscalars=1");
+
   hydro_pkg->AddParam<>("agn_feedback", *this);
 }
 
@@ -281,6 +288,7 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData<parthenon::Real> *md,
   const Real vceil2 = SQR(vceil);
   const Real eceil = eceil_;
   const Real gm1 = (hydro_pkg->Param<Real>("AdiabaticIndex") - 1.0);
+  const auto enable_tracer = enable_tracer_;
   ////////////////////////////////////////////////////////////////////////////////
 
   const parthenon::Real time = tm.time;
@@ -348,6 +356,14 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData<parthenon::Real> *md,
             cons(IM3, k, j, i) += jet_momentum * sign_jet * jet_axis_z;
             cons(IEN, k, j, i) += jet_feedback;
 
+            // Reset tracer to one for the entire material in the jet launching region as
+            // we cannot distinguish between original material in a cell and new jet
+            // material in the evolution of the jet. Eventually, we're just interested in
+            // stuff that came from here.
+            if (enable_tracer) {
+              cons(nhydro, k, j, i) = 1.0 * cons(IDN, k, j, i);
+            }
+
             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.));

src/pgen/cluster/agn_feedback.hpp

@@ -41,6 +41,9 @@ class AGNFeedback {
   const parthenon::Real kinetic_jet_radius_, kinetic_jet_thickness_, kinetic_jet_offset_;
   parthenon::Real kinetic_jet_velocity_, kinetic_jet_temperature_, kinetic_jet_e_;
 
+  // enable passive scalar to trace AGN material
+  const bool enable_tracer_;
+
   const bool disabled_;
 
   AGNFeedback(parthenon::ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg);