Merge pull request #82 from parthenon-hpc-lab/pgrete/donut

Cluster improvements and donut model

docs/cluster.md

@@ -183,14 +183,16 @@ Pertubations are controlled by the following parameters:
 <problem/cluster/init_perturb>
 # for the velocity field
 sigma_v = 0.0      # volume weighted RMS |v| in code velocity; default: 0.0 (disabled)
-k_peak_v = ???     # wavenumber in normalized units where the velocity spectrum peaks. No default value.
+l_peak_v = ???     # lengthscale (in code length units) where the velocity spectrum peaks. No default value.
+k_peak_v = ???     # (exclusive alternative to l_peak_v): wavenumber in normalized units where the velocity spectrum peaks. No default value.
 num_modes_v = 40   # (optional) number of wavemodes in spectral space; default: 40
 sol_weight_v = 1.0 # (optional) power in solenoidal (rotational) modes of the perturbation. Range between 0 (fully compressive) and 1.0 (default, fully solenoidal).
 rseed_v = 1        # (optional) integer seed for RNG for wavenumbers and amplitudes
 
 # for the magnetic field
 sigma_b = 0.0      # volume weighted RMS |B| in code magnetic; default: 0.0 (disabled)
-k_peak_b = ???     # wavenumber in normalized units where the magnetic field spectrum peaks. No default value.
+l_peak_b = ???     # lengthscale (in code length units) where the magnetic field spectrum peaks. No default value.
+k_peak_b = ???     # (exclusive alternative to l_peak_b): wavenumber in normalized units where the magnetic field spectrum peaks. No default value.
 num_modes_b = 40   # (optional) number of wavemodes in spectral space; default: 40
 rseed_b = 2        # (optional) integer seed for RNG for wavenumbers and amplitudes
 ```
@@ -306,6 +308,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 +322,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,7 +387,44 @@ 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.
 
-Magnetic feedback is injected following  ([Hui 2006](doi.org/10.1086/501499))
+#### Tracing jet material
+
+Material launched by the kinetic jet component can be traced by setting
+```
+<problem/cluster/agn_feedback>
+enable_tracer = true   # disabled by default
+```
+
+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
+
+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 
 
 $$
@@ -405,7 +448,8 @@ $$
 The parameters $\alpha$ and $\ell$ can be changed with
 ```
 <problem/cluster/magnetic_tower>
-alpha = 20
+potential_type = li
+li_alpha = 20
 l_scale = 0.001
 ```
 When injected as a fraction of 
@@ -429,9 +473,50 @@ so that the total mass injected matches the accreted mass propotioned to magneti
 l_mass_scale = 0.001
 ```
 
-A magnetic tower can also be inserted at runtime and injected at a fixed
+Mass injection by the tower is enabled by default.
+It can be disabled by setting
+```
+<problem/cluster/agn_feedback>
+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.
+
+
+```
+<problem/cluster/magnetic_tower>
+potential_type = donut
+l_scale = 0.0005         # in code length
+donut_offset = 0.001     # in code length
+donut_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
 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
@@ -454,3 +539,27 @@ where `power_per_bcg_mass` and `mass_rate_per_bcg_mass` is the power and mass
 per time respectively injected per BCG mass at a given radius. This SNIA
 feedback is otherwise fixed in time, spherically symmetric, and dependant on
 the BCG specified in `<problem/cluster/gravity>`.
+
+## Stellar feedback
+
+Cold, dense, potentially magnetically supported gas may accumulate in a small
+disk-like structure around the AGN depending on the specific setup.
+In reality, this gas would form stars.
+
+In absence of star particles (or similar) in the current setup, we use a simple
+formulation to convert some fraction of this gas to thermal energy.
+More specifically, gas within a given radius, above a critical number density
+and below a temperature threshold, will be converted to thermal energy with a
+given efficiency.
+The parameters can be set as follows (numerical values here just for illustration purpose -- by default they are all 0, i.e., stellar feedback is disabled).
+
+```
+<problem/cluster/stellar_feedback>
+stellar_radius = 0.025 # in code length
+efficiency = 5e-6
+number_density_threshold = 2.93799894e+74 # in code_length**-3
+temperature_threshold = 2e4 # in K
+```
+
+Note that all parameters need to be specified explicitly for the feedback to work
+(i.e., no hidden default values).

