Skip to content

Commit

Permalink
computes max density instead of integral
Browse files Browse the repository at this point in the history
  • Loading branch information
@mari2895
mari2895 committed Mar 11, 2024
1 parent 3a72c26 commit 691f6c1
Show file tree
Hide file tree
Showing 4 changed files with 94 additions and 99 deletions.
140 changes: 70 additions & 70 deletions src/analysis/history.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 {

Expand Down Expand Up @@ -236,85 +236,85 @@ Real ReduceMagneticFluxPhi(MeshData<Real> *md) {
return 0.5 * result; // 0.5 \int detg B^r dx2 dx3
} // Phi

// SN analysis

void ReduceCentralDensitySN(MeshData<Real> *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<p::density, diag::central_density>(resolved_pkgs.get());
auto v = desc.GetPack(md);
const int nblocks = v.GetNBlocks();
auto geom = Geometry::GetCoordinateSystem(md);

//SN analysis

void ReduceCentralDensitySN(MeshData<Real> *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<p::density, diag::central_density>(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<Real>("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<Real> *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<diag::localization_function>(resolved_pkgs.get());
auto v = desc.GetPack(md);
const int nblocks = v.GetNBlocks();
auto geom = Geometry::GetCoordinateSystem(md);

}
void ReduceLocalizationFunction(MeshData<Real> *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<diag::localization_function>(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<Real>("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
35 changes: 18 additions & 17 deletions src/analysis/history.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,13 +17,13 @@
#include <string>
#include <vector>

#include <kokkos_abstraction.hpp>
#include <parthenon/package.hpp>
#include <parthenon/driver.hpp>
#include <utils/error_checking.hpp>
#include "geometry/geometry.hpp"
#include "geometry/geometry_utils.hpp"
#include "phoebus_utils/variables.hpp"
#include <kokkos_abstraction.hpp>
#include <parthenon/driver.hpp>
#include <parthenon/package.hpp>
#include <utils/error_checking.hpp>

using namespace parthenon::package::prelude;

Expand All @@ -46,7 +46,6 @@ Real ReduceMagneticFluxPhi(MeshData<Real> *md);
void ReduceCentralDensitySN(MeshData<Real> *md);
void ReduceLocalizationFunction(MeshData<Real> *md);


template <typename Reducer_t>
Real ReduceOneVar(MeshData<Real> *md, const std::string &varname, int idx = 0) {
const auto ib = md->GetBoundsI(IndexDomain::interior);
Expand All @@ -67,8 +66,10 @@ Real ReduceOneVar(MeshData<Real> *md, const std::string &varname, int idx = 0) {
Reducer_t reducer(result);

const bool volume_weighting =
std::is_same<Reducer_t, Kokkos::Sum<Real, HostExecSpace>>::value;

std::is_same<Reducer_t, Kokkos::Sum<Real, HostExecSpace>>::value ||
std::is_same<Reducer_t, Kokkos::Sum<Real, DevExecSpace>>::value ||
std::is_same<Reducer_t, Kokkos::Sum<Real, Kokkos::DefaultExecutionSpace>>::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,
Expand All @@ -91,11 +92,13 @@ Real ReduceInGain(MeshData<Real> *md, const std::string &varname, int idx = 0) {
const auto kb = md->GetBoundsK(IndexDomain::interior);

PackIndexMap imap;
std::vector<std::string> vars = {varname};
namespace iv = internal_variables;
std::vector<std::string> 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");

Expand All @@ -106,28 +109,26 @@ Real ReduceInGain(MeshData<Real> *md, const std::string &varname, int idx = 0) {
Reducer_t reducer(result);

const bool volume_weighting =
std::is_same<Reducer_t, Kokkos::Sum<Real, HostExecSpace>>::value;
parthenon::AllReduce<bool> *pdo_gain_reducer = rad->MutableParam<parthenon::AllReduce<bool>>("do_gain_reducer");
const bool do_gain = pdo_gain_reducer->val;
std::cout<<"do_gain="<<do_gain<<std::endl;

std::is_same<Reducer_t, Kokkos::Sum<Real, HostExecSpace>>::value ||
std::is_same<Reducer_t, Kokkos::Sum<Real, DevExecSpace>>::value ||
std::is_same<Reducer_t, Kokkos::Sum<Real, Kokkos::DefaultExecutionSpace>>::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
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_
16 changes: 5 additions & 11 deletions src/fluid/fluid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -311,11 +311,11 @@ std::shared_ptr<StateDescriptor> 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<Real> *md) {
return ReduceOneVar<Kokkos::Sum<Real>>(md, fluid_cons::density::name(), 0);
};
Expand All @@ -324,20 +324,14 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
};
auto CentralDensitySN = [](MeshData<Real> *md) {
History::ReduceCentralDensitySN(md);
return ReduceOneVar<Kokkos::Sum<Real, HostExecSpace>>(md, diag::central_density::name(), 0);

};
auto norm = [](MeshData<Real> *md) {
History::ReduceLocalizationFunction(md);
return ReduceOneVar<Kokkos::Sum<Real, HostExecSpace>>(md, diag::localization_function::name(), 0);
return ReduceOneVar<Kokkos::Max<Real>>(md, diag::central_density::name(), 0);
};
auto Mgain = [](MeshData<Real> *md) {
return ReduceInGain<Kokkos::Sum<Real, HostExecSpace>>(md, fluid_prim::density::name(), 0);
return ReduceInGain<Kokkos::Sum<Real, HostExecSpace>>(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"));
Expand Down
2 changes: 1 addition & 1 deletion src/phoebus_utils/variables.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down

0 comments on commit 691f6c1

Please sign in to comment.