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

Jump to bottom

calculating Mass accretion rate for SN post-processing #211

Merged
merged 7 commits into from
Mar 21, 2024
Merged

calculating Mass accretion rate for SN post-processing #211

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
4f761b2
calculating Mdot for SN post-processing
mari2895 Mar 19, 2024
2b604be
comments applied
mari2895 Mar 20, 2024
d4073c5
moved thresholds from analysis pkg to progenitor pkg
mari2895 Mar 20, 2024
994e2cc
Mdot can now take a set or thresholds
mari2895 Mar 20, 2024
7f3bc17
few more details chainged to address suggestions
mari2895 Mar 20, 2024
4658fa0
typo
mari2895 Mar 21, 2024
c6f72c5
print statements deleted
mari2895 Mar 21, 2024
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
7 changes: 7 additions & 0 deletions src/analysis/analysis.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,10 +6,17 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
auto analysis_pkg = std::make_shared<StateDescriptor>("analysis");
Params &params = analysis_pkg->AllParams();
Real sigma = pin->GetOrAddReal("analysis", "sigma", 2e-3);
Real outside_pns_threshold =
pin->GetOrAddReal("analysis", "outside_pns_threshold",
2.42e-5); // corresponds to entropy > 3 kb/baryon
Real inside_pns_threshold = pin->GetOrAddReal("analysis", "inside_pns_threshold",
0.008); // corresponds to r < 80 km
if (sigma < 0) {
PARTHENON_THROW("sigma must be greater than 0");
}
params.Add("sigma", sigma);
params.Add("outside_pns_threshold", outside_pns_threshold);
params.Add("inside_pns_threshold", inside_pns_threshold);

return analysis_pkg;
} // initialize
Expand Down
97 changes: 97 additions & 0 deletions src/analysis/history.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include "history_utils.hpp"
#include "phoebus_utils/relativity_utils.hpp"
#include <interface/sparse_pack.hpp>
#include <parthenon/package.hpp>

namespace History {

Expand Down Expand Up @@ -279,4 +280,100 @@ void ReduceLocalizationFunction(MeshData<Real> *md) {

} // exp

using namespace parthenon::package::prelude;
template <typename F>
KOKKOS_INLINE_FUNCTION Real ComputeDivInPillbox(int ndim, int b, int k, int j, int i,
const parthenon::Coordinates_t &coords,
const F &f) {
Real div_mass_flux_integral;
div_mass_flux_integral =
-(f(b, 1, k, j, i + 1) - f(b, 1, k, j, i)) * coords.FaceArea<X1DIR>(k, j, i);
Copy link
Collaborator

Choose a reason for hiding this comment

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

Remove the minus signs here, and add them to the value in mdot as the divergence doesn't have a minus sign in its definition but mdot does, due to the orientation of surface normals.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

done


if (ndim >= 2) {
div_mass_flux_integral +=
-(f(b, 2, k, j + 1, i) - f(b, 2, k, j, i)) * coords.FaceArea<X2DIR>(k, j, i);
}
if (ndim >= 3) {
div_mass_flux_integral +=
-(f(b, 3, k + 1, j, i) - f(b, 3, k, j, i)) * coords.FaceArea<X3DIR>(k, j, i);
}
return div_mass_flux_integral;
}
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved

// This function calculates mass accretion rate which I defined as
// Mdot=Int(dV*d/dx^i{detg*rho*U^i}) where detg is the determinant of four metric, U is
// four-velocity, and dV=d^3x

Real CalculateMdot(MeshData<Real> *md, Real rc, bool gain) {
const auto ib = md->GetBoundsI(IndexDomain::interior);
const auto jb = md->GetBoundsJ(IndexDomain::interior);
const auto kb = md->GetBoundsK(IndexDomain::interior);
namespace c = fluid_cons;
using parthenon::MakePackDescriptor;
auto *pmb = md->GetParentPointer();
Mesh *pmesh = md->GetMeshPointer();
auto &resolved_pkgs = pmesh->resolved_packages;
const int ndim = pmesh->ndim;
using parthenon::PDOpt;
static auto desc = MakePackDescriptor<c::density, internal_variables::GcovHeat,
internal_variables::GcovCool>(
resolved_pkgs.get(), {}, {PDOpt::WithFluxes});
auto v = desc.GetPack(md);

const int nblocks = v.GetNBlocks();
auto geom = Geometry::GetCoordinateSystem(md);
Real result = 0.0;
const bool is_monopole_cart =
(typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::MonopoleCart));
const bool is_monopole_sph =
(typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::MonopoleSph));
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved

