diff --git a/docs/cluster.md b/docs/cluster.md index 6787d651..b4d83983 100644 --- a/docs/cluster.md +++ b/docs/cluster.md @@ -183,14 +183,16 @@ Pertubations are controlled by the following parameters: # 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 ``` @@ -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 +``` + +enable_tracer = true # disabled by default +``` + +At the moment, this also requires to enable a single passive scalar by setting +``` + +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 +``` + +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 ``` -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 +``` + +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. + + +``` + +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. ``` 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 ``. + +## 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). + +``` + +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). \ No newline at end of file diff --git a/external/parthenon b/external/parthenon index 105c7d0e..5b6bb619 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit 105c7d0eba70849f3496f3e0630a6b61bf8ce4c1 +Subproject commit 5b6bb61906f7c278f9724ee9f38e79dee8707098 diff --git a/inputs/cluster/cluster.in b/inputs/cluster/cluster.in index 61d7bc93..e3192524 100644 --- a/inputs/cluster/cluster.in +++ b/inputs/cluster/cluster.in @@ -181,8 +181,8 @@ kinetic_jet_thickness = 0.05 kinetic_jet_offset = 0.05 - -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 diff --git a/inputs/cluster/magnetic_tower.in b/inputs/cluster/magnetic_tower.in index 66564dca..44ba58a0 100644 --- a/inputs/cluster/magnetic_tower.in +++ b/inputs/cluster/magnetic_tower.in @@ -109,7 +109,8 @@ kinetic_fraction = 0 thermal_fraction = 0 -alpha = 20 +potential_type = li +li_alpha = 20 l_scale = 0.01 initial_field = 0.12431560000204142 fixed_field_rate = 12.431560000204144 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 076fbb9f..5d8fd375 100644 --- a/src/CMakeLists.txt +++ b/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 diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index 425f6365..6d0e087c 100644 --- a/src/hydro/hydro.cpp +++ b/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 *md, const SimTime &tm, const Real b if (hydro_pkg->Param("fluid") == Fluid::glmmhd) { hydro_pkg->Param("glmmhd_source")(md, beta_dt); } + const auto &enable_cooling = hydro_pkg->Param("enable_cooling"); + + if (enable_cooling == Cooling::tabular) { + const TabularCooling &tabular_cooling = + hydro_pkg->Param("tabular_cooling"); + + tabular_cooling.SrcTerm(md, beta_dt); + } if (ProblemSourceUnsplit != nullptr) { ProblemSourceUnsplit(md, tm, beta_dt); } @@ -152,14 +161,6 @@ TaskStatus AddUnsplitSources(MeshData *md, const SimTime &tm, const Real b TaskStatus AddSplitSourcesFirstOrder(MeshData *md, const SimTime &tm) { auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); - const auto &enable_cooling = hydro_pkg->Param("enable_cooling"); - - if (enable_cooling == Cooling::tabular) { - const TabularCooling &tabular_cooling = - hydro_pkg->Param("tabular_cooling"); - - tabular_cooling.SrcTerm(md, tm.dt); - } if (ProblemSourceFirstOrder != nullptr) { ProblemSourceFirstOrder(md, tm, tm.dt); } @@ -167,17 +168,6 @@ TaskStatus AddSplitSourcesFirstOrder(MeshData *md, const SimTime &tm) { } TaskStatus AddSplitSourcesStrang(MeshData *md, const SimTime &tm) { - // auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); - - // const auto &enable_cooling = hydro_pkg->Param("enable_cooling"); - - // if (enable_cooling == Cooling::tabular) { - // const TabularCooling &tabular_cooling = - // hydro_pkg->Param("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 Initialize(ParameterInput *pin) { Metadata m( {Metadata::Cell, Metadata::Independent, Metadata::FillGhost, Metadata::WithFluxes}, std::vector({nhydro + nscalars}), cons_labels); + m.RegisterRefinementOps(); pkg->AddField("cons", m); m = Metadata({Metadata::Cell, Metadata::Derived}, std::vector({nhydro + nscalars}), diff --git a/src/hydro/prolongation/custom_ops.hpp b/src/hydro/prolongation/custom_ops.hpp new file mode 100644 index 00000000..02d2b15f --- /dev/null +++ b/src/hydro/prolongation/custom_ops.hpp @@ -0,0 +1,184 @@ +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2022, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the BSD 3-Clause License (the "LICENSE"). +//======================================================================================== +// Parthenon performance portable AMR framework +// Copyright(C) 2020-2022 The Parthenon collaboration +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +// (C) (or copyright) 2021. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National Nuclear +// Security Administration. All rights in the program are reserved by Triad +// National Security, LLC, and the U.S. Department of Energy/National Nuclear +// Security Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +// in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do +// so. +//======================================================================================== + +#ifndef HYDRO_PROLONGATION_CUSTOM_OPS_HPP_ +#define HYDRO_PROLONGATION_CUSTOM_OPS_HPP_ + +#include +#include + +#include "basic_types.hpp" +#include "coordinates/coordinates.hpp" // for coordinates +#include "kokkos_abstraction.hpp" // ParArray +#include "mesh/domain.hpp" // for IndesShape +#include "prolong_restrict/pr_ops.hpp" + +namespace Hydro { +namespace refinement_ops { + +using parthenon::Coordinates_t; +using parthenon::ParArray6D; +using parthenon::TE; +using parthenon::TopologicalElement; + +// Multi-dimensional, limited prolongation: +// Multi-dim stencil corresponds to eq (5) in Stone et al. (2020). +// Limiting based on implementation in AMReX (see +// https://github.com/AMReX-Codes/amrex/blob/735c3513153f1d06f783e64f455816be85fb3602/Src/AmrCore/AMReX_MFInterp_3D_C.H#L89) +// to preserve extrema. +struct ProlongateCellMinModMultiD { + static constexpr bool OperationRequired(TopologicalElement fel, + TopologicalElement cel) { + return fel == cel; + } + template + KOKKOS_FORCEINLINE_FUNCTION static void + Do(const int l, const int m, const int n, const int k, const int j, const int i, + const IndexRange &ckb, const IndexRange &cjb, const IndexRange &cib, + const IndexRange &kb, const IndexRange &jb, const IndexRange &ib, + const Coordinates_t &coords, const Coordinates_t &coarse_coords, + const ParArrayND *pcoarse, + const ParArrayND *pfine) { + using namespace parthenon::refinement_ops::util; + auto &coarse = *pcoarse; + auto &fine = *pfine; + + constexpr int element_idx = static_cast(el) % 3; + + const int fi = (DIM > 0) ? (i - cib.s) * 2 + ib.s : ib.s; + const int fj = (DIM > 1) ? (j - cjb.s) * 2 + jb.s : jb.s; + const int fk = (DIM > 2) ? (k - ckb.s) * 2 + kb.s : kb.s; + + constexpr bool INCLUDE_X1 = + (DIM > 0) && (el == TE::CC || el == TE::F2 || el == TE::F3 || el == TE::E1); + constexpr bool INCLUDE_X2 = + (DIM > 1) && (el == TE::CC || el == TE::F3 || el == TE::F1 || el == TE::E2); + constexpr bool INCLUDE_X3 = + (DIM > 2) && (el == TE::CC || el == TE::F1 || el == TE::F2 || el == TE::E3); + + const Real fc = coarse(element_idx, l, m, n, k, j, i); + + Real dx1fm = 0; + Real dx1fp = 0; + Real gx1c = 0; + if constexpr (INCLUDE_X1) { + Real dx1m, dx1p; + GetGridSpacings<1, el>(coords, coarse_coords, cib, ib, i, fi, &dx1m, &dx1p, &dx1fm, + &dx1fp); + gx1c = GradMinMod(fc, coarse(element_idx, l, m, n, k, j, i - 1), + coarse(element_idx, l, m, n, k, j, i + 1), dx1m, dx1p); + } + + Real dx2fm = 0; + Real dx2fp = 0; + Real gx2c = 0; + if constexpr (INCLUDE_X2) { + Real dx2m, dx2p; + GetGridSpacings<2, el>(coords, coarse_coords, cjb, jb, j, fj, &dx2m, &dx2p, &dx2fm, + &dx2fp); + gx2c = GradMinMod(fc, coarse(element_idx, l, m, n, k, j - 1, i), + coarse(element_idx, l, m, n, k, j + 1, i), dx2m, dx2p); + } + Real dx3fm = 0; + Real dx3fp = 0; + Real gx3c = 0; + if constexpr (INCLUDE_X3) { + Real dx3m, dx3p; + GetGridSpacings<3, el>(coords, coarse_coords, ckb, kb, k, fk, &dx3m, &dx3p, &dx3fm, + &dx3fp); + gx3c = GradMinMod(fc, coarse(element_idx, l, m, n, k - 1, j, i), + coarse(element_idx, l, m, n, k + 1, j, i), dx3m, dx3p); + } + + // Max. expected total difference. (dx#fm/p are positive by construction) + Real dqmax = std::abs(gx1c) * std::max(dx1fm, dx1fp); + int jlim = 0; + int klim = 0; + if constexpr (DIM > 1) { + dqmax += std::abs(gx2c) * std::max(dx2fm, dx2fp); + jlim = 1; + } + if constexpr (DIM > 2) { + dqmax += std::abs(gx3c) * std::max(dx3fm, dx3fp); + klim = 1; + } + // Min/max values of all coarse cells used here + Real qmin = fc; + Real qmax = fc; + for (int koff = -klim; koff <= klim; koff++) { + for (int joff = -jlim; joff <= jlim; joff++) { + for (int ioff = -1; ioff <= 1; ioff++) { + qmin = + std::min(qmin, coarse(element_idx, l, m, n, k + koff, j + joff, i + ioff)); + qmax = + std::max(qmax, coarse(element_idx, l, m, n, k + koff, j + joff, i + ioff)); + } + } + } + + // Scaling factor to limit all slopes simultaneously + Real alpha = 1.0; + if (dqmax * alpha > (qmax - fc)) { + alpha = (qmax - fc) / dqmax; + } + if (dqmax * alpha > (fc - qmin)) { + alpha = (fc - qmin) / dqmax; + } + + // Ensure no new extrema are introduced in multi-D + gx1c *= alpha; + gx2c *= alpha; + gx3c *= alpha; + + // KGF: add the off-centered quantities first to preserve FP symmetry + // JMM: Extraneous quantities are zero + fine(element_idx, l, m, n, fk, fj, fi) = + fc - (gx1c * dx1fm + gx2c * dx2fm + gx3c * dx3fm); + if constexpr (INCLUDE_X1) + fine(element_idx, l, m, n, fk, fj, fi + 1) = + fc + (gx1c * dx1fp - gx2c * dx2fm - gx3c * dx3fm); + if constexpr (INCLUDE_X2) + fine(element_idx, l, m, n, fk, fj + 1, fi) = + fc - (gx1c * dx1fm - gx2c * dx2fp + gx3c * dx3fm); + if constexpr (INCLUDE_X2 && INCLUDE_X1) + fine(element_idx, l, m, n, fk, fj + 1, fi + 1) = + fc + (gx1c * dx1fp + gx2c * dx2fp - gx3c * dx3fm); + if constexpr (INCLUDE_X3) + fine(element_idx, l, m, n, fk + 1, fj, fi) = + fc - (gx1c * dx1fm + gx2c * dx2fm - gx3c * dx3fp); + if constexpr (INCLUDE_X3 && INCLUDE_X1) + fine(element_idx, l, m, n, fk + 1, fj, fi + 1) = + fc + (gx1c * dx1fp - gx2c * dx2fm + gx3c * dx3fp); + if constexpr (INCLUDE_X3 && INCLUDE_X2) + fine(element_idx, l, m, n, fk + 1, fj + 1, fi) = + fc - (gx1c * dx1fm - gx2c * dx2fp - gx3c * dx3fp); + if constexpr (INCLUDE_X3 && INCLUDE_X2 && INCLUDE_X1) + fine(element_idx, l, m, n, fk + 1, fj + 1, fi + 1) = + fc + (gx1c * dx1fp + gx2c * dx2fp + gx3c * dx3fp); + } +}; +} // namespace refinement_ops +} // namespace Hydro + +#endif // HYDRO_PROLONGATION_CUSTOM_OPS_HPP_ diff --git a/src/main.cpp b/src/main.cpp index 8b796610..87d82021 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -88,7 +88,8 @@ int main(int argc, char *argv[]) { pman.app_input->MeshProblemGenerator = cluster::ProblemGenerator; pman.app_input->MeshBlockUserWorkBeforeOutput = cluster::UserWorkBeforeOutput; Hydro::ProblemInitPackageData = cluster::ProblemInitPackageData; - Hydro::ProblemSourceUnsplit = cluster::ClusterSrcTerm; + Hydro::ProblemSourceUnsplit = cluster::ClusterUnsplitSrcTerm; + Hydro::ProblemSourceFirstOrder = cluster::ClusterSplitSrcTerm; Hydro::ProblemEstimateTimestep = cluster::ClusterEstimateTimestep; } else if (problem == "sod") { pman.app_input->ProblemGenerator = sod::ProblemGenerator; diff --git a/src/pgen/CMakeLists.txt b/src/pgen/CMakeLists.txt index 443c4465..382b5b7b 100644 --- a/src/pgen/CMakeLists.txt +++ b/src/pgen/CMakeLists.txt @@ -8,9 +8,12 @@ target_sources(athenaPK PRIVATE cluster.cpp cluster/agn_feedback.cpp cluster/agn_triggering.cpp + cluster/cluster_clips.cpp + cluster/cluster_reductions.cpp cluster/hydrostatic_equilibrium_sphere.cpp cluster/magnetic_tower.cpp cluster/snia_feedback.cpp + cluster/stellar_feedback.cpp cpaw.cpp diffusion.cpp field_loop.cpp diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index 1e21f534..a2b10f7e 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -9,7 +9,7 @@ // Setups up an idealized galaxy cluster with an ACCEPT-like entropy profile in // hydrostatic equilbrium with an NFW+BCG+SMBH gravitational profile, // optionally with an initial magnetic tower field. Includes AGN feedback, AGN -// triggering via cold gas, simple SNIA Feedback(TODO) +// triggering via cold gas, simple SNIA Feedback, and simple stellar feedback //======================================================================================== // C headers @@ -28,6 +28,8 @@ #include "kokkos_abstraction.hpp" #include "mesh/domain.hpp" #include "mesh/mesh.hpp" +#include "parthenon_array_generic.hpp" +#include "utils/error_checking.hpp" #include #include @@ -43,137 +45,22 @@ // Cluster headers #include "cluster/agn_feedback.hpp" #include "cluster/agn_triggering.hpp" +#include "cluster/cluster_clips.hpp" #include "cluster/cluster_gravity.hpp" +#include "cluster/cluster_reductions.hpp" #include "cluster/entropy_profiles.hpp" #include "cluster/hydrostatic_equilibrium_sphere.hpp" #include "cluster/magnetic_tower.hpp" #include "cluster/snia_feedback.hpp" -#include "parthenon_array_generic.hpp" -#include "utils/error_checking.hpp" +#include "cluster/stellar_feedback.hpp" namespace cluster { using namespace parthenon::driver::prelude; using namespace parthenon::package::prelude; using utils::few_modes_ft::FewModesFT; -template -void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, - const Real beta_dt, const EOS eos) { - - auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); - - // Apply clips -- ceilings on temperature, velocity, alfven velocity, and - // density floor -- within a radius of the AGN - const auto &dfloor = hydro_pkg->Param("cluster_dfloor"); - const auto &eceil = hydro_pkg->Param("cluster_eceil"); - const auto &vceil = hydro_pkg->Param("cluster_vceil"); - const auto &vAceil = hydro_pkg->Param("cluster_vAceil"); - const auto &clip_r = hydro_pkg->Param("cluster_clip_r"); - - if (clip_r > 0 && (dfloor > 0 || eceil < std::numeric_limits::infinity() || - vceil < std::numeric_limits::infinity() || - vAceil < std::numeric_limits::infinity())) { - // Grab some necessary variables - const auto &prim_pack = md->PackVariables(std::vector{"prim"}); - const auto &cons_pack = md->PackVariables(std::vector{"cons"}); - IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior); - IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior); - IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior); - const auto nhydro = hydro_pkg->Param("nhydro"); - const auto nscalars = hydro_pkg->Param("nscalars"); - - const Real clip_r2 = SQR(clip_r); - const Real vceil2 = SQR(vceil); - const Real vAceil2 = SQR(vAceil); - const Real gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); - - parthenon::par_for( - DEFAULT_LOOP_PATTERN, "Cluster::ApplyClusterClips", parthenon::DevExecSpace(), 0, - cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { - auto &cons = cons_pack(b); - auto &prim = prim_pack(b); - const auto &coords = cons_pack.GetCoords(b); - - const Real r2 = - SQR(coords.Xc<1>(i)) + SQR(coords.Xc<2>(j)) + SQR(coords.Xc<3>(k)); - - if (r2 < clip_r2) { - // Cell falls within clipping radius - eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); - - if (dfloor > 0) { - const Real rho = prim(IDN, k, j, i); - if (rho < dfloor) { - cons(IDN, k, j, i) = dfloor; - prim(IDN, k, j, i) = dfloor; - } - } - - if (vceil < std::numeric_limits::infinity()) { - // Apply velocity ceiling - const Real v2 = SQR(prim(IV1, k, j, i)) + SQR(prim(IV2, k, j, i)) + - SQR(prim(IV3, k, j, i)); - if (v2 > vceil2) { - // Fix the velocity to the velocity ceiling - const Real v = sqrt(v2); - cons(IM1, k, j, i) *= vceil / v; - cons(IM2, k, j, i) *= vceil / v; - cons(IM3, k, j, i) *= vceil / v; - prim(IV1, k, j, i) *= vceil / v; - prim(IV2, k, j, i) *= vceil / v; - prim(IV3, k, j, i) *= vceil / v; - - // Remove kinetic energy - cons(IEN, k, j, i) -= 0.5 * prim(IDN, k, j, i) * (v2 - vceil2); - } - } - - if (vAceil2 < std::numeric_limits::infinity()) { - // Apply Alfven velocity ceiling by raising density - const Real rho = prim(IDN, k, j, i); - const Real B2 = (SQR(prim(IB1, k, j, i)) + SQR(prim(IB2, k, j, i)) + - SQR(prim(IB3, k, j, i))); - - // compute Alfven mach number - const Real va2 = (B2 / rho); - - if (va2 > vAceil2) { - // Increase the density to match the alfven velocity ceiling - const Real rho_new = std::sqrt(B2 / vAceil2); - cons(IDN, k, j, i) = rho_new; - prim(IDN, k, j, i) = rho_new; - } - } - - if (eceil < std::numeric_limits::infinity()) { - // Apply internal energy ceiling as a pressure ceiling - const Real internal_e = prim(IPR, k, j, i) / (gm1 * prim(IDN, k, j, i)); - if (internal_e > eceil) { - cons(IEN, k, j, i) -= prim(IDN, k, j, i) * (internal_e - eceil); - prim(IPR, k, j, i) = gm1 * prim(IDN, k, j, i) * eceil; - } - } - } - }); - } -} - -void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, - const Real beta_dt) { - auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); - auto fluid = hydro_pkg->Param("fluid"); - if (fluid == Fluid::euler) { - ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param("eos")); - } else if (fluid == Fluid::glmmhd) { - ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param("eos")); - } else { - PARTHENON_FAIL("Cluster::ApplyClusterClips: Unknown EOS"); - } -} - -void ClusterSrcTerm(MeshData *md, const parthenon::SimTime &tm, - const Real beta_dt) { +void ClusterUnsplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, + const Real beta_dt) { auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); const bool &gravity_srcterm = hydro_pkg->Param("gravity_srcterm"); @@ -193,9 +80,16 @@ void ClusterSrcTerm(MeshData *md, const parthenon::SimTime &tm, const auto &snia_feedback = hydro_pkg->Param("snia_feedback"); snia_feedback.FeedbackSrcTerm(md, beta_dt, tm); - - ApplyClusterClips(md, tm, beta_dt); }; +void ClusterSplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, + const Real dt) { + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + + const auto &stellar_feedback = hydro_pkg->Param("stellar_feedback"); + stellar_feedback.FeedbackSrcTerm(md, dt, tm); + + ApplyClusterClips(md, tm, dt); +} Real ClusterEstimateTimestep(MeshData *md) { Real min_dt = std::numeric_limits::max(); @@ -341,6 +235,12 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd SNIAFeedback snia_feedback(pin, hydro_pkg); + /************************************************************ + * Read Stellar Feedback + ************************************************************/ + + StellarFeedback stellar_feedback(pin, hydro_pkg); + /************************************************************ * Read Clips (ceilings and floors) ************************************************************/ @@ -375,6 +275,63 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd hydro_pkg->AddParam("cluster_vAceil", vAceil); hydro_pkg->AddParam("cluster_clip_r", clip_r); + /************************************************************ + * Start running reductions into history outputs for clips, stellar mass, cold + * gas, and AGN extent + ************************************************************/ + + /* FIXME(forrestglines) This implementation with a reduction into Params might + be broken in several ways. + 1. Each reduction in params is Rank local. Multiple meshblocks packs per + rank adding to these params is not thread-safe + 2. These Params are not carried over between restarts. If a restart dump is + made and a history output is not, then the mass/energy between the last + history output and the restart dump is lost + */ + std::string reduction_strs[] = {"stellar_mass", "added_dfloor_mass", + "removed_eceil_energy", "removed_vceil_energy", + "added_vAceil_mass"}; + + // Add a param for each reduction, then add it as a summation reduction for + // history outputs + auto hst_vars = hydro_pkg->Param(parthenon::hist_param_key); + + for (auto reduction_str : reduction_strs) { + hydro_pkg->AddParam(reduction_str, 0.0, true); + hst_vars.emplace_back(parthenon::HistoryOutputVar( + parthenon::UserHistoryOperation::sum, + [reduction_str](MeshData *md) { + auto pmb = md->GetBlockData(0)->GetBlockPointer(); + auto hydro_pkg = pmb->packages.Get("Hydro"); + const Real reduction = hydro_pkg->Param(reduction_str); + // Reset the running count for this reduction between history outputs + hydro_pkg->UpdateParam(reduction_str, 0.0); + return reduction; + }, + reduction_str)); + } + + // Add history reduction for total cold gas using stellar mass threshold + const Real cold_thresh = + pin->GetOrAddReal("problem/cluster/reductions", "cold_temp_thresh", 0.0); + if (cold_thresh > 0) { + hydro_pkg->AddParam("reduction_cold_threshold", cold_thresh); + hst_vars.emplace_back(parthenon::HistoryOutputVar( + parthenon::UserHistoryOperation::sum, LocalReduceColdGas, "cold_mass")); + } + const Real agn_tracer_thresh = + pin->GetOrAddReal("problem/cluster/reductions", "agn_tracer_thresh", -1.0); + if (agn_tracer_thresh >= 0) { + PARTHENON_REQUIRE( + pin->GetOrAddBoolean("problem/cluster/agn_feedback", "enable_tracer", false), + "AGN Tracer must be enabled to reduce AGN tracer extent"); + hydro_pkg->AddParam("reduction_agn_tracer_threshold", agn_tracer_thresh); + hst_vars.emplace_back(parthenon::HistoryOutputVar( + parthenon::UserHistoryOperation::max, LocalReduceAGNExtent, "agn_extent")); + } + + hydro_pkg->UpdateParam(parthenon::hist_param_key, hst_vars); + /************************************************************ * Add derived fields * NOTE: these must be filled in UserWorkBeforeOutput @@ -411,7 +368,20 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd const auto sigma_v = pin->GetOrAddReal("problem/cluster/init_perturb", "sigma_v", 0.0); if (sigma_v != 0.0) { // peak of init vel perturb - auto k_peak_v = pin->GetReal("problem/cluster/init_perturb", "k_peak_v"); + auto l_peak_v = pin->GetOrAddReal("problem/cluster/init_perturb", "l_peak_v", -1.0); + auto k_peak_v = pin->GetOrAddReal("problem/cluster/init_perturb", "k_peak_v", -1.0); + + PARTHENON_REQUIRE_THROWS((l_peak_v > 0.0 && k_peak_v <= 0.0) || + (k_peak_v > 0.0 && l_peak_v <= 0.0), + "Setting initial velocity perturbation requires a single " + "length scale by either setting l_peak_v or k_peak_v."); + // Set peak wavemode as required by few_modes_fft when not directly given + if (l_peak_v > 0) { + const auto Lx = pin->GetReal("parthenon/mesh", "x1max") - + pin->GetReal("parthenon/mesh", "x1min"); + // Note that this assumes a cubic box + k_peak_v = Lx / l_peak_v; + } auto num_modes_v = pin->GetOrAddInteger("problem/cluster/init_perturb", "num_modes_v", 40); auto sol_weight_v = @@ -438,8 +408,20 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd PARTHENON_REQUIRE_THROWS(hydro_pkg->Param("fluid") == Fluid::glmmhd, "Requested initial magnetic field perturbations but not " "solving the MHD equations.") - // peak of init vel perturb - auto k_peak_b = pin->GetReal("problem/cluster/init_perturb", "k_peak_b"); + // peak of init magnetic field perturb + auto l_peak_b = pin->GetOrAddReal("problem/cluster/init_perturb", "l_peak_b", -1.0); + auto k_peak_b = pin->GetOrAddReal("problem/cluster/init_perturb", "k_peak_b", -1.0); + PARTHENON_REQUIRE_THROWS((l_peak_b > 0.0 && k_peak_b <= 0.0) || + (k_peak_b > 0.0 && l_peak_b <= 0.0), + "Setting initial B perturbation requires a single " + "length scale by either setting l_peak_b or k_peak_b."); + // Set peak wavemode as required by few_modes_fft when not directly given + if (l_peak_b > 0) { + const auto Lx = pin->GetReal("parthenon/mesh", "x1max") - + pin->GetReal("parthenon/mesh", "x1min"); + // Note that this assumes a cubic box + k_peak_b = Lx / l_peak_b; + } auto num_modes_b = pin->GetOrAddInteger("problem/cluster/init_perturb", "num_modes_b", 40); uint32_t rseed_b = pin->GetOrAddInteger("problem/cluster/init_perturb", "rseed_b", 2); diff --git a/src/pgen/cluster/agn_feedback.cpp b/src/pgen/cluster/agn_feedback.cpp index 1c76e913..bcf63c24 100644 --- a/src/pgen/cluster/agn_feedback.cpp +++ b/src/pgen/cluster/agn_feedback.cpp @@ -26,6 +26,7 @@ #include "agn_triggering.hpp" #include "cluster_utils.hpp" #include "magnetic_tower.hpp" +#include "utils/error_checking.hpp" namespace cluster { using namespace parthenon; @@ -50,7 +51,11 @@ AGNFeedback::AGNFeedback(parthenon::ParameterInput *pin, "kinetic_jet_thickness", 0.02)), kinetic_jet_offset_( pin->GetOrAddReal("problem/cluster/agn_feedback", "kinetic_jet_offset", 0.02)), - disabled_(pin->GetOrAddBoolean("problem/cluster/agn_feedback", "disabled", false)) { + enable_tracer_( + pin->GetOrAddBoolean("problem/cluster/agn_feedback", "enable_tracer", false)), + disabled_(pin->GetOrAddBoolean("problem/cluster/agn_feedback", "disabled", false)), + enable_magnetic_tower_mass_injection_(pin->GetOrAddBoolean( + "problem/cluster/agn_feedback", "enable_magnetic_tower_mass_injection", true)) { // Normalize the thermal, kinetic, and magnetic fractions to sum to 1.0 const Real total_frac = thermal_fraction_ + kinetic_fraction_ + magnetic_fraction_; @@ -64,6 +69,18 @@ AGNFeedback::AGNFeedback(parthenon::ParameterInput *pin, magnetic_fraction_ >= 0, "AGN feedback energy fractions must be non-negative."); + // Normalize the thermal, kinetic, and magnetic mass fractions to sum to 1.0 + if (enable_magnetic_tower_mass_injection_) { + thermal_mass_fraction_ = thermal_fraction_; + kinetic_mass_fraction_ = kinetic_fraction_; + magnetic_mass_fraction_ = magnetic_fraction_; + } else { + const auto total_mass_frac = thermal_fraction_ + kinetic_fraction_; + thermal_mass_fraction_ = thermal_fraction_ / total_mass_frac; + kinetic_mass_fraction_ = kinetic_fraction_ / total_mass_frac; + magnetic_mass_fraction_ = 0.0; + } + ///////////////////////////////////////////////////// // Read in or calculate jet velocity and temperature. Either and or both can // be defined but they must satify @@ -160,6 +177,10 @@ AGNFeedback::AGNFeedback(parthenon::ParameterInput *pin, } hydro_pkg->UpdateParam(parthenon::hist_param_key, hst_vars); + // Double check that tracers are also enabled in fluid solver + PARTHENON_REQUIRE_THROWS(!enable_tracer_ || hydro_pkg->Param("nscalars") == 1, + "Enabling tracer for AGN feedback requires hydro/nscalars=1"); + hydro_pkg->AddParam<>("agn_feedback", *this); } @@ -247,7 +268,7 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData *md, thermal_fraction_ * power * thermal_scaling_factor * beta_dt; // Amount of density to dump in each cell const Real thermal_density = - thermal_fraction_ * mass_rate * thermal_scaling_factor * beta_dt; + thermal_mass_fraction_ * mass_rate * thermal_scaling_factor * beta_dt; //////////////////////////////////////////////////////////////////////////////// // Kinetic Jet Quantities @@ -265,7 +286,7 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData *md, // Amount of density to dump in each cell const Real jet_density = - kinetic_fraction_ * mass_rate * kinetic_scaling_factor * beta_dt; + kinetic_mass_fraction_ * mass_rate * kinetic_scaling_factor * beta_dt; // Velocity of added gas const Real jet_velocity = kinetic_jet_velocity_; @@ -281,6 +302,7 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData *md, const Real vceil2 = SQR(vceil); const Real eceil = eceil_; const Real gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); + const auto enable_tracer = enable_tracer_; //////////////////////////////////////////////////////////////////////////////// const parthenon::Real time = tm.time; @@ -348,6 +370,14 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData *md, cons(IM3, k, j, i) += jet_momentum * sign_jet * jet_axis_z; cons(IEN, k, j, i) += jet_feedback; + // Reset tracer to one for the entire material in the jet launching region as + // we cannot distinguish between original material in a cell and new jet + // material in the evolution of the jet. Eventually, we're just interested in + // stuff that came from here. + if (enable_tracer) { + cons(nhydro, k, j, i) = 1.0 * cons(IDN, k, j, i); + } + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); const Real new_specific_internal_e = prim(IPR, k, j, i) / (prim(IDN, k, j, i) * (eos.GetGamma() - 1.)); @@ -360,7 +390,7 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData *md, // Apply velocity ceiling const Real v2 = SQR(prim(IV1, k, j, i)) + SQR(prim(IV2, k, j, i)) + SQR(prim(IV3, k, j, i)); - if (v2 > vceil2) { + if (vceil2 > 0 && v2 > vceil2) { // Fix the velocity to the velocity ceiling const Real v = sqrt(v2); cons(IM1, k, j, i) *= vceil / v; @@ -376,7 +406,7 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData *md, // Apply internal energy ceiling as a pressure ceiling const Real internal_e = prim(IPR, k, j, i) / (gm1 * prim(IDN, k, j, i)); - if (internal_e > eceil) { + if (eceil > 0 && internal_e > eceil) { cons(IEN, k, j, i) -= prim(IDN, k, j, i) * (internal_e - eceil); prim(IPR, k, j, i) = gm1 * prim(IDN, k, j, i) * eceil; } @@ -391,7 +421,7 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData *md, const auto &magnetic_tower = hydro_pkg->Param("magnetic_tower"); const Real magnetic_power = power * magnetic_fraction_; - const Real magnetic_mass_rate = mass_rate * magnetic_fraction_; + const Real magnetic_mass_rate = mass_rate * magnetic_mass_fraction_; magnetic_tower.PowerSrcTerm(magnetic_power, magnetic_mass_rate, md, beta_dt, tm); } diff --git a/src/pgen/cluster/agn_feedback.hpp b/src/pgen/cluster/agn_feedback.hpp index de516949..3bec1609 100644 --- a/src/pgen/cluster/agn_feedback.hpp +++ b/src/pgen/cluster/agn_feedback.hpp @@ -27,6 +27,7 @@ class AGNFeedback { public: const parthenon::Real fixed_power_; parthenon::Real thermal_fraction_, kinetic_fraction_, magnetic_fraction_; + parthenon::Real thermal_mass_fraction_, kinetic_mass_fraction_, magnetic_mass_fraction_; // Efficiency converting mass to energy const parthenon::Real efficiency_; @@ -41,8 +42,13 @@ class AGNFeedback { const parthenon::Real kinetic_jet_radius_, kinetic_jet_thickness_, kinetic_jet_offset_; parthenon::Real kinetic_jet_velocity_, kinetic_jet_temperature_, kinetic_jet_e_; + // enable passive scalar to trace AGN material + const bool enable_tracer_; + const bool disabled_; + const bool enable_magnetic_tower_mass_injection_; + AGNFeedback(parthenon::ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg); parthenon::Real GetFeedbackPower(parthenon::StateDescriptor *hydro_pkg) const; diff --git a/src/pgen/cluster/cluster_clips.cpp b/src/pgen/cluster/cluster_clips.cpp new file mode 100644 index 00000000..5aa86eea --- /dev/null +++ b/src/pgen/cluster/cluster_clips.cpp @@ -0,0 +1,174 @@ + +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2021-2023, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//! \file agn_triggering.cpp +// \brief Class for computing AGN triggering from Bondi-like and cold gas accretion + +// Parthenon headers +#include "kokkos_abstraction.hpp" +#include "mesh/domain.hpp" +#include "mesh/mesh.hpp" +#include "parthenon_array_generic.hpp" +#include "utils/error_checking.hpp" +#include + +// AthenaPK headers +#include "../../eos/adiabatic_glmmhd.hpp" +#include "../../eos/adiabatic_hydro.hpp" + +namespace cluster { +using namespace parthenon; + +template +void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, + const Real beta_dt, const EOS eos); + +void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, + const Real beta_dt) { + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + auto fluid = hydro_pkg->Param("fluid"); + if (fluid == Fluid::euler) { + ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param("eos")); + } else if (fluid == Fluid::glmmhd) { + ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param("eos")); + } else { + PARTHENON_FAIL("Cluster::ApplyClusterClips: Unknown EOS"); + } +} + +template +void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, + const Real beta_dt, const EOS eos) { + + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + + // Apply clips -- ceilings on temperature, velocity, alfven velocity, and + // density floor -- within a radius of the AGN + const auto &dfloor = hydro_pkg->Param("cluster_dfloor"); + const auto &eceil = hydro_pkg->Param("cluster_eceil"); + const auto &vceil = hydro_pkg->Param("cluster_vceil"); + const auto &vAceil = hydro_pkg->Param("cluster_vAceil"); + const auto &clip_r = hydro_pkg->Param("cluster_clip_r"); + + if (clip_r > 0 && (dfloor > 0 || eceil < std::numeric_limits::infinity() || + vceil < std::numeric_limits::infinity() || + vAceil < std::numeric_limits::infinity())) { + // Grab some necessary variables + const auto &prim_pack = md->PackVariables(std::vector{"prim"}); + const auto &cons_pack = md->PackVariables(std::vector{"cons"}); + IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior); + IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior); + IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior); + const auto nhydro = hydro_pkg->Param("nhydro"); + const auto nscalars = hydro_pkg->Param("nscalars"); + + const Real clip_r2 = SQR(clip_r); + const Real vceil2 = SQR(vceil); + const Real vAceil2 = SQR(vAceil); + const Real gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); + + Real added_dfloor_mass = 0.0, removed_vceil_energy = 0.0, added_vAceil_mass = 0.0, + removed_eceil_energy = 0.0; + + Kokkos::parallel_reduce( + "Cluster::ApplyClusterClips", + Kokkos::MDRangePolicy>( + DevExecSpace(), {0, kb.s, jb.s, ib.s}, + {prim_pack.GetDim(5), kb.e + 1, jb.e + 1, ib.e + 1}, + {1, 1, 1, ib.e + 1 - ib.s}), + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i, + Real &added_dfloor_mass_team, Real &removed_vceil_energy_team, + Real &added_vAceil_mass_team, Real &removed_eceil_energy_team) { + auto &cons = cons_pack(b); + auto &prim = prim_pack(b); + const auto &coords = cons_pack.GetCoords(b); + + const Real r2 = + SQR(coords.Xc<1>(i)) + SQR(coords.Xc<2>(j)) + SQR(coords.Xc<3>(k)); + + if (r2 < clip_r2) { + // Cell falls within clipping radius + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); + + if (dfloor > 0) { + const Real rho = prim(IDN, k, j, i); + if (rho < dfloor) { + added_dfloor_mass_team += (dfloor - rho) * coords.CellVolume(k, j, i); + cons(IDN, k, j, i) = dfloor; + prim(IDN, k, j, i) = dfloor; + } + } + + if (vceil < std::numeric_limits::infinity()) { + // Apply velocity ceiling + const Real v2 = SQR(prim(IV1, k, j, i)) + SQR(prim(IV2, k, j, i)) + + SQR(prim(IV3, k, j, i)); + if (v2 > vceil2) { + // Fix the velocity to the velocity ceiling + const Real v = sqrt(v2); + cons(IM1, k, j, i) *= vceil / v; + cons(IM2, k, j, i) *= vceil / v; + cons(IM3, k, j, i) *= vceil / v; + prim(IV1, k, j, i) *= vceil / v; + prim(IV2, k, j, i) *= vceil / v; + prim(IV3, k, j, i) *= vceil / v; + + // Remove kinetic energy + const Real removed_energy = 0.5 * prim(IDN, k, j, i) * (v2 - vceil2); + removed_vceil_energy_team += removed_energy * coords.CellVolume(k, j, i); + cons(IEN, k, j, i) -= removed_energy; + } + } + + if (vAceil2 < std::numeric_limits::infinity()) { + // Apply Alfven velocity ceiling by raising density + const Real rho = prim(IDN, k, j, i); + const Real B2 = (SQR(prim(IB1, k, j, i)) + SQR(prim(IB2, k, j, i)) + + SQR(prim(IB3, k, j, i))); + + // compute Alfven mach number + const Real va2 = (B2 / rho); + + if (va2 > vAceil2) { + // Increase the density to match the alfven velocity ceiling + const Real rho_new = std::sqrt(B2 / vAceil2); + added_vAceil_mass_team += (rho_new - rho) * coords.CellVolume(k, j, i); + cons(IDN, k, j, i) = rho_new; + prim(IDN, k, j, i) = rho_new; + } + } + + if (eceil < std::numeric_limits::infinity()) { + // Apply internal energy ceiling as a pressure ceiling + const Real internal_e = prim(IPR, k, j, i) / (gm1 * prim(IDN, k, j, i)); + if (internal_e > eceil) { + const Real removed_energy = prim(IDN, k, j, i) * (internal_e - eceil); + removed_eceil_energy_team += removed_energy * coords.CellVolume(k, j, i); + cons(IEN, k, j, i) -= removed_energy; + prim(IPR, k, j, i) = gm1 * prim(IDN, k, j, i) * eceil; + } + } + } + }, + added_dfloor_mass, removed_vceil_energy, added_vAceil_mass, removed_eceil_energy); + + // Add the freshly added mass/removed energy to running totals + hydro_pkg->UpdateParam("added_dfloor_mass", + added_dfloor_mass + + hydro_pkg->Param("added_dfloor_mass")); + hydro_pkg->UpdateParam("removed_vceil_energy", + removed_vceil_energy + + hydro_pkg->Param("removed_vceil_energy")); + hydro_pkg->UpdateParam("added_vAceil_mass", + added_vAceil_mass + + hydro_pkg->Param("added_vAceil_mass")); + hydro_pkg->UpdateParam("removed_eceil_energy", + removed_eceil_energy + + hydro_pkg->Param("removed_eceil_energy")); + } +} + +} // namespace cluster \ No newline at end of file diff --git a/src/pgen/cluster/cluster_clips.hpp b/src/pgen/cluster/cluster_clips.hpp new file mode 100644 index 00000000..74351f0a --- /dev/null +++ b/src/pgen/cluster/cluster_clips.hpp @@ -0,0 +1,27 @@ + + +#ifndef CLUSTER_CLUSTER_CLIPS_HPP_ +#define CLUSTER_CLUSTER_CLIPS_HPP_ +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2021-2023, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//! \file cluster_clips.hpp +// \brief Class for applying floors and ceils and reducing removed/added mass/energy + +// C++ headers +#include // sqrt() + +// Parthenon headers +#include "basic_types.hpp" +#include "mesh/mesh.hpp" + +namespace cluster { + +void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, + const Real beta_dt); + +} + +#endif // CLUSTER_AGN_TRIGGERING_HPP_ \ No newline at end of file diff --git a/src/pgen/cluster/cluster_reductions.cpp b/src/pgen/cluster/cluster_reductions.cpp new file mode 100644 index 00000000..9e9edb4b --- /dev/null +++ b/src/pgen/cluster/cluster_reductions.cpp @@ -0,0 +1,101 @@ + +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2021-2023, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//! \file cluster_reductions.cpp +// \brief Cluster-specific reductions to compute the total cold gas and maximum radius +// of AGN feedback + +// Parthenon headers +#include "kokkos_abstraction.hpp" +#include "mesh/domain.hpp" +#include "mesh/mesh.hpp" +#include "parthenon_array_generic.hpp" +#include "utils/error_checking.hpp" +#include + +// AthenaPK headers +#include "../../eos/adiabatic_glmmhd.hpp" +#include "../../eos/adiabatic_hydro.hpp" + +namespace cluster { +using namespace parthenon; + +parthenon::Real LocalReduceColdGas(parthenon::MeshData *md) { + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + + const auto &cold_thresh = hydro_pkg->Param("reduction_cold_threshold"); + auto mbar_over_kb = hydro_pkg->Param("mbar_over_kb"); + const auto e_thresh = + cold_thresh / mbar_over_kb / (hydro_pkg->Param("AdiabaticIndex") - 1.0); + + // Grab some necessary variables + const auto &prim_pack = md->PackVariables(std::vector{"prim"}); + IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior); + IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior); + IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior); + const auto nhydro = hydro_pkg->Param("nhydro"); + + const Real gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); + + Real cold_gas = 0.0; + + Kokkos::parallel_reduce( + "LocalReduceColdGas", + Kokkos::MDRangePolicy>( + DevExecSpace(), {0, kb.s, jb.s, ib.s}, + {prim_pack.GetDim(5), kb.e + 1, jb.e + 1, ib.e + 1}, + {1, 1, 1, ib.e + 1 - ib.s}), + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i, + Real &cold_gas_team) { + auto &prim = prim_pack(b); + const auto &coords = prim_pack.GetCoords(b); + + const Real internal_e = prim(IPR, k, j, i) / (gm1 * prim(IDN, k, j, i)); + if (internal_e < e_thresh) { + cold_gas_team += prim(IDN, k, j, i) * coords.CellVolume(k, j, i); + } + }, + cold_gas); + return cold_gas; +} + +parthenon::Real LocalReduceAGNExtent(parthenon::MeshData *md) { + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + + const auto &tracer_thresh = hydro_pkg->Param("reduction_agn_tracer_threshold"); + + // Grab some necessary variables + const auto &cons_pack = md->PackVariables(std::vector{"cons"}); + IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior); + IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior); + IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior); + const auto nhydro = hydro_pkg->Param("nhydro"); + + Real max_r2 = 0.0; + + Kokkos::parallel_reduce( + "LocalReduceAGNExtent", + Kokkos::MDRangePolicy>( + DevExecSpace(), {0, kb.s, jb.s, ib.s}, + {cons_pack.GetDim(5), kb.e + 1, jb.e + 1, ib.e + 1}, + {1, 1, 1, ib.e + 1 - ib.s}), + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i, + Real &max_r2_team) { + auto &cons = cons_pack(b); + const auto &coords = cons_pack.GetCoords(b); + + const auto r2 = SQR(coords.Xc<1>(k, j, i)) + SQR(coords.Xc<2>(k, j, i)) + + SQR(coords.Xc<3>(k, j, i)); + if (cons(nhydro, k, j, i) / cons(IDN, k, j, i) > tracer_thresh && r2 > max_r2) { + max_r2_team = r2; + } + }, + Kokkos::Max(max_r2)); + + return std::sqrt(max_r2); +} + +} // namespace cluster \ No newline at end of file diff --git a/src/pgen/cluster/cluster_reductions.hpp b/src/pgen/cluster/cluster_reductions.hpp new file mode 100644 index 00000000..4fd98bc2 --- /dev/null +++ b/src/pgen/cluster/cluster_reductions.hpp @@ -0,0 +1,24 @@ +#ifndef CLUSTER_CLUSTER_REDUCTIONS_HPP_ +#define CLUSTER_CLUSTER_REDUCTIONS_HPP_ +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2021-2023, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//! \file cluster_reductions.hpp +// \brief Cluster-specific reductions to compute the total cold gas and maximum radius +// of AGN feedback + +// parthenon headers +#include +#include + +namespace cluster { + +parthenon::Real LocalReduceColdGas(parthenon::MeshData *md); + +parthenon::Real LocalReduceAGNExtent(parthenon::MeshData *md); + +} // namespace cluster + +#endif // CLUSTER_CLUSTER_REDUCTIONS_HPP_ diff --git a/src/pgen/cluster/magnetic_tower.cpp b/src/pgen/cluster/magnetic_tower.cpp index 124f120d..1679ffd7 100644 --- a/src/pgen/cluster/magnetic_tower.cpp +++ b/src/pgen/cluster/magnetic_tower.cpp @@ -21,6 +21,7 @@ #include "../../units.hpp" #include "cluster_utils.hpp" #include "magnetic_tower.hpp" +#include "utils/error_checking.hpp" namespace cluster { using namespace parthenon; @@ -44,35 +45,30 @@ void MagneticTower::AddSrcTerm(parthenon::Real field_to_add, parthenon::Real mas // Grab some necessary variables const auto &prim_pack = md->PackVariables(std::vector{"prim"}); const auto &cons_pack = md->PackVariables(std::vector{"cons"}); + const auto &A_pack = md->PackVariables(std::vector{"magnetic_tower_A"}); + IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior); IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior); IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior); // Scale density_to_add to match mass_to_add when integrated over all space - const Real density_to_add = mass_to_add / (pow(l_mass_scale_, 3) * pow(M_PI, 3. / 2.)); + PARTHENON_REQUIRE_THROWS( + mass_to_add == 0.0 || l_mass_scale_ > 0.0, + "Trying to inject mass with the tower model but 0 mass lengthscale. Either disable " + "tower mass injection or set positive lengthscale."); + const auto density_to_add = + mass_to_add > 0.0 ? mass_to_add / (pow(l_mass_scale_, 3) * pow(M_PI, 3. / 2.)) + : 0.0; 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_, offset_, thickness_, + density_to_add, l_mass_scale_, jet_coords, potential_); const auto &eos = hydro_pkg->Param("eos"); // Construct magnetic vector potential then compute magnetic fields - - // Currently reallocates this vector potential everytime step and constructs - // the potential in a separate kernel. There are two solutions: - // 1. Allocate a dependant variable in the hydro package for scratch - // variables, use to store this potential. Would save time in allocations - // but would still require more DRAM memory and two kernel launches - // 2. Compute the potential (12 needed in all) in the same kernel, - // constructing the derivative without storing the potential (more - // arithmetically intensive, maybe faster) - ParArray5D A( - "magnetic_tower_A", 3, cons_pack.GetDim(5), - md->GetBlockData(0)->GetBlockPointer()->cellbounds.ncellsk(IndexDomain::entire), - md->GetBlockData(0)->GetBlockPointer()->cellbounds.ncellsj(IndexDomain::entire), - md->GetBlockData(0)->GetBlockPointer()->cellbounds.ncellsi(IndexDomain::entire)); IndexRange a_ib = ib; a_ib.s -= 1; a_ib.e += 1; @@ -90,15 +86,16 @@ void MagneticTower::AddSrcTerm(parthenon::Real field_to_add, parthenon::Real mas a_jb.e, a_ib.s, a_ib.e, KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { // Compute and apply potential + auto &A = A_pack(b); const auto &coords = cons_pack.GetCoords(b); Real a_x_, a_y_, a_z_; mt.PotentialInSimCart(coords.Xc<1>(i), coords.Xc<2>(j), coords.Xc<3>(k), a_x_, a_y_, a_z_); - A(0, b, k, j, i) = a_x_; - A(1, b, k, j, i) = a_y_; - A(2, b, k, j, i) = a_z_; + A(0, k, j, i) = a_x_; + A(1, k, j, i) = a_y_; + A(2, k, j, i) = a_z_; }); // Take the curl of the potential and apply the new magnetic field @@ -108,18 +105,19 @@ void MagneticTower::AddSrcTerm(parthenon::Real field_to_add, parthenon::Real mas ib.e, KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { auto &cons = cons_pack(b); auto &prim = prim_pack(b); + auto &A = A_pack(b); const auto &coords = cons_pack.GetCoords(b); // Take the curl of a to compute the magnetic field const Real b_x = - (A(2, b, k, j + 1, i) - A(2, b, k, j - 1, i)) / coords.Dxc<2>(j) / 2.0 - - (A(1, b, k + 1, j, i) - A(1, b, k - 1, j, i)) / coords.Dxc<3>(k) / 2.0; + (A(2, k, j + 1, i) - A(2, k, j - 1, i)) / coords.Dxc<2>(j) / 2.0 - + (A(1, k + 1, j, i) - A(1, k - 1, j, i)) / coords.Dxc<3>(k) / 2.0; const Real b_y = - (A(0, b, k + 1, j, i) - A(0, b, k - 1, j, i)) / coords.Dxc<3>(k) / 2.0 - - (A(2, b, k, j, i + 1) - A(2, b, k, j, i - 1)) / coords.Dxc<1>(i) / 2.0; + (A(0, k + 1, j, i) - A(0, k - 1, j, i)) / coords.Dxc<3>(k) / 2.0 - + (A(2, k, j, i + 1) - A(2, k, j, i - 1)) / coords.Dxc<1>(i) / 2.0; const Real b_z = - (A(1, b, k, j, i + 1) - A(1, b, k, j, i - 1)) / coords.Dxc<1>(i) / 2.0 - - (A(0, b, k, j + 1, i) - A(0, b, k, j - 1, i)) / coords.Dxc<2>(j) / 2.0; + (A(1, k, j, i + 1) - A(1, k, j, i - 1)) / coords.Dxc<1>(i) / 2.0 - + (A(0, k, j + 1, i) - A(0, k, j - 1, i)) / coords.Dxc<2>(j) / 2.0; // Add the magnetic field to the conserved variables cons(IB1, k, j, i) += b_x; @@ -133,8 +131,10 @@ void MagneticTower::AddSrcTerm(parthenon::Real field_to_add, parthenon::Real mas 0.5 * (b_x * b_x + b_y * b_y + b_z * b_z); // Add density - const Real cell_delta_rho = - mt.DensityFromSimCart(coords.Xc<1>(i), coords.Xc<2>(j), coords.Xc<3>(k)); + const auto cell_delta_rho = + density_to_add > 0.0 + ? mt.DensityFromSimCart(coords.Xc<1>(i), coords.Xc<2>(j), coords.Xc<3>(k)) + : 0.0; cons(IDN, k, j, i) += cell_delta_rho; }); } @@ -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); + 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); + 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 d6a987d7..ebcc1d44 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 @@ -28,21 +31,30 @@ 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_; 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) - : field_(field), alpha_(alpha), l_scale_(l_scale), density_(density), - l_mass_scale2_(SQR(l_mass_scale)), jet_coords_(jet_coords) { + 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), 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(l_mass_scale > 0, - "Magnetic Tower Mass Length scale must be strictly 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,11 +62,25 @@ 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) - 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; + 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) <= offset_ + thickness_) { + 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 { + PARTHENON_FAIL("Unknown magnetic tower potential."); + } } // Compute Magnetic Potential in simulation Cartesian coordinates @@ -80,12 +106,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) - 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; + 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) <= offset_ + thickness_) { + 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 { + PARTHENON_FAIL("Unknown magnetic tower potential."); + } } // Compute Magnetic field in Simulation Cartesian coordinates @@ -126,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_; @@ -133,17 +173,49 @@ 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)), + : alpha_(pin->GetOrAddReal(block, "li_alpha", 0)), l_scale_(pin->GetOrAddReal(block, "l_scale", 0)), + offset_(pin->GetOrAddReal(block, "donut_offset", 0)), + thickness_(pin->GetOrAddReal(block, "donut_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)), - 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; + PARTHENON_REQUIRE_THROWS(offset_ >= 0.0 && thickness_ > 0.0, + "Incompatible combination of donut_offset and " + "donut_thickness for magnetic donut feedback.") + PARTHENON_REQUIRE_THROWS(alpha_ == 0.0, + "Please disable (set to zero) tower li_alpha " + "for the donut model"); + } 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"); + } + + // 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 diff --git a/src/pgen/cluster/snia_feedback.hpp b/src/pgen/cluster/snia_feedback.hpp index 307f972e..41fda0ce 100644 --- a/src/pgen/cluster/snia_feedback.hpp +++ b/src/pgen/cluster/snia_feedback.hpp @@ -47,4 +47,4 @@ class SNIAFeedback { } // namespace cluster -#endif // CLUSTER_SNIAE_FEEDBACK_HPP_ +#endif // CLUSTER_SNIA_FEEDBACK_HPP_ diff --git a/src/pgen/cluster/stellar_feedback.cpp b/src/pgen/cluster/stellar_feedback.cpp new file mode 100644 index 00000000..aa130d65 --- /dev/null +++ b/src/pgen/cluster/stellar_feedback.cpp @@ -0,0 +1,176 @@ +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2023, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//! \file stellar_feedback.cpp +// \brief Class for magic heating modeling star formation + +#include + +// Parthenon headers +#include +#include +#include +#include +#include +#include + +// Athena headers +#include "../../eos/adiabatic_glmmhd.hpp" +#include "../../eos/adiabatic_hydro.hpp" +#include "../../main.hpp" +#include "../../units.hpp" +#include "cluster_gravity.hpp" +#include "cluster_utils.hpp" +#include "stellar_feedback.hpp" +#include "utils/error_checking.hpp" + +namespace cluster { +using namespace parthenon; + +StellarFeedback::StellarFeedback(parthenon::ParameterInput *pin, + parthenon::StateDescriptor *hydro_pkg) + : stellar_radius_( + pin->GetOrAddReal("problem/cluster/stellar_feedback", "stellar_radius", 0.0)), + exclusion_radius_( + pin->GetOrAddReal("problem/cluster/stellar_feedback", "exclusion_radius", 0.0)), + efficiency_( + pin->GetOrAddReal("problem/cluster/stellar_feedback", "efficiency", 0.0)), + number_density_threshold_(pin->GetOrAddReal("problem/cluster/stellar_feedback", + "number_density_threshold", 0.0)), + temperatue_threshold_(pin->GetOrAddReal("problem/cluster/stellar_feedback", + "temperature_threshold", 0.0)) { + if (stellar_radius_ == 0.0 && exclusion_radius_ == 0.0 && efficiency_ == 0.0 && + number_density_threshold_ == 0.0 && temperatue_threshold_ == 0.0) { + disabled_ = true; + } else { + disabled_ = false; + } + + if (!disabled_ && exclusion_radius_ == 0.0) { + // If exclusion_radius_ is not specified, use AGN triggering accretion radius. + // If both are zero, the PARTHENON_REQUIRE will fail + exclusion_radius_ = + pin->GetOrAddReal("problem/cluster/agn_triggering", "accretion_radius", 0); + } + + PARTHENON_REQUIRE(disabled_ || + (stellar_radius_ != 0.0 && exclusion_radius_ != 0.0 && + efficiency_ != 0.00 && number_density_threshold_ != 0.0 && + temperatue_threshold_ != 0.0), + "Enabling stellar feedback requires setting all parameters."); + + hydro_pkg->AddParam("stellar_feedback", *this); +} + +void StellarFeedback::FeedbackSrcTerm(parthenon::MeshData *md, + const parthenon::Real beta_dt, + const parthenon::SimTime &tm) const { + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + auto fluid = hydro_pkg->Param("fluid"); + if (fluid == Fluid::euler) { + FeedbackSrcTerm(md, beta_dt, tm, hydro_pkg->Param("eos")); + } else if (fluid == Fluid::glmmhd) { + FeedbackSrcTerm(md, beta_dt, tm, hydro_pkg->Param("eos")); + } else { + PARTHENON_FAIL("StellarFeedback::FeedbackSrcTerm: Unknown EOS"); + } +} +template +void StellarFeedback::FeedbackSrcTerm(parthenon::MeshData *md, + const parthenon::Real beta_dt, + const parthenon::SimTime &tm, + const EOS &eos_) const { + using parthenon::IndexDomain; + using parthenon::IndexRange; + using parthenon::Real; + + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + + if (disabled_) { + // No stellar feedback, return + return; + } + + // Grab some necessary variables + const auto &prim_pack = md->PackVariables(std::vector{"prim"}); + const auto &cons_pack = md->PackVariables(std::vector{"cons"}); + IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior); + IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior); + IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior); + const auto nhydro = hydro_pkg->Param("nhydro"); + const auto nscalars = hydro_pkg->Param("nscalars"); + + // const auto gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); + const auto units = hydro_pkg->Param("units"); + const auto mbar = hydro_pkg->Param("mu") * units.mh(); + const auto mbar_over_kb = hydro_pkg->Param("mbar_over_kb"); + + const auto mass_to_energy = efficiency_ * SQR(units.speed_of_light()); + const auto stellar_radius = stellar_radius_; + const auto exclusion_radius = exclusion_radius_; + const auto temperature_threshold = temperatue_threshold_; + const auto number_density_threshold = number_density_threshold_; + + const auto eos = eos_; + + //////////////////////////////////////////////////////////////////////////////// + + Real stellar_mass = 0.0; + + // Constant volumetric heating, reduce mass removed + Kokkos::parallel_reduce( + "StellarFeedback::FeedbackSrcTerm", + Kokkos::MDRangePolicy>( + DevExecSpace(), {0, kb.s, jb.s, ib.s}, + {prim_pack.GetDim(5), kb.e + 1, jb.e + 1, ib.e + 1}, + {1, 1, 1, ib.e + 1 - ib.s}), + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i, + Real &stellar_mass_team) { + auto &cons = cons_pack(b); + auto &prim = prim_pack(b); + const auto &coords = cons_pack.GetCoords(b); + + const auto x = coords.Xc<1>(i); + const auto y = coords.Xc<2>(j); + const auto z = coords.Xc<3>(k); + + const auto r = sqrt(x * x + y * y + z * z); + if (r > stellar_radius || r <= exclusion_radius) { + return; + } + + auto number_density = prim(IDN, k, j, i) / mbar; + if (number_density < number_density_threshold) { + return; + } + + auto temp = mbar_over_kb * prim(IPR, k, j, i) / prim(IDN, k, j, i); + if (temp > temperature_threshold) { + return; + } + + // All conditions to convert mass to energy are met + const auto cell_delta_rho = number_density_threshold * mbar - prim(IDN, k, j, i); + stellar_mass_team -= cell_delta_rho * coords.CellVolume(k, j, i); + + // First remove density at fixed temperature + AddDensityToConsAtFixedVelTemp(cell_delta_rho, cons, prim, eos.GetGamma(), k, j, + i); + // Then add thermal energy + const auto cell_delta_energy_density = -mass_to_energy * cell_delta_rho; + PARTHENON_REQUIRE( + cell_delta_energy_density > 0.0, + "Sanity check failed. Added thermal energy should be positive."); + cons(IEN, k, j, i) += cell_delta_energy_density; + + // Update prims + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); + }, + stellar_mass); + hydro_pkg->UpdateParam( + "stellar_mass", stellar_mass + hydro_pkg->Param("stellar_mass")); +} + +} // namespace cluster diff --git a/src/pgen/cluster/stellar_feedback.hpp b/src/pgen/cluster/stellar_feedback.hpp new file mode 100644 index 00000000..e3d4b20c --- /dev/null +++ b/src/pgen/cluster/stellar_feedback.hpp @@ -0,0 +1,51 @@ +#ifndef CLUSTER_STELLAR_FEEDBACK_HPP_ +#define CLUSTER_STELLAR_FEEDBACK_HPP_ +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2023, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//! \file stellar_feedback.hpp +// \brief Class for injecting Stellar feedback following BCG density + +// parthenon headers +#include +#include +#include +#include +#include + +#include "jet_coords.hpp" + +namespace cluster { + +/************************************************************ + * StellarFeedback + ************************************************************/ +class StellarFeedback { + private: + // feedback parameters in code units + const parthenon::Real stellar_radius_; // length + parthenon::Real exclusion_radius_; // length + const parthenon::Real efficiency_; // dimless + const parthenon::Real number_density_threshold_; // 1/(length^3) + const parthenon::Real temperatue_threshold_; // K + + bool disabled_; + + public: + StellarFeedback(parthenon::ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg); + + void FeedbackSrcTerm(parthenon::MeshData *md, + const parthenon::Real beta_dt, const parthenon::SimTime &tm) const; + + // Apply stellar feedback following cold gas density above a density threshold + template + void FeedbackSrcTerm(parthenon::MeshData *md, + const parthenon::Real beta_dt, const parthenon::SimTime &tm, + const EOS &eos) const; +}; + +} // namespace cluster + +#endif // CLUSTER_STELLAR_FEEDBACK_HPP_ diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 3050cbb5..182d8497 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -100,7 +100,10 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg void InitUserMeshData(ParameterInput *pin); void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md); void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin); -void ClusterSrcTerm(MeshData *md, const parthenon::SimTime &tm, const Real beta_dt); +void ClusterUnsplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, + const Real beta_dt); +void ClusterSplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, + const Real beta_dt); parthenon::Real ClusterEstimateTimestep(MeshData *md); } // namespace cluster diff --git a/src/utils/few_modes_ft.cpp b/src/utils/few_modes_ft.cpp index f80fd103..f477b211 100644 --- a/src/utils/few_modes_ft.cpp +++ b/src/utils/few_modes_ft.cpp @@ -115,9 +115,9 @@ void FewModesFT::SetPhases(MeshBlock *pmb, ParameterInput *pin) { // below take the logical grid size into account. For example, the local phases at level // 1 should be calculated assuming a grid that is twice as large as the root grid. const auto root_level = pm->GetRootLevel(); - auto gnx1 = pm->mesh_size.nx1 * std::pow(2, pmb->loc.level - root_level); - auto gnx2 = pm->mesh_size.nx2 * std::pow(2, pmb->loc.level - root_level); - auto gnx3 = pm->mesh_size.nx3 * std::pow(2, pmb->loc.level - root_level); + auto gnx1 = pm->mesh_size.nx1 * std::pow(2, pmb->loc.level() - root_level); + auto gnx2 = pm->mesh_size.nx2 * std::pow(2, pmb->loc.level() - root_level); + auto gnx3 = pm->mesh_size.nx3 * std::pow(2, pmb->loc.level() - root_level); // Restriction should also be easily fixed, just need to double check transforms and // volume weighting everywhere @@ -130,9 +130,9 @@ void FewModesFT::SetPhases(MeshBlock *pmb, ParameterInput *pin) { const auto nx2 = pmb->block_size.nx2; const auto nx3 = pmb->block_size.nx3; - const auto gis = pmb->loc.lx1 * pmb->block_size.nx1; - const auto gjs = pmb->loc.lx2 * pmb->block_size.nx2; - const auto gks = pmb->loc.lx3 * pmb->block_size.nx3; + const auto gis = pmb->loc.lx1() * pmb->block_size.nx1; + const auto gjs = pmb->loc.lx2() * pmb->block_size.nx2; + const auto gks = pmb->loc.lx3() * pmb->block_size.nx3; // make local ref to capure in lambda const auto num_modes = num_modes_; diff --git a/tst/regression/test_suites/cluster_magnetic_tower/cluster_magnetic_tower.py b/tst/regression/test_suites/cluster_magnetic_tower/cluster_magnetic_tower.py index 4da1eeb1..c80ef354 100644 --- a/tst/regression/test_suites/cluster_magnetic_tower/cluster_magnetic_tower.py +++ b/tst/regression/test_suites/cluster_magnetic_tower/cluster_magnetic_tower.py @@ -275,7 +275,7 @@ def Prepare(self, parameters, step): f"problem/cluster/agn_feedback/magnetic_fraction=1", f"problem/cluster/agn_feedback/kinetic_fraction=0", f"problem/cluster/agn_feedback/thermal_fraction=0", - f"problem/cluster/magnetic_tower/alpha={self.magnetic_tower_alpha}", + f"problem/cluster/magnetic_tower/li_alpha={self.magnetic_tower_alpha}", f"problem/cluster/magnetic_tower/l_scale={self.magnetic_tower_l_scale.in_units('code_length').v}", f"problem/cluster/magnetic_tower/initial_field={self.initial_magnetic_tower_field.in_units('sqrt(code_mass)/sqrt(code_length)/code_time').v}", f"problem/cluster/magnetic_tower/fixed_field_rate={fixed_field_rate}",