Skip to content

Commit

Permalink
comments applied
Browse files Browse the repository at this point in the history
  • Loading branch information
@mari2895
mari2895 committed Mar 20, 2024
1 parent 4f761b2 commit 2b604be
Show file tree
Hide file tree
Showing 4 changed files with 42 additions and 45 deletions.
3 changes: 3 additions & 0 deletions src/analysis/analysis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,12 +11,15 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
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
Real radius_cutoff_mdot =
pin->GetOrAddReal("analysis", "radius_cutoff_mdot", 0.04); // default 400km
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);
params.Add("radius_cutoff_mdot", radius_cutoff_mdot);

return analysis_pkg;
} // initialize
Expand Down
47 changes: 11 additions & 36 deletions src/analysis/history.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -280,26 +280,6 @@ 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);

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;
}

// 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
Expand All @@ -323,10 +303,12 @@ Real CalculateMdot(MeshData<Real> *md, Real rc, bool gain) {
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));
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;
auto analysis = pmb->packages.Get("analysis").get();
const Real inside_pns_threshold = analysis->Param<Real>("inside_pns_threshold");

parthenon::par_reduce(
parthenon::LoopPatternMDRange(), "Calculates mass accretion rate (SN)",
Expand All @@ -340,20 +322,13 @@ Real CalculateMdot(MeshData<Real> *md, Real rc, bool gain) {
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;
bool is_netheat = (v(b, internal_variables::GcovHeat(), k, j, i) -
v(b, internal_variables::GcovCool(), k, j, i) >
1.e-8); // checks that in the gain region
auto analysis = pmb->packages.Get("analysis").get();
const Real inside_pns_threshold =
analysis->Param<Real>("inside_pns_threshold");
bool is_inside_pns = (r < inside_pns_threshold); // checks that inside PNS
if (do_gain && (is_inside_pns || is_netheat)) {

lresult += ComputeDivInPillbox(
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);
});
Expand All @@ -362,10 +337,10 @@ Real CalculateMdot(MeshData<Real> *md, Real rc, bool gain) {
}

} 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);
});
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 {
Expand Down
24 changes: 22 additions & 2 deletions src/analysis/history.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,8 @@ Real ReduceInGain(MeshData<Real> *md, bool is_conserved, int idx = 0) {
Real result = 0.0;

auto geom = Geometry::GetCoordinateSystem(md);
auto analysis = pmb->packages.Get("analysis").get();
const Real outside_pns_threshold = analysis->Param<Real>("outside_pns_threshold");

parthenon::par_reduce(
parthenon::LoopPatternMDRange(),
Expand All @@ -117,8 +119,6 @@ Real ReduceInGain(MeshData<Real> *md, bool is_conserved, int idx = 0) {
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) >
1.e-8); // checks that in the gain region
auto analysis = pmb->packages.Get("analysis").get();
const Real outside_pns_threshold = analysis->Param<Real>("outside_pns_threshold");
bool is_outside_pns = (v(b, fluid_prim::entropy(), k, j, i) >
outside_pns_threshold); // checks that outside PNS
const auto &coords = v.GetCoordinates(b);
Expand All @@ -134,6 +134,26 @@ Real ReduceInGain(MeshData<Real> *md, bool is_conserved, int idx = 0) {
return result;
}

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);

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;
}

} // namespace History

#endif // ANALYSIS_HISTORY_HPP_
13 changes: 6 additions & 7 deletions src/progenitor/progenitordata.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#include <memory>
#include <vector>

#include "analysis/analysis.hpp"
#include "analysis/history.hpp"
#include "ascii_reader.hpp"
#include "geometry/geometry.hpp"
Expand Down Expand Up @@ -109,23 +110,21 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
params.Add("Srr_adm_dev", Srr_adm_dev);

// Reductions
const Real rc = pin->GetReal("analysis", "radius_cutoff_mdot");
auto HstSum = parthenon::UserHistoryOperation::sum;
auto HstMax = parthenon::UserHistoryOperation::max;
using History::ReduceInGain;
using History::ReduceOneVar;
using parthenon::HistoryOutputVar;
parthenon::HstVar_list hst_vars = {};
auto Mgain = [](MeshData<Real> *md) {
return ReduceInGain<class fluid_cons::density>(md, 1, 0);
return ReduceInGain<fluid_cons::density>(md, 1, 0);
};
auto Qgain = [](MeshData<Real> *md) {
return ReduceInGain<class internal_variables::GcovHeat>(md, 0, 0) -
ReduceInGain<class internal_variables::GcovCool>(md, 0, 0);
};
auto Mdot400 = [](MeshData<Real> *md) {
Real rc = 0.04;
return History::CalculateMdot(md, rc, 0);
return ReduceInGain<internal_variables::GcovHeat>(md, 0, 0) -
ReduceInGain<internal_variables::GcovCool>(md, 0, 0);
};
auto Mdot400 = [rc](MeshData<Real> *md) { return History::CalculateMdot(md, rc, 0); };

Real x1max = pin->GetReal("parthenon/mesh", "x1max");
auto Mdot_gain = [x1max](MeshData<Real> *md) {
Expand Down

0 comments on commit 2b604be

Please sign in to comment.