parthenon::par_reduce(
parthenon::LoopPatternMDRange(), "Calculates mass accretion rate (SN)",
DevExecSpace(), 0, nblocks - 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 &lresult) {
const auto &coords = v.GetCoordinates(b);
const Real vol = coords.CellVolume(k, j, i);

Real C[NDFULL];
geom.Coords(CellLocation::Cent, b, k, j, i, C);
Real r = std::sqrt(C[1] * C[1] + C[2] * C[2] + C[3] * C[3]);
if (r <= rc) {
if (gain) {
auto rad = pmb->packages.Get("radiation").get();
const parthenon::AllReduce<bool> *pdo_gain_reducer =
rad->MutableParam<parthenon::AllReduce<bool>>("do_gain_reducer");
const bool do_gain = pdo_gain_reducer->val;
Copy link
Collaborator

Choose a reason for hiding this comment

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

This can't be here---you can't access packages from within a par for on GPU. Move this above.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

done

bool is_netheat = (v(b, internal_variables::GcovHeat(), k, j, i) -
v(b, internal_variables::GcovCool(), k, j, i) >
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved
1.e-8); // checks that in the gain region
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved
auto analysis = pmb->packages.Get("analysis").get();
const Real inside_pns_threshold =
analysis->Param<Real>("inside_pns_threshold");
Copy link
Collaborator

Choose a reason for hiding this comment

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

Can't be here. Move this out of the kernel.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

done

bool is_inside_pns = (r < inside_pns_threshold); // checks that inside PNS
if (do_gain && (is_inside_pns || is_netheat)) {

lresult += ComputeDivInPillbox(
ndim, b, k, j, i, coords, [=](int b, int dir, int k, int j, int i) {
return v.flux(b, dir, c::density(), k, j, i);
});
} else {
lresult += 0;
}

} else {
lresult += ComputeDivInPillbox(ndim, b, k, j, i, coords,
[=](int b, int dir, int k, int j, int i) {
return v.flux(b, dir, c::density(), k, j, i);
});
}

} else {
lresult += 0.0;
}
},
Kokkos::Sum<Real>(result));
return result;
} // mdot

} // namespace History
22 changes: 17 additions & 5 deletions src/analysis/history.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ Real ReduceJetEnergyFlux(MeshData<Real> *md);
Real ReduceJetMomentumFlux(MeshData<Real> *md);
Real ReduceMagneticFluxPhi(MeshData<Real> *md);
void ReduceLocalizationFunction(MeshData<Real> *md);
Real CalculateMdot(MeshData<Real> *md, Real rc, bool gain);

