From 3a72c2601fd97fa7ee649ffc971453c43e367e1d Mon Sep 17 00:00:00 2001 From: mari2895 Date: Sun, 10 Mar 2024 10:52:38 -0400 Subject: [PATCH 1/7] gain reducer added, rhoc added, Mgain added --- src/analysis/history.cpp | 82 +++++++++++++++++++++++++++++++++ src/analysis/history.hpp | 54 +++++++++++++++++++++- src/fluid/fluid.cpp | 22 ++++++++- src/phoebus_utils/variables.hpp | 2 + 4 files changed, 157 insertions(+), 3 deletions(-) diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp index 5b044939..dacc9a94 100644 --- a/src/analysis/history.cpp +++ b/src/analysis/history.cpp @@ -16,6 +16,7 @@ #include "geometry/geometry_utils.hpp" #include "history_utils.hpp" #include "phoebus_utils/relativity_utils.hpp" +#include "analysis/analysis.hpp" namespace History { @@ -235,4 +236,85 @@ Real ReduceMagneticFluxPhi(MeshData *md) { return 0.5 * result; // 0.5 \int detg B^r dx2 dx3 } // Phi + + //SN analysis + + void ReduceCentralDensitySN(MeshData *md) { + const auto ib = md->GetBoundsI(IndexDomain::interior); + const auto jb = md->GetBoundsJ(IndexDomain::interior); + const auto kb = md->GetBoundsK(IndexDomain::interior); + namespace p = fluid_prim; + namespace diag = diagnostic_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(resolved_pkgs.get()); + auto v = desc.GetPack(md); + const int nblocks = v.GetNBlocks(); + auto geom = Geometry::GetCoordinateSystem(md); + + + parthenon::par_for( + parthenon::LoopPatternMDRange(), "Central Density for 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) { + auto &coords = v.GetCoordinates(b); + auto analysis = pmb->packages.Get("analysis").get(); + const Real x[3] = {coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i), coords.Xc<3>(k, j, i)}; + const Real sigma = analysis->Param("sigma"); + Real gam[3][3]; + Real r2 = 0; + geom.Metric(CellLocation::Cent, 0, k, j, i, gam); + for (int n = 0; n < 3; ++n) { + for (int m = 0; m < 3; ++m) { + r2 += gam[n][m] * x[n] * x[m]; + } + } + const Real rhoc = v(b, p::density(), k, j, i) *std::exp(-r2 / sigma / sigma); + v(b, diag::central_density(), k, j, i) = rhoc; + }); + + } + void ReduceLocalizationFunction(MeshData *md) { + const auto ib = md->GetBoundsI(IndexDomain::interior); + const auto jb = md->GetBoundsJ(IndexDomain::interior); + const auto kb = md->GetBoundsK(IndexDomain::interior); + namespace diag = diagnostic_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(resolved_pkgs.get()); + auto v = desc.GetPack(md); + const int nblocks = v.GetNBlocks(); + auto geom = Geometry::GetCoordinateSystem(md); + + + parthenon::par_for( + parthenon::LoopPatternMDRange(), "Central Density for 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) { + auto &coords = v.GetCoordinates(b); + auto analysis = pmb->packages.Get("analysis").get(); + const Real x[3] = {coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i), coords.Xc<3>(k, j, i)}; + const Real sigma = analysis->Param("sigma"); + Real gam[3][3]; + Real r2 = 0; + geom.Metric(CellLocation::Cent, 0, k, j, i, gam); + for (int n = 0; n < 3; ++n) { + for (int m = 0; m < 3; ++m) { + r2 += gam[n][m] * x[n] * x[m]; + } + } + v(b, diag::localization_function(), k, j, i) = std::exp(-r2 / sigma / sigma); + }); + + } // exp (this function returns normalization function that is used for localizing quantities at the center, or at some particular case. For example SN diagnostics oftec computes quantities at 400 km.) + + } // namespace History diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp index 04a61d60..6d39e903 100644 --- a/src/analysis/history.hpp +++ b/src/analysis/history.hpp @@ -19,8 +19,10 @@ #include #include +#include #include - +#include "geometry/geometry.hpp" +#include "geometry/geometry_utils.hpp" #include "phoebus_utils/variables.hpp" using namespace parthenon::package::prelude; @@ -41,7 +43,10 @@ Real ReduceMassAccretionRate(MeshData *md); Real ReduceJetEnergyFlux(MeshData *md); Real ReduceJetMomentumFlux(MeshData *md); Real ReduceMagneticFluxPhi(MeshData *md); +void ReduceCentralDensitySN(MeshData *md); +void ReduceLocalizationFunction(MeshData *md); + template Real ReduceOneVar(MeshData *md, const std::string &varname, int idx = 0) { const auto ib = md->GetBoundsI(IndexDomain::interior); @@ -58,11 +63,54 @@ Real ReduceOneVar(MeshData *md, const std::string &varname, int idx = 0) { // We choose to apply volume weighting when using the sum reduction. // This assumes that for "Sum" variables, the variable is densitized, meaning // it already contains a factor of the measure: sqrt(det(gamma)) + Real result = 0.0; + Reducer_t reducer(result); + const bool volume_weighting = std::is_same>::value; + + parthenon::par_reduce( + parthenon::LoopPatternMDRange(), "Phoebus History for " + varname, DevExecSpace(), + 0, 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 &lresult) { + // join is a Kokkos construct + // that automatically does the + // reduction operation locally + const auto &coords = pack.GetCoords(b); + const Real vol = volume_weighting ? coords.CellVolume(k, j, i) : 1.0; + reducer.join(lresult, pack(b, ivar.first + idx, k, j, i) * vol); + }, + reducer); + return result; +} +template +Real ReduceInGain(MeshData *md, const std::string &varname, int idx = 0) { + const auto ib = md->GetBoundsI(IndexDomain::interior); + const auto jb = md->GetBoundsJ(IndexDomain::interior); + const auto kb = md->GetBoundsK(IndexDomain::interior); + + PackIndexMap imap; + std::vector vars = {varname}; + const auto pack = md->PackVariables(vars, imap); + const auto ivar = imap[varname]; + auto *pmb = md->GetParentPointer(); + auto rad = pmb->packages.Get("radiation").get(); + PARTHENON_REQUIRE(ivar.first >= 0, "Var must exist"); + PARTHENON_REQUIRE(ivar.second >= ivar.first + idx, "Var must exist"); + + // We choose to apply volume weighting when using the sum reduction. + // This assumes that for "Sum" variables, the variable is densitized, meaning + // it already contains a factor of the measure: sqrt(det(gamma)) Real result = 0.0; Reducer_t reducer(result); + + const bool volume_weighting = + std::is_same>::value; + parthenon::AllReduce *pdo_gain_reducer = rad->MutableParam>("do_gain_reducer"); + const bool do_gain = pdo_gain_reducer->val; + std::cout<<"do_gain="< *md, const std::string &varname, int idx = 0) { // reduction operation locally const auto &coords = pack.GetCoords(b); const Real vol = volume_weighting ? coords.CellVolume(k, j, i) : 1.0; - reducer.join(lresult, pack(b, ivar.first + idx, k, j, i) * vol); + reducer.join(lresult, pack(b, ivar.first + idx, k, j, i) * vol * do_gain); }, reducer); return result; } + + } // namespace History #endif // ANALYSIS_HISTORY_HPP_ diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp index 66a14580..36e84703 100644 --- a/src/fluid/fluid.cpp +++ b/src/fluid/fluid.cpp @@ -207,6 +207,8 @@ std::shared_ptr Initialize(ParameterInput *pin) { physics->AddField(p::entropy::name(), mprim_scalar); physics->AddField(p::cs::name(), mprim_scalar); physics->AddField(diag::ratio_divv_cs::name(), mprim_scalar); + physics->AddField(diag::central_density::name(), mprim_scalar); + physics->AddField(diag::localization_function::name(), mprim_scalar); physics->AddField(diag::entropy_z_0::name(), mprim_scalar); physics->AddField(p::gamma1::name(), mprim_scalar); if (ye) { @@ -310,15 +312,33 @@ std::shared_ptr Initialize(ParameterInput *pin) { // By default compute integrated value of scalar conserved vars auto HstSum = parthenon::UserHistoryOperation::sum; using History::ReduceOneVar; + using History::ReduceInGain; using parthenon::HistoryOutputVar; parthenon::HstVar_list hst_vars = {}; - + auto ReduceMass = [](MeshData *md) { return ReduceOneVar>(md, fluid_cons::density::name(), 0); }; auto ReduceEn = [](MeshData *md) { return ReduceOneVar>(md, fluid_cons::energy::name(), 0); }; + auto CentralDensitySN = [](MeshData *md) { + History::ReduceCentralDensitySN(md); + return ReduceOneVar>(md, diag::central_density::name(), 0); + + }; + auto norm = [](MeshData *md) { + History::ReduceLocalizationFunction(md); + return ReduceOneVar>(md, diag::localization_function::name(), 0); + }; + auto Mgain = [](MeshData *md) { + return ReduceInGain>(md, fluid_prim::density::name(), 0); + }; + + + hst_vars.emplace_back(HistoryOutputVar(HstSum, CentralDensitySN, "central density SN")); + hst_vars.emplace_back(HistoryOutputVar(HstSum, norm, "Normalization")); + hst_vars.emplace_back(HistoryOutputVar(HstSum, Mgain, "Mgain")); hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceMass, "total baryon number")); hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceEn, "total conserved energy tau")); diff --git a/src/phoebus_utils/variables.hpp b/src/phoebus_utils/variables.hpp index 5f675402..acf66920 100644 --- a/src/phoebus_utils/variables.hpp +++ b/src/phoebus_utils/variables.hpp @@ -122,6 +122,8 @@ namespace diagnostic_variables { VARIABLE_NONS(divb); VARIABLE_NONS(ratio_divv_cs); VARIABLE_NONS(entropy_z_0); +VARIABLE_NONS(central_density); + VARIABLE_NONS(localization_function); VARIABLE_CUSTOM(divf, flux_divergence); VARIABLE_NONS(src_terms); VARIABLE_CUSTOM(r_divf, r.flux_divergence); From 691f6c122959e05d231e8450a4db93be249e2159 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Mon, 11 Mar 2024 12:37:57 -0400 Subject: [PATCH 2/7] computes max density instead of integral --- src/analysis/history.cpp | 140 ++++++++++++++++---------------- src/analysis/history.hpp | 35 ++++---- src/fluid/fluid.cpp | 16 ++-- src/phoebus_utils/variables.hpp | 2 +- 4 files changed, 94 insertions(+), 99 deletions(-) diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp index dacc9a94..fd32fe9b 100644 --- a/src/analysis/history.cpp +++ b/src/analysis/history.cpp @@ -12,11 +12,11 @@ // publicly, and to permit others to do so. #include "history.hpp" +#include "analysis/analysis.hpp" #include "geometry/geometry.hpp" #include "geometry/geometry_utils.hpp" #include "history_utils.hpp" #include "phoebus_utils/relativity_utils.hpp" -#include "analysis/analysis.hpp" namespace History { @@ -236,85 +236,85 @@ Real ReduceMagneticFluxPhi(MeshData *md) { return 0.5 * result; // 0.5 \int detg B^r dx2 dx3 } // Phi +// SN analysis + +void ReduceCentralDensitySN(MeshData *md) { + const auto ib = md->GetBoundsI(IndexDomain::interior); + const auto jb = md->GetBoundsJ(IndexDomain::interior); + const auto kb = md->GetBoundsK(IndexDomain::interior); + namespace p = fluid_prim; + namespace diag = diagnostic_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(resolved_pkgs.get()); + auto v = desc.GetPack(md); + const int nblocks = v.GetNBlocks(); + auto geom = Geometry::GetCoordinateSystem(md); - //SN analysis - - void ReduceCentralDensitySN(MeshData *md) { - const auto ib = md->GetBoundsI(IndexDomain::interior); - const auto jb = md->GetBoundsJ(IndexDomain::interior); - const auto kb = md->GetBoundsK(IndexDomain::interior); - namespace p = fluid_prim; - namespace diag = diagnostic_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(resolved_pkgs.get()); - auto v = desc.GetPack(md); - const int nblocks = v.GetNBlocks(); - auto geom = Geometry::GetCoordinateSystem(md); - - - parthenon::par_for( - parthenon::LoopPatternMDRange(), "Central Density for SN", - DevExecSpace(), 0,nblocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + parthenon::par_for( + parthenon::LoopPatternMDRange(), "Central Density for 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) { - auto &coords = v.GetCoordinates(b); - auto analysis = pmb->packages.Get("analysis").get(); - const Real x[3] = {coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i), coords.Xc<3>(k, j, i)}; + auto &coords = v.GetCoordinates(b); + auto analysis = pmb->packages.Get("analysis").get(); + const Real x[3] = {coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i), + coords.Xc<3>(k, j, i)}; const Real sigma = analysis->Param("sigma"); - Real gam[3][3]; - Real r2 = 0; - geom.Metric(CellLocation::Cent, 0, k, j, i, gam); - for (int n = 0; n < 3; ++n) { - for (int m = 0; m < 3; ++m) { + Real gam[3][3]; + Real r2 = 0; + geom.Metric(CellLocation::Cent, 0, k, j, i, gam); + for (int n = 0; n < 3; ++n) { + for (int m = 0; m < 3; ++m) { r2 += gam[n][m] * x[n] * x[m]; - } - } - const Real rhoc = v(b, p::density(), k, j, i) *std::exp(-r2 / sigma / sigma); - v(b, diag::central_density(), k, j, i) = rhoc; + } + } + const Real rhoc = v(b, p::density(), k, j, i) * std::exp(-r2 / sigma / sigma); + v(b, diag::central_density(), k, j, i) = rhoc; }); +} +void ReduceLocalizationFunction(MeshData *md) { + const auto ib = md->GetBoundsI(IndexDomain::interior); + const auto jb = md->GetBoundsJ(IndexDomain::interior); + const auto kb = md->GetBoundsK(IndexDomain::interior); + namespace diag = diagnostic_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(resolved_pkgs.get()); + auto v = desc.GetPack(md); + const int nblocks = v.GetNBlocks(); + auto geom = Geometry::GetCoordinateSystem(md); - } - void ReduceLocalizationFunction(MeshData *md) { - const auto ib = md->GetBoundsI(IndexDomain::interior); - const auto jb = md->GetBoundsJ(IndexDomain::interior); - const auto kb = md->GetBoundsK(IndexDomain::interior); - namespace diag = diagnostic_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(resolved_pkgs.get()); - auto v = desc.GetPack(md); - const int nblocks = v.GetNBlocks(); - auto geom = Geometry::GetCoordinateSystem(md); - - - parthenon::par_for( - parthenon::LoopPatternMDRange(), "Central Density for SN", - DevExecSpace(), 0,nblocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + parthenon::par_for( + parthenon::LoopPatternMDRange(), "Central Density for 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) { - auto &coords = v.GetCoordinates(b); - auto analysis = pmb->packages.Get("analysis").get(); - const Real x[3] = {coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i), coords.Xc<3>(k, j, i)}; + auto &coords = v.GetCoordinates(b); + auto analysis = pmb->packages.Get("analysis").get(); + const Real x[3] = {coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i), + coords.Xc<3>(k, j, i)}; const Real sigma = analysis->Param("sigma"); - Real gam[3][3]; - Real r2 = 0; - geom.Metric(CellLocation::Cent, 0, k, j, i, gam); - for (int n = 0; n < 3; ++n) { - for (int m = 0; m < 3; ++m) { + Real gam[3][3]; + Real r2 = 0; + geom.Metric(CellLocation::Cent, 0, k, j, i, gam); + for (int n = 0; n < 3; ++n) { + for (int m = 0; m < 3; ++m) { r2 += gam[n][m] * x[n] * x[m]; - } - } - v(b, diag::localization_function(), k, j, i) = std::exp(-r2 / sigma / sigma); + } + } + v(b, diag::localization_function(), k, j, i) = std::exp(-r2 / sigma / sigma); }); - } // exp (this function returns normalization function that is used for localizing quantities at the center, or at some particular case. For example SN diagnostics oftec computes quantities at 400 km.) +} // exp (this function returns normalization function that is used for localizing + // quantities at the center, or at some particular case. For example SN diagnostics + // oftec computes quantities at 400 km.) - } // namespace History diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp index 6d39e903..645dcda1 100644 --- a/src/analysis/history.hpp +++ b/src/analysis/history.hpp @@ -17,13 +17,13 @@ #include #include -#include -#include -#include -#include #include "geometry/geometry.hpp" #include "geometry/geometry_utils.hpp" #include "phoebus_utils/variables.hpp" +#include +#include +#include +#include using namespace parthenon::package::prelude; @@ -46,7 +46,6 @@ Real ReduceMagneticFluxPhi(MeshData *md); void ReduceCentralDensitySN(MeshData *md); void ReduceLocalizationFunction(MeshData *md); - template Real ReduceOneVar(MeshData *md, const std::string &varname, int idx = 0) { const auto ib = md->GetBoundsI(IndexDomain::interior); @@ -67,8 +66,10 @@ Real ReduceOneVar(MeshData *md, const std::string &varname, int idx = 0) { Reducer_t reducer(result); const bool volume_weighting = - std::is_same>::value; - + std::is_same>::value || + std::is_same>::value || + std::is_same>::value; + parthenon::par_reduce( parthenon::LoopPatternMDRange(), "Phoebus History for " + varname, DevExecSpace(), 0, pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, @@ -91,11 +92,13 @@ Real ReduceInGain(MeshData *md, const std::string &varname, int idx = 0) { const auto kb = md->GetBoundsK(IndexDomain::interior); PackIndexMap imap; - std::vector vars = {varname}; + namespace iv = internal_variables; + std::vector vars = {varname, iv::GcovHeat::name(), iv::GcovCool::name()}; const auto pack = md->PackVariables(vars, imap); const auto ivar = imap[varname]; + const auto iheat = imap[iv::GcovHeat::name()].first; + const auto icool = imap[iv::GcovCool::name()].first; auto *pmb = md->GetParentPointer(); - auto rad = pmb->packages.Get("radiation").get(); PARTHENON_REQUIRE(ivar.first >= 0, "Var must exist"); PARTHENON_REQUIRE(ivar.second >= ivar.first + idx, "Var must exist"); @@ -106,11 +109,10 @@ Real ReduceInGain(MeshData *md, const std::string &varname, int idx = 0) { Reducer_t reducer(result); const bool volume_weighting = - std::is_same>::value; - parthenon::AllReduce *pdo_gain_reducer = rad->MutableParam>("do_gain_reducer"); - const bool do_gain = pdo_gain_reducer->val; - std::cout<<"do_gain="<>::value || + std::is_same>::value || + std::is_same>::value; + parthenon::par_reduce( parthenon::LoopPatternMDRange(), "Phoebus History for " + varname, DevExecSpace(), 0, pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, @@ -118,16 +120,15 @@ Real ReduceInGain(MeshData *md, const std::string &varname, int idx = 0) { // join is a Kokkos construct // that automatically does the // reduction operation locally + bool is_netheat = (pack(b, iheat, k, j, i) - pack(b, icool, k, j, i) > 0); const auto &coords = pack.GetCoords(b); const Real vol = volume_weighting ? coords.CellVolume(k, j, i) : 1.0; - reducer.join(lresult, pack(b, ivar.first + idx, k, j, i) * vol * do_gain); + reducer.join(lresult, pack(b, ivar.first + idx, k, j, i) * vol * is_netheat); }, reducer); return result; } - - } // namespace History #endif // ANALYSIS_HISTORY_HPP_ diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp index 36e84703..83462166 100644 --- a/src/fluid/fluid.cpp +++ b/src/fluid/fluid.cpp @@ -311,11 +311,11 @@ std::shared_ptr Initialize(ParameterInput *pin) { // Reductions // By default compute integrated value of scalar conserved vars auto HstSum = parthenon::UserHistoryOperation::sum; - using History::ReduceOneVar; using History::ReduceInGain; + using History::ReduceOneVar; using parthenon::HistoryOutputVar; parthenon::HstVar_list hst_vars = {}; - + auto ReduceMass = [](MeshData *md) { return ReduceOneVar>(md, fluid_cons::density::name(), 0); }; @@ -324,20 +324,14 @@ std::shared_ptr Initialize(ParameterInput *pin) { }; auto CentralDensitySN = [](MeshData *md) { History::ReduceCentralDensitySN(md); - return ReduceOneVar>(md, diag::central_density::name(), 0); - - }; - auto norm = [](MeshData *md) { - History::ReduceLocalizationFunction(md); - return ReduceOneVar>(md, diag::localization_function::name(), 0); + return ReduceOneVar>(md, diag::central_density::name(), 0); }; auto Mgain = [](MeshData *md) { - return ReduceInGain>(md, fluid_prim::density::name(), 0); + return ReduceInGain>(md, fluid_prim::density::name(), + 0); }; - hst_vars.emplace_back(HistoryOutputVar(HstSum, CentralDensitySN, "central density SN")); - hst_vars.emplace_back(HistoryOutputVar(HstSum, norm, "Normalization")); hst_vars.emplace_back(HistoryOutputVar(HstSum, Mgain, "Mgain")); hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceMass, "total baryon number")); hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceEn, "total conserved energy tau")); diff --git a/src/phoebus_utils/variables.hpp b/src/phoebus_utils/variables.hpp index acf66920..14768a83 100644 --- a/src/phoebus_utils/variables.hpp +++ b/src/phoebus_utils/variables.hpp @@ -123,7 +123,7 @@ VARIABLE_NONS(divb); VARIABLE_NONS(ratio_divv_cs); VARIABLE_NONS(entropy_z_0); VARIABLE_NONS(central_density); - VARIABLE_NONS(localization_function); +VARIABLE_NONS(localization_function); VARIABLE_CUSTOM(divf, flux_divergence); VARIABLE_NONS(src_terms); VARIABLE_CUSTOM(r_divf, r.flux_divergence); From 039fdf5e227692d18447e12a8ee1d7d09d30a8a7 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Mon, 11 Mar 2024 18:06:01 -0400 Subject: [PATCH 3/7] added total net heating in history output (SN) --- src/analysis/history.hpp | 5 ++++- src/fluid/fluid.cpp | 12 ++++++++++-- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp index 645dcda1..62c2a610 100644 --- a/src/analysis/history.hpp +++ b/src/analysis/history.hpp @@ -112,6 +112,7 @@ Real ReduceInGain(MeshData *md, const std::string &varname, int idx = 0) { std::is_same>::value || std::is_same>::value || std::is_same>::value; + auto geom = Geometry::GetCoordinateSystem(md); parthenon::par_reduce( parthenon::LoopPatternMDRange(), "Phoebus History for " + varname, DevExecSpace(), @@ -120,10 +121,12 @@ Real ReduceInGain(MeshData *md, const std::string &varname, int idx = 0) { // join is a Kokkos construct // that automatically does the // reduction operation locally + Real gdet = geom.DetGamma(CellLocation::Cent, 0, k, j, i); bool is_netheat = (pack(b, iheat, k, j, i) - pack(b, icool, k, j, i) > 0); const auto &coords = pack.GetCoords(b); const Real vol = volume_weighting ? coords.CellVolume(k, j, i) : 1.0; - reducer.join(lresult, pack(b, ivar.first + idx, k, j, i) * vol * is_netheat); + reducer.join(lresult, + pack(b, ivar.first + idx, k, j, i) * gdet * vol * is_netheat); }, reducer); return result; diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp index 83462166..6fd709aa 100644 --- a/src/fluid/fluid.cpp +++ b/src/fluid/fluid.cpp @@ -311,6 +311,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { // Reductions // By default compute integrated value of scalar conserved vars auto HstSum = parthenon::UserHistoryOperation::sum; + auto HstMax = parthenon::UserHistoryOperation::max; using History::ReduceInGain; using History::ReduceOneVar; using parthenon::HistoryOutputVar; @@ -322,7 +323,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { auto ReduceEn = [](MeshData *md) { return ReduceOneVar>(md, fluid_cons::energy::name(), 0); }; - auto CentralDensitySN = [](MeshData *md) { + auto MaxDensitySN = [](MeshData *md) { History::ReduceCentralDensitySN(md); return ReduceOneVar>(md, diag::central_density::name(), 0); }; @@ -330,9 +331,16 @@ std::shared_ptr Initialize(ParameterInput *pin) { return ReduceInGain>(md, fluid_prim::density::name(), 0); }; + auto Qgain = [](MeshData *md) { + return ReduceInGain>( + md, internal_variables::GcovHeat::name(), 0) - + ReduceInGain>( + md, internal_variables::GcovCool::name(), 0); + }; - hst_vars.emplace_back(HistoryOutputVar(HstSum, CentralDensitySN, "central density SN")); + hst_vars.emplace_back(HistoryOutputVar(HstMax, MaxDensitySN, "maximum density SN")); hst_vars.emplace_back(HistoryOutputVar(HstSum, Mgain, "Mgain")); + hst_vars.emplace_back(HistoryOutputVar(HstSum, Qgain, "total net heat")); hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceMass, "total baryon number")); hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceEn, "total conserved energy tau")); From 7f64429354acabd8bd8d6725e1a0b0bb60395d9d Mon Sep 17 00:00:00 2001 From: mari2895 Date: Mon, 11 Mar 2024 18:48:44 -0400 Subject: [PATCH 4/7] cleanup --- src/analysis/history.cpp | 39 --------------------------------- src/analysis/history.hpp | 1 - src/fluid/fluid.cpp | 4 +--- src/phoebus_utils/variables.hpp | 1 - 4 files changed, 1 insertion(+), 44 deletions(-) diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp index fd32fe9b..41694a26 100644 --- a/src/analysis/history.cpp +++ b/src/analysis/history.cpp @@ -238,45 +238,6 @@ Real ReduceMagneticFluxPhi(MeshData *md) { // SN analysis -void ReduceCentralDensitySN(MeshData *md) { - const auto ib = md->GetBoundsI(IndexDomain::interior); - const auto jb = md->GetBoundsJ(IndexDomain::interior); - const auto kb = md->GetBoundsK(IndexDomain::interior); - namespace p = fluid_prim; - namespace diag = diagnostic_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(resolved_pkgs.get()); - auto v = desc.GetPack(md); - const int nblocks = v.GetNBlocks(); - auto geom = Geometry::GetCoordinateSystem(md); - - parthenon::par_for( - parthenon::LoopPatternMDRange(), "Central Density for 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) { - auto &coords = v.GetCoordinates(b); - auto analysis = pmb->packages.Get("analysis").get(); - const Real x[3] = {coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i), - coords.Xc<3>(k, j, i)}; - const Real sigma = analysis->Param("sigma"); - Real gam[3][3]; - Real r2 = 0; - geom.Metric(CellLocation::Cent, 0, k, j, i, gam); - for (int n = 0; n < 3; ++n) { - for (int m = 0; m < 3; ++m) { - r2 += gam[n][m] * x[n] * x[m]; - } - } - const Real rhoc = v(b, p::density(), k, j, i) * std::exp(-r2 / sigma / sigma); - v(b, diag::central_density(), k, j, i) = rhoc; - }); -} void ReduceLocalizationFunction(MeshData *md) { const auto ib = md->GetBoundsI(IndexDomain::interior); const auto jb = md->GetBoundsJ(IndexDomain::interior); diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp index 62c2a610..2cd97a93 100644 --- a/src/analysis/history.hpp +++ b/src/analysis/history.hpp @@ -43,7 +43,6 @@ Real ReduceMassAccretionRate(MeshData *md); Real ReduceJetEnergyFlux(MeshData *md); Real ReduceJetMomentumFlux(MeshData *md); Real ReduceMagneticFluxPhi(MeshData *md); -void ReduceCentralDensitySN(MeshData *md); void ReduceLocalizationFunction(MeshData *md); template diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp index 6fd709aa..fe096581 100644 --- a/src/fluid/fluid.cpp +++ b/src/fluid/fluid.cpp @@ -207,7 +207,6 @@ std::shared_ptr Initialize(ParameterInput *pin) { physics->AddField(p::entropy::name(), mprim_scalar); physics->AddField(p::cs::name(), mprim_scalar); physics->AddField(diag::ratio_divv_cs::name(), mprim_scalar); - physics->AddField(diag::central_density::name(), mprim_scalar); physics->AddField(diag::localization_function::name(), mprim_scalar); physics->AddField(diag::entropy_z_0::name(), mprim_scalar); physics->AddField(p::gamma1::name(), mprim_scalar); @@ -324,8 +323,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { return ReduceOneVar>(md, fluid_cons::energy::name(), 0); }; auto MaxDensitySN = [](MeshData *md) { - History::ReduceCentralDensitySN(md); - return ReduceOneVar>(md, diag::central_density::name(), 0); + return ReduceOneVar>(md, fluid_prim::density::name(), 0); }; auto Mgain = [](MeshData *md) { return ReduceInGain>(md, fluid_prim::density::name(), diff --git a/src/phoebus_utils/variables.hpp b/src/phoebus_utils/variables.hpp index 14768a83..2c127ea2 100644 --- a/src/phoebus_utils/variables.hpp +++ b/src/phoebus_utils/variables.hpp @@ -122,7 +122,6 @@ namespace diagnostic_variables { VARIABLE_NONS(divb); VARIABLE_NONS(ratio_divv_cs); VARIABLE_NONS(entropy_z_0); -VARIABLE_NONS(central_density); VARIABLE_NONS(localization_function); VARIABLE_CUSTOM(divf, flux_divergence); VARIABLE_NONS(src_terms); From be008b1a7ec1b9992611dd223272005274efc172 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Tue, 12 Mar 2024 17:43:50 -0400 Subject: [PATCH 5/7] Update and cleanup --- src/analysis/history.cpp | 9 +++--- src/analysis/history.hpp | 49 ++++++++++++------------------ src/fluid/fluid.cpp | 17 ++--------- src/pgen/progenitor.cpp | 4 ++- src/progenitor/progenitordata.cpp | 19 ++++++++++++ src/radiation/cooling_function.cpp | 6 ++-- 6 files changed, 51 insertions(+), 53 deletions(-) diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp index 41694a26..d88e3d65 100644 --- a/src/analysis/history.cpp +++ b/src/analysis/history.cpp @@ -17,6 +17,7 @@ #include "geometry/geometry_utils.hpp" #include "history_utils.hpp" #include "phoebus_utils/relativity_utils.hpp" +#include namespace History { @@ -237,7 +238,9 @@ Real ReduceMagneticFluxPhi(MeshData *md) { } // Phi // SN analysis - +// ReduceLocalizationFunction is not used currently. However this function returns +// normalization function that is used for localizing quantities at the center, or at some +// particular case. For example SN diagnostics oftec computes quantities at 400 km. void ReduceLocalizationFunction(MeshData *md) { const auto ib = md->GetBoundsI(IndexDomain::interior); const auto jb = md->GetBoundsJ(IndexDomain::interior); @@ -274,8 +277,6 @@ void ReduceLocalizationFunction(MeshData *md) { v(b, diag::localization_function(), k, j, i) = std::exp(-r2 / sigma / sigma); }); -} // exp (this function returns normalization function that is used for localizing - // quantities at the center, or at some particular case. For example SN diagnostics - // oftec computes quantities at 400 km.) +} // exp } // namespace History diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp index 2cd97a93..87c45654 100644 --- a/src/analysis/history.hpp +++ b/src/analysis/history.hpp @@ -84,50 +84,39 @@ Real ReduceOneVar(MeshData *md, const std::string &varname, int idx = 0) { return result; } -template -Real ReduceInGain(MeshData *md, const std::string &varname, int idx = 0) { +template +Real ReduceInGain(MeshData *md, int idx = 0) { const auto ib = md->GetBoundsI(IndexDomain::interior); const auto jb = md->GetBoundsJ(IndexDomain::interior); const auto kb = md->GetBoundsK(IndexDomain::interior); - PackIndexMap imap; namespace iv = internal_variables; - std::vector vars = {varname, iv::GcovHeat::name(), iv::GcovCool::name()}; - const auto pack = md->PackVariables(vars, imap); - const auto ivar = imap[varname]; - const auto iheat = imap[iv::GcovHeat::name()].first; - const auto icool = imap[iv::GcovCool::name()].first; - auto *pmb = md->GetParentPointer(); - PARTHENON_REQUIRE(ivar.first >= 0, "Var must exist"); - PARTHENON_REQUIRE(ivar.second >= ivar.first + idx, "Var must exist"); + using parthenon::MakePackDescriptor; + Mesh *pmesh = md->GetMeshPointer(); + auto &resolved_pkgs = pmesh->resolved_packages; + const int ndim = pmesh->ndim; + static auto desc = + MakePackDescriptor(resolved_pkgs.get()); + auto v = desc.GetPack(md); + const int nblocks = v.GetNBlocks(); - // We choose to apply volume weighting when using the sum reduction. - // This assumes that for "Sum" variables, the variable is densitized, meaning - // it already contains a factor of the measure: sqrt(det(gamma)) Real result = 0.0; - Reducer_t reducer(result); - const bool volume_weighting = - std::is_same>::value || - std::is_same>::value || - std::is_same>::value; auto geom = Geometry::GetCoordinateSystem(md); parthenon::par_reduce( - parthenon::LoopPatternMDRange(), "Phoebus History for " + varname, DevExecSpace(), - 0, pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + parthenon::LoopPatternMDRange(), + "Phoebus History for integrals in gain region (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) { - // join is a Kokkos construct - // that automatically does the - // reduction operation locally Real gdet = geom.DetGamma(CellLocation::Cent, 0, k, j, i); - bool is_netheat = (pack(b, iheat, k, j, i) - pack(b, icool, k, j, i) > 0); - const auto &coords = pack.GetCoords(b); - const Real vol = volume_weighting ? coords.CellVolume(k, j, i) : 1.0; - reducer.join(lresult, - pack(b, ivar.first + idx, k, j, i) * gdet * vol * is_netheat); + bool is_netheat = + (v(b, iv::GcovHeat(), k, j, i) - v(b, iv::GcovCool(), k, j, i) > 0); + 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); }, - reducer); + Kokkos::Sum(result)); return result; } diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp index fe096581..8fdaedc7 100644 --- a/src/fluid/fluid.cpp +++ b/src/fluid/fluid.cpp @@ -13,6 +13,7 @@ #include "fluid.hpp" +#include "analysis/analysis.hpp" #include "analysis/history.hpp" #include "con2prim.hpp" #include "con2prim_robust.hpp" @@ -322,23 +323,11 @@ std::shared_ptr Initialize(ParameterInput *pin) { auto ReduceEn = [](MeshData *md) { return ReduceOneVar>(md, fluid_cons::energy::name(), 0); }; - auto MaxDensitySN = [](MeshData *md) { + auto MaxDensity = [](MeshData *md) { return ReduceOneVar>(md, fluid_prim::density::name(), 0); }; - auto Mgain = [](MeshData *md) { - return ReduceInGain>(md, fluid_prim::density::name(), - 0); - }; - auto Qgain = [](MeshData *md) { - return ReduceInGain>( - md, internal_variables::GcovHeat::name(), 0) - - ReduceInGain>( - md, internal_variables::GcovCool::name(), 0); - }; - hst_vars.emplace_back(HistoryOutputVar(HstMax, MaxDensitySN, "maximum density SN")); - hst_vars.emplace_back(HistoryOutputVar(HstSum, Mgain, "Mgain")); - hst_vars.emplace_back(HistoryOutputVar(HstSum, Qgain, "total net heat")); + hst_vars.emplace_back(HistoryOutputVar(HstMax, MaxDensity, "maximum density")); hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceMass, "total baryon number")); hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceEn, "total conserved energy tau")); diff --git a/src/pgen/progenitor.cpp b/src/pgen/progenitor.cpp index 093343d5..48f076fe 100644 --- a/src/pgen/progenitor.cpp +++ b/src/pgen/progenitor.cpp @@ -16,18 +16,20 @@ #include #include +#include "analysis/history.hpp" #include "geometry/geometry.hpp" #include "microphysics/eos_phoebus/eos_phoebus.hpp" #include "monopole_gr/monopole_gr.hpp" #include "pgen/pgen.hpp" #include "progenitor/progenitordata.hpp" +#include +#include using DataBox = Spiner::DataBox; // Homologously collapsing star. namespace homologous { void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { - std::cout << "Entered into ProbGen" << std::endl; const bool is_monopole_cart = (typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::MonopoleCart)); const bool is_monopole_sph = diff --git a/src/progenitor/progenitordata.cpp b/src/progenitor/progenitordata.cpp index 3fdac85f..92043be8 100644 --- a/src/progenitor/progenitordata.cpp +++ b/src/progenitor/progenitordata.cpp @@ -1,6 +1,7 @@ #include #include +#include "analysis/history.hpp" #include "ascii_reader.hpp" #include "geometry/geometry.hpp" #include "microphysics/eos_phoebus/eos_phoebus.hpp" @@ -107,6 +108,24 @@ std::shared_ptr Initialize(ParameterInput *pin) { params.Add("S_adm_dev", S_adm_dev); params.Add("Srr_adm_dev", Srr_adm_dev); + // Reductions + 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 *md) { + return ReduceInGain(md, 0); + }; + auto Qgain = [](MeshData *md) { + return ReduceInGain(md, 0) - + ReduceInGain(md, 0); + }; + hst_vars.emplace_back(HistoryOutputVar(HstSum, Mgain, "Mgain")); + hst_vars.emplace_back(HistoryOutputVar(HstSum, Qgain, "total net heat")); + params.Add(parthenon::hist_param_key, hst_vars); + return progenitor_pkg; } // Initialize diff --git a/src/radiation/cooling_function.cpp b/src/radiation/cooling_function.cpp index b6171ec6..06b39a26 100644 --- a/src/radiation/cooling_function.cpp +++ b/src/radiation/cooling_function.cpp @@ -279,10 +279,8 @@ TaskStatus CoolingFunctionCalculateFourForce(MeshData *rc, const double dt // detg included above Kokkos::atomic_add(&(v(b, iv::Gcov(mu), k, j, i)), -Gcov_coord[mu]); } - v(b, iv::GcovHeat(), k, j, i) = - v(b, p::density(), k, j, i) * density_conversion_factor * heat; - v(b, iv::GcovCool(), k, j, i) = - v(b, p::density(), k, j, i) * density_conversion_factor * cool; + v(b, iv::GcovHeat(), k, j, i) = cdensity * H; + v(b, iv::GcovCool(), k, j, i) = cdensity * C; Kokkos::atomic_add(&(v(b, iv::Gye(), k, j, i)), Jye); }); #else From 9337f9e92f804cc1a6c9f0bdf636cf812fd79aa3 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Wed, 13 Mar 2024 13:05:38 -0400 Subject: [PATCH 6/7] commenting localization function for now --- src/fluid/fluid.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp index 8fdaedc7..2e26da5e 100644 --- a/src/fluid/fluid.cpp +++ b/src/fluid/fluid.cpp @@ -208,7 +208,8 @@ std::shared_ptr Initialize(ParameterInput *pin) { physics->AddField(p::entropy::name(), mprim_scalar); physics->AddField(p::cs::name(), mprim_scalar); physics->AddField(diag::ratio_divv_cs::name(), mprim_scalar); - physics->AddField(diag::localization_function::name(), mprim_scalar); + // physics->AddField(diag::localization_function::name(), mprim_scalar); // (MG) + // currently not in use, turn on when needed. physics->AddField(diag::entropy_z_0::name(), mprim_scalar); physics->AddField(p::gamma1::name(), mprim_scalar); if (ye) { From 9146f83c6cc37e56a945636c9c801e913c058420 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Wed, 13 Mar 2024 13:16:44 -0400 Subject: [PATCH 7/7] check added --- src/analysis/history.hpp | 2 ++ src/fluid/fluid.cpp | 4 ++-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp index 87c45654..8dbfb164 100644 --- a/src/analysis/history.hpp +++ b/src/analysis/history.hpp @@ -98,6 +98,8 @@ Real ReduceInGain(MeshData *md, int idx = 0) { static auto desc = MakePackDescriptor(resolved_pkgs.get()); auto v = desc.GetPack(md); + PARTHENON_REQUIRE_THROWS(v.ContainsHost(0, iv::GcovHeat(), iv::GcovCool()), + "Must be doing SN simulation"); const int nblocks = v.GetNBlocks(); Real result = 0.0; diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp index 67cc8df0..cddff16f 100644 --- a/src/fluid/fluid.cpp +++ b/src/fluid/fluid.cpp @@ -208,8 +208,8 @@ std::shared_ptr Initialize(ParameterInput *pin) { physics->AddField(p::entropy::name(), mprim_scalar); physics->AddField(p::cs::name(), mprim_scalar); physics->AddField(diag::ratio_divv_cs::name(), mprim_scalar); - // physics->AddField(diag::localization_function::name(), mprim_scalar); // (MG) - // currently not in use, turn on when needed. + // physics->AddField(diag::localization_function::name(), mprim_scalar); + // (MG) currently not in use, turn on when needed. physics->AddField(diag::entropy_z_0::name(), mprim_scalar); physics->AddField(p::gamma1::name(), mprim_scalar); if (ye) {