From f15c4ac1d4de9c32aa86df402e598a7dfcb8b5dd Mon Sep 17 00:00:00 2001 From: epanini Date: Wed, 18 Dec 2024 19:42:46 +0100 Subject: [PATCH 1/6] agn_triggering.cpp compute Inflow rate for total mass and cold gas mass --- src/pgen/cluster/agn_triggering.cpp | 129 +++++++++++++++++++++------- 1 file changed, 98 insertions(+), 31 deletions(-) diff --git a/src/pgen/cluster/agn_triggering.cpp b/src/pgen/cluster/agn_triggering.cpp index 242de4fa..013df6d9 100644 --- a/src/pgen/cluster/agn_triggering.cpp +++ b/src/pgen/cluster/agn_triggering.cpp @@ -66,9 +66,12 @@ AGNTriggering::AGNTriggering(parthenon::ParameterInput *pin, accretion_cfl_(pin->GetOrAddReal(block, "accretion_cfl", 1e-1)), remove_accreted_mass_(pin->GetOrAddBoolean(block, "removed_accreted_mass", true)), write_to_file_(pin->GetOrAddBoolean(block, "write_to_file", false)), + inflow_cold_(0), + inflow_tot_(0), triggering_filename_( pin->GetOrAddString(block, "triggering_filename", "agn_triggering.dat")) { + const auto units = hydro_pkg->Param("units"); const parthenon::Real He_mass_fraction = pin->GetReal("hydro", "He_mass_fraction"); const parthenon::Real H_mass_fraction = 1.0 - He_mass_fraction; @@ -85,6 +88,15 @@ AGNTriggering::AGNTriggering(parthenon::ParameterInput *pin, switch (triggering_mode_) { case AGNTriggeringMode::COLD_GAS: { hydro_pkg->AddParam("agn_triggering_cold_mass", 0, Params::Mutability::Restart); + + hydro_pkg->AddParam("agn_triggering_total_mass", 0, + Params::Mutability::Restart); + hydro_pkg->AddParam("agn_triggering_prev_cold_mass", 0, + Params::Mutability::Restart); + hydro_pkg->AddParam("agn_triggering_prev_total_mass", 0, + Params::Mutability::Restart); + + break; } case AGNTriggeringMode::BOOSTED_BONDI: @@ -122,11 +134,14 @@ AGNTriggering::AGNTriggering(parthenon::ParameterInput *pin, // and simultaneously remove cold gas (updating conserveds and primitives) template void AGNTriggering::ReduceColdMass(parthenon::Real &cold_mass, + parthenon::Real &total_mass, parthenon::MeshData *md, - const parthenon::Real dt, const EOS eos) const { + const parthenon::Real dt, const EOS eos) { using parthenon::IndexDomain; using parthenon::IndexRange; using parthenon::Real; + Real md_cold_mass = 0; + Real md_total_mass = 0; auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); @@ -153,14 +168,14 @@ void AGNTriggering::ReduceColdMass(parthenon::Real &cold_mass, const bool remove_accreted_mass = remove_accreted_mass_; - Real md_cold_mass = 0; + parthenon::par_reduce( parthenon::loop_pattern_mdrange_tag, "AGNTriggering::ReduceColdGas", 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 &team_cold_mass) { + Real &team_cold_mass, Real &team_total_mass) { auto &cons = cons_pack(b); auto &prim = prim_pack(b); const auto &coords = cons_pack.GetCoords(b); @@ -168,10 +183,11 @@ void AGNTriggering::ReduceColdMass(parthenon::Real &cold_mass, const parthenon::Real r2 = pow(coords.Xc<1>(i), 2) + pow(coords.Xc<2>(j), 2) + pow(coords.Xc<3>(k), 2); if (r2 < accretion_radius2) { - + const Real cell_total_mass = prim(IDN, k, j, i) * coords.CellVolume(k, j, i); + team_total_mass += cell_total_mass; const Real temp = mean_molecular_mass_by_kb * prim(IPR, k, j, i) / prim(IDN, k, j, i); - + if (temp <= cold_temp_thresh) { const Real cell_cold_mass = prim(IDN, k, j, i) * coords.CellVolume(k, j, i); @@ -193,8 +209,11 @@ void AGNTriggering::ReduceColdMass(parthenon::Real &cold_mass, } } }, - Kokkos::Sum(md_cold_mass)); + + Kokkos::Sum(md_cold_mass), Kokkos::Sum(md_total_mass)); cold_mass += md_cold_mass; + total_mass += md_total_mass; + } // Compute Mass-weighted total density, velocity, and sound speed and total mass @@ -234,25 +253,27 @@ void AGNTriggering::ReduceBondiTriggeringQuantities( KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i, Real <otal_mass_red, Real &lmass_weighted_density_red, Real &lmass_weighted_velocity_red, Real &lmass_weighted_cs_red) { - auto &prim = prim_pack(b); - const auto &coords = prim_pack.GetCoords(b); - const parthenon::Real r2 = + if (k >= kb.s && k <= kb.e && j >= jb.s && j <= jb.e && i >= ib.s && i <= ib.e) { + auto &prim = prim_pack(b); + const auto &coords = prim_pack.GetCoords(b); + const parthenon::Real r2 = pow(coords.Xc<1>(i), 2) + pow(coords.Xc<2>(j), 2) + pow(coords.Xc<3>(k), 2); - if (r2 < accretion_radius2) { - const Real cell_mass = prim(IDN, k, j, i) * coords.CellVolume(k, j, i); + if (r2 < accretion_radius2) { + const Real cell_mass = prim(IDN, k, j, i) * coords.CellVolume(k, j, i); - const Real cell_mass_weighted_density = cell_mass * prim(IDN, k, j, i); - const Real cell_mass_weighted_velocity = + const Real cell_mass_weighted_density = cell_mass * prim(IDN, k, j, i); + const Real cell_mass_weighted_velocity = cell_mass * sqrt(pow(prim(IV1, k, j, i), 2) + pow(prim(IV2, k, j, i), 2) + pow(prim(IV3, k, j, i), 2)); - const Real cell_mass_weighted_cs = + const Real cell_mass_weighted_cs = cell_mass * sqrt(gamma * prim(IPR, k, j, i) / prim(IDN, k, j, i)); - ltotal_mass_red += cell_mass; - lmass_weighted_density_red += cell_mass_weighted_density; - lmass_weighted_velocity_red += cell_mass_weighted_velocity; - lmass_weighted_cs_red += cell_mass_weighted_cs; - } + ltotal_mass_red += cell_mass; + lmass_weighted_density_red += cell_mass_weighted_density; + lmass_weighted_velocity_red += cell_mass_weighted_velocity; + lmass_weighted_cs_red += cell_mass_weighted_cs; + } + } }, total_mass_red, mass_weighted_density_red, mass_weighted_velocity_red, mass_weighted_cs_red); @@ -274,7 +295,7 @@ void AGNTriggering::RemoveBondiAccretedGas(parthenon::MeshData using parthenon::Real; auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); - + // Grab some necessary variables // FIXME(forrestglines) When reductions are called, is `prim` up to date? const auto &prim_pack = md->PackVariables(std::vector{"prim"}); @@ -361,16 +382,20 @@ AGNTriggering::GetAccretionRate(parthenon::StateDescriptor *hydro_pkg) const { return 0; } } - return 0; + //return 0; } parthenon::TaskStatus AGNTriggeringResetTriggering(parthenon::StateDescriptor *hydro_pkg) { - const auto &agn_triggering = hydro_pkg->Param("agn_triggering"); - + auto &agn_triggering = hydro_pkg->Param("agn_triggering"); + hydro_pkg->UpdateParam("agn_triggering_prev_cold_mass", hydro_pkg->Param("agn_triggering_cold_mass")); + hydro_pkg->UpdateParam("agn_triggering_prev_total_mass", hydro_pkg->Param("agn_triggering_total_mass")); switch (agn_triggering.triggering_mode_) { case AGNTriggeringMode::COLD_GAS: { hydro_pkg->UpdateParam("agn_triggering_cold_mass", 0); + hydro_pkg->UpdateParam("agn_triggering_total_mass", 0); + + break; } case AGNTriggeringMode::BOOSTED_BONDI: @@ -393,24 +418,28 @@ AGNTriggeringReduceTriggering(parthenon::MeshData *md, const parthenon::Real dt) { auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); - const auto &agn_triggering = hydro_pkg->Param("agn_triggering"); - + auto &agn_triggering = hydro_pkg->Param("agn_triggering"); + const parthenon::Real pre_accretion_cold_mass = hydro_pkg->Param("agn_triggering_cold_mass"); + const parthenon::Real pre_accretion_total_mass = hydro_pkg->Param("agn_triggering_total_mass"); switch (agn_triggering.triggering_mode_) { case AGNTriggeringMode::COLD_GAS: { Real cold_mass = hydro_pkg->Param("agn_triggering_cold_mass"); + Real total_mass = hydro_pkg->Param("agn_triggering_total_mass"); + auto fluid = hydro_pkg->Param("fluid"); if (fluid == Fluid::euler) { - agn_triggering.ReduceColdMass(cold_mass, md, dt, + agn_triggering.ReduceColdMass(cold_mass, total_mass,md, dt, hydro_pkg->Param("eos")); } else if (fluid == Fluid::glmmhd) { - agn_triggering.ReduceColdMass(cold_mass, md, dt, + agn_triggering.ReduceColdMass(cold_mass, total_mass,md, dt, hydro_pkg->Param("eos")); } else { PARTHENON_FAIL("AGNTriggeringReduceTriggeringQuantities: Unknown EOS"); } hydro_pkg->UpdateParam("agn_triggering_cold_mass", cold_mass); + hydro_pkg->UpdateParam("agn_triggering_total_mass", total_mass); break; } case AGNTriggeringMode::BOOSTED_BONDI: @@ -437,13 +466,22 @@ AGNTriggeringReduceTriggering(parthenon::MeshData *md, break; } } + if (agn_triggering.triggering_mode_ == AGNTriggeringMode::COLD_GAS) { + auto &agn_triggering = hydro_pkg->Param("agn_triggering"); // Make it non-const + const parthenon::Real cold_mass = hydro_pkg->Param("agn_triggering_cold_mass"); + const parthenon::Real total_mass = hydro_pkg->Param("agn_triggering_total_mass"); + + hydro_pkg->UpdateParam("agn_triggering_prev_cold_mass", cold_mass); //Move from AGNTriggeringFinalizeTriggering + hydro_pkg->UpdateParam("agn_triggering_prev_total_mass", total_mass); + } + hydro_pkg->UpdateParam("agn_triggering", agn_triggering); return TaskStatus::complete; } parthenon::TaskStatus AGNTriggeringMPIReduceTriggering(parthenon::StateDescriptor *hydro_pkg) { #ifdef MPI_PARALLEL - const auto &agn_triggering = hydro_pkg->Param("agn_triggering"); + auto &agn_triggering = hydro_pkg->Param("agn_triggering"); switch (agn_triggering.triggering_mode_) { case AGNTriggeringMode::COLD_GAS: { @@ -451,6 +489,11 @@ AGNTriggeringMPIReduceTriggering(parthenon::StateDescriptor *hydro_pkg) { PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &accretion_rate, 1, MPI_PARTHENON_REAL, MPI_SUM, MPI_COMM_WORLD)); hydro_pkg->UpdateParam("agn_triggering_cold_mass", accretion_rate); + + Real total_mass = hydro_pkg->Param("agn_triggering_total_mass"); + PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &total_mass, 1, + MPI_PARTHENON_REAL, MPI_SUM, MPI_COMM_WORLD)); + hydro_pkg->UpdateParam("agn_triggering_total_mass", total_mass); break; } case AGNTriggeringMode::BOOSTED_BONDI: @@ -485,15 +528,34 @@ parthenon::TaskStatus AGNTriggeringFinalizeTriggering(parthenon::MeshData *md, const parthenon::SimTime &tm) { auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); - const auto &agn_triggering = hydro_pkg->Param("agn_triggering"); + auto &agn_triggering = hydro_pkg->Param("agn_triggering"); + parthenon::Real cold_mass = 0.0; // Declare outside the if block + parthenon::Real total_mass = 0.0; // Declare outside the if block + if (agn_triggering.triggering_mode_ == AGNTriggeringMode::COLD_GAS) { + cold_mass = hydro_pkg->Param("agn_triggering_cold_mass"); // Assign here + total_mass = hydro_pkg->Param("agn_triggering_total_mass"); + + + agn_triggering.inflow_cold_ = (cold_mass - hydro_pkg->Param("agn_triggering_prev_cold_mass")) / tm.dt; //Now this is the only place to calculate inflow rates + agn_triggering.inflow_tot_ = (total_mass - hydro_pkg->Param("agn_triggering_prev_total_mass")) / tm.dt; //Now this is the only place to calculate inflow rates + + } + hydro_pkg->UpdateParam("agn_triggering", agn_triggering); + + hydro_pkg->UpdateParam("agn_triggering_prev_cold_mass", cold_mass); + hydro_pkg->UpdateParam("agn_triggering_prev_total_mass", total_mass); // Append quantities to file if applicable if (agn_triggering.write_to_file_ && parthenon::Globals::my_rank == 0) { std::ofstream triggering_file; triggering_file.open(agn_triggering.triggering_filename_, std::ofstream::app); - triggering_file << tm.time << " " << tm.dt << " " - << agn_triggering.GetAccretionRate(hydro_pkg.get()) << " "; + triggering_file << tm.time << " | " << tm.dt << " | " + << agn_triggering.GetAccretionRate(hydro_pkg.get()) << " | " + << cold_mass << " | " + << total_mass << " | " + << agn_triggering.inflow_cold_ << " | " + << agn_triggering.inflow_tot_ << " "; switch (agn_triggering.triggering_mode_) { case AGNTriggeringMode::COLD_GAS: { @@ -580,6 +642,11 @@ AGNTriggering::EstimateTimeStep(parthenon::MeshData *md) const return std::numeric_limits::max(); } } + //return std::numeric_limits::max(); +} + +} // namespace cluster + return std::numeric_limits::max(); } From 993b93693496612d2acf57514ba60f80751f45c5 Mon Sep 17 00:00:00 2001 From: epanini Date: Wed, 18 Dec 2024 19:44:13 +0100 Subject: [PATCH 2/6] agn_triggering.hpp modified initialization for inflow rate --- src/pgen/cluster/agn_triggering.hpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/pgen/cluster/agn_triggering.hpp b/src/pgen/cluster/agn_triggering.hpp index 8236a291..4e6c711e 100644 --- a/src/pgen/cluster/agn_triggering.hpp +++ b/src/pgen/cluster/agn_triggering.hpp @@ -15,7 +15,7 @@ #include #include #include - +#include // AthenaPK headers #include "../../units.hpp" #include "jet_coords.hpp" @@ -34,9 +34,12 @@ class AGNTriggering { private: const parthenon::Real gamma_; parthenon::Real mean_molecular_mass_; - + public: const AGNTriggeringMode triggering_mode_; + + parthenon::Real inflow_cold_=0.0; // Cold mass inflow rate + parthenon::Real inflow_tot_=0.0; // Total gas mass inflow rate const parthenon::Real accretion_radius_; @@ -71,9 +74,9 @@ class AGNTriggering { // Compute Cold gas accretion rate within the accretion radius for cold gas triggering // and simultaneously remove cold gas (updating conserveds and primitives) template - void ReduceColdMass(parthenon::Real &cold_mass, + void ReduceColdMass(parthenon::Real &cold_mass,parthenon::Real &total_mass, parthenon::MeshData *md, const parthenon::Real dt, - const EOS eos) const; + const EOS eos) ; // Compute Mass-weighted total density, velocity, and sound speed and total mass // for Bondi accretion From 9ad21139251ecc44baac1dcf5bdf58030ff8228d Mon Sep 17 00:00:00 2001 From: epanini Date: Mon, 30 Dec 2024 18:29:04 +0100 Subject: [PATCH 3/6] Update agn_triggering.hpp --- src/pgen/cluster/agn_triggering.hpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/pgen/cluster/agn_triggering.hpp b/src/pgen/cluster/agn_triggering.hpp index 4e6c711e..51c1d3db 100644 --- a/src/pgen/cluster/agn_triggering.hpp +++ b/src/pgen/cluster/agn_triggering.hpp @@ -37,9 +37,7 @@ class AGNTriggering { public: const AGNTriggeringMode triggering_mode_; - - parthenon::Real inflow_cold_=0.0; // Cold mass inflow rate - parthenon::Real inflow_tot_=0.0; // Total gas mass inflow rate + const parthenon::Real accretion_radius_; @@ -74,9 +72,9 @@ class AGNTriggering { // Compute Cold gas accretion rate within the accretion radius for cold gas triggering // and simultaneously remove cold gas (updating conserveds and primitives) template - void ReduceColdMass(parthenon::Real &cold_mass,parthenon::Real &total_mass, + void ReduceColdMass(parthenon::Real &cold_mass, parthenon::Real &total_mass, parthenon::MeshData *md, const parthenon::Real dt, - const EOS eos) ; + const EOS eos) const; // Compute Mass-weighted total density, velocity, and sound speed and total mass // for Bondi accretion From 74a8f8f901a705a2737ee0ef83646efd6a774d87 Mon Sep 17 00:00:00 2001 From: epanini Date: Mon, 30 Dec 2024 18:30:33 +0100 Subject: [PATCH 4/6] Update agn_triggering.cpp --- src/pgen/cluster/agn_triggering.cpp | 63 +++++++++++++---------------- 1 file changed, 27 insertions(+), 36 deletions(-) diff --git a/src/pgen/cluster/agn_triggering.cpp b/src/pgen/cluster/agn_triggering.cpp index 013df6d9..881805dd 100644 --- a/src/pgen/cluster/agn_triggering.cpp +++ b/src/pgen/cluster/agn_triggering.cpp @@ -66,8 +66,6 @@ AGNTriggering::AGNTriggering(parthenon::ParameterInput *pin, accretion_cfl_(pin->GetOrAddReal(block, "accretion_cfl", 1e-1)), remove_accreted_mass_(pin->GetOrAddBoolean(block, "removed_accreted_mass", true)), write_to_file_(pin->GetOrAddBoolean(block, "write_to_file", false)), - inflow_cold_(0), - inflow_tot_(0), triggering_filename_( pin->GetOrAddString(block, "triggering_filename", "agn_triggering.dat")) { @@ -136,7 +134,7 @@ template void AGNTriggering::ReduceColdMass(parthenon::Real &cold_mass, parthenon::Real &total_mass, parthenon::MeshData *md, - const parthenon::Real dt, const EOS eos) { + const parthenon::Real dt, const EOS eos) const { using parthenon::IndexDomain; using parthenon::IndexRange; using parthenon::Real; @@ -184,7 +182,12 @@ void AGNTriggering::ReduceColdMass(parthenon::Real &cold_mass, pow(coords.Xc<1>(i), 2) + pow(coords.Xc<2>(j), 2) + pow(coords.Xc<3>(k), 2); if (r2 < accretion_radius2) { const Real cell_total_mass = prim(IDN, k, j, i) * coords.CellVolume(k, j, i); - team_total_mass += cell_total_mass; + if (k >= int_kb.s && k <= int_kb.e && j >= int_jb.s && j <= int_jb.e && + i >= int_ib.s && i <= int_ib.e) { + // Only reduce the cold gas that exists on the interior grid + team_total_mass += cell_total_mass; + } + const Real temp = mean_molecular_mass_by_kb * prim(IPR, k, j, i) / prim(IDN, k, j, i); @@ -382,7 +385,7 @@ AGNTriggering::GetAccretionRate(parthenon::StateDescriptor *hydro_pkg) const { return 0; } } - //return 0; + return 0; } parthenon::TaskStatus @@ -418,7 +421,7 @@ AGNTriggeringReduceTriggering(parthenon::MeshData *md, const parthenon::Real dt) { auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); - auto &agn_triggering = hydro_pkg->Param("agn_triggering"); + const auto &agn_triggering = hydro_pkg->Param("agn_triggering"); const parthenon::Real pre_accretion_cold_mass = hydro_pkg->Param("agn_triggering_cold_mass"); const parthenon::Real pre_accretion_total_mass = hydro_pkg->Param("agn_triggering_total_mass"); switch (agn_triggering.triggering_mode_) { @@ -466,15 +469,7 @@ AGNTriggeringReduceTriggering(parthenon::MeshData *md, break; } } - if (agn_triggering.triggering_mode_ == AGNTriggeringMode::COLD_GAS) { - auto &agn_triggering = hydro_pkg->Param("agn_triggering"); // Make it non-const - const parthenon::Real cold_mass = hydro_pkg->Param("agn_triggering_cold_mass"); - const parthenon::Real total_mass = hydro_pkg->Param("agn_triggering_total_mass"); - - hydro_pkg->UpdateParam("agn_triggering_prev_cold_mass", cold_mass); //Move from AGNTriggeringFinalizeTriggering - hydro_pkg->UpdateParam("agn_triggering_prev_total_mass", total_mass); - } - hydro_pkg->UpdateParam("agn_triggering", agn_triggering); + return TaskStatus::complete; } @@ -484,16 +479,15 @@ AGNTriggeringMPIReduceTriggering(parthenon::StateDescriptor *hydro_pkg) { auto &agn_triggering = hydro_pkg->Param("agn_triggering"); switch (agn_triggering.triggering_mode_) { case AGNTriggeringMode::COLD_GAS: { - - Real accretion_rate = hydro_pkg->Param("agn_triggering_cold_mass"); - PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &accretion_rate, 1, + Real quantities[] ={ + hydro_pkg->Param("agn_triggering_cold_mass"), + hydro_pkg->Param("agn_triggering_total_mass"), + }; + + PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &quantities, 2, MPI_PARTHENON_REAL, MPI_SUM, MPI_COMM_WORLD)); - hydro_pkg->UpdateParam("agn_triggering_cold_mass", accretion_rate); - - Real total_mass = hydro_pkg->Param("agn_triggering_total_mass"); - PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &total_mass, 1, - MPI_PARTHENON_REAL, MPI_SUM, MPI_COMM_WORLD)); - hydro_pkg->UpdateParam("agn_triggering_total_mass", total_mass); + hydro_pkg->UpdateParam("agn_triggering_cold_mass", quantities[0]); + hydro_pkg->UpdateParam("agn_triggering_total_mass", quantities[1]); break; } case AGNTriggeringMode::BOOSTED_BONDI: @@ -529,15 +523,17 @@ AGNTriggeringFinalizeTriggering(parthenon::MeshData *md, const parthenon::SimTime &tm) { auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); auto &agn_triggering = hydro_pkg->Param("agn_triggering"); - parthenon::Real cold_mass = 0.0; // Declare outside the if block - parthenon::Real total_mass = 0.0; // Declare outside the if block + parthenon::Real cold_mass = 0.0; + parthenon::Real total_mass = 0.0; + parthenon::Real inflow_cold = 0.0; + parthenon::Real inflow_tot = 0.0; if (agn_triggering.triggering_mode_ == AGNTriggeringMode::COLD_GAS) { - cold_mass = hydro_pkg->Param("agn_triggering_cold_mass"); // Assign here + cold_mass = hydro_pkg->Param("agn_triggering_cold_mass"); total_mass = hydro_pkg->Param("agn_triggering_total_mass"); - agn_triggering.inflow_cold_ = (cold_mass - hydro_pkg->Param("agn_triggering_prev_cold_mass")) / tm.dt; //Now this is the only place to calculate inflow rates - agn_triggering.inflow_tot_ = (total_mass - hydro_pkg->Param("agn_triggering_prev_total_mass")) / tm.dt; //Now this is the only place to calculate inflow rates + inflow_cold = (cold_mass - hydro_pkg->Param("agn_triggering_prev_cold_mass")) / tm.dt; + inflow_tot = (total_mass - hydro_pkg->Param("agn_triggering_prev_total_mass")) / tm.dt; } hydro_pkg->UpdateParam("agn_triggering", agn_triggering); @@ -554,8 +550,8 @@ AGNTriggeringFinalizeTriggering(parthenon::MeshData *md, << agn_triggering.GetAccretionRate(hydro_pkg.get()) << " | " << cold_mass << " | " << total_mass << " | " - << agn_triggering.inflow_cold_ << " | " - << agn_triggering.inflow_tot_ << " "; + << inflow_cold << " | " + << inflow_tot << " "; switch (agn_triggering.triggering_mode_) { case AGNTriggeringMode::COLD_GAS: { @@ -642,11 +638,6 @@ AGNTriggering::EstimateTimeStep(parthenon::MeshData *md) const return std::numeric_limits::max(); } } - //return std::numeric_limits::max(); -} - -} // namespace cluster - return std::numeric_limits::max(); } From a198d1ca89439c6d98b3be2a4ee18319efc2c1a3 Mon Sep 17 00:00:00 2001 From: epanini Date: Tue, 7 Jan 2025 10:18:34 +0100 Subject: [PATCH 5/6] Update agn_triggering.cpp Updated function Reducecoldgas with Kokkos parralel_reduce intead of Parthenon --- src/pgen/cluster/agn_triggering.cpp | 30 +++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/src/pgen/cluster/agn_triggering.cpp b/src/pgen/cluster/agn_triggering.cpp index 881805dd..769d3e3e 100644 --- a/src/pgen/cluster/agn_triggering.cpp +++ b/src/pgen/cluster/agn_triggering.cpp @@ -168,24 +168,26 @@ void AGNTriggering::ReduceColdMass(parthenon::Real &cold_mass, - parthenon::par_reduce( - parthenon::loop_pattern_mdrange_tag, "AGNTriggering::ReduceColdGas", - 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 &team_cold_mass, Real &team_total_mass) { - auto &cons = cons_pack(b); - auto &prim = prim_pack(b); - const auto &coords = cons_pack.GetCoords(b); + Kokkos::parallel_reduce( + "AGNTriggering::ReduceColdGas", + 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 &team_cold_mass, Real &team_total_mass) { + auto &cons = cons_pack(b); + auto &prim = prim_pack(b); + const auto &coords = cons_pack.GetCoords(b); - const parthenon::Real r2 = + const parthenon::Real r2 = pow(coords.Xc<1>(i), 2) + pow(coords.Xc<2>(j), 2) + pow(coords.Xc<3>(k), 2); - if (r2 < accretion_radius2) { - const Real cell_total_mass = prim(IDN, k, j, i) * coords.CellVolume(k, j, i); - if (k >= int_kb.s && k <= int_kb.e && j >= int_jb.s && j <= int_jb.e && + if (r2 < accretion_radius2) { + const Real cell_total_mass = prim(IDN, k, j, i) * coords.CellVolume(k, j, i); + if (k >= int_kb.s && k <= int_kb.e && j >= int_jb.s && j <= int_jb.e && i >= int_ib.s && i <= int_ib.e) { // Only reduce the cold gas that exists on the interior grid - team_total_mass += cell_total_mass; + team_total_mass += cell_total_mass; } const Real temp = From c9ef64e6136d09a61c252f56427293123696080f Mon Sep 17 00:00:00 2001 From: epanini Date: Tue, 7 Jan 2025 10:21:33 +0100 Subject: [PATCH 6/6] Update agn_triggering.cpp deleted the updating of agn_triggering --- src/pgen/cluster/agn_triggering.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pgen/cluster/agn_triggering.cpp b/src/pgen/cluster/agn_triggering.cpp index 769d3e3e..a1562565 100644 --- a/src/pgen/cluster/agn_triggering.cpp +++ b/src/pgen/cluster/agn_triggering.cpp @@ -538,7 +538,7 @@ AGNTriggeringFinalizeTriggering(parthenon::MeshData *md, inflow_tot = (total_mass - hydro_pkg->Param("agn_triggering_prev_total_mass")) / tm.dt; } - hydro_pkg->UpdateParam("agn_triggering", agn_triggering); + hydro_pkg->UpdateParam("agn_triggering_prev_cold_mass", cold_mass); hydro_pkg->UpdateParam("agn_triggering_prev_total_mass", total_mass);