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/docs/pgen.md b/docs/pgen.md index f2ff8f53..0b58b9a2 100644 --- a/docs/pgen.md +++ b/docs/pgen.md @@ -1,11 +1,73 @@ -# Problem specific callbacks +# Problem generators -## Source functions +Different simulation setups in AthenaPK are controlled via the so-called problem generators. +New problem generators can easily be added and we are happy to accept and merge contibuted problem +generators by any user via pull requests. + +## Addding a new problem generator + +In general, four small steps are required: + +### 1. Add a new source file to the `src/pgen` folder + +The file shoud includes at least the `void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md)` +function, which is used to initialize the data at the beginning of the simulation. +Alternatively, the `MeshBlock` version (`void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin)`) can be used that +operates on a single block at a time rather than on a collection of blocks -- though this is not recommended from a performance point of view. + +The function (and all other functions in that file) should be encapsulated in their own namespace (that, ideally, is named +identical to the file itself) so that the functions name are unique across different problem generators. + +Tip: Do not write the problem generator file from scratch but mix and match from existing ones, +e.g., start with a simple one like [orszag_tang.cpp](../src/pgen/orszag_tang.cpp). + +### 2. Add new function(s) to `pgen.hpp` + +All callback functions, i.e., at least the `ProblemGenerator` plus additional optional ones (see step 5 below), +need to be added to [pgen.hpp](../src/pgen/pgen.hpp). +Again, just follow the existing pattern in the file and add the new function declarations with the appropriate namespace. + +### 3. Add callbacks to `main.cpp` + +All problem specific callback functions need to be enrolled in the [`main.cpp`](../src/main.cpp) file. +The selection (via the input file) is controlled by the `problem_id` string in the `` block. +Again, for consistency it is recommended to pick a string that matches the namespace and problem generator file. + +### 4. Ensure new problem generator is compiled + +Add the new source file to [src/pgen/CMakeLists.txt](../src/pgen/CMakeLists.txt) so that it will be compiled +along all other problem generators. + +### 5. (Optional) Add more additional callback + +In addition to the `ProblemGenerator` that initializes the data, other callback functions exists +that allow to modify the behavior of AthenaPK on a problem specific basis. +See [Callback functions](#Callback-functions)] below for available options. + +### 6. (Optional but most likely required) Write an input file + +In theory, one can hardcode all paramters in the source file (like in the +[orszag_tang.cpp](../src/pgen/orszag_tang.cpp) problem generator) but it +prevents modification of the problem setup once the binary is compiled. + +The more common usecase is to create an input file that contains a problem specific +input block. +The convention here is to have a block named `` where `NAME` is the name +of the problem generator (or namespace used). +For example, the Sod shock tube problem generator processes the input file with lines like +``` + Real rho_l = pin->GetOrAddReal("problem/sod", "rho_l", 1.0); +``` +to set the density on the left hand side of the domain. + +## Callback functions + +### Source functions Additional (physical) source terms (e.g., the ones typically on the right hand side of system of equations) can be added in various ways from an algorithmic point of view. -### Unsplit +#### Unsplit Unsplit sources are added at each stage of the (multi-stage) integration after the (conserved) fluxes are calculated. @@ -23,7 +85,7 @@ by implementing a function with that signature and assigning it to the Note, there is no requirement to call the function `MyUnsplitSource`. -### Split first order (generally not recommended) +#### Split first order (generally not recommended) If for special circumstances sources need to be added in a fully operator split way, i.e., the source function is only called *after* the full hydro (or MHD) integration, @@ -38,12 +100,27 @@ Note, as the name suggests, this method is only first order accurate (while ther is no requirement to call the function `MyFirstOrderSource`). -### Split second order (Strang splitting) +#### Split second order (Strang splitting) -Not implemented yet. +Strang splitting achieves second order by interleaving the main hydro/MHD update +with a source update. +In practice, AthenaPK calls Strang split source with `0.5*dt` before the first stage of +each integrator and with `0.5*dt` after the last stage of the integrator. +Note that for consistency, a Strang split source terms should update both conserved +and primitive variables in all zones (i.e., also in the ghost zones as those +are currently not updated after calling a source function). + +Users can enroll a custom source function +```c++ +void MyStrangSplitSource(MeshData *md, const Real beta_dt); +``` +by implementing a function with that signature and assigning it to the +`Hydro::ProblemSourceStrangSplit` callback (currently in `main.cpp`). + +Note, there is no requirement to call the function `MyStrangSplitSource`. -## Timestep restrictions +### Timestep restrictions If additional problem specific physics are implemented (e.g., through the source functions above) that require a custom timestep, it can be added via the @@ -55,4 +132,60 @@ Real ProblemEstimateTimestep(MeshData *md); The return value is expected to be the minimum value over all blocks in the contained in the `MeshData` container, cf., the hydro/mhd `EstimateTimestep` function. Note, that the hyperbolic CFL value is currently applied after the function call, i.e., -it is also applied to the problem specific function. \ No newline at end of file +it is also applied to the problem specific function. + +### Additional initialization on startup (adding variables/fields, custom history output, ...) + +Sometimes a problem generator requires a more general processing of the input file +and/or needs to make certain global adjustments (e.g., adding a custom callback function +for the history output or adding a new field). +This can be achieved via modifying the AthenaPK "package" (currently called +`parthenon::StateDescriptor`) through the following function. +```c++ +void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg) +``` + +For example, the [src/pgen/turbulence.cpp](../[src/pgen/turbulence.cpp]) problem generator +add an additional field (here for an acceleration field) by calling +```c++ + Metadata m({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, + std::vector({3})); + pkg->AddField("acc", m); +``` +in the `ProblemInitPackageData`. + +### Additional initialization on mesh creation/remeshing/load balancing + +For some problem generators it is required to initialize data on "new" blocks. +These new blocks can, for example, be created during mesh refinement +(or derefinement) and this data is typically not active data (like +conserved or primitive variables as those are handled automatically) +but more general data. +Once example is the phase information in the turbulence driver that +does not vary over time but spatially (and therefore at a per block level). + +The appropriate callback to enroll is +```c++ +void InitMeshBlockUserData(MeshBlock *pmb, ParameterInput *pin) +``` + +### UserWorkAfterLoop + +If additional work is required once the main loop finishes (i.e., once the +simulation reaches its final time) computations can be done inside the +```c++ +void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm) +``` +callback function. +This is, for example, done in the linear wave problem generator to calculate the +error norms for convergence tests. + +### MeshBlockUserWorkBeforeOutput + +Sometimes it is desirable to further process data before an output file is written. +For this the +```c++ +void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) +``` +callback is available, that is called once every time right before a data output +(hdf5 or restart) is being written. 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 06ee3df4..386aff0d 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 38c149ea..ca1916e4 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 @@ -31,6 +31,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 @@ -47,13 +49,14 @@ // 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; @@ -61,6 +64,7 @@ using namespace parthenon::package::prelude; using utils::few_modes_ft::FewModesFT; using utils::few_modes_ft_log::FewModesFTLog; +<<<<<<< HEAD template void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, const Real beta_dt, const EOS eos) { @@ -110,7 +114,7 @@ void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, const int gjs = 0; const int gis = 0; - eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i, gks, gjs, gis); // Three last parameters are just passive here + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); // Three last parameters are just passive here if (dfloor > 0) { const Real rho = prim(IDN, k, j, i); @@ -182,8 +186,12 @@ void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, } } -void ClusterSrcTerm(MeshData *md, const parthenon::SimTime &tm, - const Real beta_dt) { +//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"); @@ -203,9 +211,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(); @@ -354,6 +369,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) ************************************************************/ @@ -388,6 +409,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 @@ -460,7 +538,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 = @@ -492,8 +583,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 f12af3e9..4ddfc27e 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; @@ -338,7 +360,7 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData *md, // momentum, and total energy added depend on the triggered power. /////////////////////////////////////////////////////////////////// - eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i,0,0,0); + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); const Real old_specific_internal_e = prim(IPR, k, j, i) / (prim(IDN, k, j, i) * (eos.GetGamma() - 1.)); @@ -348,7 +370,16 @@ 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; - eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i,0,0,0); + // 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.)); PARTHENON_REQUIRE( @@ -360,7 +391,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,13 +407,13 @@ 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; } } - eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i,0,0,0); + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); PARTHENON_REQUIRE(prim(IPR, k, j, i) > 0, "Kinetic injection leads to negative pressure"); }); @@ -391,7 +422,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}",