Allow runtime def of tower potential

docs/cluster.md

@@ -412,7 +412,19 @@ material and mass added through the kinetic jet feedback mechanism.
 
 ### Magnetic feedback
 
-Magnetic feedback is injected following  ([Hui 2006](doi.org/10.1086/501499))
+Multiple options exists to inject magnetic fields:
+- a simple field loop (donut) configuration (better suited for kinetically dominated jets at larger scales)
+- a more complex pinching tower model (particularly suited for magnetically dominated jets at small scales)
+
+The option is controlled by the following parameter
+```
+<problem/cluster/magnetic_tower>
+potential_type = li # alternative: "donut"
+```
+
+#### Pinching magnetic tower
+
+Magnetic feedback is injected following  ([Li 2006](doi.org/10.1086/501499))
 where the injected magnetic field follows 
 
 $$
@@ -469,10 +481,29 @@ enable_magnetic_tower_mass_injection = false
 In this case, the injected mass through kinetic and thermal feedback
 according to their ratio.
 
+#### Simple field loop (donut) feedback
+
+Magnetic energy is injected according to the following simple potential
+$$
+A_h(r, \theta, h) = B_0 L \exp^\left ( -r^2/L^2 \right)$ for $h_\mathrm{offset} \leq |h| \leq h_\mathrm{offset} + h_\mathrm{thickness}
+$$
+resultig in a magnetic field configuration of
+$$
+B_\theta(r, \theta, h) = 2 B_0 r /L \exp^\left ( -r^2/L^2 \right)$ for $h_\mathrm{offset} \leq |h| \leq h_\mathrm{offset} + h_\mathrm{thickness}
+$$
+with all other components being zero.
+
+It is recommended to match the donut thickness to the thickness of the kinetic jet launching region
+and to choose a lengthscale that is half the jet launching radius.
+This results in almost all injected fields being confined to the launching region (and thus being
+carried along with a dominant jet).
+
+#### Fixed magnetic feedback
 
-A magnetic tower can also be inserted at runtime and injected at a fixed
+Magnetic feedback can also be inserted at runtime and injected at a fixed
 increase in magnetic field, and additional mass can be injected at a fixed
 rate.