template <typename Reducer_t>
Real ReduceOneVar(MeshData<Real> *md, const std::string &varname, int idx = 0) {
Expand Down Expand Up @@ -85,18 +86,20 @@ Real ReduceOneVar(MeshData<Real> *md, const std::string &varname, int idx = 0) {
}

template <typename Varname>
Real ReduceInGain(MeshData<Real> *md, int idx = 0) {
Real ReduceInGain(MeshData<Real> *md, bool is_conserved, int idx = 0) {
const auto ib = md->GetBoundsI(IndexDomain::interior);
const auto jb = md->GetBoundsJ(IndexDomain::interior);
const auto kb = md->GetBoundsK(IndexDomain::interior);

namespace iv = internal_variables;
using parthenon::MakePackDescriptor;
auto *pmb = md->GetParentPointer();
Mesh *pmesh = md->GetMeshPointer();
auto &resolved_pkgs = pmesh->resolved_packages;
const int ndim = pmesh->ndim;
static auto desc =
MakePackDescriptor<Varname, iv::GcovHeat, iv::GcovCool>(resolved_pkgs.get());
MakePackDescriptor<Varname, fluid_prim::entropy, iv::GcovHeat, iv::GcovCool>(
resolved_pkgs.get());
auto v = desc.GetPack(md);
PARTHENON_REQUIRE_THROWS(v.ContainsHost(0, iv::GcovHeat(), iv::GcovCool()),
"Must be doing SN simulation");
Expand All @@ -112,11 +115,20 @@ Real ReduceInGain(MeshData<Real> *md, int idx = 0) {
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 &lresult) {
Real gdet = geom.DetGamma(CellLocation::Cent, 0, k, j, i);
bool is_netheat =
(v(b, iv::GcovHeat(), k, j, i) - v(b, iv::GcovCool(), k, j, i) > 0);
bool is_netheat = (v(b, iv::GcovHeat(), k, j, i) - v(b, iv::GcovCool(), k, j, i) >
1.e-8); // checks that in the gain region
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved
auto analysis = pmb->packages.Get("analysis").get();
const Real outside_pns_threshold = analysis->Param<Real>("outside_pns_threshold");
Copy link
Collaborator

Choose a reason for hiding this comment

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

This needs to be moved outside of the kernel. It can't be called inside the par for.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

done

bool is_outside_pns = (v(b, fluid_prim::entropy(), k, j, i) >
outside_pns_threshold); // checks that outside PNS
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved
const auto &coords = v.GetCoordinates(b);
const Real vol = coords.CellVolume(k, j, i);
lresult += gdet * vol * is_netheat * v(b, Varname(idx), k, j, i);
if (is_conserved) {
lresult += vol * is_netheat * is_outside_pns * v(b, Varname(idx), k, j, i);
} else {
lresult +=
gdet * vol * is_netheat * is_outside_pns * v(b, Varname(idx), k, j, i);
}
},
Kokkos::Sum<Real>(result));
return result;
Expand Down
17 changes: 14 additions & 3 deletions src/progenitor/progenitordata.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,14 +116,25 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
using parthenon::HistoryOutputVar;
parthenon::HstVar_list hst_vars = {};
auto Mgain = [](MeshData<Real> *md) {
return ReduceInGain<class fluid_prim::density>(md, 0);
return ReduceInGain<class fluid_cons::density>(md, 1, 0);
};
auto Qgain = [](MeshData<Real> *md) {
return ReduceInGain<class internal_variables::GcovHeat>(md, 0) -
ReduceInGain<class internal_variables::GcovCool>(md, 0);
return ReduceInGain<class internal_variables::GcovHeat>(md, 0, 0) -
ReduceInGain<class internal_variables::GcovCool>(md, 0, 0);
Copy link
Collaborator

Choose a reason for hiding this comment

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

Is the word class here actually required?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

deleted and still compiles and runs so is not required.

};
auto Mdot400 = [](MeshData<Real> *md) {
Real rc = 0.04;
Copy link
Collaborator

Choose a reason for hiding this comment

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

Maybe make this a runtime set-able parameter. Might even be good to support an arbitrary number of these at different radii.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I am trying to imagine how is the second part going to work

Copy link
Collaborator

Choose a reason for hiding this comment

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

ParameterInput can read a std:::vector so you could do something like this:

auto mdot_radii = pin->GetOrAddVector("analysis", "mdot_radii", std::vector<Real>{0.04});
for (auto rc : mdot_radii) {
  auto mdot = [=](MeshData<Real> *md) { return History::CalculateMdot(md, rc, 0); };
  hst_vars.emplace_back(HistoryOutputVar(HstSum, mdot, "Mdot at r = " + std::to_string(rc)));
}

return History::CalculateMdot(md, rc, 0);
};

Real x1max = pin->GetReal("parthenon/mesh", "x1max");
auto Mdot_gain = [x1max](MeshData<Real> *md) {
return History::CalculateMdot(md, x1max, 1);
};
hst_vars.emplace_back(HistoryOutputVar(HstSum, Mgain, "Mgain"));
hst_vars.emplace_back(HistoryOutputVar(HstSum, Qgain, "total net heat"));
hst_vars.emplace_back(HistoryOutputVar(HstSum, Mdot400, "Mdot at 400km"));
hst_vars.emplace_back(HistoryOutputVar(HstSum, Mdot_gain, "Mdot gain"));
params.Add(parthenon::hist_param_key, hst_vars);

return progenitor_pkg;
Expand Down
Loading