Merge add segmentedField pull request #202 see #202

ENH: add segmentedField to represent vector of vectors

CHANGELOG.md

@@ -1,5 +1,6 @@
 # Version 0.2.0 (unreleased)
 ## Features
+- add segmentedField to represent vector of vectors [#202](https://github.com/exasim-project/NeoFOAM/pull/202)
 - Adds a minimal implementation linear algebra functionality [#219](https://github.com/exasim-project/NeoFOAM/pull/219)
 ## Fixes
 # Version 0.1.0

doc/basics/index.rst

@@ -10,6 +10,7 @@ Basics
 
     executor.rst
     fields.rst
+    segmentedField.rst
     algorithms.rst
     first_kernel.rst
     registerclass.rst

doc/basics/segmentedField.rst

@@ -0,0 +1,55 @@
+.. _fvcc_segmentedFields:
+
+SegmentedField
+^^^^^^^^^^^^^^
+
+SegmentedField is a template class that represents a field divided into multiple segments and can represent vector of vector of a defined ValueType.
+It also allows the definition of an IndexType, so each segment of the vector can be addressed.
+It can be used to represent cell to cell stencil.
+
+.. code-block:: cpp
+
+    NeoFOAM::Field<NeoFOAM::label> values(exec, {0, 1, 2, 3, 4, 5, 6, 7, 8, 9});
+    NeoFOAM::Field<NeoFOAM::localIdx> segments(exec, {0, 2, 4, 6, 8, 10});
+
+    NeoFOAM::SegmentedField<NeoFOAM::label, NeoFOAM::localIdx> segField(values, segments);
+    auto [valueSpan, segment] = segField.spans();
+    auto segView = segField.view();
+    NeoFOAM::Field<NeoFOAM::label> result(exec, 5);
+
+    NeoFOAM::fill(result, 0);
+    auto resultSpan = result.span();
+
+    parallelFor(
+        exec,
+        {0, segField.numSegments()},
+        KOKKOS_LAMBDA(const size_t segI) {
+            // check if it works with bounds
+            auto [bStart, bEnd] = segView.bounds(segI);
+            auto bVals = valueSpan.subspan(bStart, bEnd - bStart);
+            for (auto& val : bVals)
+            {
+                resultSpan[segI] += val;
+            }
+
+            // check if it works with range
+            auto [rStart, rLength] = segView.range(segI);
+            auto rVals = valueSpan.subspan(rStart, rLength);
+            for (auto& val : rVals)
+            {
+                resultSpan[segI] += val;
+            }
+
+            // check with subspan
+            auto vals = segView.span(segI);
+            for (auto& val : vals)
+            {
+                resultSpan[segI] += val;
+            }
+        }
+    );
+
+In this example, each of the five segments would have a size of two.
+This data allows the representation of stencils in a continuous memory layout, which can be beneficial for performance optimization in numerical simulations especially on GPUs.
+
+The spans method return the value and segment span and it is also possible to return a view that can also be called on a device

include/NeoFOAM/core/parallelAlgorithms.hpp

@@ -90,7 +90,6 @@ void parallelFor(Field<ValueType>& field, Kernel kernel)
     std::visit([&](const auto& e) { parallelFor(e, field, kernel); }, field.exec());
 }
 
-
 template<typename Executor, typename Kernel, typename T>
 void parallelReduce(
     [[maybe_unused]] const Executor& exec, std::pair<size_t, size_t> range, Kernel kernel, T& value
@@ -165,5 +164,49 @@ void parallelReduce(Field<ValueType>& field, Kernel kernel, T& value)
     );
 }
 
