diff --git a/docs/cluster.md b/docs/cluster.md index d31f1005..e1444e3b 100644 --- a/docs/cluster.md +++ b/docs/cluster.md @@ -448,6 +448,7 @@ $$ The parameters $\alpha$ and $\ell$ can be changed with ``` +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. + +``` + +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 index 07264d3c..1679ffd7 100644 --- 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("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("eos"); @@ -164,8 +164,8 @@ void MagneticTower::ReducePowerContribs(parthenon::Real &linear_contrib, hydro_pkg->Param("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("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 index e8e9a687..c32ede67 100644 --- 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")