Skip to content

Commit

Permalink
Output/recording feeding/feedback variables
Browse files Browse the repository at this point in the history 
Rearranged output file in agntriggering and implemented new output file for agnfeedback recording variables
  • Loading branch information
@epanini
epanini authored Jan 25, 2025
1 parent c9ef64e commit 74df575
Show file tree
Hide file tree
Showing 3 changed files with 192 additions and 82 deletions.
68 changes: 66 additions & 2 deletions src/pgen/cluster/agn_feedback.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,8 @@ namespace cluster {
using namespace parthenon;

AGNFeedback::AGNFeedback(parthenon::ParameterInput *pin,
parthenon::StateDescriptor *hydro_pkg)
parthenon::StateDescriptor *hydro_pkg,
const std::string &block)
: fixed_power_(pin->GetOrAddReal("problem/cluster/agn_feedback", "fixed_power", 0.0)),
vceil_(pin->GetOrAddReal("problem/cluster/agn_feedback", "vceil",
std::numeric_limits<Real>::infinity())),
Expand All @@ -55,7 +56,10 @@ AGNFeedback::AGNFeedback(parthenon::ParameterInput *pin,
pin->GetOrAddBoolean("problem/cluster/agn_feedback", "enable_tracer", false)),
disabled_(pin->GetOrAddBoolean("problem/cluster/agn_feedback", "disabled", false)),
enable_magnetic_tower_mass_injection_(pin->GetOrAddBoolean(
"problem/cluster/agn_feedback", "enable_magnetic_tower_mass_injection", true)) {
"problem/cluster/agn_feedback", "enable_magnetic_tower_mass_injection", true)),
write_to_file_(pin->GetOrAddBoolean(block, "write_to_file", false)),
feedback_filename_(
pin->GetOrAddString(block, "feedback_filename", "agn_feedback.dat")) {

// Normalize the thermal, kinetic, and magnetic fractions to sum to 1.0
const Real total_frac = thermal_fraction_ + kinetic_fraction_ + magnetic_fraction_;
Expand Down Expand Up @@ -181,7 +185,18 @@ AGNFeedback::AGNFeedback(parthenon::ParameterInput *pin,
PARTHENON_REQUIRE_THROWS(!enable_tracer_ || hydro_pkg->Param<int>("nscalars") == 1,
"Enabling tracer for AGN feedback requires hydro/nscalars=1");



if (write_to_file_ && parthenon::Globals::my_rank == 0) {

std::ofstream feedback_file;
feedback_file.open(feedback_filename_, std::ofstream::out | std::ofstream::trunc);
feedback_file.close();
}

hydro_pkg->AddParam<>("agn_feedback", *this);


}

parthenon::Real AGNFeedback::GetFeedbackPower(StateDescriptor *hydro_pkg) const {
Expand Down Expand Up @@ -415,6 +430,55 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData<parthenon::Real> *md,
const Real magnetic_power = power * magnetic_fraction_;
const Real magnetic_mass_rate = mass_rate * magnetic_mass_fraction_;
magnetic_tower.PowerSrcTerm(magnetic_power, magnetic_mass_rate, md, beta_dt, tm);
AGNFeedbackFinalizeFeedback(md, tm, power, mass_rate, magnetic_power,
magnetic_mass_rate, thermal_feedback, thermal_density,
jet_density, jet_feedback, jet_momentum);

}


void AGNFeedback::AGNFeedbackFinalizeFeedback(parthenon::MeshData<parthenon::Real> *md,
const parthenon::SimTime &tm, const Real power,
const Real mass_rate,
const Real magnetic_power,
const Real magnetic_mass_rate,
const Real thermal_feedback,
const Real thermal_density,
const Real jet_density,
const Real jet_feedback,
const Real jet_momentum) {

using parthenon::Real;
auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
auto &agn_feedback = hydro_pkg->Param<AGNFeedback>("agn_feedback");


// Append quantities to file
if (agn_feedback.write_to_file_ && parthenon::Globals::my_rank == 0) {
std::ofstream feedback_file;

// Open file in append mode
feedback_file.open(agn_feedback.feedback_filename_, std::ofstream::app);

// Check if the file is empty and write headers
if (feedback_file.tellp() == 0) {
feedback_file << "Time | TimeStep | Power | MassRate | MagneticPower | MagneticMassRate | "
<< "ThermalEnergy | ThermalDensity | KineticDensity | KineticEnergy | JetMomentum";
feedback_file << std::endl; // End of headers
}

// Write data
feedback_file << tm.time << " | " << tm.dt << " | "
<< power << " | "
<< mass_rate << " | "
<< magnetic_power << " | " << magnetic_mass_rate << " | "
<< thermal_feedback << " | " << thermal_density << " | "
<< jet_density << " | " << jet_feedback << " | "
<< jet_momentum;

feedback_file << std::endl;
feedback_file.close();
}
}

} // namespace cluster
23 changes: 22 additions & 1 deletion src/pgen/cluster/agn_feedback.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,17 @@ namespace cluster {
* AGNFeedback
************************************************************/
class AGNFeedback {
private:
parthenon::Real power;
parthenon::Real mass_rate;

parthenon::Real magnetic_power;
parthenon::Real magnetic_mass_rate;
parthenon::Real thermal_feedback;
parthenon::Real thermal_density;
parthenon::Real jet_density;
parthenon::Real jet_feedback;
parthenon::Real jet_momentum;
public:
const parthenon::Real fixed_power_;
parthenon::Real thermal_fraction_, kinetic_fraction_, magnetic_fraction_;
Expand All @@ -48,8 +59,11 @@ class AGNFeedback {
const bool disabled_;

const bool enable_magnetic_tower_mass_injection_;
const bool write_to_file_;
const std::string feedback_filename_;

AGNFeedback(parthenon::ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg);

AGNFeedback(parthenon::ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg, const std::string &block = "problem/cluster/agn_feedback");

parthenon::Real GetFeedbackPower(parthenon::StateDescriptor *hydro_pkg) const;
parthenon::Real GetFeedbackMassRate(parthenon::StateDescriptor *hydro_pkg) const;
Expand All @@ -62,8 +76,15 @@ class AGNFeedback {
void FeedbackSrcTerm(parthenon::MeshData<parthenon::Real> *md,
const parthenon::Real beta_dt, const parthenon::SimTime &tm,
const EOS &eos) const;
friend parthenon::TaskStatus
void AGNFeedbackFinalizeFeedback(parthenon::MeshData<parthenon::Real> *md,
const parthenon::SimTime &tm );
};

parthenon::TaskStatus
void AGNFeedbackFinalizeFeedback(parthenon::MeshData<parthenon::Real> *md,
const parthenon::SimTime &tm );

} // namespace cluster

#endif // CLUSTER_AGN_FEEDBACK_HPP_
183 changes: 104 additions & 79 deletions src/pgen/cluster/agn_triggering.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,7 +168,7 @@ void AGNTriggering::ReduceColdMass(parthenon::Real &cold_mass,



Kokkos::parallel_reduce(
Kokkos::parallel_reduce(
"AGNTriggering::ReduceColdGas",
Kokkos::MDRangePolicy<Kokkos::Rank<4>>(
DevExecSpace(), {0, kb.s, jb.s, ib.s},
Expand All @@ -180,45 +180,47 @@ void AGNTriggering::ReduceColdMass(parthenon::Real &cold_mass,
auto &prim = prim_pack(b);
const auto &coords = cons_pack.GetCoords(b);

// Compute radius squared
const parthenon::Real r2 =
pow(coords.Xc<1>(i), 2) + pow(coords.Xc<2>(j), 2) + pow(coords.Xc<3>(k), 2);

// Consider only cells within the accretion radius
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);

if (k >= int_kb.s && k <= int_kb.e && j >= int_jb.s && j <= int_jb.e &&
// Reduce total mass for interior grid only
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_cold_mass += cell_cold_mass;
team_total_mass += cell_total_mass;
}

const Real cell_delta_rho = -prim(IDN, k, j, i) / cold_t_acc * dt;

if (remove_accreted_mass) {
AddDensityToConsAtFixedVelTemp(cell_delta_rho, cons, prim, eos.GetGamma(),
k, j, i);
// Update the Primitives
eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
// Compute temperature and filter based on threshold
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);

// Reduce cold gas mass for interior grid only
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) {
team_cold_mass += cell_cold_mass;
}

// Remove cold gas if configured to do so
if (remove_accreted_mass) {
const Real cell_delta_rho = -prim(IDN, k, j, i) / cold_t_acc * dt;
AddDensityToConsAtFixedVelTemp(cell_delta_rho, cons, prim, eos.GetGamma(), k, j, i);
eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
}
}
}
}
},
Kokkos::Sum<Real>(md_cold_mass), Kokkos::Sum<Real>(md_total_mass));
},
Kokkos::Sum<Real>(md_cold_mass), Kokkos::Sum<Real>(md_total_mass));