+This works for all magnetic vector potentials.
 ```
 <problem/cluster/magnetic_tower>
 initial_field = 0.12431560000204142

inputs/cluster/cluster.in

@@ -181,7 +181,7 @@ kinetic_jet_thickness  = 0.05
 kinetic_jet_offset = 0.05
 
 <problem/cluster/magnetic_tower>
-
+potential_type = li
 alpha = 20
 l_scale = 0.001
 initial_field = 0.12431560000204142

inputs/cluster/magnetic_tower.in

@@ -109,6 +109,7 @@ kinetic_fraction = 0
 thermal_fraction = 0
 
 <problem/cluster/magnetic_tower>
+potential_type = li
 alpha = 20
 l_scale = 0.01
 initial_field = 0.12431560000204142

src/pgen/cluster/magnetic_tower.cpp

@@ -62,8 +62,9 @@ void MagneticTower::AddSrcTerm(parthenon::Real field_to_add, parthenon::Real mas
 
   const JetCoords jet_coords =
       hydro_pkg->Param<JetCoordsFactory>("jet_coords_factory").CreateJetCoords(tm.time);
-  const MagneticTowerObj mt = MagneticTowerObj(field_to_add, alpha_, l_scale_,
-                                               density_to_add, l_mass_scale_, jet_coords);
+  const MagneticTowerObj mt =
+      MagneticTowerObj(field_to_add, alpha_, l_scale_, density_to_add, l_mass_scale_,
+                       jet_coords, potential_);
 
   const auto &eos = hydro_pkg->Param<AdiabaticGLMMHDEOS>("eos");
 
@@ -164,7 +165,7 @@ void MagneticTower::ReducePowerContribs(parthenon::Real &linear_contrib,
 
   // Make a construct a copy of this with field strength 1 to send to the device
   const MagneticTowerObj mt =
-      MagneticTowerObj(1, alpha_, l_scale_, 0, l_mass_scale_, jet_coords);
+      MagneticTowerObj(1, alpha_, l_scale_, 0, l_mass_scale_, jet_coords, potential_);
 
   // Get the reduction of the linear and quadratic contributions ready
   Real linear_contrib_red, quadratic_contrib_red;
@@ -217,7 +218,7 @@ void MagneticTower::AddInitialFieldToPotential(parthenon::MeshBlock *pmb,
   const JetCoords jet_coords =
       hydro_pkg->Param<JetCoordsFactory>("jet_coords_factory").CreateJetCoords(0.0);
   const MagneticTowerObj mt(initial_field_, alpha_, l_scale_, 0, l_mass_scale_,
-                            jet_coords);
+                            jet_coords, potential_);
 
   parthenon::par_for(
       DEFAULT_LOOP_PATTERN, "MagneticTower::AddInitialFieldToPotential",

src/pgen/cluster/magnetic_tower.hpp

@@ -17,8 +17,11 @@
 #include <parthenon/package.hpp>
 
 #include "jet_coords.hpp"
+#include "utils/error_checking.hpp"
 
 namespace cluster {
+
+enum class MagneticTowerPotential { undefined, li, donut };
 /************************************************************
  *  Magnetic Tower Object, for computing magnetic field, vector potential at a
  *  fixed time with a fixed field
@@ -32,32 +35,49 @@ class MagneticTowerObj {
   const parthenon::Real density_, l_mass_scale2_;
 
   JetCoords jet_coords_;
+  // Note that this eventually might better be a template parameter, but while the number
+  // of potentials implemented is limited (and similarly complex) this should currently
+  // not be a performance concern.
+  const MagneticTowerPotential potential_;
 
  public:
   MagneticTowerObj(const parthenon::Real field, const parthenon::Real alpha,
                    const parthenon::Real l_scale, const parthenon::Real density,
-                   const parthenon::Real l_mass_scale, const JetCoords jet_coords)
+                   const parthenon::Real l_mass_scale, const JetCoords jet_coords,
+                   const MagneticTowerPotential potential)
       : field_(field), alpha_(alpha), l_scale_(l_scale), density_(density),
-        l_mass_scale2_(SQR(l_mass_scale)), jet_coords_(jet_coords) {
+        l_mass_scale2_(SQR(l_mass_scale)), jet_coords_(jet_coords),
+        potential_(potential) {
     PARTHENON_REQUIRE(l_scale > 0,
                       "Magnetic Tower Length scale must be strictly postitive");
-    PARTHENON_REQUIRE(l_mass_scale >= 0,
-                      "Magnetic Tower Mass Length scale must be zero (disabled) or postitive");
+    PARTHENON_REQUIRE(
+        l_mass_scale >= 0,
+        "Magnetic Tower Mass Length scale must be zero (disabled) or postitive");
   }
 
   // Compute Jet Potential in jet cylindrical coordinates
   KOKKOS_INLINE_FUNCTION void
   PotentialInJetCyl(const parthenon::Real r, const parthenon::Real h,
                     parthenon::Real &a_r, parthenon::Real &a_theta,
                     parthenon::Real &a_h) const __attribute__((always_inline)) {
-    const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2));
-    // Compute the potential in jet cylindrical coordinates
-    a_r = 0.0;
-    a_theta = 0.0;
-    if (fabs(h) >= 0.001 && fabs(h) <= 0.001 + 0.000390625) {
-      a_h = field_ * l_scale_ * exp_r2_h2;
+    if (potential_ == MagneticTowerPotential::donut) {
+      const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2));
+      // Compute the potential in jet cylindrical coordinates
+      a_r = 0.0;
+      a_theta = 0.0;
+      if (fabs(h) >= 0.001 && fabs(h) <= 0.001 + 0.000390625) {
+        a_h = field_ * l_scale_ * exp_r2_h2;
+      } else {
+        a_h = 0.0;
+      }
+    } else if (potential_ == MagneticTowerPotential::li) {
+      const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2) - pow(h / l_scale_, 2));
+      // Compute the potential in jet cylindrical coordinates
+      a_r = 0.0;
+      a_theta = field_ * l_scale_ * (r / l_scale_) * exp_r2_h2;
+      a_h = field_ * l_scale_ * alpha_ / 2.0 * exp_r2_h2;
     } else {
-      a_h = 0.0;
+      PARTHENON_FAIL("Unknown magnetic tower potential.");
     }
   }
 
@@ -84,16 +104,25 @@ class MagneticTowerObj {
   FieldInJetCyl(const parthenon::Real r, const parthenon::Real h, parthenon::Real &b_r,
                 parthenon::Real &b_theta, parthenon::Real &b_h) const
       __attribute__((always_inline)) {
-
-    const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2));
-    // Compute the field in jet cylindrical coordinates
-    b_r = 0.0;
-    if (fabs(h) >= 0.001 && fabs(h) <= 0.001 + 0.000390625) {
-      b_theta = 2.0 * field_ * r / l_scale_ * exp_r2_h2;
+    if (potential_ == MagneticTowerPotential::donut) {
+      const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2));
+      // Compute the field in jet cylindrical coordinates
+      b_r = 0.0;
+      if (fabs(h) >= 0.001 && fabs(h) <= 0.001 + 0.000390625) {
+        b_theta = 2.0 * field_ * r / l_scale_ * exp_r2_h2;
+      } else {
+        b_theta = 0.0;
+      }
+      b_h = 0.0;
+    } else if (potential_ == MagneticTowerPotential::li) {
+      const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2) - pow(h / l_scale_, 2));
+      // Compute the field in jet cylindrical coordinates
+      b_r = field_ * 2 * (h / l_scale_) * (r / l_scale_) * exp_r2_h2;
+      b_theta = field_ * alpha_ * (r / l_scale_) * exp_r2_h2;
+      b_h = field_ * 2 * (1 - pow(r / l_scale_, 2)) * exp_r2_h2;
     } else {
-      b_theta = 0.0;
+      PARTHENON_FAIL("Unknown magnetic tower potential.");
     }
-    b_h = 0.0;
   }
 
   // Compute Magnetic field in Simulation Cartesian coordinates
@@ -141,24 +170,41 @@ class MagneticTower {
   const parthenon::Real fixed_mass_rate_;
   const parthenon::Real l_mass_scale_;
 
+  MagneticTowerPotential potential_;
+
   MagneticTower(parthenon::ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg,
                 const std::string &block = "problem/cluster/magnetic_tower")
       : alpha_(pin->GetOrAddReal(block, "alpha", 0)),
         l_scale_(pin->GetOrAddReal(block, "l_scale", 0)),
         initial_field_(pin->GetOrAddReal(block, "initial_field", 0)),
         fixed_field_rate_(pin->GetOrAddReal(block, "fixed_field_rate", 0)),
         fixed_mass_rate_(pin->GetOrAddReal(block, "fixed_mass_rate", 0)),
-        l_mass_scale_(pin->GetOrAddReal(block, "l_mass_scale", 0)) {
-    hydro_pkg->AddParam<>("magnetic_tower", *this);
+        l_mass_scale_(pin->GetOrAddReal(block, "l_mass_scale", 0)),
+        potential_(MagneticTowerPotential::undefined) {
     hydro_pkg->AddParam<parthenon::Real>("magnetic_tower_linear_contrib", 0.0, true);
     hydro_pkg->AddParam<parthenon::Real>("magnetic_tower_quadratic_contrib", 0.0, true);
 
+    const auto potential_str = pin->GetOrAddString(block, "potential_type", "undefined");
+
+    if (potential_str == "donut") {
+      potential_ = MagneticTowerPotential::donut;
+    } else if (potential_str == "li") {
+      potential_ = MagneticTowerPotential::li;
+    } else {
+      PARTHENON_FAIL(
+          "Unknown potential for magnetic tower. Current options are: donut, li")
+    }
+
     // Vector potential is only locally used, so no need to
     // communicate/restrict/prolongate/fluxes/etc
     parthenon::Metadata m({parthenon::Metadata::Cell, parthenon::Metadata::Derived,
                            parthenon::Metadata::OneCopy},
                           std::vector<int>({3}));
     hydro_pkg->AddField("magnetic_tower_A", m);
+
+    // Finally, add object to params (should be done last as otherwise modification within
+    // this function would not survive).
+    hydro_pkg->AddParam<>("magnetic_tower", *this);
   }
 
   // Add initial magnetic field to provided potential with a single meshblock