From 14ea977bde791c5005f710621061c996f024d8d7 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 23 Aug 2023 15:51:14 -0400 Subject: [PATCH 01/19] 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 02/19] 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 03/19] 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 04/19] 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 05/19] 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 06/19] 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 07/19] 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 08/19] 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 6853b2db1828839acc3d92b3ecdd0899e2b063e7 Mon Sep 17 00:00:00 2001 From: Forrest Wolfgang Glines Date: Tue, 29 Aug 2023 16:59:43 -0600 Subject: [PATCH 09/19] Updated Parthenon to PR930 --- external/parthenon | 2 +- src/hydro/hydro.cpp | 6 ++-- src/pgen/cluster.cpp | 64 ++++++++++++++++++++---------------- src/pgen/cpaw.cpp | 6 ++-- src/pgen/field_loop.cpp | 6 ++-- src/pgen/linear_wave.cpp | 11 ++++--- src/pgen/linear_wave_mhd.cpp | 11 ++++--- src/pgen/turbulence.cpp | 14 ++++---- src/utils/few_modes_ft.cpp | 24 +++++++------- 9 files changed, 78 insertions(+), 66 deletions(-) diff --git a/external/parthenon b/external/parthenon index 5b6bb619..60d54786 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit 5b6bb61906f7c278f9724ee9f38e79dee8707098 +Subproject commit 60d5478693e26fa9b4836beb60ba7aeedfc608a4 diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index d22dddef..846adb9d 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -794,8 +794,8 @@ TaskStatus CalculateFluxes(std::shared_ptr> &md) { int il, iu, jl, ju, kl, ku; jl = jb.s, ju = jb.e, kl = kb.s, ku = kb.e; // TODO(pgrete): are these looop limits are likely too large for 2nd order - if (pmb->block_size.nx2 > 1) { - if (pmb->block_size.nx3 == 1) // 2D + if (pmb->block_size.nx(X2DIR) > 1) { + if (pmb->block_size.nx(X3DIR)== 1) // 2D jl = jb.s - 1, ju = jb.e + 1, kl = kb.s, ku = kb.e; else // 3D jl = jb.s - 1, ju = jb.e + 1, kl = kb.s - 1, ku = kb.e + 1; @@ -867,7 +867,7 @@ TaskStatus CalculateFluxes(std::shared_ptr> &md) { parthenon::ScratchPad2D::shmem_size(num_scratch_vars, nx1) * 3; // set the loop limits il = ib.s - 1, iu = ib.e + 1, kl = kb.s, ku = kb.e; - if (pmb->block_size.nx3 == 1) // 2D + if (pmb->block_size.nx(X3DIR) == 1) // 2D kl = kb.s, ku = kb.e; else // 3D kl = kb.s - 1, ku = kb.e + 1; diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index 96a781ba..de489177 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -59,8 +59,9 @@ 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"); @@ -78,11 +79,17 @@ void ClusterUnsplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, const auto &magnetic_tower = hydro_pkg->Param("magnetic_tower"); magnetic_tower.FixedFieldSrcTerm(md, beta_dt, tm); - const auto &snia_feedback = hydro_pkg->Param("snia_feedback"); - snia_feedback.FeedbackSrcTerm(md, beta_dt, tm); + //const auto &snia_feedback = hydro_pkg->Param("snia_feedback"); + //snia_feedback.FeedbackSrcTerm(md, beta_dt, tm); + + //ApplyClusterClips(md, tm, beta_dt); + + const auto &stellar_feedback = hydro_pkg->Param("stellar_feedback"); + stellar_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"); @@ -288,49 +295,50 @@ 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); /************************************************************ @@ -706,9 +714,9 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { MPI_SUM, MPI_COMM_WORLD)); #endif // MPI_PARALLEL - const auto Lx = pmesh->mesh_size.x1max - pmesh->mesh_size.x1min; - const auto Ly = pmesh->mesh_size.x2max - pmesh->mesh_size.x2min; - const auto Lz = pmesh->mesh_size.x3max - pmesh->mesh_size.x3min; + const auto Lx = pmesh->mesh_size.xmax(X1DIR) - pmesh->mesh_size.xmin(X1DIR); + const auto Ly = pmesh->mesh_size.xmax(X2DIR) - pmesh->mesh_size.xmin(X2DIR); + const auto Lz = pmesh->mesh_size.xmax(X3DIR) - pmesh->mesh_size.xmin(X3DIR); auto v_norm = std::sqrt(v2_sum / (Lx * Ly * Lz) / (SQR(sigma_v))); pmb->par_for( @@ -784,9 +792,9 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { MPI_SUM, MPI_COMM_WORLD)); #endif // MPI_PARALLEL - const auto Lx = pmesh->mesh_size.x1max - pmesh->mesh_size.x1min; - const auto Ly = pmesh->mesh_size.x2max - pmesh->mesh_size.x2min; - const auto Lz = pmesh->mesh_size.x3max - pmesh->mesh_size.x3min; + const auto Lx = pmesh->mesh_size.xmax(X1DIR) - pmesh->mesh_size.xmin(X1DIR); + const auto Ly = pmesh->mesh_size.xmax(X2DIR) - pmesh->mesh_size.xmin(X2DIR); + const auto Lz = pmesh->mesh_size.xmax(X3DIR) - pmesh->mesh_size.xmin(X3DIR); auto b_norm = std::sqrt(b2_sum / (Lx * Ly * Lz) / (SQR(sigma_b))); pmb->par_for( diff --git a/src/pgen/cpaw.cpp b/src/pgen/cpaw.cpp index c269c378..0dcefacb 100644 --- a/src/pgen/cpaw.cpp +++ b/src/pgen/cpaw.cpp @@ -34,6 +34,8 @@ namespace cpaw { using namespace parthenon::driver::prelude; +using namespace parthenon::package::prelude; + // Parameters which define initial solution -- made global so that they can be shared // with functions A1,2,3 which compute vector potentials Real den, pres, gm1, b_par, b_perp, v_perp, v_par; @@ -213,8 +215,8 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm) } // write errors - std::fprintf(pfile, "%d %d", mesh->mesh_size.nx1, mesh->mesh_size.nx2); - std::fprintf(pfile, " %d %d %e", mesh->mesh_size.nx3, tm.ncycle, rms_err); + std::fprintf(pfile, "%d %d", mesh->mesh_size.nx(X1DIR), mesh->mesh_size.nx(X2DIR)); + std::fprintf(pfile, " %d %d %e", mesh->mesh_size.nx(X3DIR), tm.ncycle, rms_err); std::fprintf(pfile, " %e %e %e %e", err[IDN], err[IM1], err[IM2], err[IM3]); std::fprintf(pfile, " %e", err[IEN]); std::fprintf(pfile, " %e %e %e", err[IB1], err[IB2], err[IB3]); diff --git a/src/pgen/field_loop.cpp b/src/pgen/field_loop.cpp index 66940d56..0ccdb69a 100644 --- a/src/pgen/field_loop.cpp +++ b/src/pgen/field_loop.cpp @@ -135,13 +135,13 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { int iprob = pin->GetInteger("problem/field_loop", "iprob"); Real ang_2, cos_a2(0.0), sin_a2(0.0), lambda(0.0); - Real x1size = pmb->pmy_mesh->mesh_size.x1max - pmb->pmy_mesh->mesh_size.x1min; - Real x2size = pmb->pmy_mesh->mesh_size.x2max - pmb->pmy_mesh->mesh_size.x2min; + Real x1size = pmb->pmy_mesh->mesh_size.xmax(X1DIR) - pmb->pmy_mesh->mesh_size.xmin(X1DIR); + Real x2size = pmb->pmy_mesh->mesh_size.xmax(X2DIR) - pmb->pmy_mesh->mesh_size.xmin(X2DIR); const bool two_d = pmb->pmy_mesh->ndim < 3; // for 2D sim set x3size to zero so that v_z is 0 below Real x3size = - two_d ? 0 : pmb->pmy_mesh->mesh_size.x3max - pmb->pmy_mesh->mesh_size.x3min; + two_d ? 0 : pmb->pmy_mesh->mesh_size.xmax(X3DIR) - pmb->pmy_mesh->mesh_size.xmin(X3DIR); // For (iprob=4) -- rotated cylinder in 3D -- set up rotation angle and wavelength if (iprob == 4) { diff --git a/src/pgen/linear_wave.cpp b/src/pgen/linear_wave.cpp index bc504b06..5480875e 100644 --- a/src/pgen/linear_wave.cpp +++ b/src/pgen/linear_wave.cpp @@ -32,6 +32,7 @@ namespace linear_wave { using namespace parthenon::driver::prelude; +using namespace parthenon::package::prelude; // TODO(pgrete) temp fix to address removal in Parthenon. Update when merging with MHD constexpr int NWAVE = 5; @@ -280,9 +281,9 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm) if (parthenon::Globals::my_rank == 0) { // normalize errors by number of cells const auto mesh_size = mesh->mesh_size; - const auto vol = (mesh_size.x1max - mesh_size.x1min) * - (mesh_size.x2max - mesh_size.x2min) * - (mesh_size.x3max - mesh_size.x3min); + const auto vol = (mesh_size.xmax(X1DIR) - mesh_size.xmin(X1DIR)) * + (mesh_size.xmax(X2DIR) - mesh_size.xmin(X2DIR)) * + (mesh_size.xmax(X3DIR) - mesh_size.xmin(X3DIR)); for (int i = 0; i < (NHYDRO + NFIELD); ++i) l1_err[i] = l1_err[i] / vol; // compute rms error @@ -320,8 +321,8 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm) } // write errors - std::fprintf(pfile, "%d %d", mesh_size.nx1, mesh_size.nx2); - std::fprintf(pfile, " %d %d", mesh_size.nx3, tm.ncycle); + std::fprintf(pfile, "%d %d", mesh_size.nx(X1DIR), mesh_size.nx(X2DIR)); + std::fprintf(pfile, " %d %d", mesh_size.nx(X2DIR), tm.ncycle); std::fprintf(pfile, " %e %e", rms_err, l1_err[IDN]); std::fprintf(pfile, " %e %e %e", l1_err[IM1], l1_err[IM2], l1_err[IM3]); std::fprintf(pfile, " %e", l1_err[IEN]); diff --git a/src/pgen/linear_wave_mhd.cpp b/src/pgen/linear_wave_mhd.cpp index 7f1b549b..67affbe7 100644 --- a/src/pgen/linear_wave_mhd.cpp +++ b/src/pgen/linear_wave_mhd.cpp @@ -32,6 +32,7 @@ namespace linear_wave_mhd { using namespace parthenon::driver::prelude; +using namespace parthenon::package::prelude; constexpr int NMHDWAVE = 7; // Parameters which define initial solution -- made global so that they can be shared @@ -300,9 +301,9 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm) if (parthenon::Globals::my_rank == 0) { // normalize errors by number of cells const auto mesh_size = mesh->mesh_size; - const auto vol = (mesh_size.x1max - mesh_size.x1min) * - (mesh_size.x2max - mesh_size.x2min) * - (mesh_size.x3max - mesh_size.x3min); + const auto vol = (mesh_size.xmax(X1DIR) - mesh_size.xmin(X1DIR)) * + (mesh_size.xmax(X2DIR) - mesh_size.xmin(X2DIR)) * + (mesh_size.xmax(X3DIR) - mesh_size.xmin(X3DIR)); for (int i = 0; i < (NGLMMHD); ++i) l1_err[i] = l1_err[i] / vol; // compute rms error @@ -342,8 +343,8 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm) } // write errors - std::fprintf(pfile, "%d %d", mesh_size.nx1, mesh_size.nx2); - std::fprintf(pfile, " %d %d", mesh_size.nx3, tm.ncycle); + std::fprintf(pfile, "%d %d", mesh_size.nx(X1DIR), mesh_size.nx(X2DIR)); + std::fprintf(pfile, " %d %d", mesh_size.nx(X3DIR), tm.ncycle); std::fprintf(pfile, " %e %e", rms_err, l1_err[IDN]); std::fprintf(pfile, " %e %e %e", l1_err[IM1], l1_err[IM2], l1_err[IM3]); std::fprintf(pfile, " %e", l1_err[IEN]); diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index 58882f52..9d0e96be 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -217,10 +217,10 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { const auto gm1 = pin->GetReal("hydro", "gamma") - 1.0; const auto p0 = pin->GetReal("problem/turbulence", "p0"); const auto rho0 = pin->GetReal("problem/turbulence", "rho0"); - const auto x3min = pmesh->mesh_size.x3min; - const auto Lx = pmesh->mesh_size.x1max - pmesh->mesh_size.x1min; - const auto Ly = pmesh->mesh_size.x2max - pmesh->mesh_size.x2min; - const auto Lz = pmesh->mesh_size.x3max - pmesh->mesh_size.x3min; + const auto x3min = pmesh->mesh_size.xmin(X3DIR); + const auto Lx = pmesh->mesh_size.xmax(X1DIR) - pmesh->mesh_size.xmin(X1DIR); + const auto Ly = pmesh->mesh_size.xmax(X2DIR) - pmesh->mesh_size.xmin(X2DIR); + const auto Lz = pmesh->mesh_size.xmax(X3DIR) - pmesh->mesh_size.xmin(X3DIR); const auto kz = 2.0 * M_PI / Lz; // already pack data here to get easy access to coords in kernels @@ -402,9 +402,9 @@ void Perturb(MeshData *md, const Real dt) { MPI_SUM, MPI_COMM_WORLD)); #endif // MPI_PARALLEL - const auto Lx = pmb->pmy_mesh->mesh_size.x1max - pmb->pmy_mesh->mesh_size.x1min; - const auto Ly = pmb->pmy_mesh->mesh_size.x2max - pmb->pmy_mesh->mesh_size.x2min; - const auto Lz = pmb->pmy_mesh->mesh_size.x3max - pmb->pmy_mesh->mesh_size.x3min; + const auto Lx = pmb->pmy_mesh->mesh_size.xmax(X1DIR) - pmb->pmy_mesh->mesh_size.xmin(X1DIR); + const auto Ly = pmb->pmy_mesh->mesh_size.xmax(X2DIR) - pmb->pmy_mesh->mesh_size.xmin(X2DIR); + const auto Lz = pmb->pmy_mesh->mesh_size.xmax(X3DIR) - pmb->pmy_mesh->mesh_size.xmin(X3DIR); const auto accel_rms = hydro_pkg->Param("turbulence/accel_rms"); auto norm = accel_rms / std::sqrt(sums[0] / (Lx * Ly * Lz)); diff --git a/src/utils/few_modes_ft.cpp b/src/utils/few_modes_ft.cpp index f477b211..e42edde2 100644 --- a/src/utils/few_modes_ft.cpp +++ b/src/utils/few_modes_ft.cpp @@ -106,18 +106,18 @@ void FewModesFT::SetPhases(MeshBlock *pmb, ParameterInput *pin) { "Few modes FT currently needs parthenon/mesh/pack_size=-1 " "to work because of global reductions.") - auto Lx1 = pm->mesh_size.x1max - pm->mesh_size.x1min; - auto Lx2 = pm->mesh_size.x2max - pm->mesh_size.x2min; - auto Lx3 = pm->mesh_size.x3max - pm->mesh_size.x3min; + const auto Lx1 = pm->mesh_size.xmax(X1DIR) - pm->mesh_size.xmin(X1DIR); + const auto Lx2 = pm->mesh_size.xmax(X2DIR) - pm->mesh_size.xmin(X2DIR); + const auto Lx3 = pm->mesh_size.xmax(X3DIR) - pm->mesh_size.xmin(X3DIR); // Adjust (logical) grid size at levels other than the root level. // This is required for simulation with mesh refinement so that the phases calculated // 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.nx(X1DIR) * std::pow(2, pmb->loc.level() - root_level); + auto gnx2 = pm->mesh_size.nx(X2DIR) * std::pow(2, pmb->loc.level() - root_level); + auto gnx3 = pm->mesh_size.nx(X3DIR) * std::pow(2, pmb->loc.level() - root_level); // Restriction should also be easily fixed, just need to double check transforms and // volume weighting everywhere @@ -126,13 +126,13 @@ void FewModesFT::SetPhases(MeshBlock *pmb, ParameterInput *pin) { "FMFT has only been tested with cubic meshes and constant " "dx/dy/dz. Remove this warning at your own risk.") - const auto nx1 = pmb->block_size.nx1; - const auto nx2 = pmb->block_size.nx2; - const auto nx3 = pmb->block_size.nx3; + const auto nx1 = pmb->block_size.nx(X1DIR); + const auto nx2 = pmb->block_size.nx(X2DIR); + const auto nx3 = pmb->block_size.nx(X3DIR); - 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.nx(X1DIR); + const auto gjs = pmb->loc.lx2() * pmb->block_size.nx(X2DIR); + const auto gks = pmb->loc.lx3() * pmb->block_size.nx(X3DIR); // make local ref to capure in lambda const auto num_modes = num_modes_; From 60551edf8a1ba07a137ccdb1f292fe8b33ea829b Mon Sep 17 00:00:00 2001 From: par-hermes Date: Tue, 29 Aug 2023 23:01:00 +0000 Subject: [PATCH 10/19] cpp-py-formatter --- src/hydro/hydro.cpp | 2 +- src/pgen/cluster.cpp | 49 +++++++++++++++++++---------------------- src/pgen/field_loop.cpp | 9 +++++--- src/pgen/turbulence.cpp | 9 +++++--- 4 files changed, 36 insertions(+), 33 deletions(-) diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index 846adb9d..b46db6b2 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -795,7 +795,7 @@ TaskStatus CalculateFluxes(std::shared_ptr> &md) { jl = jb.s, ju = jb.e, kl = kb.s, ku = kb.e; // TODO(pgrete): are these looop limits are likely too large for 2nd order if (pmb->block_size.nx(X2DIR) > 1) { - if (pmb->block_size.nx(X3DIR)== 1) // 2D + if (pmb->block_size.nx(X3DIR) == 1) // 2D jl = jb.s - 1, ju = jb.e + 1, kl = kb.s, ku = kb.e; else // 3D jl = jb.s - 1, ju = jb.e + 1, kl = kb.s - 1, ku = kb.e + 1; diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index de489177..d0f981d2 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"); @@ -79,17 +78,16 @@ void ClusterUnsplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, const auto &magnetic_tower = hydro_pkg->Param("magnetic_tower"); magnetic_tower.FixedFieldSrcTerm(md, beta_dt, tm); - //const auto &snia_feedback = hydro_pkg->Param("snia_feedback"); - //snia_feedback.FeedbackSrcTerm(md, beta_dt, tm); + // const auto &snia_feedback = hydro_pkg->Param("snia_feedback"); + // snia_feedback.FeedbackSrcTerm(md, beta_dt, tm); - //ApplyClusterClips(md, tm, beta_dt); + // ApplyClusterClips(md, tm, beta_dt); const auto &stellar_feedback = hydro_pkg->Param("stellar_feedback"); stellar_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"); @@ -295,50 +293,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/field_loop.cpp b/src/pgen/field_loop.cpp index 0ccdb69a..45eeb60e 100644 --- a/src/pgen/field_loop.cpp +++ b/src/pgen/field_loop.cpp @@ -135,13 +135,16 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { int iprob = pin->GetInteger("problem/field_loop", "iprob"); Real ang_2, cos_a2(0.0), sin_a2(0.0), lambda(0.0); - Real x1size = pmb->pmy_mesh->mesh_size.xmax(X1DIR) - pmb->pmy_mesh->mesh_size.xmin(X1DIR); - Real x2size = pmb->pmy_mesh->mesh_size.xmax(X2DIR) - pmb->pmy_mesh->mesh_size.xmin(X2DIR); + Real x1size = + pmb->pmy_mesh->mesh_size.xmax(X1DIR) - pmb->pmy_mesh->mesh_size.xmin(X1DIR); + Real x2size = + pmb->pmy_mesh->mesh_size.xmax(X2DIR) - pmb->pmy_mesh->mesh_size.xmin(X2DIR); const bool two_d = pmb->pmy_mesh->ndim < 3; // for 2D sim set x3size to zero so that v_z is 0 below Real x3size = - two_d ? 0 : pmb->pmy_mesh->mesh_size.xmax(X3DIR) - pmb->pmy_mesh->mesh_size.xmin(X3DIR); + two_d ? 0 + : pmb->pmy_mesh->mesh_size.xmax(X3DIR) - pmb->pmy_mesh->mesh_size.xmin(X3DIR); // For (iprob=4) -- rotated cylinder in 3D -- set up rotation angle and wavelength if (iprob == 4) { diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index 9d0e96be..2c9ff9a7 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -402,9 +402,12 @@ void Perturb(MeshData *md, const Real dt) { MPI_SUM, MPI_COMM_WORLD)); #endif // MPI_PARALLEL - const auto Lx = pmb->pmy_mesh->mesh_size.xmax(X1DIR) - pmb->pmy_mesh->mesh_size.xmin(X1DIR); - const auto Ly = pmb->pmy_mesh->mesh_size.xmax(X2DIR) - pmb->pmy_mesh->mesh_size.xmin(X2DIR); - const auto Lz = pmb->pmy_mesh->mesh_size.xmax(X3DIR) - pmb->pmy_mesh->mesh_size.xmin(X3DIR); + const auto Lx = + pmb->pmy_mesh->mesh_size.xmax(X1DIR) - pmb->pmy_mesh->mesh_size.xmin(X1DIR); + const auto Ly = + pmb->pmy_mesh->mesh_size.xmax(X2DIR) - pmb->pmy_mesh->mesh_size.xmin(X2DIR); + const auto Lz = + pmb->pmy_mesh->mesh_size.xmax(X3DIR) - pmb->pmy_mesh->mesh_size.xmin(X3DIR); const auto accel_rms = hydro_pkg->Param("turbulence/accel_rms"); auto norm = accel_rms / std::sqrt(sums[0] / (Lx * Ly * Lz)); From 2bdf1f8cd6aa0d57fa745c15603c51392aea96bf Mon Sep 17 00:00:00 2001 From: Forrest Wolfgang Glines Date: Sun, 10 Sep 2023 11:54:48 -0600 Subject: [PATCH 11/19] 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 12/19] 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 13/19] 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 14/19] 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 6004203365cf1b719d451746c5f7a74a8186a2d1 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Tue, 26 Sep 2023 11:48:27 +0200 Subject: [PATCH 15/19] Bump parth to current dev --- external/parthenon | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/parthenon b/external/parthenon index 60d54786..39ad758a 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit 60d5478693e26fa9b4836beb60ba7aeedfc608a4 +Subproject commit 39ad758a2c94ac5552200809aa39076a11bb8810 From f334d206f78ee2e8997f0c93390e361da99ed4f1 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 27 Sep 2023 18:35:02 +0200 Subject: [PATCH 16/19] 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. From fe532fd13e554c4044e90e5071174d40e3ec0821 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Tue, 17 Oct 2023 15:34:25 +0200 Subject: [PATCH 17/19] Bump parth to include analysis outputs --- external/parthenon | 2 +- src/pgen/cluster.cpp | 8 ++------ 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/external/parthenon b/external/parthenon index 39ad758a..8af38e42 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit 39ad758a2c94ac5552200809aa39076a11bb8810 +Subproject commit 8af38e420c56baf9c86c4c68f0787f10e8824d9d diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index a26b1da6..e1017415 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -78,13 +78,9 @@ void ClusterUnsplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, const auto &magnetic_tower = hydro_pkg->Param("magnetic_tower"); magnetic_tower.FixedFieldSrcTerm(md, beta_dt, tm); - // const auto &snia_feedback = hydro_pkg->Param("snia_feedback"); - // snia_feedback.FeedbackSrcTerm(md, beta_dt, tm); + const auto &snia_feedback = hydro_pkg->Param("snia_feedback"); + snia_feedback.FeedbackSrcTerm(md, beta_dt, tm); - // ApplyClusterClips(md, tm, beta_dt); - - const auto &stellar_feedback = hydro_pkg->Param("stellar_feedback"); - stellar_feedback.FeedbackSrcTerm(md, beta_dt, tm); }; void ClusterSplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, const Real dt) { From d83ef781d62cc49ba7de5c19648542d08a842acd Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Tue, 17 Oct 2023 15:37:16 +0200 Subject: [PATCH 18/19] Add bmag calc for cluster --- src/pgen/cluster.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index e1017415..3a0cfb9a 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -25,6 +25,7 @@ #include // c_str() // Parthenon headers +#include "Kokkos_MathematicalFunctions.hpp" #include "kokkos_abstraction.hpp" #include "mesh/domain.hpp" #include "mesh/mesh.hpp" @@ -80,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) { @@ -360,6 +360,9 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd // plasma beta hydro_pkg->AddField("plasma_beta", m); + + // plasma beta + hydro_pkg->AddField("B_mag", m); } /************************************************************ @@ -885,6 +888,7 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { if (pkg->Param("fluid") == Fluid::glmmhd) { auto &plasma_beta = data->Get("plasma_beta").data; auto &mach_alfven = data->Get("mach_alfven").data; + auto &b_mag = data->Get("B_mag").data; pmb->par_for( "Cluster::UserWorkBeforeOutput::MHD", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, @@ -897,6 +901,8 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { const Real Bz = prim(IB3, k, j, i); const Real B2 = (SQR(Bx) + SQR(By) + SQR(Bz)); + b_mag(k, j, i) = Kokkos::sqrt(B2); + // compute Alfven mach number const Real v_A = std::sqrt(B2 / rho); const Real c_s = std::sqrt(gam * P / rho); // ideal gas EOS From 2eb01991893cc3a52c68a292678fa371c4f4b2c6 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Fri, 20 Oct 2023 09:31:25 +0200 Subject: [PATCH 19/19] Bump hist parthenon ver --- external/parthenon | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/parthenon b/external/parthenon index 8af38e42..725deed2 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit 8af38e420c56baf9c86c4c68f0787f10e8824d9d +Subproject commit 725deed20ca2dff61d59f3a18d41772d8c4a7644