+template<typename Executor, typename Kernel>
+void parallelScan(
+    [[maybe_unused]] const Executor& exec, std::pair<size_t, size_t> range, Kernel kernel
+)
+{
+    auto [start, end] = range;
+    using runOn = typename Executor::exec;
+    Kokkos::parallel_scan("parallelScan", Kokkos::RangePolicy<runOn>(start, end), kernel);
+}
+
+template<typename Kernel>
+void parallelScan(const NeoFOAM::Executor& exec, std::pair<size_t, size_t> range, Kernel kernel)
+{
+    std::visit([&](const auto& e) { parallelScan(e, range, kernel); }, exec);
+}
+
+template<typename Executor, typename Kernel, typename ReturnType>
+void parallelScan(
+    [[maybe_unused]] const Executor& exec,
+    std::pair<size_t, size_t> range,
+    Kernel kernel,
+    ReturnType& returnValue
+)
+{
+    auto [start, end] = range;
+    using runOn = typename Executor::exec;
+    Kokkos::parallel_scan(
+        "parallelScan", Kokkos::RangePolicy<runOn>(start, end), kernel, returnValue
+    );
+}
+
+template<typename Kernel, typename ReturnType>
+void parallelScan(
+    const NeoFOAM::Executor& exec,
+    std::pair<size_t, size_t> range,
+    Kernel kernel,
+    ReturnType& returnValue
+)
+{
+    return std::visit(
+        [&](const auto& e) { return parallelScan(e, range, kernel, returnValue); }, exec
+    );
+}
+
 
 } // namespace NeoFOAM

include/NeoFOAM/fields/segmentedField.hpp