// Update global variables with the computed mass values
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 @@ -521,89 +523,112 @@ AGNTriggeringMPIReduceTriggering(parthenon::StateDescriptor *hydro_pkg) {
}

parthenon::TaskStatus
AGNTriggeringFinalizeTriggering(parthenon::MeshData<parthenon::Real> *md,
const parthenon::SimTime &tm) {
void AGNTriggeringFinalizeTriggering(parthenon::MeshData<parthenon::Real> *md,
const parthenon::SimTime &tm) {
auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
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");

// Declare variables at the function level
parthenon::Real cold_mass = 0.0;
parthenon::Real total_mass = 0.0;
parthenon::Real inflow_cold = 0.0;
parthenon::Real inflow_tot = 0.0;

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;
// Calculate masses and inflow rates for COLD_GAS mode
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;
}


// Update previous mass parameters
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
// Append quantities to file
if (agn_triggering.write_to_file_ && parthenon::Globals::my_rank == 0) {
std::ofstream triggering_file;

// Open file in append mode
triggering_file.open(agn_triggering.triggering_filename_, std::ofstream::app);

// Check if the file is empty and write headers
if (triggering_file.tellp() == 0) {
triggering_file << "Time | TimeStep | AccretionRate";
switch (agn_triggering.triggering_mode_) {
case AGNTriggeringMode::COLD_GAS: {
triggering_file << " | InflowCold | InflowTot | ColdMass | TotalMass";
break;
}
case AGNTriggeringMode::BOOSTED_BONDI:
case AGNTriggeringMode::BOOTH_SCHAYE: {
triggering_file << " | AvgDensity | AvgVelocity | AvgCS";
break;
}
case AGNTriggeringMode::NONE: {
break;
}
}
triggering_file << std::endl; // End of headers
}

// Write data
triggering_file << tm.time << " | " << tm.dt << " | "
<< agn_triggering.GetAccretionRate(hydro_pkg.get()) << " | "
<< cold_mass << " | "
<< total_mass << " | "
<< inflow_cold << " | "
<< inflow_tot << " ";
<< agn_triggering.GetAccretionRate(hydro_pkg.get());

switch (agn_triggering.triggering_mode_) {
case AGNTriggeringMode::COLD_GAS: {

triggering_file << hydro_pkg->Param<Real>("agn_triggering_cold_mass");
break;
}
case AGNTriggeringMode::BOOSTED_BONDI:
case AGNTriggeringMode::BOOTH_SCHAYE: {
const auto &total_mass = hydro_pkg->Param<Real>("agn_triggering_total_mass");
const auto &avg_density =
case AGNTriggeringMode::COLD_GAS: {
triggering_file << " | " << inflow_cold << " | " << inflow_tot << " | "
<< cold_mass << " | " << total_mass;
break;
}
case AGNTriggeringMode::BOOSTED_BONDI:
case AGNTriggeringMode::BOOTH_SCHAYE: {
const auto &total_mass = hydro_pkg->Param<Real>("agn_triggering_total_mass");
const auto &avg_density =
hydro_pkg->Param<Real>("agn_triggering_mass_weighted_density") / total_mass;
const auto &avg_velocity =
const auto &avg_velocity =
hydro_pkg->Param<Real>("agn_triggering_mass_weighted_velocity") / total_mass;
const auto &avg_cs =
const auto &avg_cs =
hydro_pkg->Param<Real>("agn_triggering_mass_weighted_cs") / total_mass;
triggering_file << total_mass << " " << avg_density << " " << avg_velocity << " "
<< avg_cs;
break;
}
case AGNTriggeringMode::NONE: {
break;
}
triggering_file << " | " << avg_density << " | " << avg_velocity << " | " << avg_cs;
break;
}
case AGNTriggeringMode::NONE: {
break;
}
}

triggering_file << std::endl;
triggering_file.close();
}


// Remove accreted gas if using a Bondi-like mode
if (agn_triggering.remove_accreted_mass_) {
switch (agn_triggering.triggering_mode_) {
case AGNTriggeringMode::BOOSTED_BONDI:
case AGNTriggeringMode::BOOTH_SCHAYE: {
auto fluid = hydro_pkg->Param<Fluid>("fluid");
if (fluid == Fluid::euler) {
agn_triggering.RemoveBondiAccretedGas(md, tm.dt,
case AGNTriggeringMode::BOOSTED_BONDI:
case AGNTriggeringMode::BOOTH_SCHAYE: {
auto fluid = hydro_pkg->Param<Fluid>("fluid");
if (fluid == Fluid::euler) {
agn_triggering.RemoveBondiAccretedGas(md, tm.dt,
hydro_pkg->Param<AdiabaticHydroEOS>("eos"));
} else if (fluid == Fluid::glmmhd) {
agn_triggering.RemoveBondiAccretedGas(
} else if (fluid == Fluid::glmmhd) {
agn_triggering.RemoveBondiAccretedGas(
md, tm.dt, hydro_pkg->Param<AdiabaticGLMMHDEOS>("eos"));
} else {
PARTHENON_FAIL("AGNTriggeringFinalizeTriggering: Unknown EOS");
} else {
PARTHENON_FAIL("AGNTriggeringFinalizeTriggering: Unknown EOS");
}
break;
}
case AGNTriggeringMode::COLD_GAS: // Already removed during reduction
case AGNTriggeringMode::NONE: {
break;
}
break;
}
case AGNTriggeringMode::COLD_GAS: // Already removed during reduction
case AGNTriggeringMode::NONE: {
break;
}
}
}

Expand Down

4 comments on commit 74df575

@epanini
Copy link
Author

@epanini epanini commented on 74df575 Jan 25, 2025 •
edited
Loading

Choose a reason for hiding this comment

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

Check only agn_feedback files, I add a new function FinalizeFeedback at the end (following the setup of agn_triggering)
I got this compilation problem @pgrete can you help me? I don't undertand what is wrong..all parameters of the new function are the same everytime I call the function and it is correclty initialized:
/scratch1/epanini/athenapk/src/pgen/cluster/agn_feedback.hpp(80): error: invalid combination of type specifiers
/scratch1/epanini/athenapk/src/pgen/cluster/agn_feedback.hpp(85): error: invalid combination of type specifiers

Moreover, I changed the variables that I want to record from const to parthenon::Real in order to use them in the FinalizeFeedback function, but it seems that they still not accessible:
2 errors detected in the compilation of "/scratch1/epanini/athenapk/src/pgen/cluster.cpp".
make[2]: *** [src/CMakeFiles/athenaPK.dir/build.make:314: src/CMakeFiles/athenaPK.dir/pgen/cluster.cpp.o] Error 2
make[2]: *** Waiting for unfinished jobs....
/scratch1/epanini/athenapk/src/pgen/cluster/agn_feedback.hpp(80): error: invalid combination of type specifiers

/scratch1/epanini/athenapk/src/pgen/cluster/agn_feedback.hpp(85): error: invalid combination of type specifiers

/scratch1/epanini/athenapk/src/pgen/cluster/agn_feedback.cpp(438): error: invalid combination of type specifiers

/scratch1/epanini/athenapk/src/pgen/cluster/agn_feedback.cpp(462): error: identifier "power" is undefined

/scratch1/epanini/athenapk/src/pgen/cluster/agn_feedback.cpp(463): error: identifier "mass_rate" is undefined

/scratch1/epanini/athenapk/src/pgen/cluster/agn_feedback.cpp(464): error: identifier "magnetic_power" is undefined

/scratch1/epanini/athenapk/src/pgen/cluster/agn_feedback.cpp(464): error: identifier "magnetic_mass_rate" is undefined

/scratch1/epanini/athenapk/src/pgen/cluster/agn_feedback.cpp(465): error: identifier "thermal_feedback" is undefined

/scratch1/epanini/athenapk/src/pgen/cluster/agn_feedback.cpp(465): error: identifier "thermal_density" is undefined

/scratch1/epanini/athenapk/src/pgen/cluster/agn_feedback.cpp(466): error: identifier "jet_density" is undefined

/scratch1/epanini/athenapk/src/pgen/cluster/agn_feedback.cpp(466): error: identifier "jet_feedback" is undefined

/scratch1/epanini/athenapk/src/pgen/cluster/agn_feedback.cpp(467): error: identifier "jet_momentum" is undefined

14 errors detected in the compilation of "/scratch1/epanini/athenapk/src/pgen/cluster/agn_feedback.cpp".
make[2]: *** [src/CMakeFiles/athenaPK.dir/build.make:328: src/CMakeFiles/athenaPK.dir/pgen/cluster/agn_feedback.cpp.o] Error 2
make[1]: *** [CMakeFiles/Makefile2:1695: src/CMakeFiles/athenaPK.dir/all] Error 2
make: *** [Makefile:146: all] Error 2

@pgrete
Copy link
Contributor

@pgrete pgrete commented on 74df575 Jan 27, 2025

Choose a reason for hiding this comment

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

I'm not quite sure where the issue is. I'll be able to have a closer look on 5 Feb (when I'm back at work) if the issue still persists.

@epanini
Copy link
Author

@epanini epanini commented on 74df575 Jan 27, 2025 via email

Choose a reason for hiding this comment

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

@epanini
Copy link
Author

@epanini epanini commented on 74df575 Jan 27, 2025 via email

Choose a reason for hiding this comment

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

Please sign in to comment.