From bc0d8048f72da64215e42ba3daace0a8274d154c Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Tue, 15 Aug 2023 20:59:10 +0200 Subject: [PATCH 01/26] Shoot donuts into space --- src/pgen/cluster/magnetic_tower.cpp | 6 +++--- src/pgen/cluster/magnetic_tower.hpp | 22 +++++++++++++++------- 2 files changed, 18 insertions(+), 10 deletions(-) diff --git a/src/pgen/cluster/magnetic_tower.cpp b/src/pgen/cluster/magnetic_tower.cpp index 124f120d..2208dca0 100644 --- a/src/pgen/cluster/magnetic_tower.cpp +++ b/src/pgen/cluster/magnetic_tower.cpp @@ -133,9 +133,9 @@ 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)); - cons(IDN, k, j, i) += cell_delta_rho; + // const Real cell_delta_rho = + // mt.DensityFromSimCart(coords.Xc<1>(i), coords.Xc<2>(j), coords.Xc<3>(k)); + // cons(IDN, k, j, i) += cell_delta_rho; }); } diff --git a/src/pgen/cluster/magnetic_tower.hpp b/src/pgen/cluster/magnetic_tower.hpp index d6a987d7..2dd2faf6 100644 --- a/src/pgen/cluster/magnetic_tower.hpp +++ b/src/pgen/cluster/magnetic_tower.hpp @@ -50,11 +50,15 @@ 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)); + 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 = field_ * l_scale_ * (r / l_scale_) * exp_r2_h2; - a_h = field_ * l_scale_ * alpha_ / 2.0 * exp_r2_h2; + a_theta = 0.0; + if (fabs(h) >= 0.001 && fabs(h) <= 0.001 + 0.000390625) { + a_h = field_ * l_scale_ * alpha_ / 2.0 * exp_r2_h2; + } else { + a_h = 0.0; + } } // Compute Magnetic Potential in simulation Cartesian coordinates @@ -81,11 +85,15 @@ class MagneticTowerObj { 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)); + const parthenon::Real exp_r2_h2 = exp(-pow(r / 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; + b_r = 0.0; + if (fabs(h) >= 0.001 && fabs(h) <= 0.001 + 0.000390625) { + b_theta = 2.0 * field_ * r / l_scale_ * exp_r2_h2; + } else { + b_theta = 0.0; + } + b_h = 0.0; } // Compute Magnetic field in Simulation Cartesian coordinates From 72a2b8ec90f5c5d7b217eb9b05450e8caccdb0c6 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 16 Aug 2023 08:59:19 +0200 Subject: [PATCH 02/26] Fix tower prefactor --- src/pgen/cluster/magnetic_tower.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pgen/cluster/magnetic_tower.hpp b/src/pgen/cluster/magnetic_tower.hpp index 2dd2faf6..d63acb55 100644 --- a/src/pgen/cluster/magnetic_tower.hpp +++ b/src/pgen/cluster/magnetic_tower.hpp @@ -55,7 +55,7 @@ class MagneticTowerObj { a_r = 0.0; a_theta = 0.0; if (fabs(h) >= 0.001 && fabs(h) <= 0.001 + 0.000390625) { - a_h = field_ * l_scale_ * alpha_ / 2.0 * exp_r2_h2; + a_h = field_ * l_scale_ * exp_r2_h2; } else { a_h = 0.0; } From ec7a620228be3e0565fefc8e32104b63c9287c24 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 16 Aug 2023 11:22:51 +0200 Subject: [PATCH 03/26] Bump parth --- external/parthenon | 2 +- src/utils/few_modes_ft.cpp | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/external/parthenon b/external/parthenon index 105c7d0e..ab9d1671 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit 105c7d0eba70849f3496f3e0630a6b61bf8ce4c1 +Subproject commit ab9d16714d31e853d000a0757f3487bcb9804023 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_; From 4f1078a7812233f3336dac1e47944932ebd0f72f Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 14 Aug 2023 16:53:04 +0200 Subject: [PATCH 04/26] Fix interface --- src/hydro/prolongation/custom_ops.hpp | 179 ++++++++++++++++++++++++++ 1 file changed, 179 insertions(+) create mode 100644 src/hydro/prolongation/custom_ops.hpp diff --git a/src/hydro/prolongation/custom_ops.hpp b/src/hydro/prolongation/custom_ops.hpp new file mode 100644 index 00000000..a86eaba8 --- /dev/null +++ b/src/hydro/prolongation/custom_ops.hpp @@ -0,0 +1,179 @@ +//======================================================================================== +// 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::TopologicalElement; +using parthenon::TE; + +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_X2) { + 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_ From 3ffecd5171416b36af1d1598fa1cedf171f48b33 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 14 Aug 2023 16:52:36 +0200 Subject: [PATCH 05/26] Add multi-d lim prolong --- src/CMakeLists.txt | 1 + src/hydro/hydro.cpp | 3 +++ src/hydro/prolongation/custom_ops.hpp | 2 +- 3 files changed, 5 insertions(+), 1 deletion(-) 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..95218748 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" @@ -560,6 +561,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 index a86eaba8..9f42b3a4 100644 --- a/src/hydro/prolongation/custom_ops.hpp +++ b/src/hydro/prolongation/custom_ops.hpp @@ -38,8 +38,8 @@ namespace refinement_ops { using parthenon::Coordinates_t; using parthenon::ParArray6D; -using parthenon::TopologicalElement; using parthenon::TE; +using parthenon::TopologicalElement; struct ProlongateCellMinModMultiD { static constexpr bool OperationRequired(TopologicalElement fel, From 8104ec2cea68e7a9b575de317b06debf6f618c85 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Thu, 17 Aug 2023 10:15:08 +0200 Subject: [PATCH 06/26] Add jet tracing --- docs/cluster.md | 29 +++++++++++++++++++++++++++++ src/pgen/cluster/agn_feedback.cpp | 16 ++++++++++++++++ src/pgen/cluster/agn_feedback.hpp | 3 +++ 3 files changed, 48 insertions(+) diff --git a/docs/cluster.md b/docs/cluster.md index 6787d651..5a4c0ffe 100644 --- a/docs/cluster.md +++ b/docs/cluster.md @@ -306,6 +306,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 +320,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,6 +385,31 @@ 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. +#### Tracing jet material + +Material launched by the kinetic jet component can be traced by setting +``` + +enable_tracer = true +``` + +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 + Magnetic feedback is injected following ([Hui 2006](doi.org/10.1086/501499)) where the injected magnetic field follows diff --git a/src/pgen/cluster/agn_feedback.cpp b/src/pgen/cluster/agn_feedback.cpp index 1c76e913..c539d2d3 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,6 +51,8 @@ AGNFeedback::AGNFeedback(parthenon::ParameterInput *pin, "kinetic_jet_thickness", 0.02)), kinetic_jet_offset_( pin->GetOrAddReal("problem/cluster/agn_feedback", "kinetic_jet_offset", 0.02)), + enable_tracer_( + pin->GetOrAddBoolean("problem/cluster/agn_feedback", "enable_tracer", false)), disabled_(pin->GetOrAddBoolean("problem/cluster/agn_feedback", "disabled", false)) { // Normalize the thermal, kinetic, and magnetic fractions to sum to 1.0 @@ -160,6 +163,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); } @@ -281,6 +288,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 +356,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.)); diff --git a/src/pgen/cluster/agn_feedback.hpp b/src/pgen/cluster/agn_feedback.hpp index de516949..64d501bd 100644 --- a/src/pgen/cluster/agn_feedback.hpp +++ b/src/pgen/cluster/agn_feedback.hpp @@ -41,6 +41,9 @@ 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_; AGNFeedback(parthenon::ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg); From d449c773115f47fe8338210d93d1219c4c456b5f Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Thu, 17 Aug 2023 11:57:17 +0200 Subject: [PATCH 07/26] Add magic heating from stellar feedback --- docs/cluster.md | 26 ++++- src/pgen/CMakeLists.txt | 1 + src/pgen/cluster.cpp | 16 ++- src/pgen/cluster/agn_feedback.cpp | 2 +- src/pgen/cluster/snia_feedback.hpp | 2 +- src/pgen/cluster/stellar_feedback.cpp | 155 ++++++++++++++++++++++++++ src/pgen/cluster/stellar_feedback.hpp | 50 +++++++++ 7 files changed, 246 insertions(+), 6 deletions(-) create mode 100644 src/pgen/cluster/stellar_feedback.cpp create mode 100644 src/pgen/cluster/stellar_feedback.hpp diff --git a/docs/cluster.md b/docs/cluster.md index 5a4c0ffe..0356171d 100644 --- a/docs/cluster.md +++ b/docs/cluster.md @@ -390,7 +390,7 @@ with kinetic jet feedback. Material launched by the kinetic jet component can be traced by setting ``` -enable_tracer = true +enable_tracer = true # disabled by default ``` At the moment, this also requires to enable a single passive scalar by setting @@ -483,3 +483,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/src/pgen/CMakeLists.txt b/src/pgen/CMakeLists.txt index 443c4465..02b51c3f 100644 --- a/src/pgen/CMakeLists.txt +++ b/src/pgen/CMakeLists.txt @@ -11,6 +11,7 @@ target_sources(athenaPK PRIVATE 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..c9de63cd 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 @@ -48,8 +50,7 @@ #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; @@ -194,6 +195,9 @@ void ClusterSrcTerm(MeshData *md, const parthenon::SimTime &tm, const auto &snia_feedback = hydro_pkg->Param("snia_feedback"); snia_feedback.FeedbackSrcTerm(md, beta_dt, tm); + const auto &stellar_feedback = hydro_pkg->Param("stellar_feedback"); + stellar_feedback.FeedbackSrcTerm(md, beta_dt, tm); + ApplyClusterClips(md, tm, beta_dt); }; @@ -341,6 +345,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) ************************************************************/ diff --git a/src/pgen/cluster/agn_feedback.cpp b/src/pgen/cluster/agn_feedback.cpp index c539d2d3..3d1569f3 100644 --- a/src/pgen/cluster/agn_feedback.cpp +++ b/src/pgen/cluster/agn_feedback.cpp @@ -164,7 +164,7 @@ 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, + 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); 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..8691cde3 --- /dev/null +++ b/src/pgen/cluster/stellar_feedback.cpp @@ -0,0 +1,155 @@ +//======================================================================================== +// 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)), + 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 && efficiency_ == 0.0 && number_density_threshold_ == 0.0 && + temperatue_threshold_ == 0.0) { + disabled_ = true; + } else { + disabled_ = false; + } + PARTHENON_REQUIRE(disabled_ || (stellar_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 temperature_threshold = temperatue_threshold_; + const auto number_density_threshold = number_density_threshold_; + + const auto eos = eos_; + + //////////////////////////////////////////////////////////////////////////////// + + // Constant volumetric heating + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "StellarFeedback::FeedbackSrcTerm", 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 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) { + 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); + + // 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_mass = -cell_delta_rho * coords.CellVolume(k, j, i); + const auto cell_delta_energy_density = + mass_to_energy * cell_delta_mass / coords.CellVolume(k, j, i); + 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); + }); +} + +} // namespace cluster diff --git a/src/pgen/cluster/stellar_feedback.hpp b/src/pgen/cluster/stellar_feedback.hpp new file mode 100644 index 00000000..d91dff6a --- /dev/null +++ b/src/pgen/cluster/stellar_feedback.hpp @@ -0,0 +1,50 @@ +#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 + 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 the feedback from stellare tied to the BCG density + 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_ From b546bdc6f81d5cb6f973b1ae590e98cca231f5d6 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Fri, 18 Aug 2023 11:32:17 +0200 Subject: [PATCH 08/26] Reuse potential in tower as persistent field. --- src/pgen/cluster/magnetic_tower.cpp | 36 +++++++++++------------------ src/pgen/cluster/magnetic_tower.hpp | 7 ++++++ 2 files changed, 20 insertions(+), 23 deletions(-) diff --git a/src/pgen/cluster/magnetic_tower.cpp b/src/pgen/cluster/magnetic_tower.cpp index 2208dca0..75b78bb8 100644 --- a/src/pgen/cluster/magnetic_tower.cpp +++ b/src/pgen/cluster/magnetic_tower.cpp @@ -44,6 +44,8 @@ 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); @@ -59,20 +61,6 @@ void MagneticTower::AddSrcTerm(parthenon::Real field_to_add, parthenon::Real mas 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 +78,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 +97,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; diff --git a/src/pgen/cluster/magnetic_tower.hpp b/src/pgen/cluster/magnetic_tower.hpp index d63acb55..2a4a06fe 100644 --- a/src/pgen/cluster/magnetic_tower.hpp +++ b/src/pgen/cluster/magnetic_tower.hpp @@ -152,6 +152,13 @@ class MagneticTower { hydro_pkg->AddParam<>("magnetic_tower", *this); hydro_pkg->AddParam("magnetic_tower_linear_contrib", 0.0, true); hydro_pkg->AddParam("magnetic_tower_quadratic_contrib", 0.0, true); + + // 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); } // Add initial magnetic field to provided potential with a single meshblock From 373e0a5914a88db852ec392fcfcc8c7ead5f5341 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Fri, 18 Aug 2023 15:02:45 +0200 Subject: [PATCH 09/26] Allow disabling of mass injection by the tower --- docs/cluster.md | 10 ++++++++++ src/pgen/cluster/agn_feedback.cpp | 22 ++++++++++++++++++---- src/pgen/cluster/agn_feedback.hpp | 3 +++ src/pgen/cluster/magnetic_tower.cpp | 17 +++++++++++++---- src/pgen/cluster/magnetic_tower.hpp | 4 ++-- 5 files changed, 46 insertions(+), 10 deletions(-) diff --git a/docs/cluster.md b/docs/cluster.md index 0356171d..e69771ce 100644 --- a/docs/cluster.md +++ b/docs/cluster.md @@ -458,6 +458,16 @@ so that the total mass injected matches the accreted mass propotioned to magneti l_mass_scale = 0.001 ``` +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. + + A magnetic tower 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. diff --git a/src/pgen/cluster/agn_feedback.cpp b/src/pgen/cluster/agn_feedback.cpp index 3d1569f3..094efdb6 100644 --- a/src/pgen/cluster/agn_feedback.cpp +++ b/src/pgen/cluster/agn_feedback.cpp @@ -53,7 +53,9 @@ AGNFeedback::AGNFeedback(parthenon::ParameterInput *pin, pin->GetOrAddReal("problem/cluster/agn_feedback", "kinetic_jet_offset", 0.02)), enable_tracer_( pin->GetOrAddBoolean("problem/cluster/agn_feedback", "enable_tracer", false)), - disabled_(pin->GetOrAddBoolean("problem/cluster/agn_feedback", "disabled", 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_; @@ -67,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 @@ -254,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 @@ -272,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_; @@ -407,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 64d501bd..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_; @@ -46,6 +47,8 @@ class AGNFeedback { 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/magnetic_tower.cpp b/src/pgen/cluster/magnetic_tower.cpp index 75b78bb8..c2cd1663 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; @@ -51,7 +52,13 @@ void MagneticTower::AddSrcTerm(parthenon::Real field_to_add, parthenon::Real mas 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); @@ -123,9 +130,11 @@ 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)); - // cons(IDN, k, j, i) += cell_delta_rho; + 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; }); } diff --git a/src/pgen/cluster/magnetic_tower.hpp b/src/pgen/cluster/magnetic_tower.hpp index 2a4a06fe..258678ac 100644 --- a/src/pgen/cluster/magnetic_tower.hpp +++ b/src/pgen/cluster/magnetic_tower.hpp @@ -41,8 +41,8 @@ class MagneticTowerObj { l_mass_scale2_(SQR(l_mass_scale)), jet_coords_(jet_coords) { 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 From 1f4c457413a0ae2de4d6c7266290be7248be3703 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Fri, 18 Aug 2023 17:00:55 +0200 Subject: [PATCH 10/26] Allow for specifying perturbation lengthscales in cluster init --- docs/cluster.md | 6 ++++-- src/pgen/cluster.cpp | 31 ++++++++++++++++++++++++++++--- 2 files changed, 32 insertions(+), 5 deletions(-) diff --git a/docs/cluster.md b/docs/cluster.md index e69771ce..e5256a7c 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 = ??? # (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 = ??? # (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 ``` diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index c9de63cd..44319061 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -421,7 +421,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 = @@ -448,8 +461,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); From ecd1f6d0de55980a6795e824879b57f65471eda7 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Fri, 18 Aug 2023 17:55:47 +0200 Subject: [PATCH 11/26] Allow runtime def of tower potential --- docs/cluster.md | 35 +++++++++++- inputs/cluster/cluster.in | 2 +- inputs/cluster/magnetic_tower.in | 1 + src/pgen/cluster/magnetic_tower.cpp | 9 +-- src/pgen/cluster/magnetic_tower.hpp | 88 ++++++++++++++++++++++------- 5 files changed, 107 insertions(+), 28 deletions(-) diff --git a/docs/cluster.md b/docs/cluster.md index e5256a7c..d31f1005 100644 --- a/docs/cluster.md +++ b/docs/cluster.md @@ -412,7 +412,19 @@ material and mass added through the kinetic jet feedback mechanism. ### Magnetic feedback -Magnetic feedback is injected following ([Hui 2006](doi.org/10.1086/501499)) +Multiple options exists to inject magnetic fields: +- a simple field loop (donut) configuration (better suited for kinetically dominated jets at larger scales) +- a more complex pinching tower model (particularly suited for magnetically dominated jets at small scales) + +The option is controlled by the following parameter +``` + +potential_type = li # alternative: "donut" +``` + +#### Pinching magnetic tower + +Magnetic feedback is injected following ([Li 2006](doi.org/10.1086/501499)) where the injected magnetic field follows $$ @@ -469,10 +481,29 @@ enable_magnetic_tower_mass_injection = false In this case, the injected mass through kinetic and thermal feedback according to their ratio. +#### Simple field loop (donut) feedback + +Magnetic energy is injected according to the following simple potential +$$ +A_h(r, \theta, h) = B_0 L \exp^\left ( -r^2/L^2 \right)$ for $h_\mathrm{offset} \leq |h| \leq h_\mathrm{offset} + h_\mathrm{thickness} +$$ +resultig in a magnetic field configuration of +$$ +B_\theta(r, \theta, h) = 2 B_0 r /L \exp^\left ( -r^2/L^2 \right)$ for $h_\mathrm{offset} \leq |h| \leq h_\mathrm{offset} + h_\mathrm{thickness} +$$ +with all other components being zero. + +It is recommended to match the donut thickness to the thickness of the kinetic jet launching region +and to choose a lengthscale that is half the jet launching radius. +This results in almost all injected fields being confined to the launching region (and thus being +carried along with a dominant jet). + +#### Fixed magnetic feedback -A magnetic tower can also be inserted at runtime and injected at a fixed +Magnetic feedback can also be inserted at runtime and injected at a fixed increase in magnetic field, and additional mass can be injected at a fixed rate. +This works for all magnetic vector potentials. ``` initial_field = 0.12431560000204142 diff --git a/inputs/cluster/cluster.in b/inputs/cluster/cluster.in index 61d7bc93..fe1ff02b 100644 --- a/inputs/cluster/cluster.in +++ b/inputs/cluster/cluster.in @@ -181,7 +181,7 @@ kinetic_jet_thickness = 0.05 kinetic_jet_offset = 0.05 - +potential_type = li alpha = 20 l_scale = 0.001 initial_field = 0.12431560000204142 diff --git a/inputs/cluster/magnetic_tower.in b/inputs/cluster/magnetic_tower.in index 66564dca..15b5647d 100644 --- a/inputs/cluster/magnetic_tower.in +++ b/inputs/cluster/magnetic_tower.in @@ -109,6 +109,7 @@ kinetic_fraction = 0 thermal_fraction = 0 +potential_type = li alpha = 20 l_scale = 0.01 initial_field = 0.12431560000204142 diff --git a/src/pgen/cluster/magnetic_tower.cpp b/src/pgen/cluster/magnetic_tower.cpp index c2cd1663..07264d3c 100644 --- a/src/pgen/cluster/magnetic_tower.cpp +++ b/src/pgen/cluster/magnetic_tower.cpp @@ -62,8 +62,9 @@ void MagneticTower::AddSrcTerm(parthenon::Real field_to_add, parthenon::Real mas const JetCoords jet_coords = hydro_pkg->Param("jet_coords_factory").CreateJetCoords(tm.time); - const MagneticTowerObj mt = MagneticTowerObj(field_to_add, alpha_, l_scale_, - density_to_add, l_mass_scale_, jet_coords); + const MagneticTowerObj mt = + MagneticTowerObj(field_to_add, alpha_, l_scale_, density_to_add, l_mass_scale_, + jet_coords, potential_); const auto &eos = hydro_pkg->Param("eos"); @@ -164,7 +165,7 @@ void MagneticTower::ReducePowerContribs(parthenon::Real &linear_contrib, // Make a construct a copy of this with field strength 1 to send to the device const MagneticTowerObj mt = - MagneticTowerObj(1, alpha_, l_scale_, 0, l_mass_scale_, jet_coords); + MagneticTowerObj(1, alpha_, l_scale_, 0, l_mass_scale_, jet_coords, potential_); // Get the reduction of the linear and quadratic contributions ready Real linear_contrib_red, quadratic_contrib_red; @@ -217,7 +218,7 @@ void MagneticTower::AddInitialFieldToPotential(parthenon::MeshBlock *pmb, const JetCoords jet_coords = hydro_pkg->Param("jet_coords_factory").CreateJetCoords(0.0); const MagneticTowerObj mt(initial_field_, alpha_, l_scale_, 0, l_mass_scale_, - jet_coords); + jet_coords, potential_); parthenon::par_for( DEFAULT_LOOP_PATTERN, "MagneticTower::AddInitialFieldToPotential", diff --git a/src/pgen/cluster/magnetic_tower.hpp b/src/pgen/cluster/magnetic_tower.hpp index 258678ac..e8e9a687 100644 --- a/src/pgen/cluster/magnetic_tower.hpp +++ b/src/pgen/cluster/magnetic_tower.hpp @@ -17,8 +17,11 @@ #include #include "jet_coords.hpp" +#include "utils/error_checking.hpp" namespace cluster { + +enum class MagneticTowerPotential { undefined, li, donut }; /************************************************************ * Magnetic Tower Object, for computing magnetic field, vector potential at a * fixed time with a fixed field @@ -32,17 +35,24 @@ class MagneticTowerObj { const parthenon::Real density_, l_mass_scale2_; JetCoords jet_coords_; + // Note that this eventually might better be a template parameter, but while the number + // of potentials implemented is limited (and similarly complex) this should currently + // not be a performance concern. + const MagneticTowerPotential potential_; public: MagneticTowerObj(const parthenon::Real field, const parthenon::Real alpha, const parthenon::Real l_scale, const parthenon::Real density, - const parthenon::Real l_mass_scale, const JetCoords jet_coords) + const parthenon::Real l_mass_scale, const JetCoords jet_coords, + const MagneticTowerPotential potential) : field_(field), alpha_(alpha), l_scale_(l_scale), density_(density), - l_mass_scale2_(SQR(l_mass_scale)), jet_coords_(jet_coords) { + l_mass_scale2_(SQR(l_mass_scale)), jet_coords_(jet_coords), + potential_(potential) { PARTHENON_REQUIRE(l_scale > 0, "Magnetic Tower Length scale must be strictly postitive"); - PARTHENON_REQUIRE(l_mass_scale >= 0, - "Magnetic Tower Mass Length scale must be zero (disabled) or postitive"); + PARTHENON_REQUIRE( + l_mass_scale >= 0, + "Magnetic Tower Mass Length scale must be zero (disabled) or postitive"); } // Compute Jet Potential in jet cylindrical coordinates @@ -50,14 +60,24 @@ class MagneticTowerObj { PotentialInJetCyl(const parthenon::Real r, const parthenon::Real h, parthenon::Real &a_r, parthenon::Real &a_theta, parthenon::Real &a_h) const __attribute__((always_inline)) { - const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2)); - // Compute the potential in jet cylindrical coordinates - a_r = 0.0; - a_theta = 0.0; - if (fabs(h) >= 0.001 && fabs(h) <= 0.001 + 0.000390625) { - a_h = field_ * l_scale_ * exp_r2_h2; + if (potential_ == MagneticTowerPotential::donut) { + const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2)); + // Compute the potential in jet cylindrical coordinates + a_r = 0.0; + a_theta = 0.0; + if (fabs(h) >= 0.001 && fabs(h) <= 0.001 + 0.000390625) { + a_h = field_ * l_scale_ * exp_r2_h2; + } else { + a_h = 0.0; + } + } else if (potential_ == MagneticTowerPotential::li) { + const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2) - pow(h / l_scale_, 2)); + // Compute the potential in jet cylindrical coordinates + a_r = 0.0; + a_theta = field_ * l_scale_ * (r / l_scale_) * exp_r2_h2; + a_h = field_ * l_scale_ * alpha_ / 2.0 * exp_r2_h2; } else { - a_h = 0.0; + PARTHENON_FAIL("Unknown magnetic tower potential."); } } @@ -84,16 +104,25 @@ class MagneticTowerObj { FieldInJetCyl(const parthenon::Real r, const parthenon::Real h, parthenon::Real &b_r, parthenon::Real &b_theta, parthenon::Real &b_h) const __attribute__((always_inline)) { - - const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2)); - // Compute the field in jet cylindrical coordinates - b_r = 0.0; - if (fabs(h) >= 0.001 && fabs(h) <= 0.001 + 0.000390625) { - b_theta = 2.0 * field_ * r / l_scale_ * exp_r2_h2; + if (potential_ == MagneticTowerPotential::donut) { + const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2)); + // Compute the field in jet cylindrical coordinates + b_r = 0.0; + if (fabs(h) >= 0.001 && fabs(h) <= 0.001 + 0.000390625) { + b_theta = 2.0 * field_ * r / l_scale_ * exp_r2_h2; + } else { + b_theta = 0.0; + } + b_h = 0.0; + } else if (potential_ == MagneticTowerPotential::li) { + const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2) - pow(h / l_scale_, 2)); + // Compute the field in jet cylindrical coordinates + b_r = field_ * 2 * (h / l_scale_) * (r / l_scale_) * exp_r2_h2; + b_theta = field_ * alpha_ * (r / l_scale_) * exp_r2_h2; + b_h = field_ * 2 * (1 - pow(r / l_scale_, 2)) * exp_r2_h2; } else { - b_theta = 0.0; + PARTHENON_FAIL("Unknown magnetic tower potential."); } - b_h = 0.0; } // Compute Magnetic field in Simulation Cartesian coordinates @@ -141,6 +170,8 @@ class MagneticTower { const parthenon::Real fixed_mass_rate_; const parthenon::Real l_mass_scale_; + MagneticTowerPotential potential_; + MagneticTower(parthenon::ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg, const std::string &block = "problem/cluster/magnetic_tower") : alpha_(pin->GetOrAddReal(block, "alpha", 0)), @@ -148,17 +179,32 @@ class MagneticTower { initial_field_(pin->GetOrAddReal(block, "initial_field", 0)), fixed_field_rate_(pin->GetOrAddReal(block, "fixed_field_rate", 0)), fixed_mass_rate_(pin->GetOrAddReal(block, "fixed_mass_rate", 0)), - l_mass_scale_(pin->GetOrAddReal(block, "l_mass_scale", 0)) { - hydro_pkg->AddParam<>("magnetic_tower", *this); + l_mass_scale_(pin->GetOrAddReal(block, "l_mass_scale", 0)), + potential_(MagneticTowerPotential::undefined) { hydro_pkg->AddParam("magnetic_tower_linear_contrib", 0.0, true); hydro_pkg->AddParam("magnetic_tower_quadratic_contrib", 0.0, true); + const auto potential_str = pin->GetOrAddString(block, "potential_type", "undefined"); + + if (potential_str == "donut") { + potential_ = MagneticTowerPotential::donut; + } else if (potential_str == "li") { + potential_ = MagneticTowerPotential::li; + } else { + PARTHENON_FAIL( + "Unknown potential for magnetic tower. Current options are: donut, li") + } + // Vector potential is only locally used, so no need to // communicate/restrict/prolongate/fluxes/etc parthenon::Metadata m({parthenon::Metadata::Cell, parthenon::Metadata::Derived, parthenon::Metadata::OneCopy}, std::vector({3})); hydro_pkg->AddField("magnetic_tower_A", m); + + // Finally, add object to params (should be done last as otherwise modification within + // this function would not survive). + hydro_pkg->AddParam<>("magnetic_tower", *this); } // Add initial magnetic field to provided potential with a single meshblock From 1b27b3ab701089d53c3a9c3a7fc74d026dd9a2b8 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Fri, 18 Aug 2023 19:03:32 +0200 Subject: [PATCH 12/26] Add offset and thickness params to donut feedback --- docs/cluster.md | 13 +++++++++++++ src/pgen/cluster/magnetic_tower.cpp | 12 ++++++------ src/pgen/cluster/magnetic_tower.hpp | 23 +++++++++++++++++------ 3 files changed, 36 insertions(+), 12 deletions(-) diff --git a/docs/cluster.md b/docs/cluster.md index d31f1005..e1444e3b 100644 --- a/docs/cluster.md +++ b/docs/cluster.md @@ -448,6 +448,7 @@ $$ The parameters $\alpha$ and $\ell$ can be changed with ``` +potential_type = li alpha = 20 l_scale = 0.001 ``` @@ -493,11 +494,23 @@ B_\theta(r, \theta, h) = 2 B_0 r /L \exp^\left ( -r^2/L^2 \right)$ for $h_\mathr $$ with all other components being zero. + +``` + +potential_type = donut +l_scale = 0.0005 # in code length +offset = 0.001 # in code length +thickness = 0.0001 # in code length +``` + It is recommended to match the donut thickness to the thickness of the kinetic jet launching region and to choose a lengthscale that is half the jet launching radius. This results in almost all injected fields being confined to the launching region (and thus being carried along with a dominant jet). +Mass feedback for the donut model is currently identical to the Li tower above +(and, thus, the parameters above pertaining to the mass injection equally apply here). + #### Fixed magnetic feedback Magnetic feedback can also be inserted at runtime and injected at a fixed diff --git a/src/pgen/cluster/magnetic_tower.cpp b/src/pgen/cluster/magnetic_tower.cpp index 07264d3c..1679ffd7 100644 --- a/src/pgen/cluster/magnetic_tower.cpp +++ b/src/pgen/cluster/magnetic_tower.cpp @@ -63,8 +63,8 @@ void MagneticTower::AddSrcTerm(parthenon::Real field_to_add, parthenon::Real mas const JetCoords jet_coords = hydro_pkg->Param("jet_coords_factory").CreateJetCoords(tm.time); const MagneticTowerObj mt = - MagneticTowerObj(field_to_add, alpha_, l_scale_, density_to_add, l_mass_scale_, - jet_coords, potential_); + MagneticTowerObj(field_to_add, alpha_, l_scale_, offset_, thickness_, + density_to_add, l_mass_scale_, jet_coords, potential_); const auto &eos = hydro_pkg->Param("eos"); @@ -164,8 +164,8 @@ void MagneticTower::ReducePowerContribs(parthenon::Real &linear_contrib, hydro_pkg->Param("jet_coords_factory").CreateJetCoords(tm.time); // Make a construct a copy of this with field strength 1 to send to the device - const MagneticTowerObj mt = - MagneticTowerObj(1, alpha_, l_scale_, 0, l_mass_scale_, jet_coords, potential_); + const MagneticTowerObj mt = MagneticTowerObj(1, alpha_, l_scale_, offset_, thickness_, + 0, l_mass_scale_, jet_coords, potential_); // Get the reduction of the linear and quadratic contributions ready Real linear_contrib_red, quadratic_contrib_red; @@ -217,8 +217,8 @@ void MagneticTower::AddInitialFieldToPotential(parthenon::MeshBlock *pmb, const JetCoords jet_coords = hydro_pkg->Param("jet_coords_factory").CreateJetCoords(0.0); - const MagneticTowerObj mt(initial_field_, alpha_, l_scale_, 0, l_mass_scale_, - jet_coords, potential_); + const MagneticTowerObj mt(initial_field_, alpha_, l_scale_, offset_, thickness_, 0, + l_mass_scale_, jet_coords, potential_); parthenon::par_for( DEFAULT_LOOP_PATTERN, "MagneticTower::AddInitialFieldToPotential", diff --git a/src/pgen/cluster/magnetic_tower.hpp b/src/pgen/cluster/magnetic_tower.hpp index e8e9a687..c32ede67 100644 --- a/src/pgen/cluster/magnetic_tower.hpp +++ b/src/pgen/cluster/magnetic_tower.hpp @@ -31,6 +31,7 @@ class MagneticTowerObj { private: const parthenon::Real field_; const parthenon::Real alpha_, l_scale_; + const parthenon::Real offset_, thickness_; const parthenon::Real density_, l_mass_scale2_; @@ -42,12 +43,13 @@ class MagneticTowerObj { public: MagneticTowerObj(const parthenon::Real field, const parthenon::Real alpha, - const parthenon::Real l_scale, const parthenon::Real density, + const parthenon::Real l_scale, const parthenon::Real offset, + const parthenon::Real thickness, const parthenon::Real density, const parthenon::Real l_mass_scale, const JetCoords jet_coords, const MagneticTowerPotential potential) - : field_(field), alpha_(alpha), l_scale_(l_scale), density_(density), - l_mass_scale2_(SQR(l_mass_scale)), jet_coords_(jet_coords), - potential_(potential) { + : field_(field), alpha_(alpha), l_scale_(l_scale), offset_(offset), + thickness_(thickness), density_(density), l_mass_scale2_(SQR(l_mass_scale)), + jet_coords_(jet_coords), potential_(potential) { PARTHENON_REQUIRE(l_scale > 0, "Magnetic Tower Length scale must be strictly postitive"); PARTHENON_REQUIRE( @@ -65,7 +67,7 @@ class MagneticTowerObj { // Compute the potential in jet cylindrical coordinates a_r = 0.0; a_theta = 0.0; - if (fabs(h) >= 0.001 && fabs(h) <= 0.001 + 0.000390625) { + if (fabs(h) >= 0.001 && fabs(h) <= offset_ + thickness_) { a_h = field_ * l_scale_ * exp_r2_h2; } else { a_h = 0.0; @@ -108,7 +110,7 @@ class MagneticTowerObj { const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2)); // Compute the field in jet cylindrical coordinates b_r = 0.0; - if (fabs(h) >= 0.001 && fabs(h) <= 0.001 + 0.000390625) { + if (fabs(h) >= 0.001 && fabs(h) <= offset_ + thickness_) { b_theta = 2.0 * field_ * r / l_scale_ * exp_r2_h2; } else { b_theta = 0.0; @@ -163,6 +165,7 @@ class MagneticTowerObj { class MagneticTower { public: const parthenon::Real alpha_, l_scale_; + const parthenon::Real offset_, thickness_; const parthenon::Real initial_field_; const parthenon::Real fixed_field_rate_; @@ -176,6 +179,8 @@ class MagneticTower { const std::string &block = "problem/cluster/magnetic_tower") : alpha_(pin->GetOrAddReal(block, "alpha", 0)), l_scale_(pin->GetOrAddReal(block, "l_scale", 0)), + offset_(pin->GetOrAddReal(block, "offset", 0)), + thickness_(pin->GetOrAddReal(block, "thickness", 0)), initial_field_(pin->GetOrAddReal(block, "initial_field", 0)), fixed_field_rate_(pin->GetOrAddReal(block, "fixed_field_rate", 0)), fixed_mass_rate_(pin->GetOrAddReal(block, "fixed_mass_rate", 0)), @@ -188,8 +193,14 @@ class MagneticTower { if (potential_str == "donut") { potential_ = MagneticTowerPotential::donut; + PARTHENON_REQUIRE_THROWS( + offset_ >= 0.0 && thickness_ > 0.0, + "Incompatible combination of offset and thickness for magnetic donut feedback.") } else if (potential_str == "li") { potential_ = MagneticTowerPotential::li; + PARTHENON_REQUIRE_THROWS(offset_ <= 0.0 && thickness_ <= 0.0, + "Please disable (set to zero) tower offset and thickness " + "for the Li tower model"); } else { PARTHENON_FAIL( "Unknown potential for magnetic tower. Current options are: donut, li") From 261fd549e6af81c7c272bf3e2154e27a75b92130 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 21 Aug 2023 11:10:59 +0200 Subject: [PATCH 13/26] Address review comments --- docs/cluster.md | 4 ++-- src/hydro/prolongation/custom_ops.hpp | 7 ++++++- src/pgen/cluster/magnetic_tower.hpp | 3 --- src/pgen/cluster/stellar_feedback.cpp | 4 +--- 4 files changed, 9 insertions(+), 9 deletions(-) diff --git a/docs/cluster.md b/docs/cluster.md index e1444e3b..7e650cfa 100644 --- a/docs/cluster.md +++ b/docs/cluster.md @@ -184,7 +184,7 @@ 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) l_peak_v = ??? # lengthscale (in code length units) where the velocity spectrum peaks. No default value. -k_peak_v = ??? # (alternative to l_peak_v): wavenumber in normalized 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 @@ -192,7 +192,7 @@ rseed_v = 1 # (optional) integer seed for RNG for wavenumbers and amplitu # for the magnetic field sigma_b = 0.0 # volume weighted RMS |B| in code magnetic; default: 0.0 (disabled) l_peak_b = ??? # lengthscale (in code length units) where the magnetic field spectrum peaks. No default value. -k_peak_b = ??? # (alternative to l_peak_b): wavenumber in normalized 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 ``` diff --git a/src/hydro/prolongation/custom_ops.hpp b/src/hydro/prolongation/custom_ops.hpp index 9f42b3a4..02d2b15f 100644 --- a/src/hydro/prolongation/custom_ops.hpp +++ b/src/hydro/prolongation/custom_ops.hpp @@ -41,6 +41,11 @@ 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) { @@ -77,7 +82,7 @@ struct ProlongateCellMinModMultiD { Real dx1fm = 0; Real dx1fp = 0; Real gx1c = 0; - if constexpr (INCLUDE_X2) { + if constexpr (INCLUDE_X1) { Real dx1m, dx1p; GetGridSpacings<1, el>(coords, coarse_coords, cib, ib, i, fi, &dx1m, &dx1p, &dx1fm, &dx1fp); diff --git a/src/pgen/cluster/magnetic_tower.hpp b/src/pgen/cluster/magnetic_tower.hpp index c32ede67..743aa58b 100644 --- a/src/pgen/cluster/magnetic_tower.hpp +++ b/src/pgen/cluster/magnetic_tower.hpp @@ -201,9 +201,6 @@ class MagneticTower { PARTHENON_REQUIRE_THROWS(offset_ <= 0.0 && thickness_ <= 0.0, "Please disable (set to zero) tower offset and thickness " "for the Li tower model"); - } else { - PARTHENON_FAIL( - "Unknown potential for magnetic tower. Current options are: donut, li") } // Vector potential is only locally used, so no need to diff --git a/src/pgen/cluster/stellar_feedback.cpp b/src/pgen/cluster/stellar_feedback.cpp index 8691cde3..6c410d71 100644 --- a/src/pgen/cluster/stellar_feedback.cpp +++ b/src/pgen/cluster/stellar_feedback.cpp @@ -140,9 +140,7 @@ void StellarFeedback::FeedbackSrcTerm(parthenon::MeshData *md, AddDensityToConsAtFixedVelTemp(cell_delta_rho, cons, prim, eos.GetGamma(), k, j, i); // Then add thermal energy - const auto cell_delta_mass = -cell_delta_rho * coords.CellVolume(k, j, i); - const auto cell_delta_energy_density = - mass_to_energy * cell_delta_mass / coords.CellVolume(k, j, i); + 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; From 14ea977bde791c5005f710621061c996f024d8d7 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 23 Aug 2023 15:51:14 -0400 Subject: [PATCH 14/26] Cool first 'unsplit' then problem source --- src/hydro/hydro.cpp | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index 95218748..d22dddef 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -143,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); } @@ -153,14 +161,14 @@ 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"); + // const auto &enable_cooling = hydro_pkg->Param("enable_cooling"); - if (enable_cooling == Cooling::tabular) { - const TabularCooling &tabular_cooling = - hydro_pkg->Param("tabular_cooling"); + // if (enable_cooling == Cooling::tabular) { + // const TabularCooling &tabular_cooling = + // hydro_pkg->Param("tabular_cooling"); - tabular_cooling.SrcTerm(md, tm.dt); - } + // tabular_cooling.SrcTerm(md, tm.dt); + // } if (ProblemSourceFirstOrder != nullptr) { ProblemSourceFirstOrder(md, tm, tm.dt); } From a3c2855b6a89d719d69db0ebd4348d0495325e27 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 23 Aug 2023 15:55:11 -0400 Subject: [PATCH 15/26] Bump parth to enable chunking by default --- external/parthenon | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/parthenon b/external/parthenon index ab9d1671..5b6bb619 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit ab9d16714d31e853d000a0757f3487bcb9804023 +Subproject commit 5b6bb61906f7c278f9724ee9f38e79dee8707098 From 1408069a66a0fe97fa4ace6ef7f1f5edf513e422 Mon Sep 17 00:00:00 2001 From: Forrest Wolfgang Glines Date: Thu, 24 Aug 2023 10:59:07 -0600 Subject: [PATCH 16/26] Added reductions for stellar mass, clips --- src/main.cpp | 3 +- src/pgen/cluster.cpp | 87 +++++++++++++++++++++++---- src/pgen/cluster/stellar_feedback.cpp | 23 ++++--- src/pgen/pgen.hpp | 3 +- 4 files changed, 96 insertions(+), 20 deletions(-) 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/cluster.cpp b/src/pgen/cluster.cpp index 44319061..ba5a7263 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -88,10 +88,18 @@ void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, 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) { + 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); @@ -106,6 +114,7 @@ void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, 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; } @@ -126,7 +135,9 @@ void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, prim(IV3, k, j, i) *= vceil / v; // Remove kinetic energy - cons(IEN, k, j, i) -= 0.5 * prim(IDN, k, j, i) * (v2 - vceil2); + 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; } } @@ -142,6 +153,7 @@ void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, 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; } @@ -151,12 +163,25 @@ void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, // 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); + 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")); } } @@ -173,7 +198,7 @@ void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, } } -void ClusterSrcTerm(MeshData *md, const parthenon::SimTime &tm, +void ClusterUnsplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, const Real beta_dt) { auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); @@ -195,11 +220,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); +}; +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, beta_dt, tm); + stellar_feedback.FeedbackSrcTerm(md, dt, tm); - ApplyClusterClips(md, tm, beta_dt); -}; + ApplyClusterClips(md, tm, dt); +} Real ClusterEstimateTimestep(MeshData *md) { Real min_dt = std::numeric_limits::max(); @@ -385,6 +415,41 @@ 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 and stellar mass + ************************************************************/ + + /* 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)); + } + hydro_pkg->UpdateParam(parthenon::hist_param_key, hst_vars); + /************************************************************ * Add derived fields * NOTE: these must be filled in UserWorkBeforeOutput diff --git a/src/pgen/cluster/stellar_feedback.cpp b/src/pgen/cluster/stellar_feedback.cpp index 6c410d71..4921ce7d 100644 --- a/src/pgen/cluster/stellar_feedback.cpp +++ b/src/pgen/cluster/stellar_feedback.cpp @@ -105,11 +105,17 @@ void StellarFeedback::FeedbackSrcTerm(parthenon::MeshData *md, //////////////////////////////////////////////////////////////////////////////// - // Constant volumetric heating - parthenon::par_for( - DEFAULT_LOOP_PATTERN, "StellarFeedback::FeedbackSrcTerm", 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) { + 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); @@ -135,6 +141,7 @@ void StellarFeedback::FeedbackSrcTerm(parthenon::MeshData *md, // 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, @@ -142,12 +149,14 @@ void StellarFeedback::FeedbackSrcTerm(parthenon::MeshData *md, // 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.") + "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/pgen.hpp b/src/pgen/pgen.hpp index 3050cbb5..11dfb101 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -100,7 +100,8 @@ 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 From 2e2f5542a245e2ca30db6000987f9cee5c7f6cf0 Mon Sep 17 00:00:00 2001 From: Forrest Wolfgang Glines Date: Thu, 24 Aug 2023 16:16:02 -0600 Subject: [PATCH 17/26] Moved Clusterclips into own src/header --- src/pgen/CMakeLists.txt | 1 + src/pgen/cluster.cpp | 141 +---------------------- src/pgen/cluster/cluster_clips.cpp | 172 ++++++++++++++++++++++++++++ src/pgen/cluster/cluster_clips.hpp | 27 +++++ src/pgen/cluster/magnetic_tower.hpp | 11 +- 5 files changed, 208 insertions(+), 144 deletions(-) create mode 100644 src/pgen/cluster/cluster_clips.cpp create mode 100644 src/pgen/cluster/cluster_clips.hpp diff --git a/src/pgen/CMakeLists.txt b/src/pgen/CMakeLists.txt index 02b51c3f..58ac4627 100644 --- a/src/pgen/CMakeLists.txt +++ b/src/pgen/CMakeLists.txt @@ -8,6 +8,7 @@ target_sources(athenaPK PRIVATE cluster.cpp cluster/agn_feedback.cpp cluster/agn_triggering.cpp + cluster/cluster_clips.cpp cluster/hydrostatic_equilibrium_sphere.cpp cluster/magnetic_tower.cpp cluster/snia_feedback.cpp diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index ba5a7263..744d6492 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -46,6 +46,7 @@ #include "cluster/agn_feedback.hpp" #include "cluster/agn_triggering.hpp" #include "cluster/cluster_gravity.hpp" +#include "cluster/cluster_clips.hpp" #include "cluster/entropy_profiles.hpp" #include "cluster/hydrostatic_equilibrium_sphere.hpp" #include "cluster/magnetic_tower.hpp" @@ -57,146 +58,6 @@ 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); - - 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")); - } -} - -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 ClusterUnsplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, const Real beta_dt) { diff --git a/src/pgen/cluster/cluster_clips.cpp b/src/pgen/cluster/cluster_clips.cpp new file mode 100644 index 00000000..0f5bb569 --- /dev/null +++ b/src/pgen/cluster/cluster_clips.cpp @@ -0,0 +1,172 @@ + +//======================================================================================== +// 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..6da56fa0 --- /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/magnetic_tower.hpp b/src/pgen/cluster/magnetic_tower.hpp index 743aa58b..3d6d536f 100644 --- a/src/pgen/cluster/magnetic_tower.hpp +++ b/src/pgen/cluster/magnetic_tower.hpp @@ -177,10 +177,10 @@ class MagneticTower { 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, "offset", 0)), - thickness_(pin->GetOrAddReal(block, "thickness", 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)), @@ -195,7 +195,10 @@ class MagneticTower { potential_ = MagneticTowerPotential::donut; PARTHENON_REQUIRE_THROWS( offset_ >= 0.0 && thickness_ > 0.0, - "Incompatible combination of offset and thickness for magnetic donut feedback.") + "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, From aa13207bb9a7bd5d5d652036781fd8cd35671afc Mon Sep 17 00:00:00 2001 From: Forrest Wolfgang Glines Date: Fri, 25 Aug 2023 13:12:39 -0600 Subject: [PATCH 18/26] Added cold gas, AGN extent reductions --- src/pgen/CMakeLists.txt | 1 + src/pgen/cluster.cpp | 29 +++++++- src/pgen/cluster/cluster_clips.cpp | 2 - src/pgen/cluster/cluster_reductions.cpp | 99 +++++++++++++++++++++++++ src/pgen/cluster/cluster_reductions.hpp | 25 +++++++ 5 files changed, 152 insertions(+), 4 deletions(-) create mode 100644 src/pgen/cluster/cluster_reductions.cpp create mode 100644 src/pgen/cluster/cluster_reductions.hpp diff --git a/src/pgen/CMakeLists.txt b/src/pgen/CMakeLists.txt index 58ac4627..382b5b7b 100644 --- a/src/pgen/CMakeLists.txt +++ b/src/pgen/CMakeLists.txt @@ -9,6 +9,7 @@ target_sources(athenaPK PRIVATE 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 diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index 744d6492..7be64315 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -45,8 +45,9 @@ // Cluster headers #include "cluster/agn_feedback.hpp" #include "cluster/agn_triggering.hpp" -#include "cluster/cluster_gravity.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" @@ -277,7 +278,8 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd hydro_pkg->AddParam("cluster_clip_r", clip_r); /************************************************************ - * Start running reductions into history outputs for clips and stellar mass + * 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 @@ -309,6 +311,29 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd }, 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){ + auto mbar_over_kb = hydro_pkg->Param("mbar_over_kb"); + 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); /************************************************************ diff --git a/src/pgen/cluster/cluster_clips.cpp b/src/pgen/cluster/cluster_clips.cpp index 0f5bb569..80b16a16 100644 --- a/src/pgen/cluster/cluster_clips.cpp +++ b/src/pgen/cluster/cluster_clips.cpp @@ -26,7 +26,6 @@ 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"); @@ -168,5 +167,4 @@ void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, } } - } // namespace cluster \ 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..05d9e2f7 --- /dev/null +++ b/src/pgen/cluster/cluster_reductions.cpp @@ -0,0 +1,99 @@ + +//======================================================================================== +// 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( + "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 &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( + "Cluster::ApplyClusterClips", + 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) > 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..57d56b66 --- /dev/null +++ b/src/pgen/cluster/cluster_reductions.hpp @@ -0,0 +1,25 @@ +#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_ From 3a8b38de86a80986e2b9b3722f3a31e4ea168683 Mon Sep 17 00:00:00 2001 From: Forrest Wolfgang Glines Date: Fri, 25 Aug 2023 15:33:25 -0600 Subject: [PATCH 19/26] Fixed kernel names --- src/pgen/cluster/cluster_reductions.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pgen/cluster/cluster_reductions.cpp b/src/pgen/cluster/cluster_reductions.cpp index 05d9e2f7..6118566a 100644 --- a/src/pgen/cluster/cluster_reductions.cpp +++ b/src/pgen/cluster/cluster_reductions.cpp @@ -42,7 +42,7 @@ LocalReduceColdGas(parthenon::MeshData *md){ Real cold_gas = 0.0; Kokkos::parallel_reduce( - "Cluster::ApplyClusterClips", + "LocalReduceColdGas", Kokkos::MDRangePolicy>( DevExecSpace(), {0, kb.s, jb.s, ib.s}, {prim_pack.GetDim(5), kb.e + 1, jb.e + 1, ib.e + 1}, @@ -76,7 +76,7 @@ LocalReduceAGNExtent(parthenon::MeshData *md){ Real max_r2 = 0.0; Kokkos::parallel_reduce( - "Cluster::ApplyClusterClips", + "LocalReduceAGNExtent", Kokkos::MDRangePolicy>( DevExecSpace(), {0, kb.s, jb.s, ib.s}, {cons_pack.GetDim(5), kb.e + 1, jb.e + 1, ib.e + 1}, From 5772f59c0b34009cd956707d5a874849278f13bd Mon Sep 17 00:00:00 2001 From: Forrest Wolfgang Glines Date: Mon, 28 Aug 2023 14:10:20 -0600 Subject: [PATCH 20/26] Non-positive vceil Tceil in AGN do nothing --- src/pgen/cluster.cpp | 1 - src/pgen/cluster/agn_feedback.cpp | 4 ++-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index 7be64315..2c2fe049 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -81,7 +81,6 @@ void ClusterUnsplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, const auto &snia_feedback = hydro_pkg->Param("snia_feedback"); snia_feedback.FeedbackSrcTerm(md, beta_dt, tm); - }; void ClusterSplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, const Real dt) { diff --git a/src/pgen/cluster/agn_feedback.cpp b/src/pgen/cluster/agn_feedback.cpp index 094efdb6..bcf63c24 100644 --- a/src/pgen/cluster/agn_feedback.cpp +++ b/src/pgen/cluster/agn_feedback.cpp @@ -390,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; @@ -406,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; } From 711973bedae296eb5737964df8a025548a16beb6 Mon Sep 17 00:00:00 2001 From: par-hermes Date: Mon, 28 Aug 2023 20:13:20 +0000 Subject: [PATCH 21/26] cpp-py-formatter --- src/pgen/cluster.cpp | 42 ++++++++++----------- src/pgen/cluster/cluster_clips.cpp | 50 +++++++++++++------------ src/pgen/cluster/cluster_clips.hpp | 4 +- src/pgen/cluster/cluster_reductions.cpp | 44 +++++++++++----------- src/pgen/cluster/cluster_reductions.hpp | 9 ++--- src/pgen/cluster/magnetic_tower.hpp | 6 +-- src/pgen/cluster/stellar_feedback.cpp | 24 ++++++------ src/pgen/pgen.hpp | 6 ++- 8 files changed, 96 insertions(+), 89 deletions(-) diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index 2c2fe049..96a781ba 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -59,9 +59,8 @@ using namespace parthenon::driver::prelude; using namespace parthenon::package::prelude; using utils::few_modes_ft::FewModesFT; - void ClusterUnsplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, - const Real beta_dt) { + const Real beta_dt) { auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); const bool &gravity_srcterm = hydro_pkg->Param("gravity_srcterm"); @@ -83,7 +82,7 @@ void ClusterUnsplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, snia_feedback.FeedbackSrcTerm(md, beta_dt, tm); }; void ClusterSplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, - const Real dt) { + const Real dt) { auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); const auto &stellar_feedback = hydro_pkg->Param("stellar_feedback"); @@ -289,50 +288,49 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd 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"}; + 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 + // 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); + 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); + // 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 + // 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){ + if (cold_thresh > 0) { hydro_pkg->AddParam("reduction_cold_threshold", cold_thresh); hst_vars.emplace_back(parthenon::HistoryOutputVar( - parthenon::UserHistoryOperation::sum, LocalReduceColdGas, - "cold_mass")); + 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){ - auto mbar_over_kb = hydro_pkg->Param("mbar_over_kb"); - PARTHENON_REQUIRE(pin->GetOrAddBoolean("problem/cluster/agn_feedback", "enable_tracer", false), - "AGN Tracer must be enabled to reduce AGN tracer extent"); + if (agn_tracer_thresh >= 0) { + auto mbar_over_kb = hydro_pkg->Param("mbar_over_kb"); + 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")); + parthenon::UserHistoryOperation::max, LocalReduceAGNExtent, "agn_extent")); } - hydro_pkg->UpdateParam(parthenon::hist_param_key, hst_vars); /************************************************************ diff --git a/src/pgen/cluster/cluster_clips.cpp b/src/pgen/cluster/cluster_clips.cpp index 80b16a16..5aa86eea 100644 --- a/src/pgen/cluster/cluster_clips.cpp +++ b/src/pgen/cluster/cluster_clips.cpp @@ -74,14 +74,14 @@ void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, 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}), + "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) { + 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); @@ -96,7 +96,7 @@ void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, 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); + added_dfloor_mass_team += (dfloor - rho) * coords.CellVolume(k, j, i); cons(IDN, k, j, i) = dfloor; prim(IDN, k, j, i) = dfloor; } @@ -118,7 +118,7 @@ void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, // 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); + removed_vceil_energy_team += removed_energy * coords.CellVolume(k, j, i); cons(IEN, k, j, i) -= removed_energy; } } @@ -135,7 +135,7 @@ void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, 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); + 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; } @@ -146,24 +146,28 @@ void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, 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); + 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")); + }, + 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")); } } diff --git a/src/pgen/cluster/cluster_clips.hpp b/src/pgen/cluster/cluster_clips.hpp index 6da56fa0..74351f0a 100644 --- a/src/pgen/cluster/cluster_clips.hpp +++ b/src/pgen/cluster/cluster_clips.hpp @@ -11,13 +11,13 @@ // \brief Class for applying floors and ceils and reducing removed/added mass/energy // C++ headers -#include // sqrt() +#include // sqrt() // Parthenon headers #include "basic_types.hpp" #include "mesh/mesh.hpp" -namespace cluster{ +namespace cluster { void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, const Real beta_dt); diff --git a/src/pgen/cluster/cluster_reductions.cpp b/src/pgen/cluster/cluster_reductions.cpp index 6118566a..0e8c88ed 100644 --- a/src/pgen/cluster/cluster_reductions.cpp +++ b/src/pgen/cluster/cluster_reductions.cpp @@ -5,7 +5,8 @@ // 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 +// \brief Cluster-specific reductions to compute the total cold gas and maximum radius +// of AGN feedback // Parthenon headers #include "kokkos_abstraction.hpp" @@ -22,13 +23,13 @@ namespace cluster { using namespace parthenon; -parthenon::Real -LocalReduceColdGas(parthenon::MeshData *md){ +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); + 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"}); @@ -42,26 +43,26 @@ LocalReduceColdGas(parthenon::MeshData *md){ 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}), + "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); + if (internal_e < e_thresh) { + cold_gas_team += prim(IDN, k, j, i) * coords.CellVolume(k, j, i); } - }, cold_gas); + }, + cold_gas); return cold_gas; } -parthenon::Real -LocalReduceAGNExtent(parthenon::MeshData *md){ +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"); @@ -76,11 +77,11 @@ LocalReduceAGNExtent(parthenon::MeshData *md){ 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}), + "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); @@ -88,10 +89,11 @@ LocalReduceAGNExtent(parthenon::MeshData *md){ 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) > tracer_thresh && r2 > max_r2 ){ + if (cons(nhydro, k, j, i) > tracer_thresh && r2 > max_r2) { max_r2_team = r2; } - }, Kokkos::Max(max_r2)); + }, + Kokkos::Max(max_r2)); return std::sqrt(max_r2); } diff --git a/src/pgen/cluster/cluster_reductions.hpp b/src/pgen/cluster/cluster_reductions.hpp index 57d56b66..4fd98bc2 100644 --- a/src/pgen/cluster/cluster_reductions.hpp +++ b/src/pgen/cluster/cluster_reductions.hpp @@ -6,7 +6,8 @@ // 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 +// \brief Cluster-specific reductions to compute the total cold gas and maximum radius +// of AGN feedback // parthenon headers #include @@ -14,11 +15,9 @@ namespace cluster { -parthenon::Real -LocalReduceColdGas(parthenon::MeshData *md); +parthenon::Real LocalReduceColdGas(parthenon::MeshData *md); -parthenon::Real -LocalReduceAGNExtent(parthenon::MeshData *md); +parthenon::Real LocalReduceAGNExtent(parthenon::MeshData *md); } // namespace cluster diff --git a/src/pgen/cluster/magnetic_tower.hpp b/src/pgen/cluster/magnetic_tower.hpp index 3d6d536f..ebcc1d44 100644 --- a/src/pgen/cluster/magnetic_tower.hpp +++ b/src/pgen/cluster/magnetic_tower.hpp @@ -193,9 +193,9 @@ class MagneticTower { 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(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"); diff --git a/src/pgen/cluster/stellar_feedback.cpp b/src/pgen/cluster/stellar_feedback.cpp index 4921ce7d..f4eb774e 100644 --- a/src/pgen/cluster/stellar_feedback.cpp +++ b/src/pgen/cluster/stellar_feedback.cpp @@ -109,11 +109,11 @@ void StellarFeedback::FeedbackSrcTerm(parthenon::MeshData *md, // 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}), + "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); @@ -141,22 +141,24 @@ void StellarFeedback::FeedbackSrcTerm(parthenon::MeshData *md, // 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); + 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."); + 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")); + }, + stellar_mass); + hydro_pkg->UpdateParam( + "stellar_mass", stellar_mass + hydro_pkg->Param("stellar_mass")); } } // namespace cluster diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 11dfb101..182d8497 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -100,8 +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 ClusterUnsplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, const Real beta_dt); -void ClusterSplitSrcTerm(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 From 2bdf1f8cd6aa0d57fa745c15603c51392aea96bf Mon Sep 17 00:00:00 2001 From: Forrest Wolfgang Glines Date: Sun, 10 Sep 2023 11:54:48 -0600 Subject: [PATCH 22/26] Added exclusion radius to stellar feedback --- src/pgen/cluster/stellar_feedback.cpp | 16 +++++++++++++--- src/pgen/cluster/stellar_feedback.hpp | 3 ++- 2 files changed, 15 insertions(+), 4 deletions(-) diff --git a/src/pgen/cluster/stellar_feedback.cpp b/src/pgen/cluster/stellar_feedback.cpp index f4eb774e..db5c6eae 100644 --- a/src/pgen/cluster/stellar_feedback.cpp +++ b/src/pgen/cluster/stellar_feedback.cpp @@ -33,19 +33,28 @@ 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 && efficiency_ == 0.0 && number_density_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; } - PARTHENON_REQUIRE(disabled_ || (stellar_radius_ != 0.0 && efficiency_ != 0.00 && + + 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."); @@ -98,6 +107,7 @@ void StellarFeedback::FeedbackSrcTerm(parthenon::MeshData *md, 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_; @@ -125,7 +135,7 @@ void StellarFeedback::FeedbackSrcTerm(parthenon::MeshData *md, const auto z = coords.Xc<3>(k); const auto r = sqrt(x * x + y * y + z * z); - if (r > stellar_radius) { + if (r > stellar_radius || r <= exclusion_radius) { return; } diff --git a/src/pgen/cluster/stellar_feedback.hpp b/src/pgen/cluster/stellar_feedback.hpp index d91dff6a..973716b4 100644 --- a/src/pgen/cluster/stellar_feedback.hpp +++ b/src/pgen/cluster/stellar_feedback.hpp @@ -26,6 +26,7 @@ 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 @@ -38,7 +39,7 @@ class StellarFeedback { void FeedbackSrcTerm(parthenon::MeshData *md, const parthenon::Real beta_dt, const parthenon::SimTime &tm) const; - // Apply the feedback from stellare tied to the BCG density + // 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, From d4a3d19edc9d2b1865370ce8cb3044b1d0541d1e Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 20 Sep 2023 11:00:46 +0200 Subject: [PATCH 23/26] Use prim scalar for AGN extent --- src/pgen/cluster.cpp | 1 - src/pgen/cluster/cluster_reductions.cpp | 2 +- src/pgen/cluster/stellar_feedback.cpp | 16 +++++++++------- src/pgen/cluster/stellar_feedback.hpp | 2 +- 4 files changed, 11 insertions(+), 10 deletions(-) diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index 96a781ba..a2b10f7e 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -322,7 +322,6 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd const Real agn_tracer_thresh = pin->GetOrAddReal("problem/cluster/reductions", "agn_tracer_thresh", -1.0); if (agn_tracer_thresh >= 0) { - auto mbar_over_kb = hydro_pkg->Param("mbar_over_kb"); PARTHENON_REQUIRE( pin->GetOrAddBoolean("problem/cluster/agn_feedback", "enable_tracer", false), "AGN Tracer must be enabled to reduce AGN tracer extent"); diff --git a/src/pgen/cluster/cluster_reductions.cpp b/src/pgen/cluster/cluster_reductions.cpp index 0e8c88ed..9e9edb4b 100644 --- a/src/pgen/cluster/cluster_reductions.cpp +++ b/src/pgen/cluster/cluster_reductions.cpp @@ -89,7 +89,7 @@ parthenon::Real LocalReduceAGNExtent(parthenon::MeshData *md) { 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) > tracer_thresh && r2 > max_r2) { + if (cons(nhydro, k, j, i) / cons(IDN, k, j, i) > tracer_thresh && r2 > max_r2) { max_r2_team = r2; } }, diff --git a/src/pgen/cluster/stellar_feedback.cpp b/src/pgen/cluster/stellar_feedback.cpp index db5c6eae..aa130d65 100644 --- a/src/pgen/cluster/stellar_feedback.cpp +++ b/src/pgen/cluster/stellar_feedback.cpp @@ -41,22 +41,24 @@ StellarFeedback::StellarFeedback(parthenon::ParameterInput *pin, "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) { + 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 (!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); + 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), + 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); diff --git a/src/pgen/cluster/stellar_feedback.hpp b/src/pgen/cluster/stellar_feedback.hpp index 973716b4..e3d4b20c 100644 --- a/src/pgen/cluster/stellar_feedback.hpp +++ b/src/pgen/cluster/stellar_feedback.hpp @@ -26,7 +26,7 @@ class StellarFeedback { private: // feedback parameters in code units const parthenon::Real stellar_radius_; // length - parthenon::Real exclusion_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 From 36d51399ea607f7b7c42994191743be794cff969 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 20 Sep 2023 12:19:04 +0200 Subject: [PATCH 24/26] Fix updated param names in docs and test inputs --- docs/cluster.md | 8 ++++---- inputs/cluster/cluster.in | 2 +- inputs/cluster/magnetic_tower.in | 2 +- .../cluster_magnetic_tower/cluster_magnetic_tower.py | 2 +- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/cluster.md b/docs/cluster.md index 7e650cfa..b4d83983 100644 --- a/docs/cluster.md +++ b/docs/cluster.md @@ -449,7 +449,7 @@ The parameters $\alpha$ and $\ell$ can be changed with ``` potential_type = li -alpha = 20 +li_alpha = 20 l_scale = 0.001 ``` When injected as a fraction of @@ -498,9 +498,9 @@ with all other components being zero. ``` potential_type = donut -l_scale = 0.0005 # in code length -offset = 0.001 # in code length -thickness = 0.0001 # in code length +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 diff --git a/inputs/cluster/cluster.in b/inputs/cluster/cluster.in index fe1ff02b..e3192524 100644 --- a/inputs/cluster/cluster.in +++ b/inputs/cluster/cluster.in @@ -182,7 +182,7 @@ kinetic_jet_offset = 0.05 potential_type = li -alpha = 20 +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 15b5647d..44ba58a0 100644 --- a/inputs/cluster/magnetic_tower.in +++ b/inputs/cluster/magnetic_tower.in @@ -110,7 +110,7 @@ thermal_fraction = 0 potential_type = li -alpha = 20 +li_alpha = 20 l_scale = 0.01 initial_field = 0.12431560000204142 fixed_field_rate = 12.431560000204144 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}", From bd73195b10cb72287b4f3eb30751bf93414f90aa Mon Sep 17 00:00:00 2001 From: Forrest Wolfgang Glines Date: Thu, 21 Sep 2023 20:42:00 -0600 Subject: [PATCH 25/26] Removed commented alt. cooling terms --- src/hydro/hydro.cpp | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index d22dddef..6d0e087c 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -161,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); } @@ -176,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); } From f334d206f78ee2e8997f0c93390e361da99ed4f1 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 27 Sep 2023 18:35:02 +0200 Subject: [PATCH 26/26] Update pgen.md --- docs/pgen.md | 149 ++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 141 insertions(+), 8 deletions(-) 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.