@@ -0,0 +1,215 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2023 NeoFOAM authors
+#pragma once
+
+#include "NeoFOAM/core/primitives/label.hpp"
+#include "NeoFOAM/fields/field.hpp"
+
+namespace NeoFOAM
+{
+
+/**
+ * @brief Compute segment offsets from an input field corresponding to lengths by computing a prefix
+ * sum.
+ *
+ * The offsets are computed by a prefix sum of the input values. So, with given
+ * input of {1, 2, 3, 4, 5} the offsets are {0, 1, 3, 6, 10, 15}.
+ * Note that the length of offSpan must be  length of valSpan + 1
+ * and are all elements of offSpan are required to be zero
+ *
+ * @param[in] in The values to compute the offsets from.
+ * @param[in,out] offsets The field to store the resulting offsets in.
+ */
+template<typename IndexType>
+IndexType segmentsFromIntervals(const Field<IndexType>& intervals, Field<IndexType>& offsets)
+{
+    IndexType finalValue = 0;
+    auto inSpan = intervals.span();
+    auto offsSpan = offsets.span();
+    NF_ASSERT_EQUAL(inSpan.size() + 1, offsSpan.size());
+    NeoFOAM::parallelScan(
+        intervals.exec(),
+        {1, offsSpan.size()},
+        KOKKOS_LAMBDA(const std::size_t i, NeoFOAM::localIdx& update, const bool final) {
+            update += inSpan[i - 1];
+            if (final)
+            {
+                offsSpan[i] = update;
+            }
+        },
+        finalValue
+    );
+    return finalValue;
+}
+
+/**
+ * @brief A class representing a segment of indices.
+ *
+ * @tparam IndexType The type of the indices.
+ */
+template<typename ValueType, typename IndexType = NeoFOAM::localIdx>
+class SegmentedFieldView
+{
+public:
+
+    /**
+     * @brief A span with the values.
+     */
+    std::span<ValueType> values;
+
+    /**
+     * @brief A span of indices representing the segments.
+     */
+    std::span<IndexType> segments;
+
+    /**
+     * @brief Get the bounds of a segment.
+     *
+     * @param segI The index of the segment.
+     * @return A pair of indices representing the start and end of the segment.
+     */
+    KOKKOS_INLINE_FUNCTION
+    Kokkos::pair<IndexType, IndexType> bounds(std::size_t segI) const
+    {
+        return Kokkos::pair<IndexType, IndexType> {segments[segI], segments[segI + 1]};
+    }
+
+    /**
+     * @brief Get the range, ie. [start,end), of a segment.
+     *
+     * @param segI The index of the segment.
+     * @return A pair of indices representing the start and length of the segment.
+     */
+    KOKKOS_INLINE_FUNCTION
+    Kokkos::pair<IndexType, IndexType> range(std::size_t segI) const
+    {
+        return Kokkos::pair<IndexType, IndexType> {
+            segments[segI], segments[segI + 1] - segments[segI]
+        };
+    }
+
+    /**
+     * @brief Get a subspan of values corresponding to a segment.
+     *
+     * @tparam ValueType The type of the values.
+     * @param segI The index of the segment.
+     * @return A subspan of values corresponding to the segment.
+     */
+    KOKKOS_INLINE_FUNCTION std::span<ValueType> span(std::size_t segI) const
+    {
+        auto [start, length] = range(segI);
+        return values.subspan(start, length);
+    }
+
+    /**
+     * @brief Access an element of the segments.
+     *
+     * @param i The index of the element.
+     * @return The value of the element at the specified index.
+     */
+    KOKKOS_INLINE_FUNCTION
+    IndexType operator[](std::size_t i) const { return segments[i]; }
+};
+
+/**
+ * @class SegmentedField
+ * @brief Data structure that stores a segmented fields or a vector of vectors
+ *
+ * @ingroup Fields
+ */
+template<typename ValueType, typename IndexType>
+class SegmentedField
+{
+public:
+
+
+    /**
+     * @brief Create a segmented field with a given size and number of segments.
+     * @param exec  Executor associated to the matrix
+     * @param size  size of the matrix
+     * @param numSegments  number of segments
+     */
+    SegmentedField(const Executor& exec, size_t size, size_t numSegments)
+        : values_(exec, size), segments_(exec, numSegments + 1)
+    {}
+
+    /*
+     * @brief Create a segmented field from intervals.
+     * @param intervals The intervals to create the segmented field from.
+     * @note The intervals are the lengths of each segment
+     */
+    SegmentedField(const Field<IndexType>& intervals)
+        : values_(intervals.exec(), 0),
+          segments_(intervals.exec(), intervals.size() + 1, IndexType(0))
+    {
+        IndexType valueSize = segmentsFromIntervals(intervals, segments_);
+        values_ = Field<ValueType>(intervals.exec(), valueSize);
+    }
+
+
+    /**
+     * @brief Constructor to create a segmentedField from values and the segments.
+     * @param values The values of the segmented field.
+     * @param segments The segments of the segmented field.
+     */
+    SegmentedField(const Field<ValueType>& values, const Field<IndexType>& segments)
+        : values_(values), segments_(segments)
+    {
+        NF_ASSERT(values.exec() == segments.exec(), "Executors are not the same.");
+    }
+
+
+    /**
+     * @brief Get the executor associated with the segmented field.
+     * @return Reference to the executor.
+     */
+    const Executor& exec() const { return values_.exec(); }
+
+    /**
+     * @brief Get the size of the segmented field.
+     * @return The size of the segmented field.
+     */
+    size_t size() const { return values_.size(); }
+
+    /**
+     * @brief Get the number of segments in the segmented field.
+     * @return The number of segments.
+     */
+    size_t numSegments() const { return segments_.size() - 1; }
+
+
+    /**
+     * @brief get a view of the segmented field
+     * @return Span of the fields
+     */
+    [[nodiscard]] SegmentedFieldView<ValueType, IndexType> view() &
+    {
+        return SegmentedFieldView<ValueType, IndexType> {values_.span(), segments_.span()};
+    }
+
+    // ensures no return a span of a temporary object --> invalid memory access
+    [[nodiscard]] SegmentedFieldView<ValueType, IndexType> view() && = delete;
+
+    /**
+     * @brief get the combined value and range spans of the segmented field
+     * @return Combined value and range spans of the fields
+     */
+    [[nodiscard]] std::pair<std::span<ValueType>, std::span<IndexType>> spans() &
+    {
+        return {values_.span(), segments_.span()};
+    }
+
+    // ensures not to return a span of a temporary object --> invalid memory access
+    [[nodiscard]] std::pair<std::span<ValueType>, std::span<IndexType>> spans() && = delete;
+
+    const Field<ValueType>& values() const { return values_; }
+
+    const Field<IndexType>& segments() const { return segments_; }
+
+private:
+
+    Field<ValueType> values_;
+    Field<IndexType> segments_; //!< stores the [start, end) of segment i at index i, i+1
+};
+
+} // namespace NeoFOAM