From 3b3bc9f3e48dbcb8377a6dce16aade48bbec33e5 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Wed, 13 May 2026 00:30:29 +0200 Subject: [PATCH 01/12] Handle strided field indexing --- src/Field/BareField.h | 9 ++ src/Field/BareField.hpp | 25 +++ src/Field/IndexedBareField.h | 170 ++++++++++++++++++++ src/Field/SparseIndexedBareField.h | 125 +++++++++++++++ src/Index/Index.h | 2 +- src/Index/Index.hpp | 32 +++- src/Index/SIndex.h | 120 +++++++++++++++ src/IpplCore.h | 1 + unit_tests/BareField/BareField.cpp | 239 +++++++++++++++++++++++++++++ unit_tests/CMakeLists.txt | 1 + unit_tests/Index/CMakeLists.txt | 4 + unit_tests/Index/Index.cpp | 46 ++++++ 12 files changed, 770 insertions(+), 4 deletions(-) create mode 100644 src/Field/IndexedBareField.h create mode 100644 src/Field/SparseIndexedBareField.h create mode 100644 src/Index/SIndex.h create mode 100644 unit_tests/Index/CMakeLists.txt create mode 100644 unit_tests/Index/Index.cpp diff --git a/src/Field/BareField.h b/src/Field/BareField.h index 33908ef53..ff2c014d0 100644 --- a/src/Field/BareField.h +++ b/src/Field/BareField.h @@ -20,6 +20,7 @@ #include "Field/HaloCells.h" #include "FieldLayout/FieldLayout.h" +#include "Index/SIndex.h" namespace ippl { class Index; @@ -166,6 +167,14 @@ namespace ippl { return dview_m(args...); } + auto operator[](const Domain_t& domain); + + auto operator[](const Index& index); + + auto operator[](int index); + + auto operator[](const SIndex& sindex); + view_type& getView() { return dview_m; } const view_type& getView() const { return dview_m; } diff --git a/src/Field/BareField.hpp b/src/Field/BareField.hpp index efcddb5f1..44a3f1ff9 100644 --- a/src/Field/BareField.hpp +++ b/src/Field/BareField.hpp @@ -14,6 +14,8 @@ #include "Utility/Inform.h" #include "Utility/IpplInfo.h" +#include "Field/IndexedBareField.h" +#include "Field/SparseIndexedBareField.h" namespace Kokkos { template struct reduction_identity> { @@ -139,6 +141,29 @@ namespace ippl { Kokkos::resize(dview_m, args...); } + template + auto BareField::operator[](const Domain_t& domain) { + return IndexedBareField>(*this, domain); + } + + template + auto BareField::operator[](const Index& index) { + Domain_t domain = getDomain(); + domain[0] = index; + return IndexedBareField, 1>(*this, domain); + } + + template + auto BareField::operator[](int index) { + static_assert(Dim == 1, "Integer overload is only available for one-dimensional fields."); + return operator[](Index(index, index)); + } + + template + auto BareField::operator[](const SIndex& sindex) { + return SparseIndexedBareField>(*this, sindex); + } + template void BareField::fillHalo() { if (layout_m->comm.size() > 1) { diff --git a/src/Field/IndexedBareField.h b/src/Field/IndexedBareField.h new file mode 100644 index 000000000..195c9a1d2 --- /dev/null +++ b/src/Field/IndexedBareField.h @@ -0,0 +1,170 @@ +// +// Class IndexedBareField +// Lightweight indexed view into a BareField. +// +#ifndef IPPL_INDEXED_BARE_FIELD_H +#define IPPL_INDEXED_BARE_FIELD_H + +#include + +#include "Expression/IpplExpressions.h" +#include "Index/NDIndex.h" +#include "Utility/ParallelDispatch.h" + +namespace ippl { + + template + class IndexedBareField + : public detail::Expression, + sizeof(Field) + sizeof(NDIndex)> { + public: + using field_type = Field; + using value_type = typename field_type::value_type; + using view_type = typename field_type::view_type; + using execution_space = typename field_type::execution_space; + using Domain_t = NDIndex; + + constexpr static unsigned dim = field_type::dim; + constexpr static size_t expression_size = sizeof(Field) + sizeof(Domain_t); + + IndexedBareField(field_type& field, const Domain_t& domain) + : view_m(field.getView()) + , owned_m(field.getOwned()) + , domain_m(domain) + , nghost_m(field.getNghost()) {} + + IndexedBareField(const IndexedBareField&) = default; + + IndexedBareField(const view_type& view, const Domain_t& owned, const Domain_t& domain, + int nghost) + : view_m(view) + , owned_m(owned) + , domain_m(domain) + , nghost_m(nghost) {} + + auto operator[](const Index& index) const { + static_assert(Brackets < dim, "Too many Index brackets for field dimension."); + Domain_t domain = domain_m; + domain[Brackets] = index; + return IndexedBareField(view_m, owned_m, domain, nghost_m); + } + + auto operator[](int index) const { + return operator[](Index(index, index)); + } + + template + KOKKOS_INLINE_FUNCTION value_type operator()(Args... args) const { + static_assert(Brackets == dim, "IndexedBareField expression requires all dimensions."); + static_assert(sizeof...(Args) == dim); + typename RangePolicy::index_array_type rel{args...}; + return apply(view_m, relativeToView(rel)); + } + + IndexedBareField& operator=(value_type value) { + static_assert(Brackets == dim, "IndexedBareField assignment requires all dimensions."); + assign(ValueAssign{value}); + return *this; + } + + template + IndexedBareField& operator=(const detail::Expression& expr) { + static_assert(Brackets == dim, "IndexedBareField assignment requires all dimensions."); + using capture_type = detail::CapturedExpression; + capture_type expr_ = reinterpret_cast(expr); + assign(ExpressionAssign{expr_}); + return *this; + } + + IndexedBareField& operator=(const IndexedBareField& rhs) { + static_assert(Brackets == dim, "IndexedBareField assignment requires all dimensions."); + const detail::Expression, expression_size>& expr = rhs; + return operator=(expr); + } + + const Domain_t& getDomain() const { return domain_m; } + + private: + view_type view_m; + Domain_t owned_m; + Domain_t domain_m; + int nghost_m; + + struct ValueAssign { + value_type value; + + template + KOKKOS_INLINE_FUNCTION value_type operator()(const Coords&) const { + return value; + } + }; + + template + struct ExpressionAssign { + detail::CapturedExpression expr; + + template + KOKKOS_INLINE_FUNCTION auto operator()(const Coords& rel) const { + return apply(expr, rel); + } + }; + + template + KOKKOS_INLINE_FUNCTION auto relativeToView(const Coords& rel) const { + typename RangePolicy::index_array_type viewCoords; + for (unsigned d = 0; d < dim; ++d) { + const auto global = domain_m[d].first() + rel[d] * domain_m[d].stride(); + viewCoords[d] = (global - owned_m[d].first()) / owned_m[d].stride() + nghost_m; + } + return viewCoords; + } + + template + KOKKOS_INLINE_FUNCTION auto globalToRelative(const Coords& global) const { + typename RangePolicy::index_array_type rel; + for (unsigned d = 0; d < dim; ++d) { + rel[d] = (global[d] - domain_m[d].first()) / domain_m[d].stride(); + } + return rel; + } + + template + KOKKOS_INLINE_FUNCTION auto globalToView(const Coords& global) const { + typename RangePolicy::index_array_type viewCoords; + for (unsigned d = 0; d < dim; ++d) { + viewCoords[d] = (global[d] - owned_m[d].first()) / owned_m[d].stride() + nghost_m; + } + return viewCoords; + } + + template + void assign(const Functor& functor) { + Domain_t local = domain_m.intersect(owned_m); + if (local.empty()) { + return; + } + + using range_policy_type = RangePolicy; + using index_type = typename range_policy_type::index_type; + Kokkos::Array begin, end; + for (unsigned d = 0; d < dim; ++d) { + begin[d] = 0; + end[d] = local[d].length(); + } + + ippl::parallel_for( + "IndexedBareField::operator=", createRangePolicy(begin, end), + KOKKOS_CLASS_LAMBDA(const typename range_policy_type::index_array_type& args) { + typename range_policy_type::index_array_type global; + for (unsigned d = 0; d < dim; ++d) { + global[d] = local[d].first() + args[d] * local[d].stride(); + } + const auto rel = globalToRelative(global); + apply(view_m, globalToView(global)) = functor(rel); + }); + } + }; + +} // namespace ippl + +#endif diff --git a/src/Field/SparseIndexedBareField.h b/src/Field/SparseIndexedBareField.h new file mode 100644 index 000000000..189f25377 --- /dev/null +++ b/src/Field/SparseIndexedBareField.h @@ -0,0 +1,125 @@ +// +// Class SparseIndexedBareField +// Sparse point view into a BareField. +// +#ifndef IPPL_SPARSE_INDEXED_BARE_FIELD_H +#define IPPL_SPARSE_INDEXED_BARE_FIELD_H + +#include + +#include "Expression/IpplExpressions.h" +#include "Index/SIndex.h" +#include "Utility/ParallelDispatch.h" + +namespace ippl { + + template + class SparseIndexedBareField + : public detail::Expression, + sizeof(Field) + sizeof(SIndex)> { + public: + using field_type = Field; + using value_type = typename field_type::value_type; + using view_type = typename field_type::view_type; + using execution_space = typename field_type::execution_space; + using sindex_type = SIndex; + using point_type = typename sindex_type::point_type; + using points_view_type = Kokkos::View; + + constexpr static unsigned dim = 1; + constexpr static unsigned field_dim = field_type::dim; + constexpr static size_t expression_size = sizeof(Field) + sizeof(sindex_type); + + SparseIndexedBareField(field_type& field, const sindex_type& sindex) + : view_m(field.getView()) + , owned_m(field.getOwned()) + , points_m(sindex.template getDevicePoints()) + , nghost_m(field.getNghost()) {} + + SparseIndexedBareField(const SparseIndexedBareField&) = default; + + template + KOKKOS_INLINE_FUNCTION value_type operator()(Idx i) const { + return apply(view_m, pointToView(points_m(i))); + } + + SparseIndexedBareField& operator=(value_type value) { + assign(ValueAssign{value}); + return *this; + } + + template + SparseIndexedBareField& operator=(const detail::Expression& expr) { + using capture_type = detail::CapturedExpression; + capture_type expr_ = reinterpret_cast(expr); + assign(ExpressionAssign{expr_}); + return *this; + } + + SparseIndexedBareField& operator=(const SparseIndexedBareField& rhs) { + const detail::Expression& expr = rhs; + return operator=(expr); + } + + private: + view_type view_m; + NDIndex owned_m; + points_view_type points_m; + int nghost_m; + + struct ValueAssign { + value_type value; + + template + KOKKOS_INLINE_FUNCTION value_type operator()(const Coords&) const { + return value; + } + }; + + template + struct ExpressionAssign { + detail::CapturedExpression expr; + + template + KOKKOS_INLINE_FUNCTION auto operator()(const Coords& rel) const { + return apply(expr, rel); + } + }; + + KOKKOS_INLINE_FUNCTION bool isLocal(const point_type& point) const { + for (unsigned d = 0; d < field_dim; ++d) { + if (!owned_m[d].contains(Index(point[d], point[d]))) { + return false; + } + } + return true; + } + + KOKKOS_INLINE_FUNCTION auto pointToView(const point_type& point) const { + typename RangePolicy::index_array_type viewCoords; + for (unsigned d = 0; d < field_dim; ++d) { + viewCoords[d] = (point[d] - owned_m[d].first()) / owned_m[d].stride() + nghost_m; + } + return viewCoords; + } + + template + void assign(const Functor& functor) { + ippl::parallel_for( + "SparseIndexedBareField::operator=", getRangePolicy(points_m), + KOKKOS_CLASS_LAMBDA( + const typename RangePolicy<1, execution_space>::index_array_type& args) { + const auto i = args[0]; + const point_type point = points_m(i); + if (isLocal(point)) { + typename RangePolicy<1, execution_space>::index_array_type rel; + rel[0] = i; + apply(view_m, pointToView(point)) = functor(rel); + } + }); + } + }; + +} // namespace ippl + +#endif diff --git a/src/Index/Index.h b/src/Index/Index.h index b390f5f6e..13c725dd5 100644 --- a/src/Index/Index.h +++ b/src/Index/Index.h @@ -206,7 +206,7 @@ namespace ippl { // Test to see if there is any overlap between two Indexes. KOKKOS_INLINE_FUNCTION bool touches(const Index& a) const; - // Test to see if one contains another (endpoints only) + // Test to see if one contains all points of another. KOKKOS_INLINE_FUNCTION bool contains(const Index& a) const; // Split one into two. KOKKOS_INLINE_FUNCTION bool split(Index& l, Index& r) const; diff --git a/src/Index/Index.hpp b/src/Index/Index.hpp index 4c930c3fc..9c669ae88 100644 --- a/src/Index/Index.hpp +++ b/src/Index/Index.hpp @@ -152,11 +152,37 @@ namespace ippl { } KOKKOS_INLINE_FUNCTION bool Index::touches(const Index& a) const { - return (min() <= a.max()) && (max() >= a.min()); + return !intersect(a).empty(); } KOKKOS_INLINE_FUNCTION bool Index::contains(const Index& a) const { - return (min() <= a.min()) && (max() >= a.max()); + if (a.empty()) { + return true; + } + const int thisMin = min(); + const int thisMax = max(); + const int thatMin = a.min(); + const int thatMax = a.max(); + int thisStride = stride(); + int thatStride = a.stride(); + if (thisStride < 0) { + thisStride = -thisStride; + } + if (thatStride < 0) { + thatStride = -thatStride; + } + + const bool endpointsContained = (thisMin <= thatMin) && (thisMax >= thatMax); + if (!endpointsContained || thisStride == 1) { + return endpointsContained; + } + + const bool endpointsOnStride = + ((thatMin - thisMin) % thisStride == 0) && ((thisMax - thatMax) % thisStride == 0); + if (a.length() == 1) { + return endpointsOnStride; + } + return endpointsOnStride && (thatStride % thisStride == 0); } KOKKOS_INLINE_FUNCTION bool Index::split(Index& l, Index& r) const { @@ -261,7 +287,7 @@ namespace ippl { KOKKOS_INLINE_FUNCTION Index Index::grow(int ncells) const { Index index; - index.first_m = this->first_m - ncells; + index.first_m = this->first_m - ncells * this->stride_m; index.length_m = this->length_m + 2 * ncells; index.stride_m = this->stride_m; diff --git a/src/Index/SIndex.h b/src/Index/SIndex.h new file mode 100644 index 000000000..d299faa83 --- /dev/null +++ b/src/Index/SIndex.h @@ -0,0 +1,120 @@ +// +// Class SIndex +// Sparse set of field index points. +// +#ifndef IPPL_SINDEX_H +#define IPPL_SINDEX_H + +#include + +#include + +#include "FieldLayout/FieldLayout.h" +#include "Index/NDIndex.h" +#include "Types/Vector.h" + +namespace ippl { + + template + class SIndex { + public: + using point_type = Vector; + + SIndex() = default; + + explicit SIndex(FieldLayout& layout) + : layout_m(&layout) + , domain_m(layout.getDomain()) {} + + void initialize(FieldLayout& layout) { + layout_m = &layout; + domain_m = layout.getDomain(); + points_m.clear(); + } + + bool needInitialize() const { return layout_m == nullptr; } + + bool addIndex(const point_type& point) { + if (!contains(domain_m, point) || hasIndex(point)) { + return false; + } + points_m.push_back(point); + return true; + } + + bool addIndex(const NDIndex& point) { + point_type values; + for (unsigned d = 0; d < Dim; ++d) { + if (point[d].length() != 1) { + return false; + } + values[d] = point[d].first(); + } + return addIndex(values); + } + + void clear() { points_m.clear(); } + + std::size_t size() const { return points_m.size(); } + + bool empty() const { return points_m.empty(); } + + const std::vector& points() const { return points_m; } + + const NDIndex& getDomain() const { return domain_m; } + + FieldLayout& getFieldLayout() const { return *layout_m; } + + bool hasIndex(const point_type& point) const { + for (const auto& existing : points_m) { + bool same = true; + for (unsigned d = 0; d < Dim; ++d) { + same = same && existing[d] == point[d]; + } + if (same) { + return true; + } + } + return false; + } + + bool hasIndex(const NDIndex& point) const { + point_type values; + for (unsigned d = 0; d < Dim; ++d) { + if (point[d].length() != 1) { + return false; + } + values[d] = point[d].first(); + } + return hasIndex(values); + } + + template + Kokkos::View getDevicePoints() const { + Kokkos::View devicePoints("SIndex::points", points_m.size()); + auto hostPoints = Kokkos::create_mirror_view(devicePoints); + for (std::size_t i = 0; i < points_m.size(); ++i) { + hostPoints(i) = points_m[i]; + } + Kokkos::deep_copy(devicePoints, hostPoints); + return devicePoints; + } + + private: + FieldLayout* layout_m = nullptr; + NDIndex domain_m; + std::vector points_m; + + static bool contains(const NDIndex& domain, const point_type& point) { + for (unsigned d = 0; d < Dim; ++d) { + if (!domain[d].contains(Index(point[d], point[d]))) { + return false; + } + } + return true; + } + }; + +} // namespace ippl + +#endif diff --git a/src/IpplCore.h b/src/IpplCore.h index 4ee5d796c..71f1e8b8b 100644 --- a/src/IpplCore.h +++ b/src/IpplCore.h @@ -18,6 +18,7 @@ // #include "Index/Index.h" // #include "Index/NDIndex.h" +#include "Index/SIndex.h" #include "FieldLayout/FieldLayout.h" #include "FieldLayout/SubFieldLayout.h" diff --git a/unit_tests/BareField/BareField.cpp b/unit_tests/BareField/BareField.cpp index fe439dc6f..2e3725867 100644 --- a/unit_tests/BareField/BareField.cpp +++ b/unit_tests/BareField/BareField.cpp @@ -106,6 +106,245 @@ TYPED_TEST(BareFieldTest, DeepCopy) { }); } +TYPED_TEST(BareFieldTest, IndexedSubdomainScalarAssignment) { + using T = typename TestFixture::value_type; + + auto& field = this->field; + *field = 0; + + ippl::NDIndex subdomain = field->getDomain(); + subdomain[0] = + ippl::Index(subdomain[0].first(), subdomain[0].last(), 2 * subdomain[0].stride()); + + (*field)[subdomain] = T(7); + + auto mirror = field->getHostMirror(); + Kokkos::deep_copy(mirror, field->getView()); + + const auto owned = field->getOwned(); + const int nghost = field->getNghost(); + nestedViewLoop(mirror, nghost, [&](const Idx... args) { + std::array pointIndices; + const std::array viewCoords{static_cast(args)...}; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + const int global = owned[d].first() + + static_cast(viewCoords[d] - nghost) * owned[d].stride(); + pointIndices[d] = ippl::Index(global, global); + } + auto point = std::make_from_tuple>(pointIndices); + const T expected = subdomain.contains(point) ? T(7) : T(0); + assertEqual(expected, mirror(args...)); + }); +} + +TYPED_TEST(BareFieldTest, IndexedSubdomainExpressionAssignment) { + using T = typename TestFixture::value_type; + + auto& lhs = this->field; + lhs->operator=(0); + + typename TestFixture::field_type rhs(this->layout); + rhs = T(5); + + ippl::NDIndex subdomain = lhs->getDomain(); + subdomain[0] = + ippl::Index(subdomain[0].first(), subdomain[0].last(), 2 * subdomain[0].stride()); + + (*lhs)[subdomain] = rhs[subdomain] + T(2); + + auto mirror = lhs->getHostMirror(); + Kokkos::deep_copy(mirror, lhs->getView()); + + const auto owned = lhs->getOwned(); + const int nghost = lhs->getNghost(); + nestedViewLoop(mirror, nghost, [&](const Idx... args) { + std::array pointIndices; + const std::array viewCoords{static_cast(args)...}; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + const int global = owned[d].first() + + static_cast(viewCoords[d] - nghost) * owned[d].stride(); + pointIndices[d] = ippl::Index(global, global); + } + auto point = std::make_from_tuple>(pointIndices); + const T expected = subdomain.contains(point) ? T(7) : T(0); + assertEqual(expected, mirror(args...)); + }); +} + +TYPED_TEST(BareFieldTest, ChainedIndexedSubdomainScalarAssignment) { + using T = typename TestFixture::value_type; + + auto& field = this->field; + *field = 0; + + ippl::NDIndex subdomain = field->getDomain(); + subdomain[0] = + ippl::Index(subdomain[0].first(), subdomain[0].last(), 2 * subdomain[0].stride()); + + if constexpr (TestFixture::dim == 1) { + (*field)[subdomain[0]] = T(7); + } else if constexpr (TestFixture::dim == 2) { + (*field)[subdomain[0]][subdomain[1]] = T(7); + } else if constexpr (TestFixture::dim == 3) { + (*field)[subdomain[0]][subdomain[1]][subdomain[2]] = T(7); + } else if constexpr (TestFixture::dim == 4) { + (*field)[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]] = T(7); + } else if constexpr (TestFixture::dim == 5) { + (*field)[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]][subdomain[4]] = T(7); + } else if constexpr (TestFixture::dim == 6) { + (*field)[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]][subdomain[4]] + [subdomain[5]] = T(7); + } + + auto mirror = field->getHostMirror(); + Kokkos::deep_copy(mirror, field->getView()); + + const auto owned = field->getOwned(); + const int nghost = field->getNghost(); + nestedViewLoop(mirror, nghost, [&](const Idx... args) { + std::array pointIndices; + const std::array viewCoords{static_cast(args)...}; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + const int global = owned[d].first() + + static_cast(viewCoords[d] - nghost) * owned[d].stride(); + pointIndices[d] = ippl::Index(global, global); + } + auto point = std::make_from_tuple>(pointIndices); + const T expected = subdomain.contains(point) ? T(7) : T(0); + assertEqual(expected, mirror(args...)); + }); +} + +TYPED_TEST(BareFieldTest, ChainedIndexedSubdomainExpressionAssignment) { + using T = typename TestFixture::value_type; + + auto& lhs = this->field; + lhs->operator=(0); + + typename TestFixture::field_type rhs(this->layout); + rhs = T(5); + + ippl::NDIndex subdomain = lhs->getDomain(); + subdomain[0] = + ippl::Index(subdomain[0].first(), subdomain[0].last(), 2 * subdomain[0].stride()); + + if constexpr (TestFixture::dim == 1) { + (*lhs)[subdomain[0]] = rhs[subdomain[0]] + T(2); + } else if constexpr (TestFixture::dim == 2) { + (*lhs)[subdomain[0]][subdomain[1]] = rhs[subdomain[0]][subdomain[1]] + T(2); + } else if constexpr (TestFixture::dim == 3) { + (*lhs)[subdomain[0]][subdomain[1]][subdomain[2]] = + rhs[subdomain[0]][subdomain[1]][subdomain[2]] + T(2); + } else if constexpr (TestFixture::dim == 4) { + (*lhs)[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]] = + rhs[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]] + T(2); + } else if constexpr (TestFixture::dim == 5) { + (*lhs)[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]][subdomain[4]] = + rhs[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]][subdomain[4]] + T(2); + } else if constexpr (TestFixture::dim == 6) { + (*lhs)[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]][subdomain[4]] + [subdomain[5]] = + rhs[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]][subdomain[4]] + [subdomain[5]] + + T(2); + } + + auto mirror = lhs->getHostMirror(); + Kokkos::deep_copy(mirror, lhs->getView()); + + const auto owned = lhs->getOwned(); + const int nghost = lhs->getNghost(); + nestedViewLoop(mirror, nghost, [&](const Idx... args) { + std::array pointIndices; + const std::array viewCoords{static_cast(args)...}; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + const int global = owned[d].first() + + static_cast(viewCoords[d] - nghost) * owned[d].stride(); + pointIndices[d] = ippl::Index(global, global); + } + auto point = std::make_from_tuple>(pointIndices); + const T expected = subdomain.contains(point) ? T(7) : T(0); + assertEqual(expected, mirror(args...)); + }); +} + +TYPED_TEST(BareFieldTest, SparseIndexedSubdomainScalarAssignment) { + using T = typename TestFixture::value_type; + + auto& field = this->field; + *field = 0; + + ippl::SIndex sindex(this->layout); + const auto domain = field->getDomain(); + + for (int i = domain[0].first(); i <= domain[0].last(); i += 2 * domain[0].stride()) { + typename ippl::SIndex::point_type point; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + point[d] = domain[d].first(); + } + point[0] = i; + sindex.addIndex(point); + } + + (*field)[sindex] = T(7); + + auto mirror = field->getHostMirror(); + Kokkos::deep_copy(mirror, field->getView()); + + const auto owned = field->getOwned(); + const int nghost = field->getNghost(); + nestedViewLoop(mirror, nghost, [&](const Idx... args) { + typename ippl::SIndex::point_type point; + const std::array viewCoords{static_cast(args)...}; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + point[d] = owned[d].first() + + static_cast(viewCoords[d] - nghost) * owned[d].stride(); + } + const T expected = sindex.hasIndex(point) ? T(7) : T(0); + assertEqual(expected, mirror(args...)); + }); +} + +TYPED_TEST(BareFieldTest, SparseIndexedSubdomainExpressionAssignment) { + using T = typename TestFixture::value_type; + + auto& lhs = this->field; + lhs->operator=(0); + + typename TestFixture::field_type rhs(this->layout); + rhs = T(5); + + ippl::SIndex sindex(this->layout); + const auto domain = lhs->getDomain(); + + for (int i = domain[0].first(); i <= domain[0].last(); i += 2 * domain[0].stride()) { + typename ippl::SIndex::point_type point; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + point[d] = domain[d].first(); + } + point[0] = i; + sindex.addIndex(point); + } + + (*lhs)[sindex] = rhs[sindex] + T(2); + + auto mirror = lhs->getHostMirror(); + Kokkos::deep_copy(mirror, lhs->getView()); + + const auto owned = lhs->getOwned(); + const int nghost = lhs->getNghost(); + nestedViewLoop(mirror, nghost, [&](const Idx... args) { + typename ippl::SIndex::point_type point; + const std::array viewCoords{static_cast(args)...}; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + point[d] = owned[d].first() + + static_cast(viewCoords[d] - nghost) * owned[d].stride(); + } + const T expected = sindex.hasIndex(point) ? T(7) : T(0); + assertEqual(expected, mirror(args...)); + }); +} + TYPED_TEST(BareFieldTest, Sum) { using T = typename TestFixture::value_type; diff --git a/unit_tests/CMakeLists.txt b/unit_tests/CMakeLists.txt index af778d67b..90f555c0b 100644 --- a/unit_tests/CMakeLists.txt +++ b/unit_tests/CMakeLists.txt @@ -30,6 +30,7 @@ add_subdirectory(FEM) add_subdirectory(Interpolation) add_subdirectory(FFT) add_subdirectory(Field) +add_subdirectory(Index) add_subdirectory(Meshes) add_subdirectory(Particle) add_subdirectory(PIC) diff --git a/unit_tests/Index/CMakeLists.txt b/unit_tests/Index/CMakeLists.txt new file mode 100644 index 000000000..54cd13576 --- /dev/null +++ b/unit_tests/Index/CMakeLists.txt @@ -0,0 +1,4 @@ +file(RELATIVE_PATH _relPath "${PROJECT_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}") +message(STATUS "Adding unit tests found in ${_relPath}") + +add_ippl_test(Index USE_GTEST_MAIN) diff --git a/unit_tests/Index/Index.cpp b/unit_tests/Index/Index.cpp new file mode 100644 index 000000000..0fc0e11aa --- /dev/null +++ b/unit_tests/Index/Index.cpp @@ -0,0 +1,46 @@ +// +// Unit test Index +// Test strided Index semantics. +// +#include "Ippl.h" + +#include "gtest/gtest.h" + +TEST(IndexTest, StridedTouchesRequiresCommonElement) { + ippl::Index even(0, 10, 2); + ippl::Index odd(1, 9, 2); + + EXPECT_FALSE(even.touches(odd)); + EXPECT_FALSE(odd.touches(even)); + EXPECT_TRUE(even.intersect(odd).empty()); +} + +TEST(IndexTest, StridedContainsRequiresContainedElements) { + ippl::Index even(0, 10, 2); + ippl::Index odd(1, 9, 2); + ippl::Index everyFourth(0, 8, 4); + + EXPECT_FALSE(even.contains(odd)); + EXPECT_TRUE(even.contains(everyFourth)); + EXPECT_FALSE(everyFourth.contains(even)); +} + +TEST(IndexTest, GrowUsesIndexStrideForCoordinateBounds) { + ippl::Index even(2, 10, 2); + ippl::Index grown = even.grow(2); + + EXPECT_EQ(grown.first(), -2); + EXPECT_EQ(grown.last(), 14); + EXPECT_EQ(grown.stride(), 2); + EXPECT_EQ(grown.length(), 9); +} + +TEST(IndexTest, NegativeStrideGrowKeepsOrdering) { + ippl::Index descending(10, 2, -2); + ippl::Index grown = descending.grow(2); + + EXPECT_EQ(grown.first(), 14); + EXPECT_EQ(grown.last(), -2); + EXPECT_EQ(grown.stride(), -2); + EXPECT_EQ(grown.length(), 9); +} From c81a954d662360f4cc1fb07e4e5e8e122d1a871f Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Wed, 13 May 2026 00:39:31 +0200 Subject: [PATCH 02/12] Add sparse index offsets --- src/Index/SIndex.h | 219 ++++++++++++++++++++++++++++- src/Index/SOffset.h | 133 ++++++++++++++++++ src/IpplCore.h | 1 + unit_tests/BareField/BareField.cpp | 125 ++++++++++++++++ 4 files changed, 472 insertions(+), 6 deletions(-) create mode 100644 src/Index/SOffset.h diff --git a/src/Index/SIndex.h b/src/Index/SIndex.h index d299faa83..165952d48 100644 --- a/src/Index/SIndex.h +++ b/src/Index/SIndex.h @@ -7,10 +7,13 @@ #include +#include +#include #include #include "FieldLayout/FieldLayout.h" #include "Index/NDIndex.h" +#include "Index/SOffset.h" #include "Types/Vector.h" namespace ippl { @@ -26,6 +29,14 @@ namespace ippl { : layout_m(&layout) , domain_m(layout.getDomain()) {} + SIndex(const SIndex&) = default; + + SIndex(const SIndex& sindex, const SOffset& offset) + : layout_m(sindex.layout_m) + , domain_m(sindex.domain_m) + , points_m(sindex.points_m) + , offset_m(sindex.offset_m + offset) {} + void initialize(FieldLayout& layout) { layout_m = &layout; domain_m = layout.getDomain(); @@ -34,6 +45,59 @@ namespace ippl { bool needInitialize() const { return layout_m == nullptr; } + SIndex& operator=(const SIndex&) = default; + + SIndex& operator=(const SOffset& offset) { + clear(); + offset_m = SOffset(); + addIndex(pointFromOffset(offset)); + return *this; + } + + SIndex& operator=(const NDIndex& domain) { + clear(); + addIndex(domain); + return *this; + } + + SIndex& operator&=(const SIndex& rhs) { + eraseIf([&](const point_type& point) { return !rhs.hasIndex(point); }); + domain_m = domain_m.intersect(rhs.domain_m); + return *this; + } + + SIndex& operator&=(const SOffset& offset) { + const point_type point = pointFromOffset(offset); + eraseIf([&](const point_type& existing) { return !samePoint(existing, point); }); + return *this; + } + + SIndex& operator&=(const NDIndex& domain) { + eraseIf([&](const point_type& point) { return !contains(domain, effectivePoint(point)); }); + domain_m = domain_m.intersect(domain); + return *this; + } + + SIndex& operator|=(const SIndex& rhs) { + domain_m = boundingUnion(domain_m, rhs.domain_m); + for (const auto& point : rhs.points_m) { + addIndex(point); + } + return *this; + } + + SIndex& operator|=(const SOffset& offset) { + domain_m = boundingUnion(domain_m, singlePointDomain(pointFromOffset(offset))); + addIndex(pointFromOffset(offset)); + return *this; + } + + SIndex& operator|=(const NDIndex& domain) { + domain_m = boundingUnion(domain_m, domain); + addIndex(domain); + return *this; + } + bool addIndex(const point_type& point) { if (!contains(domain_m, point) || hasIndex(point)) { return false; @@ -43,14 +107,33 @@ namespace ippl { } bool addIndex(const NDIndex& point) { + if (point.empty()) { + return false; + } + + bool addedAny = false; point_type values; + addIndex(point, values, 0, addedAny); + return addedAny; + } + + bool removeIndex(const point_type& point) { + const auto originalSize = points_m.size(); + eraseIf([&](const point_type& existing) { return samePoint(existing, point); }); + return points_m.size() != originalSize; + } + + bool removeIndex(const SOffset& offset) { + return removeIndex(pointFromOffset(offset)); + } + + bool removeIndex(const NDIndex& point) { for (unsigned d = 0; d < Dim; ++d) { if (point[d].length() != 1) { return false; } - values[d] = point[d].first(); } - return addIndex(values); + return removeIndex(pointFromNDIndex(point)); } void clear() { points_m.clear(); } @@ -63,8 +146,14 @@ namespace ippl { const NDIndex& getDomain() const { return domain_m; } + void setDomain(const NDIndex& domain) { domain_m = domain; } + FieldLayout& getFieldLayout() const { return *layout_m; } + SOffset& getOffset() { return offset_m; } + + const SOffset& getOffset() const { return offset_m; } + bool hasIndex(const point_type& point) const { for (const auto& existing : points_m) { bool same = true; @@ -79,14 +168,44 @@ namespace ippl { } bool hasIndex(const NDIndex& point) const { - point_type values; for (unsigned d = 0; d < Dim; ++d) { if (point[d].length() != 1) { return false; } - values[d] = point[d].first(); } - return hasIndex(values); + return hasIndex(pointFromNDIndex(point)); + } + + SIndex operator()(const SOffset& offset) const { return SIndex(*this, offset); } + + SIndex operator()(int i0) const { + static_assert(Dim == 1); + return operator()(SOffset(i0)); + } + + SIndex operator()(int i0, int i1) const { + static_assert(Dim == 2); + return operator()(SOffset(i0, i1)); + } + + SIndex operator()(int i0, int i1, int i2) const { + static_assert(Dim == 3); + return operator()(SOffset(i0, i1, i2)); + } + + SIndex operator()(int i0, int i1, int i2, int i3) const { + static_assert(Dim == 4); + return operator()(SOffset(i0, i1, i2, i3)); + } + + SIndex operator()(int i0, int i1, int i2, int i3, int i4) const { + static_assert(Dim == 5); + return operator()(SOffset(i0, i1, i2, i3, i4)); + } + + SIndex operator()(int i0, int i1, int i2, int i3, int i4, int i5) const { + static_assert(Dim == 6); + return operator()(SOffset(i0, i1, i2, i3, i4, i5)); } template @@ -94,16 +213,29 @@ namespace ippl { Kokkos::View devicePoints("SIndex::points", points_m.size()); auto hostPoints = Kokkos::create_mirror_view(devicePoints); for (std::size_t i = 0; i < points_m.size(); ++i) { - hostPoints(i) = points_m[i]; + hostPoints(i) = effectivePoint(points_m[i]); } Kokkos::deep_copy(devicePoints, hostPoints); return devicePoints; } + friend SIndex operator+(const SIndex& sindex, const SOffset& offset) { + return SIndex(sindex, offset); + } + + friend SIndex operator+(const SOffset& offset, const SIndex& sindex) { + return SIndex(sindex, offset); + } + + friend SIndex operator-(const SIndex& sindex, const SOffset& offset) { + return SIndex(sindex, -offset); + } + private: FieldLayout* layout_m = nullptr; NDIndex domain_m; std::vector points_m; + SOffset offset_m; static bool contains(const NDIndex& domain, const point_type& point) { for (unsigned d = 0; d < Dim; ++d) { @@ -113,6 +245,81 @@ namespace ippl { } return true; } + + static bool samePoint(const point_type& lhs, const point_type& rhs) { + for (unsigned d = 0; d < Dim; ++d) { + if (lhs[d] != rhs[d]) { + return false; + } + } + return true; + } + + static point_type pointFromOffset(const SOffset& offset) { + point_type point; + for (unsigned d = 0; d < Dim; ++d) { + point[d] = offset[d]; + } + return point; + } + + static point_type pointFromNDIndex(const NDIndex& domain) { + point_type point; + for (unsigned d = 0; d < Dim; ++d) { + point[d] = domain[d].first(); + } + return point; + } + + point_type effectivePoint(const point_type& point) const { + point_type shifted = point; + for (unsigned d = 0; d < Dim; ++d) { + shifted[d] += offset_m[d]; + } + return shifted; + } + + template + void eraseIf(const Predicate& predicate) { + auto out = points_m.begin(); + for (auto in = points_m.begin(); in != points_m.end(); ++in) { + if (!predicate(*in)) { + *out = *in; + ++out; + } + } + points_m.erase(out, points_m.end()); + } + + void addIndex(const NDIndex& domain, point_type& values, unsigned d, bool& addedAny) { + if (d == Dim) { + addedAny = addIndex(values) || addedAny; + return; + } + + for (int value = domain[d].first(); value <= domain[d].last(); + value += domain[d].stride()) { + values[d] = value; + addIndex(domain, values, d + 1, addedAny); + } + } + + static NDIndex boundingUnion(const NDIndex& lhs, const NDIndex& rhs) { + NDIndex result; + for (unsigned d = 0; d < Dim; ++d) { + result[d] = Index(std::min(lhs[d].first(), rhs[d].first()), + std::max(lhs[d].last(), rhs[d].last()), lhs[d].stride()); + } + return result; + } + + static NDIndex singlePointDomain(const point_type& point) { + NDIndex result; + for (unsigned d = 0; d < Dim; ++d) { + result[d] = Index(point[d], point[d]); + } + return result; + } }; } // namespace ippl diff --git a/src/Index/SOffset.h b/src/Index/SOffset.h new file mode 100644 index 000000000..2a871ff65 --- /dev/null +++ b/src/Index/SOffset.h @@ -0,0 +1,133 @@ +// +// Class SOffset +// Integer offset for sparse field index points. +// +#ifndef IPPL_SOFFSET_H +#define IPPL_SOFFSET_H + +#include +#include + +#include "Index/NDIndex.h" +#include "Types/Vector.h" + +namespace ippl { + + template + class SOffset { + public: + using point_type = Vector; + + SOffset() = default; + + explicit SOffset(int value) + : values_m(value) {} + + template ::type = true> + explicit SOffset(Args... args) + : values_m(static_cast(args)...) {} + + explicit SOffset(const point_type& values) + : values_m(values) {} + + int& operator[](unsigned d) { return values_m[d]; } + + int operator[](unsigned d) const { return values_m[d]; } + + const point_type& values() const { return values_m; } + + SOffset& operator+=(const SOffset& rhs) { + for (unsigned d = 0; d < Dim; ++d) { + values_m[d] += rhs[d]; + } + return *this; + } + + SOffset& operator-=(const SOffset& rhs) { + for (unsigned d = 0; d < Dim; ++d) { + values_m[d] -= rhs[d]; + } + return *this; + } + + bool operator==(const SOffset& rhs) const { + for (unsigned d = 0; d < Dim; ++d) { + if (values_m[d] != rhs[d]) { + return false; + } + } + return true; + } + + bool operator!=(const SOffset& rhs) const { return !(*this == rhs); } + + bool inside(const NDIndex& domain) const { + for (unsigned d = 0; d < Dim; ++d) { + if (!domain[d].contains(Index(values_m[d], values_m[d]))) { + return false; + } + } + return true; + } + + private: + point_type values_m; + }; + + template + SOffset operator+(SOffset lhs, const SOffset& rhs) { + lhs += rhs; + return lhs; + } + + template + SOffset operator-(SOffset lhs, const SOffset& rhs) { + lhs -= rhs; + return lhs; + } + + template + SOffset operator-(const SOffset& value) { + SOffset result; + for (unsigned d = 0; d < Dim; ++d) { + result[d] = -value[d]; + } + return result; + } + + template + NDIndex operator+(NDIndex domain, const SOffset& offset) { + for (unsigned d = 0; d < Dim; ++d) { + domain[d] = Index(domain[d].first() + offset[d], domain[d].last() + offset[d], + domain[d].stride()); + } + return domain; + } + + template + NDIndex operator+(const SOffset& offset, const NDIndex& domain) { + return domain + offset; + } + + template + NDIndex operator-(const NDIndex& domain, const SOffset& offset) { + return domain + (-offset); + } + + template + std::ostream& operator<<(std::ostream& out, const SOffset& offset) { + out << "["; + for (unsigned d = 0; d < Dim; ++d) { + out << offset[d]; + if (d + 1 < Dim) { + out << ","; + } + } + out << "]"; + return out; + } + +} // namespace ippl + +#endif diff --git a/src/IpplCore.h b/src/IpplCore.h index 71f1e8b8b..51fffe3e2 100644 --- a/src/IpplCore.h +++ b/src/IpplCore.h @@ -18,6 +18,7 @@ // #include "Index/Index.h" // #include "Index/NDIndex.h" +#include "Index/SOffset.h" #include "Index/SIndex.h" #include "FieldLayout/FieldLayout.h" diff --git a/unit_tests/BareField/BareField.cpp b/unit_tests/BareField/BareField.cpp index 2e3725867..450920a6d 100644 --- a/unit_tests/BareField/BareField.cpp +++ b/unit_tests/BareField/BareField.cpp @@ -345,6 +345,131 @@ TYPED_TEST(BareFieldTest, SparseIndexedSubdomainExpressionAssignment) { }); } +TYPED_TEST(BareFieldTest, SparseIndexedOffsetExpressionAssignment) { + using T = typename TestFixture::value_type; + + auto& lhs = this->field; + lhs->operator=(0); + + typename TestFixture::field_type rhs(this->layout); + rhs = T(5); + rhs.fillHalo(); + + ippl::SIndex sindex(this->layout); + const auto domain = lhs->getDomain(); + const int offsetValue = domain[0].stride(); + + for (int i = domain[0].first(); i <= domain[0].last() - offsetValue; + i += 2 * domain[0].stride()) { + typename ippl::SIndex::point_type point; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + point[d] = domain[d].first(); + } + point[0] = i; + sindex.addIndex(point); + } + + ippl::SOffset offset; + offset[0] = offsetValue; + + (*lhs)[sindex] = rhs[sindex + offset] + T(2); + + auto mirror = lhs->getHostMirror(); + Kokkos::deep_copy(mirror, lhs->getView()); + + const auto owned = lhs->getOwned(); + const int nghost = lhs->getNghost(); + nestedViewLoop(mirror, nghost, [&](const Idx... args) { + typename ippl::SIndex::point_type point; + const std::array viewCoords{static_cast(args)...}; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + point[d] = owned[d].first() + + static_cast(viewCoords[d] - nghost) * owned[d].stride(); + } + const T expected = sindex.hasIndex(point) ? T(7) : T(0); + assertEqual(expected, mirror(args...)); + }); +} + +TYPED_TEST(BareFieldTest, SparseIndexedCallOffsetExpressionAssignment) { + using T = typename TestFixture::value_type; + + auto& lhs = this->field; + lhs->operator=(0); + + typename TestFixture::field_type rhs(this->layout); + rhs = T(5); + rhs.fillHalo(); + + ippl::SIndex sindex(this->layout); + const auto domain = lhs->getDomain(); + const int offsetValue = domain[0].stride(); + + for (int i = domain[0].first(); i <= domain[0].last() - offsetValue; + i += 2 * domain[0].stride()) { + typename ippl::SIndex::point_type point; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + point[d] = domain[d].first(); + } + point[0] = i; + sindex.addIndex(point); + } + + ippl::SOffset offset; + offset[0] = offsetValue; + + (*lhs)[sindex] = rhs[sindex(offset)] + T(2); + + auto mirror = lhs->getHostMirror(); + Kokkos::deep_copy(mirror, lhs->getView()); + + const auto owned = lhs->getOwned(); + const int nghost = lhs->getNghost(); + nestedViewLoop(mirror, nghost, [&](const Idx... args) { + typename ippl::SIndex::point_type point; + const std::array viewCoords{static_cast(args)...}; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + point[d] = owned[d].first() + + static_cast(viewCoords[d] - nghost) * owned[d].stride(); + } + const T expected = sindex.hasIndex(point) ? T(7) : T(0); + assertEqual(expected, mirror(args...)); + }); +} + +TYPED_TEST(BareFieldTest, SparseIndexSetOperations) { + ippl::SIndex sindex(this->layout); + const auto domain = this->field->getDomain(); + + typename ippl::SIndex::point_type firstPoint; + typename ippl::SIndex::point_type secondPoint; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + firstPoint[d] = domain[d].first(); + secondPoint[d] = domain[d].first(); + } + secondPoint[0] = domain[0].first() + domain[0].stride(); + + EXPECT_TRUE(sindex.addIndex(firstPoint)); + EXPECT_TRUE(sindex.addIndex(secondPoint)); + EXPECT_FALSE(sindex.addIndex(firstPoint)); + + ippl::NDIndex subset = domain; + subset[0] = ippl::Index(firstPoint[0], firstPoint[0], domain[0].stride()); + sindex &= subset; + + EXPECT_TRUE(sindex.hasIndex(firstPoint)); + EXPECT_FALSE(sindex.hasIndex(secondPoint)); + + ippl::SOffset secondOffset; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + secondOffset[d] = secondPoint[d]; + } + sindex |= secondOffset; + + EXPECT_TRUE(sindex.hasIndex(secondPoint)); + EXPECT_EQ(sindex.size(), 2u); +} + TYPED_TEST(BareFieldTest, Sum) { using T = typename TestFixture::value_type; From ee58cf0b86cb484a8c7ccc9130584eee760a525e Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Wed, 13 May 2026 00:44:58 +0200 Subject: [PATCH 03/12] Build sparse indexes from predicates --- src/Index/SIndex.h | 139 +++++++++++++++++++++++++++++ unit_tests/BareField/BareField.cpp | 74 +++++++++++++++ 2 files changed, 213 insertions(+) diff --git a/src/Index/SIndex.h b/src/Index/SIndex.h index 165952d48..6e495e4ff 100644 --- a/src/Index/SIndex.h +++ b/src/Index/SIndex.h @@ -8,16 +8,22 @@ #include #include +#include #include #include #include "FieldLayout/FieldLayout.h" +#include "Expression/IpplExpressions.h" #include "Index/NDIndex.h" #include "Index/SOffset.h" #include "Types/Vector.h" +#include "Utility/ParallelDispatch.h" namespace ippl { + template + class IndexedSIndex; + template class SIndex { public: @@ -60,6 +66,12 @@ namespace ippl { return *this; } + template + SIndex& operator=(const detail::Expression& expr) { + assignExpression(expr, domain_m, EvaluationMode::ViewCoordinates); + return *this; + } + SIndex& operator&=(const SIndex& rhs) { eraseIf([&](const point_type& point) { return !rhs.hasIndex(point); }); domain_m = domain_m.intersect(rhs.domain_m); @@ -178,6 +190,16 @@ namespace ippl { SIndex operator()(const SOffset& offset) const { return SIndex(*this, offset); } + IndexedSIndex operator[](const NDIndex& domain) { + return IndexedSIndex(*this, domain); + } + + IndexedSIndex operator[](const Index& index) { + NDIndex domain = domain_m; + domain[0] = index; + return operator[](domain); + } + SIndex operator()(int i0) const { static_assert(Dim == 1); return operator()(SOffset(i0)); @@ -237,6 +259,8 @@ namespace ippl { std::vector points_m; SOffset offset_m; + enum class EvaluationMode { ViewCoordinates, RelativeCoordinates }; + static bool contains(const NDIndex& domain, const point_type& point) { for (unsigned d = 0; d < Dim; ++d) { if (!domain[d].contains(Index(point[d], point[d]))) { @@ -320,8 +344,123 @@ namespace ippl { } return result; } + + template + KOKKOS_INLINE_FUNCTION static typename RangePolicy::index_type flatIndex( + const Coords& coords, const Kokkos::Array< + typename RangePolicy:: + index_type, + Dim>& extents) { + typename RangePolicy::index_type flat = 0; + for (unsigned d = 0; d < Dim; ++d) { + flat = flat * extents[d] + coords[d]; + } + return flat; + } + + template + void assignExpression(const detail::Expression& expr, const NDIndex& domain, + EvaluationMode mode) { + clear(); + setDomain(domain); + + NDIndex local = domain.intersect(getFieldLayout().getLocalNDIndex()); + if (local.empty()) { + return; + } + + using execution_space = Kokkos::DefaultExecutionSpace; + using range_policy_type = RangePolicy; + using index_type = typename range_policy_type::index_type; + using capture_type = detail::CapturedExpression; + + Kokkos::Array begin, end, extents; + index_type size = 1; + for (unsigned d = 0; d < Dim; ++d) { + begin[d] = 0; + end[d] = local[d].length(); + extents[d] = local[d].length(); + size *= extents[d]; + } + + Kokkos::View selected("SIndex::selected", size); + capture_type expr_ = reinterpret_cast(expr); + const int nghost = 1; + auto localCopy = local; + auto domainCopy = domain; + + ippl::parallel_for( + "SIndex::operator=", createRangePolicy(begin, end), + KOKKOS_LAMBDA(const typename range_policy_type::index_array_type& args) { + typename range_policy_type::index_array_type exprCoords; + for (unsigned d = 0; d < Dim; ++d) { + const auto global = localCopy[d].first() + args[d] * localCopy[d].stride(); + if (mode == EvaluationMode::ViewCoordinates) { + exprCoords[d] = + (global - localCopy[d].first()) / localCopy[d].stride() + nghost; + } else { + exprCoords[d] = + (global - domainCopy[d].first()) / domainCopy[d].stride(); + } + } + selected(flatIndex(args, extents)) = apply(expr_, exprCoords) ? 1 : 0; + }); + + auto hostSelected = Kokkos::create_mirror_view(selected); + Kokkos::deep_copy(hostSelected, selected); + + point_type global; + addSelected(local, extents, hostSelected, global, 0, 0); + } + + template + void addSelected(const NDIndex& local, const Kokkos::Array< + typename RangePolicy< + Dim, Kokkos::DefaultExecutionSpace>:: + index_type, + Dim>& extents, + const HostView& selected, point_type& global, unsigned d, + typename RangePolicy::index_type flat) { + if (d == Dim) { + if (selected(flat)) { + addIndex(global); + } + return; + } + + for (typename RangePolicy::index_type i = 0; + i < extents[d]; ++i) { + global[d] = local[d].first() + static_cast(i) * local[d].stride(); + addSelected(local, extents, selected, global, d + 1, flat * extents[d] + i); + } + } + + friend class IndexedSIndex; }; + template + class IndexedSIndex { + public: + IndexedSIndex(SIndex& sindex, const NDIndex& domain) + : sindex_m(sindex) + , domain_m(domain) {} + + template + IndexedSIndex& operator=(const detail::Expression& expr) { + sindex_m.assignExpression(expr, domain_m, SIndex::EvaluationMode::RelativeCoordinates); + return *this; + } + + private: + SIndex& sindex_m; + NDIndex domain_m; + }; + + template >> + auto gt(const detail::Expression& expr, const T& value) { + return expr > value; + } + } // namespace ippl #endif diff --git a/unit_tests/BareField/BareField.cpp b/unit_tests/BareField/BareField.cpp index 450920a6d..ba1974808 100644 --- a/unit_tests/BareField/BareField.cpp +++ b/unit_tests/BareField/BareField.cpp @@ -470,6 +470,80 @@ TYPED_TEST(BareFieldTest, SparseIndexSetOperations) { EXPECT_EQ(sindex.size(), 2u); } +TYPED_TEST(BareFieldTest, SparseIndexPredicateAssignment) { + using T = typename TestFixture::value_type; + + auto& source = this->field; + auto view = source->getView(); + const auto lDom = source->getLayout().getLocalNDIndex(); + + Kokkos::parallel_for("Set field", source->getFieldRangePolicy(), + FieldVal{view, lDom}); + Kokkos::fence(); + + typename TestFixture::field_type selected(this->layout); + selected = T(0); + + ippl::SIndex sindex(this->layout); + sindex = *source > T(2); + selected[sindex] = (*source)[sindex]; + + auto sourceMirror = source->getHostMirror(); + auto selectedMirror = selected.getHostMirror(); + Kokkos::deep_copy(sourceMirror, source->getView()); + Kokkos::deep_copy(selectedMirror, selected.getView()); + + const int nghost = source->getNghost(); + nestedViewLoop(sourceMirror, nghost, [&](const Idx... args) { + const T value = sourceMirror(args...); + const T expected = value > T(2) ? value : T(0); + assertEqual(expected, selectedMirror(args...)); + }); +} + +TYPED_TEST(BareFieldTest, SparseIndexSubsetPredicateAssignment) { + using T = typename TestFixture::value_type; + + auto& source = this->field; + auto view = source->getView(); + const auto lDom = source->getLayout().getLocalNDIndex(); + + Kokkos::parallel_for("Set field", source->getFieldRangePolicy(), + FieldVal{view, lDom}); + Kokkos::fence(); + + typename TestFixture::field_type selected(this->layout); + selected = T(0); + + ippl::NDIndex subset = source->getDomain(); + subset[0] = ippl::Index(subset[0].first(), subset[0].last(), 2 * subset[0].stride()); + + ippl::SIndex sindex(this->layout); + sindex[subset] = (*source)[subset] > T(2); + selected[sindex] = (*source)[sindex]; + + auto sourceMirror = source->getHostMirror(); + auto selectedMirror = selected.getHostMirror(); + Kokkos::deep_copy(sourceMirror, source->getView()); + Kokkos::deep_copy(selectedMirror, selected.getView()); + + const auto owned = source->getOwned(); + const int nghost = source->getNghost(); + nestedViewLoop(sourceMirror, nghost, [&](const Idx... args) { + std::array pointIndices; + const std::array viewCoords{static_cast(args)...}; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + const int global = owned[d].first() + + static_cast(viewCoords[d] - nghost) * owned[d].stride(); + pointIndices[d] = ippl::Index(global, global); + } + auto point = std::make_from_tuple>(pointIndices); + const T value = sourceMirror(args...); + const T expected = subset.contains(point) && value > T(2) ? value : T(0); + assertEqual(expected, selectedMirror(args...)); + }); +} + TYPED_TEST(BareFieldTest, Sum) { using T = typename TestFixture::value_type; From b060c7f35cf9825266b0f7fada85c95a80e61bb6 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Wed, 13 May 2026 00:48:58 +0200 Subject: [PATCH 04/12] Support sparse compound field assignment --- src/Field/SparseIndexedBareField.h | 145 +++++++++++++++++++++++++++++ unit_tests/BareField/BareField.cpp | 32 +++++++ 2 files changed, 177 insertions(+) diff --git a/src/Field/SparseIndexedBareField.h b/src/Field/SparseIndexedBareField.h index 189f25377..264b89cc3 100644 --- a/src/Field/SparseIndexedBareField.h +++ b/src/Field/SparseIndexedBareField.h @@ -61,6 +61,58 @@ namespace ippl { return operator=(expr); } + SparseIndexedBareField& operator+=(value_type value) { + update(PlusValueAssign{value}); + return *this; + } + + SparseIndexedBareField& operator-=(value_type value) { + update(MinusValueAssign{value}); + return *this; + } + + SparseIndexedBareField& operator*=(value_type value) { + update(MultiplyValueAssign{value}); + return *this; + } + + SparseIndexedBareField& operator/=(value_type value) { + update(DivideValueAssign{value}); + return *this; + } + + template + SparseIndexedBareField& operator+=(const detail::Expression& expr) { + using capture_type = detail::CapturedExpression; + capture_type expr_ = reinterpret_cast(expr); + update(PlusExpressionAssign{expr_}); + return *this; + } + + template + SparseIndexedBareField& operator-=(const detail::Expression& expr) { + using capture_type = detail::CapturedExpression; + capture_type expr_ = reinterpret_cast(expr); + update(MinusExpressionAssign{expr_}); + return *this; + } + + template + SparseIndexedBareField& operator*=(const detail::Expression& expr) { + using capture_type = detail::CapturedExpression; + capture_type expr_ = reinterpret_cast(expr); + update(MultiplyExpressionAssign{expr_}); + return *this; + } + + template + SparseIndexedBareField& operator/=(const detail::Expression& expr) { + using capture_type = detail::CapturedExpression; + capture_type expr_ = reinterpret_cast(expr); + update(DivideExpressionAssign{expr_}); + return *this; + } + private: view_type view_m; NDIndex owned_m; @@ -76,6 +128,42 @@ namespace ippl { } }; + struct PlusValueAssign { + value_type value; + + template + KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords&) const { + return current + value; + } + }; + + struct MinusValueAssign { + value_type value; + + template + KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords&) const { + return current - value; + } + }; + + struct MultiplyValueAssign { + value_type value; + + template + KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords&) const { + return current * value; + } + }; + + struct DivideValueAssign { + value_type value; + + template + KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords&) const { + return current / value; + } + }; + template struct ExpressionAssign { detail::CapturedExpression expr; @@ -86,6 +174,46 @@ namespace ippl { } }; + template + struct PlusExpressionAssign { + detail::CapturedExpression expr; + + template + KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords& rel) const { + return current + apply(expr, rel); + } + }; + + template + struct MinusExpressionAssign { + detail::CapturedExpression expr; + + template + KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords& rel) const { + return current - apply(expr, rel); + } + }; + + template + struct MultiplyExpressionAssign { + detail::CapturedExpression expr; + + template + KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords& rel) const { + return current * apply(expr, rel); + } + }; + + template + struct DivideExpressionAssign { + detail::CapturedExpression expr; + + template + KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords& rel) const { + return current / apply(expr, rel); + } + }; + KOKKOS_INLINE_FUNCTION bool isLocal(const point_type& point) const { for (unsigned d = 0; d < field_dim; ++d) { if (!owned_m[d].contains(Index(point[d], point[d]))) { @@ -118,6 +246,23 @@ namespace ippl { } }); } + + template + void update(const Functor& functor) { + ippl::parallel_for( + "SparseIndexedBareField::operator op=", getRangePolicy(points_m), + KOKKOS_CLASS_LAMBDA( + const typename RangePolicy<1, execution_space>::index_array_type& args) { + const auto i = args[0]; + const point_type point = points_m(i); + if (isLocal(point)) { + typename RangePolicy<1, execution_space>::index_array_type rel; + rel[0] = i; + const auto viewCoords = pointToView(point); + apply(view_m, viewCoords) = functor(apply(view_m, viewCoords), rel); + } + }); + } }; } // namespace ippl diff --git a/unit_tests/BareField/BareField.cpp b/unit_tests/BareField/BareField.cpp index ba1974808..ba6b25b93 100644 --- a/unit_tests/BareField/BareField.cpp +++ b/unit_tests/BareField/BareField.cpp @@ -544,6 +544,38 @@ TYPED_TEST(BareFieldTest, SparseIndexSubsetPredicateAssignment) { }); } +TYPED_TEST(BareFieldTest, SparseIndexedCompoundAssignment) { + using T = typename TestFixture::value_type; + + auto& source = this->field; + auto view = source->getView(); + const auto lDom = source->getLayout().getLocalNDIndex(); + + Kokkos::parallel_for("Set field", source->getFieldRangePolicy(), + FieldVal{view, lDom}); + Kokkos::fence(); + + typename TestFixture::field_type values(this->layout); + values = T(1); + + ippl::SIndex sindex(this->layout); + sindex = *source > T(2); + values[sindex] = (*source)[sindex]; + values[sindex] *= T(1.01); + + auto sourceMirror = source->getHostMirror(); + auto valuesMirror = values.getHostMirror(); + Kokkos::deep_copy(sourceMirror, source->getView()); + Kokkos::deep_copy(valuesMirror, values.getView()); + + const int nghost = source->getNghost(); + nestedViewLoop(sourceMirror, nghost, [&](const Idx... args) { + const T value = sourceMirror(args...); + const T expected = value > T(2) ? value * T(1.01) : T(1); + assertEqual(expected, valuesMirror(args...)); + }); +} + TYPED_TEST(BareFieldTest, Sum) { using T = typename TestFixture::value_type; From d4e27ca8aab61fdfd511d4e5e20ba172e1f0641e Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Wed, 13 May 2026 07:27:02 +0200 Subject: [PATCH 05/12] Fix CUDA extended lambda access --- src/Field/IndexedBareField.h | 1 + src/Field/SparseIndexedBareField.h | 1 + src/Index/SIndex.h | 1 + 3 files changed, 3 insertions(+) diff --git a/src/Field/IndexedBareField.h b/src/Field/IndexedBareField.h index 195c9a1d2..0786fc99d 100644 --- a/src/Field/IndexedBareField.h +++ b/src/Field/IndexedBareField.h @@ -137,6 +137,7 @@ namespace ippl { return viewCoords; } + public: template void assign(const Functor& functor) { Domain_t local = domain_m.intersect(owned_m); diff --git a/src/Field/SparseIndexedBareField.h b/src/Field/SparseIndexedBareField.h index 264b89cc3..16104734b 100644 --- a/src/Field/SparseIndexedBareField.h +++ b/src/Field/SparseIndexedBareField.h @@ -231,6 +231,7 @@ namespace ippl { return viewCoords; } + public: template void assign(const Functor& functor) { ippl::parallel_for( diff --git a/src/Index/SIndex.h b/src/Index/SIndex.h index 6e495e4ff..3e53bee4b 100644 --- a/src/Index/SIndex.h +++ b/src/Index/SIndex.h @@ -358,6 +358,7 @@ namespace ippl { return flat; } + public: template void assignExpression(const detail::Expression& expr, const NDIndex& domain, EvaluationMode mode) { From c7e612045fabd655b8c75ab9f0a5a5621c9c181d Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Wed, 13 May 2026 07:44:53 +0200 Subject: [PATCH 06/12] Expose CUDA lambda argument types --- src/Field/IndexedBareField.h | 2 ++ src/Field/SparseIndexedBareField.h | 2 ++ src/Index/SIndex.h | 2 ++ 3 files changed, 6 insertions(+) diff --git a/src/Field/IndexedBareField.h b/src/Field/IndexedBareField.h index 0786fc99d..67e74abe5 100644 --- a/src/Field/IndexedBareField.h +++ b/src/Field/IndexedBareField.h @@ -90,6 +90,7 @@ namespace ippl { Domain_t domain_m; int nghost_m; + public: struct ValueAssign { value_type value; @@ -109,6 +110,7 @@ namespace ippl { } }; + private: template KOKKOS_INLINE_FUNCTION auto relativeToView(const Coords& rel) const { typename RangePolicy::index_array_type viewCoords; diff --git a/src/Field/SparseIndexedBareField.h b/src/Field/SparseIndexedBareField.h index 16104734b..5f102e3c7 100644 --- a/src/Field/SparseIndexedBareField.h +++ b/src/Field/SparseIndexedBareField.h @@ -119,6 +119,7 @@ namespace ippl { points_view_type points_m; int nghost_m; + public: struct ValueAssign { value_type value; @@ -214,6 +215,7 @@ namespace ippl { } }; + private: KOKKOS_INLINE_FUNCTION bool isLocal(const point_type& point) const { for (unsigned d = 0; d < field_dim; ++d) { if (!owned_m[d].contains(Index(point[d], point[d]))) { diff --git a/src/Index/SIndex.h b/src/Index/SIndex.h index 3e53bee4b..449900484 100644 --- a/src/Index/SIndex.h +++ b/src/Index/SIndex.h @@ -259,8 +259,10 @@ namespace ippl { std::vector points_m; SOffset offset_m; + public: enum class EvaluationMode { ViewCoordinates, RelativeCoordinates }; + private: static bool contains(const NDIndex& domain, const point_type& point) { for (unsigned d = 0; d < Dim; ++d) { if (!domain[d].contains(Index(point[d], point[d]))) { From 1149dd70ebcedca55f94fd5d70db893e59c8bbe8 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Wed, 13 May 2026 08:18:55 +0200 Subject: [PATCH 07/12] Evaluate sparse predicates in expression space --- src/Expression/IpplExpressions.h | 18 ++++++++++++++++++ src/Expression/IpplOperations.h | 2 ++ src/Index/SIndex.h | 2 +- 3 files changed, 21 insertions(+), 1 deletion(-) diff --git a/src/Expression/IpplExpressions.h b/src/Expression/IpplExpressions.h index 632238eb2..1e76b5c4e 100644 --- a/src/Expression/IpplExpressions.h +++ b/src/Expression/IpplExpressions.h @@ -5,6 +5,8 @@ #ifndef IPPL_EXPRESSIONS_H #define IPPL_EXPRESSIONS_H +#include + #include namespace ippl { @@ -94,6 +96,22 @@ namespace ippl { template struct isExpression> : std::true_type {}; + template + struct ExecutionSpaceOf { + using type = Kokkos::DefaultExecutionSpace; + }; + + template + struct ExecutionSpaceOf> { + using type = typename E::execution_space; + }; + + template + struct BinaryExecutionSpace { + using type = std::conditional_t<(E1::dim != 0), typename ExecutionSpaceOf::type, + typename ExecutionSpaceOf::type>; + }; + } // namespace detail } // namespace ippl diff --git a/src/Expression/IpplOperations.h b/src/Expression/IpplOperations.h index 1e02a6b8a..a63d07887 100644 --- a/src/Expression/IpplOperations.h +++ b/src/Expression/IpplOperations.h @@ -71,6 +71,7 @@ namespace ippl { struct fun : public detail::Expression, sizeof(E)> { \ constexpr static unsigned dim = E::dim; \ using value_type = typename E::value_type; \ + using execution_space = typename detail::ExecutionSpaceOf::type; \ \ KOKKOS_FUNCTION \ fun(const E& u) \ @@ -132,6 +133,7 @@ namespace ippl { struct fun : public detail::Expression, sizeof(E1) + sizeof(E2)> { \ constexpr static unsigned dim = std::max(E1::dim, E2::dim); \ using value_type = typename E1::value_type; \ + using execution_space = typename detail::BinaryExecutionSpace::type; \ \ KOKKOS_FUNCTION \ fun(const E1& u, const E2& v) \ diff --git a/src/Index/SIndex.h b/src/Index/SIndex.h index 449900484..8ae57df24 100644 --- a/src/Index/SIndex.h +++ b/src/Index/SIndex.h @@ -372,7 +372,7 @@ namespace ippl { return; } - using execution_space = Kokkos::DefaultExecutionSpace; + using execution_space = typename detail::ExecutionSpaceOf::type; using range_policy_type = RangePolicy; using index_type = typename range_policy_type::index_type; using capture_type = detail::CapturedExpression; From 266b39e13434d17c23d56b744a5a78e70b78d9a0 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Wed, 13 May 2026 08:29:29 +0200 Subject: [PATCH 08/12] Use predicate index type in SIndex helpers --- src/Index/SIndex.h | 27 +++++++++------------------ 1 file changed, 9 insertions(+), 18 deletions(-) diff --git a/src/Index/SIndex.h b/src/Index/SIndex.h index 8ae57df24..2f8fa50ac 100644 --- a/src/Index/SIndex.h +++ b/src/Index/SIndex.h @@ -347,13 +347,10 @@ namespace ippl { return result; } - template - KOKKOS_INLINE_FUNCTION static typename RangePolicy::index_type flatIndex( - const Coords& coords, const Kokkos::Array< - typename RangePolicy:: - index_type, - Dim>& extents) { - typename RangePolicy::index_type flat = 0; + template + KOKKOS_INLINE_FUNCTION static IndexType flatIndex( + const Coords& coords, const Kokkos::Array& extents) { + IndexType flat = 0; for (unsigned d = 0; d < Dim; ++d) { flat = flat * extents[d] + coords[d]; } @@ -413,17 +410,12 @@ namespace ippl { Kokkos::deep_copy(hostSelected, selected); point_type global; - addSelected(local, extents, hostSelected, global, 0, 0); + addSelected(local, extents, hostSelected, global, 0, index_type{0}); } - template - void addSelected(const NDIndex& local, const Kokkos::Array< - typename RangePolicy< - Dim, Kokkos::DefaultExecutionSpace>:: - index_type, - Dim>& extents, - const HostView& selected, point_type& global, unsigned d, - typename RangePolicy::index_type flat) { + template + void addSelected(const NDIndex& local, const Kokkos::Array& extents, + const HostView& selected, point_type& global, unsigned d, IndexType flat) { if (d == Dim) { if (selected(flat)) { addIndex(global); @@ -431,8 +423,7 @@ namespace ippl { return; } - for (typename RangePolicy::index_type i = 0; - i < extents[d]; ++i) { + for (IndexType i = 0; i < extents[d]; ++i) { global[d] = local[d].first() + static_cast(i) * local[d].stride(); addSelected(local, extents, selected, global, d + 1, flat * extents[d] + i); } From 31e338f8ad5d819e2180f559c5143e90fb3f1545 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Sun, 17 May 2026 22:17:51 +0200 Subject: [PATCH 09/12] Speed up predicate-built sparse indexes Avoid routing SIndex predicate assignment through addIndex(), which performs a linear duplicate scan for every selected point. Predicate assignment clears the index and visits each local candidate point once, so selected points can be appended directly after reserving capacity. This removes the quadratic host-side construction cost that made the high-rank SparseIndexSubsetPredicateAssignment tests look like a hang, especially under MPI/GPU test runs. --- src/Index/SIndex.h | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/src/Index/SIndex.h b/src/Index/SIndex.h index 2f8fa50ac..5abbefb05 100644 --- a/src/Index/SIndex.h +++ b/src/Index/SIndex.h @@ -8,6 +8,7 @@ #include #include +#include #include #include #include @@ -409,23 +410,30 @@ namespace ippl { auto hostSelected = Kokkos::create_mirror_view(selected); Kokkos::deep_copy(hostSelected, selected); + std::size_t selectedCount = 0; + for (index_type flat = 0; flat < size; ++flat) { + selectedCount += hostSelected(flat) != 0; + } + points_m.reserve(selectedCount); + point_type global; - addSelected(local, extents, hostSelected, global, 0, index_type{0}); + appendSelected(local, extents, hostSelected, global, 0, index_type{0}); } template - void addSelected(const NDIndex& local, const Kokkos::Array& extents, - const HostView& selected, point_type& global, unsigned d, IndexType flat) { + void appendSelected(const NDIndex& local, const Kokkos::Array& extents, + const HostView& selected, point_type& global, unsigned d, + IndexType flat) { if (d == Dim) { if (selected(flat)) { - addIndex(global); + points_m.push_back(global); } return; } for (IndexType i = 0; i < extents[d]; ++i) { global[d] = local[d].first() + static_cast(i) * local[d].stride(); - addSelected(local, extents, selected, global, d + 1, flat * extents[d] + i); + appendSelected(local, extents, selected, global, d + 1, flat * extents[d] + i); } } From 7ed04e555a88ce62da01f64e885c485dedd6a7c1 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Mon, 18 May 2026 06:09:44 +0200 Subject: [PATCH 10/12] Align captured expression storage for device kernels Keep the existing reinterpret-cast capture model, but align CapturedExpression and its byte buffer to the captured expression type. This avoids reinterpreting byte-aligned storage as expression objects that may contain stricter-aligned members such as Kokkos views. Also speed up predicate-built SIndex construction by reserving selected capacity and appending selected points directly. Predicate assignment clears the sparse index and visits each candidate point once, so the public addIndex duplicate scan is unnecessary there and caused quadratic host-side cost in high-rank sparse-index tests. --- src/Expression/IpplExpressions.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Expression/IpplExpressions.h b/src/Expression/IpplExpressions.h index 1e76b5c4e..fff85df4a 100644 --- a/src/Expression/IpplExpressions.h +++ b/src/Expression/IpplExpressions.h @@ -43,7 +43,7 @@ namespace ippl { * of the BareField class. */ template - struct CapturedExpression { + struct alignas(E) CapturedExpression { constexpr static unsigned dim = E::dim; template @@ -52,7 +52,7 @@ namespace ippl { return reinterpret_cast(*this)(args...); } - char buffer[N]; + alignas(E) char buffer[N]; }; /*! From ab30b464c28865b2238fc6274c9957a2824444bf Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Mon, 18 May 2026 14:28:52 +0200 Subject: [PATCH 11/12] Add 3D subfield sum regression test Add a focused Field regression test for power-of-two 3D subdomains. The test initializes a cubic field to one, verifies the full-field sum, then repeatedly assigns smaller N/2, N/4, ..., 2^3 subdomains into a zeroed target field and checks that the summed value matches the expected subdomain volume. The target keeps the same global layout as the source so the test is stable under MPI-backed ctest runs while still exercising the Field subdomain access and assignment path. --- unit_tests/Field/Field.cpp | 75 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) diff --git a/unit_tests/Field/Field.cpp b/unit_tests/Field/Field.cpp index 505832b5d..f34f9c970 100644 --- a/unit_tests/Field/Field.cpp +++ b/unit_tests/Field/Field.cpp @@ -14,6 +14,9 @@ template class FieldTest; +template +class Field3DRegressionTest; + template class FieldTest>> : public ::testing::Test { public: @@ -63,6 +66,48 @@ class FieldTest>> : public ::testing::Test { std::array domain; }; +template +class Field3DRegressionTest>> : public ::testing::Test { +public: + using value_type = T; + using exec_space = ExecSpace; + constexpr static unsigned dim = 3; + constexpr static size_t initial_size = 16; + + using mesh_type = ippl::UniformCartesian; + using centering_type = typename mesh_type::DefaultCentering; + using field_type = ippl::Field; + using layout_type = ippl::FieldLayout; + + struct FieldBundle { + std::shared_ptr layout; + std::shared_ptr mesh; + std::shared_ptr field; + }; + + FieldBundle makeField(size_t n) const { + std::array indices; + indices.fill(ippl::Index(n)); + auto owned = std::make_from_tuple>(indices); + + ippl::Vector hx; + ippl::Vector origin; + for (unsigned d = 0; d < dim; ++d) { + hx[d] = T(1) / static_cast(n); + origin[d] = T(0); + } + + std::array isParallel; + isParallel.fill(true); + + FieldBundle bundle; + bundle.layout = std::make_shared(MPI_COMM_WORLD, owned, isParallel); + bundle.mesh = std::make_shared(owned, hx, origin); + bundle.field = std::make_shared(*bundle.mesh, *bundle.layout); + return bundle; + } +}; + template struct VFieldVal { using vfield_view_type = typename FieldTest::vfield_type::view_type; @@ -151,6 +196,9 @@ struct FieldVal { using Tests = TestParams::tests<1, 2, 3, 4, 5, 6>; TYPED_TEST_SUITE(FieldTest, Tests); +using Field3DRegressionTests = TestParams::tests<3>; +TYPED_TEST_SUITE(Field3DRegressionTest, Field3DRegressionTests); + TYPED_TEST(FieldTest, DeepCopy) { auto& field = this->field; @@ -184,6 +232,33 @@ TYPED_TEST(FieldTest, Sum) { assertEqual(expected, sum); } +TYPED_TEST(Field3DRegressionTest, PowerOfTwoSubfieldSums) { + using T = typename TestFixture::value_type; + constexpr size_t N = TestFixture::initial_size; + + auto source = this->makeField(N); + *source.field = T(1); + + assertEqual(static_cast(N * N * N), source.field->sum()); + + for (size_t n = N / 2; n >= 2; n /= 2) { + auto target = this->makeField(N); + *target.field = T(0); + + std::array indices; + indices.fill(ippl::Index(n)); + const auto subdomain = std::make_from_tuple>(indices); + + (*target.field)[subdomain] = (*source.field)[subdomain]; + + assertEqual(static_cast(n * n * n), target.field->sum()); + + if (n == 2) { + break; + } + } +} + TYPED_TEST(FieldTest, Norm1) { using T = typename TestFixture::value_type; From 0abfa7a1f3336da14ff35b62375bef5051a04107 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Tue, 19 May 2026 21:33:54 +0200 Subject: [PATCH 12/12] Test differently-strided BareField subset assignment Add a BareField regression for assigning between compatible rectangular subsets with different physical strides: A[a] = B[b] + scalar The test constructs local owned subdomains where the left-hand side uses a larger stride than the right-hand side while both subsets have the same logical extent. It verifies that IndexedBareField evaluates each side in its own subdomain coordinates instead of assuming identical strides. --- unit_tests/BareField/BareField.cpp | 54 ++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/unit_tests/BareField/BareField.cpp b/unit_tests/BareField/BareField.cpp index ba6b25b93..4bf8b0ecc 100644 --- a/unit_tests/BareField/BareField.cpp +++ b/unit_tests/BareField/BareField.cpp @@ -171,6 +171,60 @@ TYPED_TEST(BareFieldTest, IndexedSubdomainExpressionAssignment) { }); } +TYPED_TEST(BareFieldTest, IndexedSubdomainDifferentStrideExpressionAssignment) { + using T = typename TestFixture::value_type; + + auto& lhs = this->field; + lhs->operator=(0); + + typename TestFixture::field_type rhs(this->layout); + auto rhsView = rhs.getView(); + const auto lDom = rhs.getLayout().getLocalNDIndex(); + Kokkos::parallel_for("Set rhs field", rhs.getFieldRangePolicy(), FieldVal{rhsView, lDom}); + Kokkos::fence(); + + ippl::NDIndex lhsSubdomain = lhs->getOwned(); + lhsSubdomain[0] = ippl::Index(lhsSubdomain[0].first(), lhsSubdomain[0].last(), + 2 * lhsSubdomain[0].stride()); + + ippl::NDIndex rhsSubdomain = rhs.getOwned(); + const int rhsLast = rhsSubdomain[0].first() + + (lhsSubdomain[0].length() - 1) * rhsSubdomain[0].stride(); + rhsSubdomain[0] = ippl::Index(rhsSubdomain[0].first(), rhsLast, rhsSubdomain[0].stride()); + + ASSERT_EQ(lhsSubdomain[0].length(), rhsSubdomain[0].length()); + ASSERT_NE(lhsSubdomain[0].stride(), rhsSubdomain[0].stride()); + + (*lhs)[lhsSubdomain] = rhs[rhsSubdomain] + T(2); + + auto lhsMirror = lhs->getHostMirror(); + auto rhsMirror = rhs.getHostMirror(); + Kokkos::deep_copy(lhsMirror, lhs->getView()); + Kokkos::deep_copy(rhsMirror, rhs.getView()); + + const auto owned = lhs->getOwned(); + const int nghost = lhs->getNghost(); + nestedViewLoop(lhsMirror, nghost, [&](const Idx... args) { + std::array pointIndices; + typename ippl::RangePolicy::index_array_type + rhsCoords; + const std::array viewCoords{static_cast(args)...}; + for (unsigned d = 0; d < TestFixture::dim; ++d) { + const int global = owned[d].first() + + static_cast(viewCoords[d] - nghost) * owned[d].stride(); + pointIndices[d] = ippl::Index(global, global); + + const int rel = (global - lhsSubdomain[d].first()) / lhsSubdomain[d].stride(); + const int rhsGlobal = rhsSubdomain[d].first() + rel * rhsSubdomain[d].stride(); + rhsCoords[d] = (rhsGlobal - owned[d].first()) / owned[d].stride() + nghost; + } + auto point = std::make_from_tuple>(pointIndices); + const T expected = lhsSubdomain.contains(point) ? ippl::apply(rhsMirror, rhsCoords) + T(2) + : T(0); + assertEqual(expected, lhsMirror(args...)); + }); +} + TYPED_TEST(BareFieldTest, ChainedIndexedSubdomainScalarAssignment) { using T = typename TestFixture::value_type;