Add offset and thickness params to donut feedback

parthenon-hpc-lab · Aug 18, 2023 · 1b27b3a · 1b27b3a
1 parent ecd1f6d
commit 1b27b3a
Show file tree

Hide file tree

Showing 3 changed files with 36 additions and 12 deletions.
diff --git a/docs/cluster.md b/docs/cluster.md
@@ -448,6 +448,7 @@ $$
 The parameters $\alpha$ and $\ell$ can be changed with
 ```
 <problem/cluster/magnetic_tower>
+potential_type = li
 alpha = 20
 l_scale = 0.001
 ```
@@ -493,11 +494,23 @@ B_\theta(r, \theta, h) = 2 B_0 r /L \exp^\left ( -r^2/L^2 \right)$ for $h_\mathr
 $$
 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
+```
+
 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).
 
+Mass feedback for the donut model is currently identical to the Li tower above
+(and, thus, the parameters above pertaining to the mass injection equally apply here).
+
 #### Fixed magnetic feedback
 
 Magnetic feedback can also be inserted at runtime and injected at a fixed

diff --git a/src/pgen/cluster/magnetic_tower.cpp b/src/pgen/cluster/magnetic_tower.cpp
@@ -63,8 +63,8 @@ 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, potential_);
+      MagneticTowerObj(field_to_add, alpha_, l_scale_, offset_, thickness_,
+                       density_to_add, l_mass_scale_, jet_coords, potential_);
 
   const auto &eos = hydro_pkg->Param<AdiabaticGLMMHDEOS>("eos");
 
@@ -164,8 +164,8 @@ void MagneticTower::ReducePowerContribs(parthenon::Real &linear_contrib,
       hydro_pkg->Param<JetCoordsFactory>("jet_coords_factory").CreateJetCoords(tm.time);
 
   // 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, potential_);
+  const MagneticTowerObj mt = MagneticTowerObj(1, alpha_, l_scale_, offset_, thickness_,
+                                               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,8 +217,8 @@ 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, potential_);
+  const MagneticTowerObj mt(initial_field_, alpha_, l_scale_, offset_, thickness_, 0,
+                            l_mass_scale_, jet_coords, potential_);
 
   parthenon::par_for(
       DEFAULT_LOOP_PATTERN, "MagneticTower::AddInitialFieldToPotential",

diff --git a/src/pgen/cluster/magnetic_tower.hpp b/src/pgen/cluster/magnetic_tower.hpp
@@ -31,6 +31,7 @@ class MagneticTowerObj {
  private:
   const parthenon::Real field_;
   const parthenon::Real alpha_, l_scale_;
+  const parthenon::Real offset_, thickness_;
 
   const parthenon::Real density_, l_mass_scale2_;
 
@@ -42,12 +43,13 @@ class MagneticTowerObj {
 
  public:
   MagneticTowerObj(const parthenon::Real field, const parthenon::Real alpha,
-                   const parthenon::Real l_scale, const parthenon::Real density,
+                   const parthenon::Real l_scale, const parthenon::Real offset,
+                   const parthenon::Real thickness, const parthenon::Real density,
                    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),
-        potential_(potential) {
+      : field_(field), alpha_(alpha), l_scale_(l_scale), offset_(offset),
+        thickness_(thickness), density_(density), 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(
@@ -65,7 +67,7 @@ class MagneticTowerObj {
       // 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) {
+      if (fabs(h) >= 0.001 && fabs(h) <= offset_ + thickness_) {
         a_h = field_ * l_scale_ * exp_r2_h2;
       } else {
         a_h = 0.0;
@@ -108,7 +110,7 @@ class MagneticTowerObj {
       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) {
+      if (fabs(h) >= 0.001 && fabs(h) <= offset_ + thickness_) {
         b_theta = 2.0 * field_ * r / l_scale_ * exp_r2_h2;
       } else {
         b_theta = 0.0;
@@ -163,6 +165,7 @@ class MagneticTowerObj {
 class MagneticTower {
  public:
   const parthenon::Real alpha_, l_scale_;
+  const parthenon::Real offset_, thickness_;
 
   const parthenon::Real initial_field_;
   const parthenon::Real fixed_field_rate_;
@@ -176,6 +179,8 @@ class MagneticTower {
                 const std::string &block = "problem/cluster/magnetic_tower")
       : alpha_(pin->GetOrAddReal(block, "alpha", 0)),
         l_scale_(pin->GetOrAddReal(block, "l_scale", 0)),
+        offset_(pin->GetOrAddReal(block, "offset", 0)),
+        thickness_(pin->GetOrAddReal(block, "thickness", 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)),
@@ -188,8 +193,14 @@ class MagneticTower {
 
     if (potential_str == "donut") {
       potential_ = MagneticTowerPotential::donut;
+      PARTHENON_REQUIRE_THROWS(
+          offset_ >= 0.0 && thickness_ > 0.0,
+          "Incompatible combination of offset and thickness for magnetic donut feedback.")
     } else if (potential_str == "li") {
       potential_ = MagneticTowerPotential::li;
+      PARTHENON_REQUIRE_THROWS(offset_ <= 0.0 && thickness_ <= 0.0,
+                               "Please disable (set to zero) tower offset and thickness "
+                               "for the Li tower model");
     } else {
       PARTHENON_FAIL(
           "Unknown potential for magnetic tower. Current options are: donut, li")