Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

AGN triggering_inflow rate #133

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open

AGN triggering_inflow rate #133

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
f15c4ac
agn_triggering.cpp
epanini Dec 18, 2024
993b936
agn_triggering.hpp
epanini Dec 18, 2024
9ad2113
Update agn_triggering.hpp
epanini Dec 30, 2024
74a8f8f
Update agn_triggering.cpp
epanini Dec 30, 2024
a198d1c
Update agn_triggering.cpp
epanini Jan 7, 2025
c9ef64e
Update agn_triggering.cpp
epanini Jan 7, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
124 changes: 91 additions & 33 deletions src/pgen/cluster/agn_triggering.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,7 @@ AGNTriggering::AGNTriggering(parthenon::ParameterInput *pin,
triggering_filename_(
pin->GetOrAddString(block, "triggering_filename", "agn_triggering.dat")) {


const auto units = hydro_pkg->Param<Units>("units");
const parthenon::Real He_mass_fraction = pin->GetReal("hydro", "He_mass_fraction");
const parthenon::Real H_mass_fraction = 1.0 - He_mass_fraction;
Expand All @@ -85,6 +86,15 @@ AGNTriggering::AGNTriggering(parthenon::ParameterInput *pin,
switch (triggering_mode_) {
case AGNTriggeringMode::COLD_GAS: {
hydro_pkg->AddParam<Real>("agn_triggering_cold_mass", 0, Params::Mutability::Restart);

hydro_pkg->AddParam<Real>("agn_triggering_total_mass", 0,
Params::Mutability::Restart);
hydro_pkg->AddParam<Real>("agn_triggering_prev_cold_mass", 0,
Params::Mutability::Restart);
hydro_pkg->AddParam<Real>("agn_triggering_prev_total_mass", 0,
Params::Mutability::Restart);


break;
}
case AGNTriggeringMode::BOOSTED_BONDI:
Expand Down Expand Up @@ -122,11 +132,14 @@ AGNTriggering::AGNTriggering(parthenon::ParameterInput *pin,
// and simultaneously remove cold gas (updating conserveds and primitives)
template <typename EOS>
void AGNTriggering::ReduceColdMass(parthenon::Real &cold_mass,
parthenon::Real &total_mass,
parthenon::MeshData<parthenon::Real> *md,
const parthenon::Real dt, const EOS eos) const {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the const can stay as the function does not modify the member variables of the AGNTriggering object.

const parthenon::Real dt, const EOS eos) const {
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");

Expand All @@ -153,25 +166,31 @@ 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) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks like the abstraction for a par reduce to two scalars is currently not available within Parthenon (though we're working on this).
In the meantime I suggest to use a raw Kokkos pattern, i.e.,

    Kokkos::parallel_reduce(
        "AGNTriggering::ReduceColdGas",
        Kokkos::MDRangePolicy<Kokkos::Rank<4>>(
            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 =
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 &&
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);

if (temp <= cold_temp_thresh) {

const Real cell_cold_mass = prim(IDN, k, j, i) * coords.CellVolume(k, j, i);
Expand All @@ -193,8 +212,11 @@ void AGNTriggering::ReduceColdMass(parthenon::Real &cold_mass,
}
}
},
Kokkos::Sum<Real>(md_cold_mass));

Kokkos::Sum<Real>(md_cold_mass), Kokkos::Sum<Real>(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
Expand Down Expand Up @@ -234,25 +256,27 @@ void AGNTriggering::ReduceBondiTriggeringQuantities(
KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i,
Real &ltotal_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;
}
}
},
Comment on lines +261 to 282
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What's going on with these changes here?

As far as I can tell the if k>=kb.s is new, but it does't change anything as the indices only go over the interior indices anyway.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this ia a bug because I didn't change this part of the code

total_mass_red, mass_weighted_density_red, mass_weighted_velocity_red,
mass_weighted_cs_red);
Expand All @@ -274,7 +298,7 @@ void AGNTriggering::RemoveBondiAccretedGas(parthenon::MeshData<parthenon::Real>
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<std::string>{"prim"});
Expand Down Expand Up @@ -366,11 +390,15 @@ AGNTriggering::GetAccretionRate(parthenon::StateDescriptor *hydro_pkg) const {

parthenon::TaskStatus
AGNTriggeringResetTriggering(parthenon::StateDescriptor *hydro_pkg) {
const auto &agn_triggering = hydro_pkg->Param<AGNTriggering>("agn_triggering");

auto &agn_triggering = hydro_pkg->Param<AGNTriggering>("agn_triggering");
hydro_pkg->UpdateParam<Real>("agn_triggering_prev_cold_mass", hydro_pkg->Param<Real>("agn_triggering_cold_mass"));
hydro_pkg->UpdateParam<Real>("agn_triggering_prev_total_mass", hydro_pkg->Param<Real>("agn_triggering_total_mass"));
switch (agn_triggering.triggering_mode_) {
case AGNTriggeringMode::COLD_GAS: {
hydro_pkg->UpdateParam<Real>("agn_triggering_cold_mass", 0);
hydro_pkg->UpdateParam<Real>("agn_triggering_total_mass", 0);


break;
}
case AGNTriggeringMode::BOOSTED_BONDI:
Expand All @@ -394,23 +422,27 @@ AGNTriggeringReduceTriggering(parthenon::MeshData<parthenon::Real> *md,

auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
const auto &agn_triggering = hydro_pkg->Param<AGNTriggering>("agn_triggering");

const parthenon::Real pre_accretion_cold_mass = hydro_pkg->Param<Real>("agn_triggering_cold_mass");
const parthenon::Real pre_accretion_total_mass = hydro_pkg->Param<Real>("agn_triggering_total_mass");
switch (agn_triggering.triggering_mode_) {
case AGNTriggeringMode::COLD_GAS: {
Real cold_mass = hydro_pkg->Param<parthenon::Real>("agn_triggering_cold_mass");
Real total_mass = hydro_pkg->Param<parthenon::Real>("agn_triggering_total_mass");


auto fluid = hydro_pkg->Param<Fluid>("fluid");
if (fluid == Fluid::euler) {
agn_triggering.ReduceColdMass(cold_mass, md, dt,
agn_triggering.ReduceColdMass(cold_mass, total_mass,md, dt,
hydro_pkg->Param<AdiabaticHydroEOS>("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<AdiabaticGLMMHDEOS>("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:
Expand All @@ -437,20 +469,25 @@ AGNTriggeringReduceTriggering(parthenon::MeshData<parthenon::Real> *md,
break;
}
}

return TaskStatus::complete;
}

parthenon::TaskStatus
AGNTriggeringMPIReduceTriggering(parthenon::StateDescriptor *hydro_pkg) {
#ifdef MPI_PARALLEL
const auto &agn_triggering = hydro_pkg->Param<AGNTriggering>("agn_triggering");
auto &agn_triggering = hydro_pkg->Param<AGNTriggering>("agn_triggering");
switch (agn_triggering.triggering_mode_) {
case AGNTriggeringMode::COLD_GAS: {

Real accretion_rate = hydro_pkg->Param<Real>("agn_triggering_cold_mass");
PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &accretion_rate, 1,
Real quantities[] ={
hydro_pkg->Param<Real>("agn_triggering_cold_mass"),
hydro_pkg->Param<Real>("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);
hydro_pkg->UpdateParam("agn_triggering_cold_mass", quantities[0]);
hydro_pkg->UpdateParam("agn_triggering_total_mass", quantities[1]);
break;
}
case AGNTriggeringMode::BOOSTED_BONDI:
Expand Down Expand Up @@ -485,15 +522,36 @@ parthenon::TaskStatus
AGNTriggeringFinalizeTriggering(parthenon::MeshData<parthenon::Real> *md,
const parthenon::SimTime &tm) {
auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
const auto &agn_triggering = hydro_pkg->Param<AGNTriggering>("agn_triggering");
auto &agn_triggering = hydro_pkg->Param<AGNTriggering>("agn_triggering");
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<Real>("agn_triggering_cold_mass");
total_mass = hydro_pkg->Param<Real>("agn_triggering_total_mass");


inflow_cold = (cold_mass - hydro_pkg->Param<Real>("agn_triggering_prev_cold_mass")) / tm.dt;
inflow_tot = (total_mass - hydro_pkg->Param<Real>("agn_triggering_prev_total_mass")) / tm.dt;

}
hydro_pkg->UpdateParam<AGNTriggering>("agn_triggering", agn_triggering);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line can probably go as it's not required (agn_triggering is not modified in the code above) and might raise an error because you're updating a non-mutable object.


hydro_pkg->UpdateParam<Real>("agn_triggering_prev_cold_mass", cold_mass);
hydro_pkg->UpdateParam<Real>("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 << " | "
<< inflow_cold << " | "
<< inflow_tot << " ";

switch (agn_triggering.triggering_mode_) {
case AGNTriggeringMode::COLD_GAS: {
Expand Down
7 changes: 4 additions & 3 deletions src/pgen/cluster/agn_triggering.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@
#include <mesh/mesh.hpp>
#include <parameter_input.hpp>
#include <parthenon/package.hpp>

#include <string>
// AthenaPK headers
#include "../../units.hpp"
#include "jet_coords.hpp"
Expand All @@ -34,10 +34,11 @@ class AGNTriggering {
private:
const parthenon::Real gamma_;
parthenon::Real mean_molecular_mass_;

public:
const AGNTriggeringMode triggering_mode_;


const parthenon::Real accretion_radius_;

// Parameters for cold-gas triggering
Expand Down Expand Up @@ -71,7 +72,7 @@ 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 <typename EOS>
void ReduceColdMass(parthenon::Real &cold_mass,
void ReduceColdMass(parthenon::Real &cold_mass, parthenon::Real &total_mass,
parthenon::MeshData<parthenon::Real> *md, const parthenon::Real dt,
const EOS eos) const;

Expand Down