inputs/cluster/cluster.in

@@ -181,8 +181,8 @@ kinetic_jet_thickness  = 0.05
 kinetic_jet_offset = 0.05
 
 <problem/cluster/magnetic_tower>
-
-alpha = 20
+potential_type = li
+li_alpha = 20
 l_scale = 0.001
 initial_field = 0.12431560000204142
 #NOTE: Change to this line to disable initial magnetic tower

inputs/cluster/magnetic_tower.in

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

src/CMakeLists.txt

@@ -12,6 +12,7 @@ add_executable(
         hydro/hydro_driver.cpp
         hydro/hydro.cpp
         hydro/glmmhd/dedner_source.cpp
+        hydro/prolongation/custom_ops.hpp
         hydro/srcterms/gravitational_field.hpp
         hydro/srcterms/tabular_cooling.hpp
         hydro/srcterms/tabular_cooling.cpp

src/hydro/hydro.cpp

@@ -30,6 +30,7 @@
 #include "glmmhd/glmmhd.hpp"
 #include "hydro.hpp"
 #include "outputs/outputs.hpp"
+#include "prolongation/custom_ops.hpp"
 #include "rsolvers/rsolvers.hpp"
 #include "srcterms/tabular_cooling.hpp"
 #include "utils/error_checking.hpp"
@@ -142,6 +143,14 @@ TaskStatus AddUnsplitSources(MeshData<Real> *md, const SimTime &tm, const Real b
   if (hydro_pkg->Param<Fluid>("fluid") == Fluid::glmmhd) {
     hydro_pkg->Param<GLMMHD::SourceFun_t>("glmmhd_source")(md, beta_dt);
   }
+  const auto &enable_cooling = hydro_pkg->Param<Cooling>("enable_cooling");
+
+  if (enable_cooling == Cooling::tabular) {
+    const TabularCooling &tabular_cooling =
+        hydro_pkg->Param<TabularCooling>("tabular_cooling");
+
+    tabular_cooling.SrcTerm(md, beta_dt);
+  }
   if (ProblemSourceUnsplit != nullptr) {
     ProblemSourceUnsplit(md, tm, beta_dt);
   }
@@ -152,32 +161,13 @@ TaskStatus AddUnsplitSources(MeshData<Real> *md, const SimTime &tm, const Real b
 TaskStatus AddSplitSourcesFirstOrder(MeshData<Real> *md, const SimTime &tm) {
   auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
 
-  const auto &enable_cooling = hydro_pkg->Param<Cooling>("enable_cooling");
-
-  if (enable_cooling == Cooling::tabular) {
-    const TabularCooling &tabular_cooling =
-        hydro_pkg->Param<TabularCooling>("tabular_cooling");
-
-    tabular_cooling.SrcTerm(md, tm.dt);
-  }
   if (ProblemSourceFirstOrder != nullptr) {
     ProblemSourceFirstOrder(md, tm, tm.dt);
   }
   return TaskStatus::complete;
 }
 
 TaskStatus AddSplitSourcesStrang(MeshData<Real> *md, const SimTime &tm) {
-  // auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
-
-  // const auto &enable_cooling = hydro_pkg->Param<Cooling>("enable_cooling");
-
-  // if (enable_cooling == Cooling::tabular) {
-  //   const TabularCooling &tabular_cooling =
-  //       hydro_pkg->Param<TabularCooling>("tabular_cooling");
-
-  //   tabular_cooling.SrcTerm(md, 0.5 * tm.dt);
-  // }
-
   if (ProblemSourceStrangSplit != nullptr) {
     ProblemSourceStrangSplit(md, tm, tm.dt);
   }
@@ -560,6 +550,8 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   Metadata m(
       {Metadata::Cell, Metadata::Independent, Metadata::FillGhost, Metadata::WithFluxes},
       std::vector<int>({nhydro + nscalars}), cons_labels);
+  m.RegisterRefinementOps<refinement_ops::ProlongateCellMinModMultiD,
+                          parthenon::refinement_ops::RestrictAverage>();
   pkg->AddField("cons", m);
 
   m = Metadata({Metadata::Cell, Metadata::Derived}, std::vector<int>({nhydro + nscalars}),