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.

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

Add flat sparse packs #877

Merged
merged 20 commits into from
Jun 1, 2023
Merged
Show file tree
Hide file tree
Changes from 14 commits
Commits
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
1 change: 1 addition & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
Date: 2023-05-26

### Added (new features/APIs/variables/...)
- [[PR 877]](https://github.com/parthenon-hpc-lab/parthenon/pull/877) Add flat sparse packs
- [[PR 830]](https://github.com/parthenon-hpc-lab/parthenon/pull/830) Add particle output
- [[PR 840]](https://github.com/parthenon-hpc-lab/parthenon/pull/840) Generalized integrators infrastructure in a backwards compatible way
- [[PR 810]](https://github.com/parthenon-hpc-lab/parthenon/pull/810) Add suport for Ascent in-situ visualization
Expand Down
10 changes: 10 additions & 0 deletions doc/sphinx/src/interface/sparse.rst
Original file line number Diff line number Diff line change
Expand Up @@ -137,6 +137,16 @@ simulation. The set of possible operations for a sparse field are:
``lroberts36/merge-sparse-with-jdolence-sparse`` that is being used
in ``Riot``. This should probably be brought into ``develop``, since
the “sparse sparse pack” access pattern is probably not desirable.
- *Flatten outer indices:* Most common packs over ``MeshData`` produce
a 5-dimensional data structure, where the slowest two moving indices
are over ``MeshBlock``\ s then tensor/vector components of
``Variable``\ s. However, it is sometimes desirable to construct a
4-dimensional data structure where the slowest moving index ranges
over all variables and all blocks. So for example if we had a system
with 9 blocks and 5 variables per block, the slowest moving index
would be of size 45. This can be enabled ``SparsePack``s with the
``GetFlat`` series of factory functions or by passing the optional
``flat`` boolean into the constructor.

In comparison to a sparse field, a dense field only requires the
operation *Access*.
Expand Down
129 changes: 113 additions & 16 deletions src/interface/sparse_pack.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -125,10 +125,10 @@ class SparsePack : public SparsePackBase {
// The pack will be created and accessible on the device
template <class T>
static SparsePack Get(T *pmd, const std::vector<MetadataFlag> &flags = {},
bool fluxes = false, bool coarse = false) {
bool fluxes = false, bool coarse = false, bool flatten = false) {
const impl::PackDescriptor desc(std::vector<std::string>{Ts::name()...},
std::vector<bool>{Ts::regex()...}, flags, fluxes,
coarse);
coarse, flatten);
return SparsePack(SparsePackBase::GetPack(pmd, desc));
}

Expand All @@ -143,9 +143,9 @@ class SparsePack : public SparsePackBase {
template <class T, class VAR_VEC>
static std::tuple<SparsePack, SparsePackIdxMap>
Get(T *pmd, const VAR_VEC &vars, const std::vector<MetadataFlag> &flags = {},
bool fluxes = false, bool coarse = false) {
bool fluxes = false, bool coarse = false, bool flatten = false) {
static_assert(sizeof...(Ts) == 0, "Cannot create a string/type hybrid pack");
impl::PackDescriptor desc(vars, flags, fluxes, coarse);
impl::PackDescriptor desc(vars, flags, fluxes, coarse, flatten);
return {SparsePack(SparsePackBase::GetPack(pmd, desc)),
SparsePackBase::GetIdxMap(desc)};
}
Expand All @@ -155,7 +155,7 @@ class SparsePack : public SparsePackBase {
static SparsePack GetWithFluxes(T *pmd, const std::vector<MetadataFlag> &flags = {}) {
const bool coarse = false;
const bool fluxes = true;
return Get(pmd, flags, fluxes, coarse);
return Get(pmd, AddFlag_(flags), fluxes, coarse);
}

template <class T, class VAR_VEC>
Expand All @@ -164,7 +164,7 @@ class SparsePack : public SparsePackBase {
const std::vector<MetadataFlag> &flags = {}) {
const bool coarse = false;
const bool fluxes = true;
Get(pmd, vars, flags, fluxes, coarse);
Get(pmd, vars, AddFlag_(flags), fluxes, coarse);
}

template <class T>
Expand All @@ -183,73 +183,133 @@ class SparsePack : public SparsePackBase {
Get(pmd, vars, flags, fluxes, coarse);
}

template <class T>
static SparsePack GetFlat(T *pmd, const std::vector<MetadataFlag> &flags = {}) {
const bool coarse = false;
const bool fluxes = false;
const bool flatten = true;
return Get(pmd, flags, fluxes, coarse, flatten);
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved
}

template <class T, class VAR_VEC>
static std::tuple<SparsePack, SparsePackIdxMap>
GetFlat(T *pmd, const VAR_VEC &vars, const std::vector<MetadataFlag> &flags = {}) {
const bool coarse = false;
const bool fluxes = false;
const bool flatten = true;
return Get(pmd, vars, flags, fluxes, coarse, flatten);
}

template <class T>
static SparsePack GetFlatWithFluxes(T *pmd,
const std::vector<MetadataFlag> &flags = {}) {
const bool coarse = false;
const bool fluxes = true;
const bool flatten = true;
return Get(pmd, AddFlag_(flags), fluxes, coarse, flatten);
}

template <class T, class VAR_VEC>
static std::tuple<SparsePack, SparsePackIdxMap>
GetFlatWithFluxes(T *pmd, const VAR_VEC &vars,
const std::vector<MetadataFlag> &flags = {}) {
const bool coarse = false;
const bool fluxes = true;
const bool flatten = true;
return Get(pmd, vars, AddFlag_(flags), fluxes, coarse, flatten);
}

template <class T>
static SparsePack GetFlatWithCoarse(T *pmd,
const std::vector<MetadataFlag> &flags = {}) {
const bool coarse = true;
const bool fluxes = false;
const bool flatten = true;
return Get(pmd, flags, fluxes, coarse, flatten);
}

template <class T, class VAR_VEC>
static std::tuple<SparsePack, SparsePackIdxMap>
GetFlatWithCoarse(T *pmd, const VAR_VEC &vars,
const std::vector<MetadataFlag> &flags = {}) {
const bool coarse = true;
const bool fluxes = false;
const bool flatten = true;
return Get(pmd, vars, flags, fluxes, coarse, flatten);
}

// Methods for getting parts of the shape of the pack
KOKKOS_FORCEINLINE_FUNCTION
int GetNBlocks() const { return nblocks_; }

KOKKOS_FORCEINLINE_FUNCTION
int GetNDim() const { return ndim_; }

KOKKOS_FORCEINLINE_FUNCTION
int GetMaxNumberOfVars() const { return pack_.extent_int(2); }

KOKKOS_INLINE_FUNCTION
const Coordinates_t &GetCoordinates(const int b) const { return coords_(b)(); }
const Coordinates_t &GetCoordinates(const int b = 0) const { return coords_(b)(); }

// Bound overloads
KOKKOS_INLINE_FUNCTION int GetLowerBound(const int b) const { return 0; }
KOKKOS_INLINE_FUNCTION int GetLowerBound(const int b = 0) const { return 0; }
Copy link
Collaborator

Choose a reason for hiding this comment

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

I think this should change for flat packs. GetLowerBound(b) with b > 0 should return bounds_(0, b, 0) and GetUpperBound(b) should return bounds_(1, b, nvar_).


KOKKOS_INLINE_FUNCTION int GetUpperBound(const int b) const {
return bounds_(1, b, nvar_);
KOKKOS_INLINE_FUNCTION int GetUpperBound(const int b = 0) const {
return bounds_(1, !flat_ * b, !flat_ * nvar_);
Copy link
Collaborator

Choose a reason for hiding this comment

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

I don't think there should be any check for if the pack is flat here. This is why I think that there actually needs to be a new function that returns the total upper bound

KOKKOS_INLINE_FUNCTION int GetUpperBound() const { 
  assert(flat_); 
  return bounds_(1, GetNBlocks() - 1, nvar_); 
}

}

KOKKOS_INLINE_FUNCTION int GetLowerBound(const int b, PackIdx idx) const {
static_assert(sizeof...(Ts) == 0, "Cannot create a string/type hybrid pack");
PARTHENON_DEBUG_REQUIRE(!flat_, "Not valid for flat packs");
return bounds_(0, b, idx.VariableIdx());
}

KOKKOS_INLINE_FUNCTION int GetUpperBound(const int b, PackIdx idx) const {
static_assert(sizeof...(Ts) == 0, "Cannot create a string/type hybrid pack");
PARTHENON_DEBUG_REQUIRE(!flat_, "Not valid for flat packs");
return bounds_(1, b, idx.VariableIdx());
}

template <class TIn, REQUIRES(IncludesType<TIn, Ts...>::value)>
KOKKOS_INLINE_FUNCTION int GetLowerBound(const int b, const TIn &) const {
PARTHENON_DEBUG_REQUIRE(!flat_, "Not valid for flat packs");
const int vidx = GetTypeIdx<TIn, Ts...>::value;
return bounds_(0, b, vidx);
}

template <class TIn, REQUIRES(IncludesType<TIn, Ts...>::value)>
KOKKOS_INLINE_FUNCTION int GetUpperBound(const int b, const TIn &) const {
PARTHENON_DEBUG_REQUIRE(!flat_, "Not valid for flat packs");
const int vidx = GetTypeIdx<TIn, Ts...>::value;
return bounds_(1, b, vidx);
}

// Host Bound overloads
KOKKOS_INLINE_FUNCTION int GetLowerBoundHost(const int b) const { return 0; }
KOKKOS_INLINE_FUNCTION int GetLowerBoundHost(const int b = 0) const { return 0; }

KOKKOS_INLINE_FUNCTION int GetUpperBoundHost(const int b) const {
return bounds_h_(1, b, nvar_);
KOKKOS_INLINE_FUNCTION int GetUpperBoundHost(const int b = 0) const {
return bounds_h_(1, !flat_ * b, !flat_ * nvar_);
Copy link
Collaborator

Choose a reason for hiding this comment

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

Also no flat check here.

}

KOKKOS_INLINE_FUNCTION int GetLowerBoundHost(const int b, PackIdx idx) const {
static_assert(sizeof...(Ts) == 0);
PARTHENON_DEBUG_REQUIRE(!flat_, "Not valid for flat packs");
return bounds_h_(0, b, idx.VariableIdx());
}

KOKKOS_INLINE_FUNCTION int GetUpperBoundHost(const int b, PackIdx idx) const {
static_assert(sizeof...(Ts) == 0);
PARTHENON_DEBUG_REQUIRE(!flat_, "Not valid for flat packs");
return bounds_h_(1, b, idx.VariableIdx());
}

template <class TIn, REQUIRES(IncludesType<TIn, Ts...>::value)>
KOKKOS_INLINE_FUNCTION int GetLowerBoundHost(const int b, const TIn &) const {
PARTHENON_DEBUG_REQUIRE(!flat_, "Not valid for flat packs");
const int vidx = GetTypeIdx<TIn, Ts...>::value;
return bounds_h_(0, b, vidx);
}

template <class TIn, REQUIRES(IncludesType<TIn, Ts...>::value)>
KOKKOS_INLINE_FUNCTION int GetUpperBoundHost(const int b, const TIn &) const {
PARTHENON_DEBUG_REQUIRE(!flat_, "Not valid for flat packs");
const int vidx = GetTypeIdx<TIn, Ts...>::value;
return bounds_h_(1, b, vidx);
}
Expand All @@ -261,55 +321,75 @@ class SparsePack : public SparsePackBase {
KOKKOS_INLINE_FUNCTION
auto &operator()(const int b, PackIdx idx) const {
static_assert(sizeof...(Ts) == 0);
PARTHENON_DEBUG_REQUIRE(!flat_, "Accessor cannot be used for flat packs");
const int n = bounds_(0, b, idx.VariableIdx()) + idx.Offset();
return pack_(0, b, n);
}

template <class TIn, REQUIRES(IncludesType<TIn, Ts...>::value)>
KOKKOS_INLINE_FUNCTION auto &operator()(const int b, const TIn &t) const {
PARTHENON_DEBUG_REQUIRE(!flat_, "Accessor cannot be used for flat packs");
const int vidx = GetLowerBound(b, t) + t.idx;
return pack_(0, b, vidx);
}

KOKKOS_INLINE_FUNCTION
Real &operator()(const int b, const int idx, const int k, const int j,
const int i) const {
PARTHENON_DEBUG_REQUIRE(!flat_, "Accessor cannot be used for flat packs");
return pack_(0, b, idx)(k, j, i);
}
KOKKOS_INLINE_FUNCTION
Real &operator()(int idx, const int k, const int j, const int i) const {
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved
PARTHENON_DEBUG_REQUIRE(flat_, "Accessor valid only for flat packs");
return pack_(0, 0, idx)(k, j, i);
}

KOKKOS_INLINE_FUNCTION
Real &operator()(const int b, PackIdx idx, const int k, const int j,
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved
const int i) const {
static_assert(sizeof...(Ts) == 0, "Cannot create a string/type hybrid pack");
PARTHENON_DEBUG_REQUIRE(!flat_, "Accessor cannot be used for flat packs");
const int n = bounds_(0, b, idx.VariableIdx()) + idx.Offset();
return pack_(0, b, n)(k, j, i);
}

template <class TIn, REQUIRES(IncludesType<TIn, Ts...>::value)>
KOKKOS_INLINE_FUNCTION Real &operator()(const int b, const TIn &t, const int k,
const int j, const int i) const {
PARTHENON_DEBUG_REQUIRE(!flat_, "Accessor cannot be used for flat packs");
const int vidx = GetLowerBound(b, t) + t.idx;
return pack_(0, b, vidx)(k, j, i);
}

// flux() overloads
KOKKOS_INLINE_FUNCTION
auto &flux(const int b, const int dir, const int idx) const {
PARTHENON_DEBUG_REQUIRE(!flat_, "Accessor cannot be used for flat packs");
PARTHENON_DEBUG_REQUIRE(dir > 0 && dir < 4 && with_fluxes_, "Bad input to flux call");
return pack_(dir, b, idx);
}

KOKKOS_INLINE_FUNCTION
Real &flux(const int b, const int dir, const int idx, const int k, const int j,
const int i) const {
PARTHENON_DEBUG_REQUIRE(!flat_, "Accessor cannot be used for flat packs");
PARTHENON_DEBUG_REQUIRE(dir > 0 && dir < 4 && with_fluxes_, "Bad input to flux call");
return pack_(dir, b, idx)(k, j, i);
}

KOKKOS_INLINE_FUNCTION
Real &flux(const int dir, const int idx, const int k, const int j, const int i) const {
PARTHENON_DEBUG_REQUIRE(flat_, "Accessor must only be used for flat packs");
PARTHENON_DEBUG_REQUIRE(dir > 0 && dir < 4 && with_fluxes_, "Bad input to flux call");
return pack_(dir, 0, idx)(k, j, i);
}

KOKKOS_INLINE_FUNCTION
Real &flux(const int b, const int dir, PackIdx idx, const int k, const int j,
const int i) const {
static_assert(sizeof...(Ts) == 0, "Cannot create a string/type hybrid pack");
PARTHENON_DEBUG_REQUIRE(!flat_, "Accessor cannot be used for flat packs");
PARTHENON_DEBUG_REQUIRE(dir > 0 && dir < 4 && with_fluxes_, "Bad input to flux call");
const int n = bounds_(0, b, idx.VariableIdx()) + idx.Offset();
return pack_(dir, b, n)(k, j, i);
Expand All @@ -318,10 +398,27 @@ class SparsePack : public SparsePackBase {
template <class TIn, REQUIRES(IncludesType<TIn, Ts...>::value)>
KOKKOS_INLINE_FUNCTION Real &flux(const int b, const int dir, const TIn &t, const int k,
const int j, const int i) const {
PARTHENON_DEBUG_REQUIRE(!flat_, "Accessor cannot be used for flat packs");
PARTHENON_DEBUG_REQUIRE(dir > 0 && dir < 4 && with_fluxes_, "Bad input to flux call");
const int vidx = GetLowerBound(b, t) + t.idx;
return pack_(dir, b, vidx)(k, j, i);
}

private:
// Must make a copy of the vector, since input vector is const,
// and it may not even be an lvalue.
static std::vector<MetadataFlag> AddFlag_(const std::vector<MetadataFlag> &flags,
MetadataFlag mf = Metadata::WithFluxes) {
lroberts36 marked this conversation as resolved.
Show resolved Hide resolved
if (std::find(flags.begin(), flags.end(), mf) == flags.end()) {
std::vector<MetadataFlag> out;
out.reserve(flags.size() + 1);
out.insert(out.begin(), flags.begin(), flags.end());
out.push_back(mf);
return out;
} else {
return flags;
}
}
};

} // namespace parthenon
Expand Down
Loading