diff --git a/src/include/serif/discretization/blocks/base.hpp b/src/include/serif/discretization/blocks/base.hpp index 6b194140..1ca25297 100644 --- a/src/include/serif/discretization/blocks/base.hpp +++ b/src/include/serif/discretization/blocks/base.hpp @@ -14,6 +14,17 @@ namespace serif::discretization::blocks { */ struct block {}; + /** + * @enum RepresentedValueKind + * @brief Enum representing the kind of value a block represents, either a physical quantity (value) or a residual (residual). + * + * Generally this is mostly used in the solver/normalization system. + */ + enum class RepresentedValueKind { + value, ///< A value representing a physical quantity (e.g., temperature, pressure). + residual ///< A residual representing the difference between the current state and the desired state in a numerical method. + }; + /** * @struct residual_block_base * @brief Base class for residual blocks, derived from block. diff --git a/src/include/serif/discretization/blocks/keys.hpp b/src/include/serif/discretization/blocks/keys.hpp new file mode 100644 index 00000000..a85ca168 --- /dev/null +++ b/src/include/serif/discretization/blocks/keys.hpp @@ -0,0 +1,16 @@ +#pragma once + +/** + * @file keys.hpp + * @author Emily M. Boudreaux + * + * @brief Keys for the shared layout of a discretized state and residual in a block-structured system. + */ + + +namespace serif::discretization::blocks::keys { + /** + * @brief Share one independently ordered state and residual layout across numerical services. + */ + struct Layout final {}; +} diff --git a/src/include/serif/discretization/blocks/selection.hpp b/src/include/serif/discretization/blocks/selection.hpp new file mode 100644 index 00000000..ad6556bb --- /dev/null +++ b/src/include/serif/discretization/blocks/selection.hpp @@ -0,0 +1,47 @@ +#pragma once + +#include +#include +#include "serif/discretization/blocks/system_layout.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::discretization::blocks { + /** + * @brief Select named packed coordinates in their requested order, independently for values and residuals. + * + * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals. + * @tparam Layout The block layout type, e.g. SystemLayout. + * @tparam Coordinates The coordinate tags to select, e.g. keys::Density, + * + * @param layout The block layout from which to select coordinates. + * @param keys A list of coordinate tags to select. Note that this is not used at runtime, but rather to select the correct overloads of layout.size and layout.offset. + * + * @returns A vector of packed indices corresponding to the requested coordinates in the order requested. + * + * Given some layout (see @ref SystemLayout) and a list of coordinate tags return the corresponding packed indices in the order requested. + * The result is suitable for indexing a contiguous vector of values or residuals. + */ + template + [[nodiscard]] std::vector coordinate_indices( + const Layout& layout, + utils::dependencies::KeyList /*keys*/ + ) { + utils::errors::serif_ensure(layout.template total_size() <= static_cast(std::numeric_limits::max()),"Named coordinate selection exceeds the backend index range; use a layout with representable packed indices."); + + std::vector result; + result.reserve((std::size_t{0} + ... + layout.template size())); + + const auto append = [&]() { + const int offset = static_cast(layout.template offset()); + const int size = static_cast(layout.template size()); + + for (int index = 0; index < size; ++index) { + result.push_back(offset + index); + } + }; + + (append.template operator()(), ...); + + return result; + } +} diff --git a/src/include/serif/discretization/blocks/system_layout.hpp b/src/include/serif/discretization/blocks/system_layout.hpp new file mode 100644 index 00000000..370d21aa --- /dev/null +++ b/src/include/serif/discretization/blocks/system_layout.hpp @@ -0,0 +1,188 @@ +#pragma once + +/** + * @file system_layout.hpp + * @author Emily M. Boudreaux + * + * @brief Describe the independently ordered physical state and residual coordinates of a discrete system. This is an + * accounting layer to keep index offsets straight in memory. + */ + +#include +#include +#include +#include + +#include "serif/discretization/blocks/base.hpp" +#include "serif/utils/dependencies/key_list_operations.hpp" + +namespace serif::discretization::blocks +{ + /** + * @brief Sum some memory span with checking for if any element is too large to fit in a size_t + * + * @param sizes A span of sizes to sum + * @return The sum of the sizes + * + * @throws std::overflow_error If an element of sizes is too large to fit in a size_t + * + * @todo Really this does not need to be specific to this module, this might more generally be named something like "safe sum" and moved to a more general location. + */ + [[nodiscard]] std::size_t checked_block_size(std::span sizes); + + /** + * @brief Describe the independently ordered physical state and residual coordinates of a discrete system. + * + * @tparam ValuesT A KeyList of the physical coordinates represented in the state vector. + * @tparam ResidualsT A KeyList of the physical coordinates represented in the residual vector. + * + */ + template + struct SystemLayout final { + using Values = std::remove_cvref_t; ///< The physical coordinates represented in the state vector. + using Residuals = std::remove_cvref_t; ///< The physical coordinates represented in the residual vector. + + /** + * @brief Select the coordinate tag at compile time based on the kind of the represented values selected (e.g. RepresentedValueKind::values or RepresentedValueKind::residuals). + * + * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals. + * + * @requires Kind == RepresentedValueKind::value || Kind == RepresentedValueKind::residual + * + * @note This requires clause is not strictly necessary right now since we have a template type argument; however, if we ever add + * more RepresentedValueKind types, this will ensure that we don't accidentally allow them to be used here. + */ + template + requires(Kind == RepresentedValueKind::value || Kind == RepresentedValueKind::residual) + using Coordinates = std::conditional_t; + + std::array> state_sizes; ///< The coefficient counts for each physical coordinate in the state vector, in the order of Values. + std::array> residual_sizes; ///< The coefficient counts for each physical coordinate in the residual vector, in the order of Residuals. + + /** + * @brief Get the sizes of the physical coordinates for the requested kind. + * + * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals. + * @returns A span of sizes for the physical coordinates of the requested kind. + * + * @requires Kind == RepresentedValueKind::value || Kind == RepresentedValueKind::residual + */ + template + requires(Kind == RepresentedValueKind::value || Kind == RepresentedValueKind::residual) + [[nodiscard]] std::span sizes() const noexcept { + if constexpr (Kind == RepresentedValueKind::value) { + return state_sizes; + } + else { + return residual_sizes; + } + } + + /** + * @brief Locate a physical identity on the requested side of the system. + * + * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals. + * @tparam Coordinate The coordinate type to locate. + * @returns The index of the coordinate in the requested kind. + * + * @requires Coordinate is contained in the KeyList of the requested kind. + */ + template + requires utils::dependencies::ContainsKey> + [[nodiscard]] static consteval std::size_t index() { + if constexpr (Kind == RepresentedValueKind::value) { + return utils::dependencies::key_index(); + } + else { + return utils::dependencies::key_index(); + } + } + + /** + * @brief Read the coefficient count for a physical coordinate. + * + * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals. + * @tparam Coordinate The coordinate type to locate. + * @returns The coefficient count for the specified coordinate. + * + * @requires Coordinate is contained in the KeyList of the requested kind. + */ + template + requires utils::dependencies::ContainsKey> + [[nodiscard]] std::size_t size() const noexcept { + return sizes()[index()]; // A note from EMB on Oct 14 2026 to a future dev. Maybe I am sorry I picked C++? We end up with stuff like the syntax on the left... But perhaps we would one way or the other. + } + + /** + * @brief Locate a coordinate in packed storage, rejecting size overflow before creating a view. Note + * that this specifically locates the first coefficient of a coordinate, or the offset to where that coordinate's coefficients begin in packed storage. + * + * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals. + * @tparam Coordinate The coordinate type to locate. + * @returns The offset of the coordinate in the requested kind. (e.g. the index of the first coefficient of that coordinate in packed storage.) + * + * @requires Coordinate is contained in the KeyList of the requested kind. + */ + template + requires utils::dependencies::ContainsKey> + [[nodiscard]] std::size_t offset() const { + return checked_block_size(sizes().first(index())); + } + + /** + * @brief Count all coefficients of one kind + * + * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals. + * @returns The total coefficient count for the requested kind. + */ + template + requires(Kind == RepresentedValueKind::value || Kind == RepresentedValueKind::residual) + [[nodiscard]] std::size_t total_size() const { + return checked_block_size(sizes()); + } + + /** + * @brief Default equality comparison + */ + [[nodiscard]] bool operator==(const SystemLayout&) const = default; + }; + + namespace detail + { + /** + * @brief Default false specialization for is_system_layout, specialized to true for SystemLayout. + * @tparam Candidate The type to check for being a SystemLayout. + */ + template + inline constexpr bool is_system_layout = false; + + /** + * @brief True specialization for is_system_layout, specialized to true for SystemLayout. + * @tparam Values The KeyList of physical coordinates represented in the state vector. + * @tparam Residuals The KeyList of physical coordinates represented in the residual vector. + */ + template + inline constexpr bool is_system_layout> = true; + } // namespace detail + + /** + * @brief Compile time check for a SystemLayout type, e.g. SystemLayout. + */ + template + concept SystemLayoutType = detail::is_system_layout>; + + /** + * @brief Build both coordinate orderings from the same physical size resolver. + * @param size A callable whose templated operator()() returns that coordinate's coefficient count. + */ + template + [[nodiscard]] constexpr std::remove_cvref_t make_system_layout(const SizeResolver& size) { + using Result = std::remove_cvref_t; + + const auto sizes = [&](utils::dependencies::KeyList) { + return std::array{static_cast(size.template operator()())...}; + }; + + return {.state_sizes = sizes(typename Result::Values{}), .residual_sizes = sizes(typename Result::Residuals{})}; + } +} // namespace serif::discretization::blocks diff --git a/src/include/serif/discretization/blocks/traits.hpp b/src/include/serif/discretization/blocks/traits.hpp index b6d4782f..ad7c1da9 100644 --- a/src/include/serif/discretization/blocks/traits.hpp +++ b/src/include/serif/discretization/blocks/traits.hpp @@ -67,4 +67,7 @@ namespace serif::discretization::blocks { static constexpr int order = FunctionalSpace::order; ///< The order of the functional space associated with the term. e.g. 1st order, 2nd order, etc. static constexpr std::size_t rank = FunctionalSpace::rank::value; ///< The rank of the functional space associated with the term. }; -} \ No newline at end of file + /** @brief A unique scalar boundary value is available for continuous scalar fields and scalar globals. */ + template + concept ScalarTraceTerm = GlobalScalarTerm || (SpatialTerm && SpatialTermTraits::rank == 0 && SpatialTermTraits::family == functional_space::Family::h1); +} diff --git a/src/include/serif/discretization/blocks/vector_view.hpp b/src/include/serif/discretization/blocks/vector_view.hpp new file mode 100644 index 00000000..0cf0ef59 --- /dev/null +++ b/src/include/serif/discretization/blocks/vector_view.hpp @@ -0,0 +1,113 @@ +#pragma once + +/** + * @file vector_view.hpp + * @author Emily M. Boudreaux + * + * @brief This file includes a class for viewing some SystemLayout as a span of coefficients without copying the underlying storage. + */ + +#include +#include +#include +#include + +#include "serif/discretization/blocks/system_layout.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::discretization::blocks +{ + /** + * @brief Build a non-owning view of a packed layout as a span. + * + * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals. + * @tparam Layout The block layout type, e.g. SystemLayout. + * @tparam Scalar The scalar type of the coefficients, e.g. double. + * + * @requires Kind == RepresentedValueKind::value || Kind == RepresentedValueKind::residual + * + * @note The layout must remain unchanged while the view is used. Probably we want to, in the future, add some way to + * validate this; however, I am at the moment at a loss for how to do so. + */ + template + requires(Kind == RepresentedValueKind::value || Kind == RepresentedValueKind::residual) + class BlockVectorView final + { + public: + using LayoutType = std::remove_cvref_t; ///< Alias for the layout type, e.g. SystemLayout. + + /** + * @brief Construct one of these views from a layout and a span of coefficients. + * + * @param layout The layout to use for the view. + * @param values The span of coefficients to view. + * + * @throws std::domain_error If the span of coefficients does not match the total size of the layout. + * + * @todo Really this could go in many places; however, I think we should move away from throwing std exceptions and introduce our own exception hierarchy. + */ + BlockVectorView( + const LayoutType& layout, + std::span values + ) : + m_layout(layout), + m_values(values) { + const std::size_t required = layout.template total_size(); + + utils::errors::serif_ensure(values.size() == required, "Packed storage has {} coefficients but the system layout requires {}. Bind the matching state or residual vector.", values.size(), required); + } + + /** + * @brief Reject temporary metadata whose offsets would dangle before the coefficients are accessed. + */ + BlockVectorView(const LayoutType&&, std::span) = delete; + + /** + * @brief Prevent implicit const-span conversion from retaining a temporary owning container. + */ + template + requires(!std::is_lvalue_reference_v && !std::ranges::borrowed_range) + BlockVectorView(const LayoutType&, Range&&) = delete; + + /** + * @brief Access the requested physical coordinate while retaining storage constness. + * + * @tparam Coordinate The coordinate tag to access + * + * @requires A layout type with kind Kind must contain the requested coordinate tag. + */ + template + requires utils::dependencies::ContainsKey> + [[nodiscard]] std::span get() const + { + const std::size_t offset = m_layout.template offset(); + const std::size_t count = m_layout.template size(); + + utils::errors::serif_ensure(offset <= m_values.size() && count <= m_values.size() - offset,"The system layout changed after this vector view was bound. Rebind the view to matching storage."); + + return m_values.subspan(offset, count); + } + + private: + const LayoutType& m_layout; ///< The layout to use for the view. The layout must remain unchanged while the view is used. + std::span m_values; ///< The span of coefficients to view. + }; + + /** + * @brief A convenient alias for a view of the physical state coefficients. + * + * @tparam Layout The block layout type, e.g. SystemLayout. + * @tparam Scalar The scalar type of the coefficients, e.g. double. + */ + template + using StateView = BlockVectorView; + + /** + * @brief A convenient alias for a view of the residual coefficients. + * + * @tparam Layout The block layout type, e.g. SystemLayout. + * @tparam Scalar The scalar type of the coefficients, e.g. double. + */ + template + using ResidualView = BlockVectorView; +} // namespace serif::discretization::blocks