diff --git a/docs/cluster.md b/docs/cluster.md index e5256a7c..d31f1005 100644 --- a/docs/cluster.md +++ b/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 +``` + +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. ``` initial_field = 0.12431560000204142 diff --git a/inputs/cluster/cluster.in b/inputs/cluster/cluster.in index 61d7bc93..fe1ff02b 100644 --- a/inputs/cluster/cluster.in +++ b/inputs/cluster/cluster.in @@ -181,7 +181,7 @@ kinetic_jet_thickness = 0.05 kinetic_jet_offset = 0.05 - +potential_type = li alpha = 20 l_scale = 0.001 initial_field = 0.12431560000204142 diff --git a/inputs/cluster/magnetic_tower.in b/inputs/cluster/magnetic_tower.in index 66564dca..15b5647d 100644 --- a/inputs/cluster/magnetic_tower.in +++ b/inputs/cluster/magnetic_tower.in @@ -109,6 +109,7 @@ kinetic_fraction = 0 thermal_fraction = 0 +potential_type = li alpha = 20 l_scale = 0.01 initial_field = 0.12431560000204142 diff --git a/src/pgen/cluster/magnetic_tower.cpp b/src/pgen/cluster/magnetic_tower.cpp index c2cd1663..07264d3c 100644 --- a/src/pgen/cluster/magnetic_tower.cpp +++ b/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("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("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("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", diff --git a/src/pgen/cluster/magnetic_tower.hpp b/src/pgen/cluster/magnetic_tower.hpp index 258678ac..e8e9a687 100644 --- a/src/pgen/cluster/magnetic_tower.hpp +++ b/src/pgen/cluster/magnetic_tower.hpp @@ -17,8 +17,11 @@ #include #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,17 +35,24 @@ 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 @@ -50,14 +60,24 @@ class MagneticTowerObj { 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,6 +170,8 @@ 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)), @@ -148,17 +179,32 @@ class MagneticTower { 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("magnetic_tower_linear_contrib", 0.0, true); hydro_pkg->AddParam("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({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