From 11079d9a7264c55470b59f92bde85b17b76897cb Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Thu, 19 Mar 2026 16:11:04 -0500 Subject: [PATCH 01/10] SolutionTransfer: do not force const interface to vector --- source/solution_transfer.template.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/solution_transfer.template.h b/source/solution_transfer.template.h index d9b1f4b09..ae622f742 100644 --- a/source/solution_transfer.template.h +++ b/source/solution_transfer.template.h @@ -727,7 +727,7 @@ namespace ryujin update_new_state_vector(); - const auto &precomputed = std::get<1>(new_state_vector); + auto &precomputed = std::get<1>(new_state_vector); /* The limiter requires valid precomputed values. Therefore, update: */ const auto update_precomputed_values = [&]() { From aa601fb8f0299eebceca052c290aeb2db3213521 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Fri, 20 Mar 2026 19:42:08 -0500 Subject: [PATCH 02/10] MultiComponentVector: split into base class and view --- source/multicomponent_vector.h | 579 +++++++++++++++++++++++++-------- 1 file changed, 437 insertions(+), 142 deletions(-) diff --git a/source/multicomponent_vector.h b/source/multicomponent_vector.h index 16a49487a..57071ab80 100644 --- a/source/multicomponent_vector.h +++ b/source/multicomponent_vector.h @@ -20,7 +20,7 @@ namespace ryujin * This function takes a scalar MPI partitioner @p scalar_partitioner as * argument and returns a shared pointer to a new "vector" multicomponent * partitioner that defines storage and MPI synchronization for a vector - * consisting of @p n_components components. The vector partitioner is + * consisting of @p n_comp components. The vector partitioner is * intended to efficiently store non-scalar vectors such as the state * vectors U. Let (U_i)_k denote the k-th component of a state vector * element U_i, we then store @@ -39,75 +39,189 @@ namespace ryujin create_vector_partitioner( const std::shared_ptr &scalar_partitioner, - const unsigned int n_components); + const unsigned int n_comp); + + + template ::size(), + typename MemorySpace = dealii::MemorySpace::Host, + bool writable = true> + class MultiComponentVectorView; /** * A wrapper around dealii::LinearAlgebra::distributed::Vector - * that stores a vector element of @p n_components components per entry + * that stores a vector element of @p n_comp components per entry * (instead of a scalar value). * * @ingroup SIMD */ template ::size()> class MultiComponentVector - : private dealii::LinearAlgebra::distributed::Vector + : public MultiComponentVectorView { public: /** - * @name Typedefs and constexpr constants + * @name Constructor, reinitialization, assignment */ //@{ + MultiComponentVector() = default; + /** - * Shorthand typedef for the underlying dealii::VectorizedArray type - * used to insert and extract SIMD packed values from the - * MultiComponentVector. + * Reinitializes the MultiComponentVector with a vector MPI partitioner + * that was created first with create_vector_partitioner(). */ - using VectorizedArray = dealii::VectorizedArray; + void reinit_with_vector_partitioner( + const std::shared_ptr + &vector_partitioner); /** - * Shorthand typedef for the underlying scalar - * dealii::LinearAlgebra::distributed::Vector used to insert - * and extract a single component of the MultiComponentVector. + * Reinitializes the MultiComponentVector with a scalar MPI partitioner. + * The function calls create_vector_partitioner() internally to + * create and store a corresponding "vector" MPI partitioner. */ - using ScalarHostVector = - dealii::LinearAlgebra::distributed::Vector; + void reinit_with_scalar_partitioner( + const std::shared_ptr + &scalar_partitioner); //@} /** - * @name Constructor and reinitialization + * Memory space access and synchronization: */ //@{ /** - * Default constructor + * Return a writable view on the sparse matrix for the selected memory + * space. */ - MultiComponentVector() = default; + template + MultiComponentVectorView + get_view(); /** - * Reinitializes the MultiComponentVector with a vector MPI partitioner - * that was created first with create_vector_partitioner(). + * Return a read-only view on the sparse matrix for the selected memory + * space. */ - void reinit_with_vector_partitioner( - const std::shared_ptr - &vector_partitioner); + template + MultiComponentVectorView + get_view() const; /** - * Reinitializes the MultiComponentVector with a scalar MPI partitioner. - * The function calls create_vector_partitioner() internally to - * create and store a corresponding "vector" MPI partitioner. + * Returns true if the templated memory space is the currently active + * memory space. */ - void reinit_with_scalar_partitioner( - const std::shared_ptr - &scalar_partitioner); + template + bool is_active_memory_space() const; + + /** + * Move internal data from the currently active memory space to the + * templated memory space. + */ + template + void move_to_memory_space(); + + //@} + + //@} + /** + * MPI synchronization. + */ + //@{ + + /** + * MPI synchronization: Zero out all ghost rows. + */ + template + void zero_out_ghost_values_on_memory_space(); + + /** + * MPI synchronization: Import all ghost values from neighboring MPI + * ranks on the templated memory space. + */ + template + void update_ghost_values_on_memory_space(); + + /** + * MPI synchronization: Copy the data that has accumulated in the + * ghost range to the owning processor. This function operates on the + * templated memory space. + */ + template + void compress_on_memory_space(dealii::VectorOperation::values operation); + + private: + //@} + /** + * @name Internal fields, methods, and friends + */ + //@{ + + dealii::LinearAlgebra::distributed::Vector + host_vector_; + + dealii::LinearAlgebra::distributed::Vector + default_vector_; + + bool host_space_active_ = true; + + template + friend class MultiComponentVectorView; + + //@} + }; + + + template + class MultiComponentVectorView + { + public: + /** + * @name Constructor and reinitialization + */ + //@{ + + MultiComponentVectorView() = default; + + MultiComponentVectorView(MultiComponentVector + &multi_component_vector) + requires(writable); + + MultiComponentVectorView( + const MultiComponentVector + &multi_component_vector) + requires(!writable); + + template + void reinit(MultiComponentVector &multi_component_vector) + requires(writable != std::is_const_v); + + //@} + /** + * @name Typedefs and constexpr constants + */ + //@{ + + /** + * Shorthand typedef for the underlying scalar + * dealii::LinearAlgebra::distributed::Vector used to + * insert and extract a single component of the MultiComponentVector. + */ + using ScalarVector = + dealii::LinearAlgebra::distributed::Vector; //@} /** - * @name Extracting and inserting a single component stored in a - * ScalarHostVector + * @name Extracting and inserting components, scaled addition */ //@{ @@ -130,7 +244,7 @@ namespace ryujin * vectors). */ template - void extract_component(ScalarHostVector &scalar_vector, + void extract_component(ScalarVector &scalar_vector, unsigned int component, const Functor &functor = std::identity{}) const; @@ -152,9 +266,10 @@ namespace ryujin * single scalar vectors by deal.II interpolation functions. */ template - void insert_component(const ScalarHostVector &scalar_vector, + void insert_component(const ScalarVector &scalar_vector, unsigned int component, - const Functor &functor = std::identity{}); + const Functor &functor = std::identity{}) + requires writable; /** * Variant of the method above that reads values out of a @@ -163,7 +278,16 @@ namespace ryujin template void insert_component(const dealii::Vector &scalar_vector, unsigned int component, - const Functor &functor = std::identity{}); + const Functor &functor = std::identity{}) + requires writable; + + /** + * Scaled addition of the given vector $U$ and the argument vector + * $V$: $U\leftarrow s\,U+a\,V$. + */ + void + sadd(const Number s, const Number a, const MultiComponentVectorView &V) + requires writable; //@} /** @@ -177,26 +301,26 @@ namespace ryujin * * If the template parameter @a Number2 is a VectorizedArray then * the function returns a SIMD vectorized dealii::Tensor populated with - * entries from the @p n_components component vectors stored at + * entries from the @p n_comp component vectors stored at * indices i, i+1, ..., i+simd_length-1. * - * @note This function is only available if `n_components` is equal to 1. + * @note This function is only available if `n_comp` is equal to 1. */ template - Number2 read_entry(const unsigned int i) const; + DEAL_II_HOST_DEVICE Number2 read_entry(const unsigned int i) const; /** * Variant of above function. * * Returns a SIMD vectorized dealii::Tensor populated with entries from - * the @p n_components component vectors stored at indices *(js), *(js+1), + * the @p n_comp component vectors stored at indices *(js), *(js+1), * ..., *(js+simd_length-1), i.e., @p js has to point to an array of * size @p simd_length containing all indices. * - * @note This function is only available if `n_components` is equal to 1. + * @note This function is only available if `n_comp` is equal to 1. */ template - Number2 read_entry(const unsigned int *js) const; + DEAL_II_HOST_DEVICE Number2 read_entry(const unsigned int *js) const; /** * Return the tensor-valued entry indexed by @p i. @@ -207,42 +331,44 @@ namespace ryujin * indices i, i+1, ..., i+simd_length-1. */ template > - Tensor read_tensor(const unsigned int i) const; + typename Tensor = dealii::Tensor<1, n_comp, Number2>> + DEAL_II_HOST_DEVICE Tensor read_tensor(const unsigned int i) const; /** * Variant of above function. * * Returns a SIMD vectorized dealii::Tensor populated with entries from - * the @p n_components component vectors stored at indices *(js), *(js+1), + * the @p n_comp component vectors stored at indices *(js), *(js+1), * ..., *(js+simd_length-1), i.e., @p js has to point to an array of * size @p simd_length containing all indices. */ template > - Tensor read_tensor(const unsigned int *js) const; + typename Tensor = dealii::Tensor<1, n_comp, Number2>> + DEAL_II_HOST_DEVICE Tensor read_tensor(const unsigned int *js) const; /** * Write a (scalar valued) @p entry to the vector at position by @p i. * * If the template parameter @a Number2 is a VectorizedArray then * the function takes a SIMD vectorized @p tensor as argument instead - * and updates the values of the @p n_components component vectors at + * and updates the values of the @p n_comp component vectors at * indices i, i+1, ..., i+simd_length_1. with the values supplied by @p * tensor. * - * @note This function is only available if `n_components` is equal to 1. + * @note This function is only available if `n_comp` is equal to 1. */ template - void write_entry(const Number2 &entry, const unsigned int i); + DEAL_II_HOST_DEVICE void write_entry(const Number2 &entry, + const unsigned int i) const + requires writable; /** - * Update the values of the @p n_components component vector at index + * Update the values of the @p n_comp component vector at index * @p i with the values supplied by @p tensor. * * If the template parameter @a Number2 is a VectorizedArray then * the function takes a SIMD vectorized @p tensor as argument instead - * and updates the values of the @p n_components component vectors at + * and updates the values of the @p n_comp component vectors at * indices i, i+1, ..., i+simd_length_1. with the values supplied by @p * tensor. * @@ -251,33 +377,37 @@ namespace ryujin * Number, and has a type trait `value_type`. */ template > - void write_tensor(const Tensor &tensor, const unsigned int i); + typename Tensor = dealii::Tensor<1, n_comp, Number2>> + DEAL_II_HOST_DEVICE void write_tensor(const Tensor &tensor, + const unsigned int i) const + requires writable; /** * Add a (scalar valued) @p entry to the vector at position @p i. - * Update the values of the @p n_components component vector at index @p i + * Update the values of the @p n_comp component vector at index @p i * by adding the values supplied by @p tensor. * * If the template parameter @a Number2 is a VectorizedArray then * the function takes a SIMD vectorized @p tensor as argument instead - * and updates the values of the @p n_components component vectors at + * and updates the values of the @p n_comp component vectors at * indices i, i+1, ..., i+simd_length_1. with the values supplied by @p * tensor. * - * @note This function is only available if `n_components` is equal to 1. + * @note This function is only available if `n_comp` is equal to 1. */ template - void add_entry(const Number2 &entry, const unsigned int i); + DEAL_II_HOST_DEVICE void add_entry(const Number2 &entry, + const unsigned int i) const + requires writable; /** - * Update the values of the @p n_components component vector at index @p i + * Update the values of the @p n_comp component vector at index @p i * by adding the values supplied by @p tensor. * * If the template parameter @a Number2 is a VectorizedArray then * the function takes a SIMD vectorized @p tensor as argument instead - * and updates the values of the @p n_components component vectors at + * and updates the values of the @p n_comp component vectors at * indices i, i+1, ..., i+simd_length_1. with the values supplied by @p * tensor. * @@ -286,27 +416,49 @@ namespace ryujin * Number, and has a type trait `value_type`. */ template > - void add_tensor(const Tensor &tensor, const unsigned int i); + typename Tensor = dealii::Tensor<1, n_comp, Number2>> + DEAL_II_HOST_DEVICE void add_tensor(const Tensor &tensor, + const unsigned int i) const + requires writable; //@} /** - * @name Vector interface + * MPI synchronization. */ //@{ - void sadd(const Number s, - const Number a, - const MultiComponentVector &V) - { - ScalarHostVector::sadd(s, a, V); - } + /** + * MPI synchronization: Zero out all ghost rows. + */ + void zero_out_ghost_values() const + requires(writable); - using ScalarHostVector::update_ghost_values; + /** + * MPI synchronization: Import all ghost values from neighboring MPI + * ranks on the templated memory space. + */ + void update_ghost_values() const + requires(writable); - using ScalarHostVector::zero_out_ghost_values; + /** + * MPI synchronization: Copy the data that has accumulated in the + * ghost range to the owning processor. This function operates on the + * templated memory space. + */ + void compress(dealii::VectorOperation::values operation) const + requires(writable); - using ScalarHostVector::compress; + private: + //@} + /** + * Internal data fields: + */ + //@{ + + MultiComponentVector + *multi_component_vector_; + + Number *data_; //@} }; @@ -320,54 +472,118 @@ namespace ryujin */ - template - void MultiComponentVector:: + template + MultiComponentVectorView:: + MultiComponentVectorView(MultiComponentVector + &multi_component_vector) + requires(writable) + { + reinit(multi_component_vector); + } + + + template + MultiComponentVectorView:: + MultiComponentVectorView( + const MultiComponentVector + &multi_component_vector) + requires(!writable) + { + reinit(multi_component_vector); + } + + + template + template + void + MultiComponentVectorView:: + reinit(MultiComponentVector &multi_component_vector) + requires(writable != std::is_const_v) + { + using HostSpace = dealii::MemorySpace::Host; + using DefaultSpace = dealii::MemorySpace::Default; + + static_assert(std::is_same_v || + std::is_same_v, + "Unexpected memory space"); + + multi_component_vector_ = &multi_component_vector; + + if constexpr (std::is_same_v) { + data_ = multi_component_vector_->data(); + } else { + data_ = multi_component_vector_->data(); + } + } + + +#if 0 + /* + * ------------------------------------------------------------------------- + * Inline function definitions + * ------------------------------------------------------------------------- + */ + + template + void MultiComponentVector:: reinit_with_vector_partitioner( const std::shared_ptr &vector_partitioner) { /* Special case of a zero component vector */ - if (n_components == 0) + if (n_comp == 0) return; ScalarHostVector::reinit(vector_partitioner); } - template - void MultiComponentVector:: + template + void MultiComponentVector:: reinit_with_scalar_partitioner( const std::shared_ptr &scalar_partitioner) { /* Special case of a zero component vector: */ - if (n_components == 0) + if (n_comp == 0) return; /* Special case of a scalar vector: */ - if (n_components == 1) + if (n_comp == 1) ScalarHostVector::reinit(scalar_partitioner); auto vector_partitioner = - create_vector_partitioner(scalar_partitioner, n_components); + create_vector_partitioner(scalar_partitioner, n_comp); ScalarHostVector::reinit(vector_partitioner); } - template + template template void - MultiComponentVector::extract_component( + MultiComponentVector::extract_component( ScalarHostVector &scalar_vector, unsigned int component, const Functor &functor) const { - Assert(n_components > 0, + Assert(n_comp > 0, dealii::ExcMessage( "Cannot extract from a vector with zero components.")); - AssertIndexRange(component, n_components); + AssertIndexRange(component, n_comp); - Assert(n_components * + Assert(n_comp * scalar_vector.get_partitioner()->locally_owned_size() == this->get_partitioner()->locally_owned_size(), dealii::ExcMessage("Called with a scalar_vector argument that has " @@ -376,25 +592,25 @@ namespace ryujin scalar_vector.get_partitioner()->locally_owned_size(); for (unsigned int i = 0; i < local_size; ++i) scalar_vector.local_element(i) = - functor(this->local_element(i * n_components + component)); + functor(this->local_element(i * n_comp + component)); scalar_vector.update_ghost_values(); } - template + template template void - MultiComponentVector::insert_component( + MultiComponentVector::insert_component( const ScalarHostVector &scalar_vector, unsigned int component, const Functor &functor) { - Assert(n_components > 0, + Assert(n_comp > 0, dealii::ExcMessage( "Cannot insert into a vector with zero components.")); - AssertIndexRange(component, n_components); + AssertIndexRange(component, n_comp); - Assert(n_components * + Assert(n_comp * scalar_vector.get_partitioner()->locally_owned_size() == this->get_partitioner()->locally_owned_size(), dealii::ExcMessage("Called with a scalar_vector argument that has " @@ -402,43 +618,43 @@ namespace ryujin const auto local_size = scalar_vector.get_partitioner()->locally_owned_size(); for (unsigned int i = 0; i < local_size; ++i) - this->local_element(i * n_components + component) = + this->local_element(i * n_comp + component) = functor(scalar_vector.local_element(i)); } - template + template template void - MultiComponentVector::insert_component( + MultiComponentVector::insert_component( const dealii::Vector &scalar_vector, unsigned int component, const Functor &functor) { - Assert(n_components > 0, + Assert(n_comp > 0, dealii::ExcMessage( "Cannot insert into a vector with zero components.")); - AssertIndexRange(component, n_components); + AssertIndexRange(component, n_comp); - Assert(n_components * scalar_vector.size() >= + Assert(n_comp * scalar_vector.size() >= this->get_partitioner()->locally_owned_size(), dealii::ExcMessage("Called with a scalar_vector argument that has " "incompatible local range.")); const auto local_size = scalar_vector.size(); for (unsigned int i = 0; i < local_size; ++i) - this->local_element(i * n_components + component) = + this->local_element(i * n_comp + component) = functor(scalar_vector[i]); } - template + template template DEAL_II_ALWAYS_INLINE inline Number2 - MultiComponentVector::read_entry( + MultiComponentVector::read_entry( const unsigned int i) const { static_assert( - n_components == 1, + n_comp == 1, "Attempted to read a scalar value from a tensor-valued vector entry"); AssertIndexRange(i, @@ -450,10 +666,10 @@ namespace ryujin } - template + template template DEAL_II_ALWAYS_INLINE inline Tensor - MultiComponentVector::read_tensor( + MultiComponentVector::read_tensor( const unsigned int i) const { static_assert(std::is_same_v, @@ -466,38 +682,38 @@ namespace ryujin Tensor tensor; /* Special case of a zero component vector */ - if constexpr (n_components == 0) + if constexpr (n_comp == 0) return tensor; if constexpr (std::is_same_v) { /* Vectorized fast access. index must be divisible by simd_length */ std::array indices; for (unsigned int k = 0; k < VectorizedArray::size(); ++k) - indices[k] = k * n_components; + indices[k] = k * n_comp; - dealii::vectorized_load_and_transpose(n_components, - this->begin() + i * n_components, + dealii::vectorized_load_and_transpose(n_comp, + this->begin() + i * n_comp, indices.data(), &tensor[0]); } else { /* Non-vectorized sequential access. */ - for (unsigned int d = 0; d < n_components; ++d) - tensor[d] = this->local_element(i * n_components + d); + for (unsigned int d = 0; d < n_comp; ++d) + tensor[d] = this->local_element(i * n_comp + d); } return tensor; } - template + template template DEAL_II_ALWAYS_INLINE inline Number2 - MultiComponentVector::read_entry( + MultiComponentVector::read_entry( const unsigned int *js) const { static_assert( - n_components == 1, + n_comp == 1, "Attempted to read a scalar value from a tensor-valued vector entry"); const auto result = read_tensor(js); @@ -505,10 +721,10 @@ namespace ryujin } - template + template template DEAL_II_ALWAYS_INLINE inline Tensor - MultiComponentVector::read_tensor( + MultiComponentVector::read_tensor( const unsigned int *js) const { static_assert(std::is_same_v, @@ -516,7 +732,7 @@ namespace ryujin Tensor tensor; /* Special case of a zero component vector */ - if constexpr (n_components == 0) + if constexpr (n_comp == 0) return tensor; if constexpr (std::is_same_v) { @@ -527,11 +743,11 @@ namespace ryujin AssertIndexRange(js[k], this->get_partitioner()->locally_owned_size() + this->get_partitioner()->n_ghost_indices()); - indices[k] = js[k] * n_components; + indices[k] = js[k] * n_comp; } dealii::vectorized_load_and_transpose( - n_components, this->begin(), indices.data(), &tensor[0]); + n_comp, this->begin(), indices.data(), &tensor[0]); } else { /* Non-vectorized sequential access. */ @@ -540,21 +756,21 @@ namespace ryujin this->get_partitioner()->locally_owned_size() + this->get_partitioner()->n_ghost_indices()); - for (unsigned int d = 0; d < n_components; ++d) - tensor[d] = this->local_element(js[0] * n_components + d); + for (unsigned int d = 0; d < n_comp; ++d) + tensor[d] = this->local_element(js[0] * n_comp + d); } return tensor; } - template + template template DEAL_II_ALWAYS_INLINE inline void - MultiComponentVector::write_entry( + MultiComponentVector::write_entry( const Number2 &entry, const unsigned int i) { - static_assert(n_components == 1, + static_assert(n_comp == 1, "Attempted to write a scalar value into a tensor-valued " "vector entry"); @@ -562,17 +778,17 @@ namespace ryujin this->get_partitioner()->locally_owned_size() + this->get_partitioner()->n_ghost_indices()); - dealii::Tensor<1, n_components, Number2> tensor; + dealii::Tensor<1, n_comp, Number2> tensor; tensor[0] = entry; write_tensor(tensor, i); } - template + template template DEAL_II_ALWAYS_INLINE inline void - MultiComponentVector::write_tensor( + MultiComponentVector::write_tensor( const Tensor &tensor, const unsigned int i) { static_assert(std::is_same_v, @@ -583,7 +799,7 @@ namespace ryujin this->get_partitioner()->n_ghost_indices()); /* Special case of a zero component vector */ - if constexpr (n_components == 0) + if constexpr (n_comp == 0) return; if constexpr (std::is_same_v) { @@ -591,31 +807,31 @@ namespace ryujin std::array indices; for (unsigned int k = 0; k < VectorizedArray::size(); ++k) - indices[k] = k * n_components; + indices[k] = k * n_comp; dealii::vectorized_transpose_and_store(/*add into*/ false, - n_components, + n_comp, &tensor[0], indices.data(), this->begin() + - i * n_components); + i * n_comp); } else { /* Non-vectorized sequential access. */ - for (unsigned int d = 0; d < n_components; ++d) - this->local_element(i * n_components + d) = tensor[d]; + for (unsigned int d = 0; d < n_comp; ++d) + this->local_element(i * n_comp + d) = tensor[d]; } } - template + template template DEAL_II_ALWAYS_INLINE inline void - MultiComponentVector::add_entry( + MultiComponentVector::add_entry( const Number2 &entry, const unsigned int i) { - static_assert(n_components == 1, + static_assert(n_comp == 1, "Attempted to write a scalar value into a tensor-valued " "matrix entry"); @@ -623,17 +839,17 @@ namespace ryujin this->get_partitioner()->locally_owned_size() + this->get_partitioner()->n_ghost_indices()); - dealii::Tensor<1, n_components, Number2> tensor; + dealii::Tensor<1, n_comp, Number2> tensor; tensor[0] = entry; add_tensor(tensor, i); } - template + template template DEAL_II_ALWAYS_INLINE inline void - MultiComponentVector::add_tensor( + MultiComponentVector::add_tensor( const Tensor &tensor, const unsigned int i) { static_assert(std::is_same_v, @@ -644,7 +860,7 @@ namespace ryujin this->get_partitioner()->n_ghost_indices()); /* Special case of a zero component vector */ - if constexpr (n_components == 0) + if constexpr (n_comp == 0) return; if constexpr (std::is_same_v) { @@ -652,23 +868,102 @@ namespace ryujin std::array indices; for (unsigned int k = 0; k < VectorizedArray::size(); ++k) - indices[k] = k * n_components; + indices[k] = k * n_comp; dealii::vectorized_transpose_and_store(/*add into*/ true, - n_components, + n_comp, &tensor[0], indices.data(), this->begin() + - i * n_components); + i * n_comp); } else { /* Non-vectorized sequential access. */ - for (unsigned int d = 0; d < n_components; ++d) - this->local_element(i * n_components + d) += tensor[d]; + for (unsigned int d = 0; d < n_comp; ++d) + this->local_element(i * n_comp + d) += tensor[d]; } } +#endif + + template + void + MultiComponentVectorView:: + zero_out_ghost_values() const + requires(writable) + { + using HostSpace = dealii::MemorySpace::Host; + using DefaultSpace = dealii::MemorySpace::Default; + + static_assert(std::is_same_v || + std::is_same_v, + "Unexpected memory space"); + + Assert(multi_component_vector_ + ->template is_active_memory_space(), + dealii::ExcMessage("The chosen memory space is not active.")); + + multi_component_vector_ + ->template zero_out_ghost_values_on_memory_space(); + } + + + template + void + MultiComponentVectorView:: + update_ghost_values() const + requires(writable) + { + using HostSpace = dealii::MemorySpace::Host; + using DefaultSpace = dealii::MemorySpace::Default; + + static_assert(std::is_same_v || + std::is_same_v, + "Unexpected memory space"); + + Assert(multi_component_vector_ + ->template is_active_memory_space(), + dealii::ExcMessage("The chosen memory space is not active.")); + + multi_component_vector_ + ->template update_ghost_values_on_memory_space(); + } + + + template + void + MultiComponentVectorView:: + compress(dealii::VectorOperation::values operation) const + requires(writable) + { + using HostSpace = dealii::MemorySpace::Host; + using DefaultSpace = dealii::MemorySpace::Default; + + static_assert(std::is_same_v || + std::is_same_v, + "Unexpected memory space"); + + Assert(multi_component_vector_ + ->template is_active_memory_space(), + dealii::ExcMessage("The chosen memory space is not active.")); + + multi_component_vector_->template compress_on_memory_space( + operation); + } + #endif } // namespace Vectors } // namespace ryujin From f45c248a1a1b5cefacc1b395865cc9107e661153 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Fri, 20 Mar 2026 20:37:06 -0500 Subject: [PATCH 03/10] MultiComponentVector: implement cons, operator=, reinit() --- source/multicomponent_vector.h | 141 +++++++++++++++++++++++---------- 1 file changed, 99 insertions(+), 42 deletions(-) diff --git a/source/multicomponent_vector.h b/source/multicomponent_vector.h index 57071ab80..4806785c6 100644 --- a/source/multicomponent_vector.h +++ b/source/multicomponent_vector.h @@ -71,6 +71,10 @@ namespace ryujin MultiComponentVector() = default; + MultiComponentVector(const MultiComponentVector &other); + + MultiComponentVector(MultiComponentVector &&other) noexcept; + /** * Reinitializes the MultiComponentVector with a vector MPI partitioner * that was created first with create_vector_partitioner(). @@ -88,6 +92,12 @@ namespace ryujin const std::shared_ptr &scalar_partitioner); + MultiComponentVector & + operator=(const MultiComponentVector &other); + + MultiComponentVector & + operator=(MultiComponentVector &&other) noexcept; + //@} /** * Memory space access and synchronization: @@ -472,6 +482,93 @@ namespace ryujin */ + template + MultiComponentVector::MultiComponentVector( + const MultiComponentVector &other) + { + *this = other; + } + + + template + MultiComponentVector::MultiComponentVector( + MultiComponentVector &&other) noexcept + { + *this = other; + } + + + template + void MultiComponentVector:: + reinit_with_vector_partitioner( + const std::shared_ptr + &vector_partitioner) + { + /* Special case of a zero component vector */ + if (n_comp == 0) + return; + + host_vector_.reinit(vector_partitioner); + + /* Reinitialize view to point to the correct vector data: */ + MultiComponentVectorView::reinit(*this); + } + + + template + void MultiComponentVector:: + reinit_with_scalar_partitioner( + const std::shared_ptr + &scalar_partitioner) + { + /* Special case of a zero component vector: */ + if (n_comp == 0) + return; + + /* Special case of a scalar vector: */ + if (n_comp == 1) + host_vector_.reinit(scalar_partitioner); + + auto vector_partitioner = + create_vector_partitioner(scalar_partitioner, n_comp); + + host_vector_.reinit(vector_partitioner); + + /* Reinitialize view to point to the correct vector data: */ + MultiComponentVectorView::reinit(*this); + } + + + template + auto MultiComponentVector::operator=( + const MultiComponentVector &other) -> MultiComponentVector & + { + host_vector_ = other.host_vector_; + default_vector_ = other.default_vector_; + host_space_active_ = other.host_space_active_; + + /* Reinitialize view to point to the correct vector data: */ + MultiComponentVectorView::reinit(*this); + + return *this; + } + + + template + auto MultiComponentVector::operator=( + MultiComponentVector &&other) noexcept -> MultiComponentVector & + { + host_vector_ = std::move(other.host_vector_); + default_vector_ = std::move(other.default_vector_); + host_space_active_ = other.host_space_active_; + + /* Reinitialize view to point to the correct vector data: */ + MultiComponentVectorView::reinit(*this); + + return *this; + } + + template ) { - data_ = multi_component_vector_->data(); + data_ = multi_component_vector_->host_vector_.begin(); } else { - data_ = multi_component_vector_->data(); + data_ = multi_component_vector_->default_vector_.begin(); } } #if 0 - /* - * ------------------------------------------------------------------------- - * Inline function definitions - * ------------------------------------------------------------------------- - */ - - template - void MultiComponentVector:: - reinit_with_vector_partitioner( - const std::shared_ptr - &vector_partitioner) - { - /* Special case of a zero component vector */ - if (n_comp == 0) - return; - - ScalarHostVector::reinit(vector_partitioner); - } - - template - void MultiComponentVector:: - reinit_with_scalar_partitioner( - const std::shared_ptr - &scalar_partitioner) - { - /* Special case of a zero component vector: */ - if (n_comp == 0) - return; - - /* Special case of a scalar vector: */ - if (n_comp == 1) - ScalarHostVector::reinit(scalar_partitioner); - - auto vector_partitioner = - create_vector_partitioner(scalar_partitioner, n_comp); - - ScalarHostVector::reinit(vector_partitioner); - } - - template template void From ff6fb3f026029739b653d557338a4bbbc79928a7 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Fri, 20 Mar 2026 22:55:34 -0500 Subject: [PATCH 04/10] MultiComponentVector: implement MemorySpace management --- source/multicomponent_vector.h | 155 +++++++++++++++++++++++++++++++-- 1 file changed, 149 insertions(+), 6 deletions(-) diff --git a/source/multicomponent_vector.h b/source/multicomponent_vector.h index 4806785c6..66322f633 100644 --- a/source/multicomponent_vector.h +++ b/source/multicomponent_vector.h @@ -92,11 +92,9 @@ namespace ryujin const std::shared_ptr &scalar_partitioner); - MultiComponentVector & - operator=(const MultiComponentVector &other); + MultiComponentVector &operator=(const MultiComponentVector &other); - MultiComponentVector & - operator=(MultiComponentVector &&other) noexcept; + MultiComponentVector &operator=(MultiComponentVector &&other) noexcept; //@} /** @@ -465,8 +463,8 @@ namespace ryujin */ //@{ - MultiComponentVector - *multi_component_vector_; + using MCV = MultiComponentVector; + std::conditional_t multi_component_vector_; Number *data_; @@ -510,6 +508,11 @@ namespace ryujin host_vector_.reinit(vector_partitioner); + using KokkosHost = dealii::MemorySpace::Host::kokkos_space; + using KokkosDefault = dealii::MemorySpace::Default::kokkos_space; + if (!std::is_same_v) + default_vector_.reinit(vector_partitioner); + /* Reinitialize view to point to the correct vector data: */ MultiComponentVectorView::reinit(*this); } @@ -534,6 +537,11 @@ namespace ryujin host_vector_.reinit(vector_partitioner); + using KokkosHost = dealii::MemorySpace::Host::kokkos_space; + using KokkosDefault = dealii::MemorySpace::Default::kokkos_space; + if (!std::is_same_v) + default_vector_.reinit(vector_partitioner); + /* Reinitialize view to point to the correct vector data: */ MultiComponentVectorView::reinit(*this); } @@ -569,6 +577,141 @@ namespace ryujin } + template + template + MultiComponentVectorView + MultiComponentVector::get_view() + { + Assert(is_active_memory_space(), + dealii::ExcMessage("The chosen memory space is not active.")); + + return MultiComponentVectorView(*this); + } + + + template + template + MultiComponentVectorView + MultiComponentVector::get_view() const + { + Assert(is_active_memory_space(), + dealii::ExcMessage("The chosen memory space is not active.")); + + return MultiComponentVectorView(*this); + } + + + template + template + bool + MultiComponentVector::is_active_memory_space() + const + { + using HostSpace = dealii::MemorySpace::Host; + using DefaultSpace = dealii::MemorySpace::Default; + static_assert(std::is_same_v || + std::is_same_v, + "Unexpected Kokkos memory space"); + + return host_space_active_ == std::is_same_v; + } + + + template + template + void + MultiComponentVector::move_to_memory_space() + { + using HostSpace = dealii::MemorySpace::Host::kokkos_space; + using DefaultSpace = dealii::MemorySpace::Default::kokkos_space; + static_assert(std::is_same_v || + std::is_same_v, + "Unexpected Kokkos memory space"); + + if (is_active_memory_space()) + return; + + if constexpr (std::is_same_v) { + host_vector_.import_elements(default_vector_, + dealii::VectorOperation::insert); + host_space_active_ = true; + + } else if constexpr (std::is_same_v) { + default_vector_.import_elements(host_vector_, + dealii::VectorOperation::insert); + host_space_active_ = false; + } + } + + + template + template + void MultiComponentVector:: + zero_out_ghost_values_on_memory_space() + { + using HostSpace = dealii::MemorySpace::Host; + using DefaultSpace = dealii::MemorySpace::Default; + + static_assert(std::is_same_v || + std::is_same_v, + "Unexpected memory space"); + + if constexpr (std::is_same_v) { + host_vector_.zero_out_ghost_values(); + } else { + default_vector_.zero_out_ghost_values(); + } + } + + + template + template + void MultiComponentVector:: + update_ghost_values_on_memory_space() + { + using HostSpace = dealii::MemorySpace::Host; + using DefaultSpace = dealii::MemorySpace::Default; + + static_assert(std::is_same_v || + std::is_same_v, + "Unexpected memory space"); + + if constexpr (std::is_same_v) { + host_vector_.update_ghost_values(); + } else { + default_vector_.update_ghost_values(); + } + } + + + template + template + void MultiComponentVector:: + compress_on_memory_space(dealii::VectorOperation::values operation) + { + using HostSpace = dealii::MemorySpace::Host; + using DefaultSpace = dealii::MemorySpace::Default; + + static_assert(std::is_same_v || + std::is_same_v, + "Unexpected memory space"); + + if constexpr (std::is_same_v) { + host_vector_.compress(operation); + } else { + default_vector_.compress(operation); + } + } + + template Date: Fri, 20 Mar 2026 23:26:57 -0500 Subject: [PATCH 05/10] MultiComponentVector: implement read_* and write_* functions --- source/multicomponent_vector.h | 235 +++++++++++++++++++++------------ 1 file changed, 149 insertions(+), 86 deletions(-) diff --git a/source/multicomponent_vector.h b/source/multicomponent_vector.h index 66322f633..ac4ad7623 100644 --- a/source/multicomponent_vector.h +++ b/source/multicomponent_vector.h @@ -467,7 +467,8 @@ namespace ryujin std::conditional_t multi_component_vector_; Number *data_; - + unsigned int n_locally_owned_; + unsigned int n_locally_relevant_; //@} }; @@ -762,9 +763,21 @@ namespace ryujin multi_component_vector_ = &multi_component_vector; if constexpr (std::is_same_v) { - data_ = multi_component_vector_->host_vector_.begin(); + auto &vector = multi_component_vector_->host_vector_; + const auto &partitioner = vector.get_partitioner(); + + data_ = vector.begin(); + n_locally_owned_ = partitioner->locally_owned_size(); + n_locally_relevant_ = n_locally_owned_ + partitioner->n_ghost_indices(); + } else { - data_ = multi_component_vector_->default_vector_.begin(); + + auto &vector = multi_component_vector_->default_vector_; + const auto &partitioner = vector.get_partitioner(); + + data_ = vector.begin(); + n_locally_owned_ = partitioner->locally_owned_size(); + n_locally_relevant_ = n_locally_owned_ + partitioner->n_ghost_indices(); } } @@ -845,39 +858,50 @@ namespace ryujin this->local_element(i * n_comp + component) = functor(scalar_vector[i]); } +#endif - template + template template - DEAL_II_ALWAYS_INLINE inline Number2 - MultiComponentVector::read_entry( - const unsigned int i) const + DEAL_II_HOST_DEVICE_ALWAYS_INLINE Number2 + MultiComponentVectorView::read_entry(const unsigned int i) const { static_assert( n_comp == 1, "Attempted to read a scalar value from a tensor-valued vector entry"); - AssertIndexRange(i, - this->get_partitioner()->locally_owned_size() + - this->get_partitioner()->n_ghost_indices()); + AssertIndexRange(i, n_locally_relevant_); const auto result = read_tensor(i); return result[0]; } - template + template template - DEAL_II_ALWAYS_INLINE inline Tensor - MultiComponentVector::read_tensor( - const unsigned int i) const + DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor + MultiComponentVectorView::read_tensor(const unsigned int i) const { static_assert(std::is_same_v, "type mismatch"); - AssertIndexRange(i, - this->get_partitioner()->locally_owned_size() + - this->get_partitioner()->n_ghost_indices()); + AssertIndexRange(i, n_locally_relevant_); Tensor tensor; @@ -885,32 +909,38 @@ namespace ryujin if constexpr (n_comp == 0) return tensor; - if constexpr (std::is_same_v) { + using VA = dealii::VectorizedArray; + if constexpr (std::is_same_v) { /* Vectorized fast access. index must be divisible by simd_length */ - std::array indices; - for (unsigned int k = 0; k < VectorizedArray::size(); ++k) + std::array indices; + for (unsigned int k = 0; k < VA::size(); ++k) indices[k] = k * n_comp; - dealii::vectorized_load_and_transpose(n_comp, - this->begin() + i * n_comp, - indices.data(), - &tensor[0]); + dealii::vectorized_load_and_transpose( + n_comp, data_ + i * n_comp, indices.data(), &tensor[0]); + } else { /* Non-vectorized sequential access. */ - for (unsigned int d = 0; d < n_comp; ++d) - tensor[d] = this->local_element(i * n_comp + d); + tensor[d] = data_[i * n_comp + d]; } return tensor; } - template + template template - DEAL_II_ALWAYS_INLINE inline Number2 - MultiComponentVector::read_entry( - const unsigned int *js) const + DEAL_II_HOST_DEVICE_ALWAYS_INLINE Number2 + MultiComponentVectorView::read_entry(const unsigned int *js) const { static_assert( n_comp == 1, @@ -921,62 +951,75 @@ namespace ryujin } - template + template template - DEAL_II_ALWAYS_INLINE inline Tensor - MultiComponentVector::read_tensor( - const unsigned int *js) const + DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor + MultiComponentVectorView::read_tensor(const unsigned int *js) + const { static_assert(std::is_same_v, "type mismatch"); + Tensor tensor; /* Special case of a zero component vector */ if constexpr (n_comp == 0) return tensor; - if constexpr (std::is_same_v) { + using VA = dealii::VectorizedArray; + if constexpr (std::is_same_v) { /* Vectorized fast access. index must be divisible by simd_length */ - std::array indices; - for (unsigned int k = 0; k < VectorizedArray::size(); ++k) { - AssertIndexRange(js[k], - this->get_partitioner()->locally_owned_size() + - this->get_partitioner()->n_ghost_indices()); + std::array indices; + for (unsigned int k = 0; k < VA::size(); ++k) { + AssertIndexRange(js[k], n_locally_relevant_); indices[k] = js[k] * n_comp; } dealii::vectorized_load_and_transpose( - n_comp, this->begin(), indices.data(), &tensor[0]); + n_comp, data_, indices.data(), &tensor[0]); } else { /* Non-vectorized sequential access. */ - AssertIndexRange(*js, - this->get_partitioner()->locally_owned_size() + - this->get_partitioner()->n_ghost_indices()); + AssertIndexRange(*js, n_locally_relevant_); for (unsigned int d = 0; d < n_comp; ++d) - tensor[d] = this->local_element(js[0] * n_comp + d); + tensor[d] = data_[js[0] * n_comp + d]; } return tensor; } - template + template template - DEAL_II_ALWAYS_INLINE inline void - MultiComponentVector::write_entry( - const Number2 &entry, const unsigned int i) + DEAL_II_HOST_DEVICE_ALWAYS_INLINE void + MultiComponentVectorView::write_entry(const Number2 &entry, + const unsigned int i) const + requires writable { static_assert(n_comp == 1, "Attempted to write a scalar value into a tensor-valued " "vector entry"); - AssertIndexRange(i, - this->get_partitioner()->locally_owned_size() + - this->get_partitioner()->n_ghost_indices()); + AssertIndexRange(i, n_locally_relevant_); dealii::Tensor<1, n_comp, Number2> tensor; tensor[0] = entry; @@ -985,59 +1028,73 @@ namespace ryujin } - template + template template - DEAL_II_ALWAYS_INLINE inline void - MultiComponentVector::write_tensor( - const Tensor &tensor, const unsigned int i) + DEAL_II_HOST_DEVICE_ALWAYS_INLINE void + MultiComponentVectorView::write_tensor(const Tensor &tensor, + const unsigned int i) const + requires writable { static_assert(std::is_same_v, "type mismatch"); - AssertIndexRange(i, - this->get_partitioner()->locally_owned_size() + - this->get_partitioner()->n_ghost_indices()); + AssertIndexRange(i, n_locally_relevant_); /* Special case of a zero component vector */ if constexpr (n_comp == 0) return; - if constexpr (std::is_same_v) { + using VA = dealii::VectorizedArray; + if constexpr (std::is_same_v) { /* Vectorized fast access. index must be divisible by simd_length */ - std::array indices; - for (unsigned int k = 0; k < VectorizedArray::size(); ++k) + std::array indices; + for (unsigned int k = 0; k < VA::size(); ++k) indices[k] = k * n_comp; dealii::vectorized_transpose_and_store(/*add into*/ false, n_comp, &tensor[0], indices.data(), - this->begin() + - i * n_comp); + data_ + i * n_comp); } else { /* Non-vectorized sequential access. */ for (unsigned int d = 0; d < n_comp; ++d) - this->local_element(i * n_comp + d) = tensor[d]; + data_[i * n_comp + d] = tensor[d]; } } - template + template template - DEAL_II_ALWAYS_INLINE inline void - MultiComponentVector::add_entry( - const Number2 &entry, const unsigned int i) + DEAL_II_HOST_DEVICE_ALWAYS_INLINE void + MultiComponentVectorView::add_entry(const Number2 &entry, + const unsigned int i) const + requires writable { static_assert(n_comp == 1, "Attempted to write a scalar value into a tensor-valued " "matrix entry"); - AssertIndexRange(i, - this->get_partitioner()->locally_owned_size() + - this->get_partitioner()->n_ghost_indices()); + AssertIndexRange(i, n_locally_relevant_); dealii::Tensor<1, n_comp, Number2> tensor; tensor[0] = entry; @@ -1046,46 +1103,52 @@ namespace ryujin } - template + template template - DEAL_II_ALWAYS_INLINE inline void - MultiComponentVector::add_tensor( - const Tensor &tensor, const unsigned int i) + DEAL_II_HOST_DEVICE_ALWAYS_INLINE void + MultiComponentVectorView::add_tensor(const Tensor &tensor, + const unsigned int i) const + requires writable { static_assert(std::is_same_v, "type mismatch"); - AssertIndexRange(i, - this->get_partitioner()->locally_owned_size() + - this->get_partitioner()->n_ghost_indices()); + AssertIndexRange(i, n_locally_relevant_); /* Special case of a zero component vector */ if constexpr (n_comp == 0) return; - if constexpr (std::is_same_v) { + using VA = dealii::VectorizedArray; + if constexpr (std::is_same_v) { /* Vectorized fast access. index must be divisible by simd_length */ - std::array indices; - for (unsigned int k = 0; k < VectorizedArray::size(); ++k) + std::array indices; + for (unsigned int k = 0; k < VA::size(); ++k) indices[k] = k * n_comp; dealii::vectorized_transpose_and_store(/*add into*/ true, n_comp, &tensor[0], indices.data(), - this->begin() + - i * n_comp); + data_ + i * n_comp); } else { /* Non-vectorized sequential access. */ for (unsigned int d = 0; d < n_comp; ++d) - this->local_element(i * n_comp + d) += tensor[d]; + data_[i * n_comp + d] += tensor[d]; } } -#endif template Date: Fri, 20 Mar 2026 23:54:30 -0500 Subject: [PATCH 06/10] MultiComponentVector: implement insert_component / extract_component --- source/multicomponent_vector.h | 113 +++++++++++++++++++++------------ 1 file changed, 72 insertions(+), 41 deletions(-) diff --git a/source/multicomponent_vector.h b/source/multicomponent_vector.h index ac4ad7623..313c0588b 100644 --- a/source/multicomponent_vector.h +++ b/source/multicomponent_vector.h @@ -276,7 +276,7 @@ namespace ryujin template void insert_component(const ScalarVector &scalar_vector, unsigned int component, - const Functor &functor = std::identity{}) + const Functor &functor = std::identity{}) const requires writable; /** @@ -286,7 +286,7 @@ namespace ryujin template void insert_component(const dealii::Vector &scalar_vector, unsigned int component, - const Functor &functor = std::identity{}) + const Functor &functor = std::identity{}) const requires writable; /** @@ -782,83 +782,114 @@ namespace ryujin } -#if 0 - template + template template - void - MultiComponentVector::extract_component( - ScalarHostVector &scalar_vector, - unsigned int component, - const Functor &functor) const + void MultiComponentVectorView< + Number, + n_comp, + simd_length, + MemorySpace, + writable>::extract_component(ScalarVector &scalar_vector, + unsigned int component, + const Functor &functor) const { + using HostSpace = dealii::MemorySpace::Host; + AssertThrow((std::is_same_v), + dealii::ExcNotImplemented()); + Assert(n_comp > 0, dealii::ExcMessage( "Cannot extract from a vector with zero components.")); AssertIndexRange(component, n_comp); - Assert(n_comp * - scalar_vector.get_partitioner()->locally_owned_size() == - this->get_partitioner()->locally_owned_size(), - dealii::ExcMessage("Called with a scalar_vector argument that has " - "incompatible local range.")); const auto local_size = scalar_vector.get_partitioner()->locally_owned_size(); + + Assert(n_comp * local_size == n_locally_owned_, + dealii::ExcMessage("Called with a scalar_vector argument that has " + "incompatible local range.")); + for (unsigned int i = 0; i < local_size; ++i) - scalar_vector.local_element(i) = - functor(this->local_element(i * n_comp + component)); + scalar_vector.local_element(i) = functor(data_[i * n_comp + component]); scalar_vector.update_ghost_values(); } - template + template template - void - MultiComponentVector::insert_component( - const ScalarHostVector &scalar_vector, - unsigned int component, - const Functor &functor) + void MultiComponentVectorView< + Number, + n_comp, + simd_length, + MemorySpace, + writable>::insert_component(const ScalarVector &scalar_vector, + unsigned int component, + const Functor &functor) const + requires writable { + using HostSpace = dealii::MemorySpace::Host; + AssertThrow((std::is_same_v), + dealii::ExcNotImplemented()); + Assert(n_comp > 0, dealii::ExcMessage( "Cannot insert into a vector with zero components.")); AssertIndexRange(component, n_comp); - Assert(n_comp * - scalar_vector.get_partitioner()->locally_owned_size() == - this->get_partitioner()->locally_owned_size(), - dealii::ExcMessage("Called with a scalar_vector argument that has " - "incompatible local range.")); const auto local_size = scalar_vector.get_partitioner()->locally_owned_size(); + + Assert(n_comp * local_size == n_locally_owned_, + dealii::ExcMessage("Called with a scalar_vector argument that has " + "incompatible local range.")); + for (unsigned int i = 0; i < local_size; ++i) - this->local_element(i * n_comp + component) = - functor(scalar_vector.local_element(i)); + data_[i * n_comp + component] = functor(scalar_vector.local_element(i)); } - template + template template - void - MultiComponentVector::insert_component( - const dealii::Vector &scalar_vector, - unsigned int component, - const Functor &functor) + void MultiComponentVectorView< + Number, + n_comp, + simd_length, + MemorySpace, + writable>::insert_component(const dealii::Vector &scalar_vector, + unsigned int component, + const Functor &functor) const + requires writable { + using HostSpace = dealii::MemorySpace::Host; + AssertThrow((std::is_same_v), + dealii::ExcInternalError()); + Assert(n_comp > 0, dealii::ExcMessage( "Cannot insert into a vector with zero components.")); AssertIndexRange(component, n_comp); - Assert(n_comp * scalar_vector.size() >= - this->get_partitioner()->locally_owned_size(), + const auto local_size = scalar_vector.size(); + + Assert(n_comp * local_size >= n_locally_owned_, dealii::ExcMessage("Called with a scalar_vector argument that has " "incompatible local range.")); - const auto local_size = scalar_vector.size(); + for (unsigned int i = 0; i < local_size; ++i) - this->local_element(i * n_comp + component) = - functor(scalar_vector[i]); + data_[i * n_comp + component] = functor(scalar_vector[i]); } -#endif template Date: Sat, 21 Mar 2026 00:14:02 -0500 Subject: [PATCH 07/10] MultiComponentVector: implement sadd() --- source/multicomponent_vector.h | 31 ++++++++++++++++++++++++++----- 1 file changed, 26 insertions(+), 5 deletions(-) diff --git a/source/multicomponent_vector.h b/source/multicomponent_vector.h index 313c0588b..0b5d0221a 100644 --- a/source/multicomponent_vector.h +++ b/source/multicomponent_vector.h @@ -293,8 +293,9 @@ namespace ryujin * Scaled addition of the given vector $U$ and the argument vector * $V$: $U\leftarrow s\,U+a\,V$. */ - void - sadd(const Number s, const Number a, const MultiComponentVectorView &V) + void sadd(const Number s, + const Number a, + const MultiComponentVectorView &v) const requires writable; //@} @@ -631,11 +632,11 @@ namespace ryujin void MultiComponentVector::move_to_memory_space() { - using HostSpace = dealii::MemorySpace::Host::kokkos_space; - using DefaultSpace = dealii::MemorySpace::Default::kokkos_space; + using HostSpace = dealii::MemorySpace::Host; + using DefaultSpace = dealii::MemorySpace::Default; static_assert(std::is_same_v || std::is_same_v, - "Unexpected Kokkos memory space"); + "Unexpected memory space"); if (is_active_memory_space()) return; @@ -892,6 +893,26 @@ namespace ryujin } + template + void MultiComponentVectorView::sadd( + const Num s, const Num a, const MultiComponentVectorView &v) const + requires writable + { + using HS = dealii::MemorySpace::Host; + using DS = dealii::MemorySpace::Default; + static_assert(std::is_same_v || std::is_same_v, + "Unexpected memory space"); + + if constexpr (std::is_same_v) { + const auto &other = v.multi_component_vector_->host_vector_; + multi_component_vector_->host_vector_.sadd(s, a, other); + } else if constexpr (std::is_same_v) { + const auto &other = v.multi_component_vector_->default_vector_; + multi_component_vector_->default_vector_.sadd(s, a, other); + } + } + + template Date: Sat, 21 Mar 2026 13:36:31 -0500 Subject: [PATCH 08/10] MultiComponentVector: silence gcc warning --- source/multicomponent_vector.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/source/multicomponent_vector.h b/source/multicomponent_vector.h index 0b5d0221a..6213eabd9 100644 --- a/source/multicomponent_vector.h +++ b/source/multicomponent_vector.h @@ -485,6 +485,7 @@ namespace ryujin template MultiComponentVector::MultiComponentVector( const MultiComponentVector &other) + : MultiComponentVectorView() { *this = other; } @@ -493,6 +494,7 @@ namespace ryujin template MultiComponentVector::MultiComponentVector( MultiComponentVector &&other) noexcept + : MultiComponentVectorView() { *this = other; } From d015324d5703dcca3bc4e073deb881a1fc848cd5 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Sat, 21 Mar 2026 19:33:58 -0500 Subject: [PATCH 09/10] MultiComponentVector: avoid using default_vector_ when host and default are the same spaces --- source/multicomponent_vector.h | 76 ++++++++++++++++++---------------- 1 file changed, 41 insertions(+), 35 deletions(-) diff --git a/source/multicomponent_vector.h b/source/multicomponent_vector.h index 6213eabd9..585e5ad38 100644 --- a/source/multicomponent_vector.h +++ b/source/multicomponent_vector.h @@ -172,6 +172,14 @@ namespace ryujin dealii::MemorySpace::Host> host_vector_; + /* + * We avoid setting up the default_vector_ if default and host happen + * to be the same memory space. + */ + static constexpr bool have_separate_memory_spaces = + !std::is_same_v; + dealii::LinearAlgebra::distributed::Vector default_vector_; @@ -511,10 +519,7 @@ namespace ryujin return; host_vector_.reinit(vector_partitioner); - - using KokkosHost = dealii::MemorySpace::Host::kokkos_space; - using KokkosDefault = dealii::MemorySpace::Default::kokkos_space; - if (!std::is_same_v) + if constexpr (have_separate_memory_spaces) default_vector_.reinit(vector_partitioner); /* Reinitialize view to point to the correct vector data: */ @@ -540,10 +545,7 @@ namespace ryujin create_vector_partitioner(scalar_partitioner, n_comp); host_vector_.reinit(vector_partitioner); - - using KokkosHost = dealii::MemorySpace::Host::kokkos_space; - using KokkosDefault = dealii::MemorySpace::Default::kokkos_space; - if (!std::is_same_v) + if constexpr (have_separate_memory_spaces) default_vector_.reinit(vector_partitioner); /* Reinitialize view to point to the correct vector data: */ @@ -556,7 +558,8 @@ namespace ryujin const MultiComponentVector &other) -> MultiComponentVector & { host_vector_ = other.host_vector_; - default_vector_ = other.default_vector_; + if constexpr (have_separate_memory_spaces) + default_vector_ = other.default_vector_; host_space_active_ = other.host_space_active_; /* Reinitialize view to point to the correct vector data: */ @@ -571,7 +574,8 @@ namespace ryujin MultiComponentVector &&other) noexcept -> MultiComponentVector & { host_vector_ = std::move(other.host_vector_); - default_vector_ = std::move(other.default_vector_); + if constexpr (have_separate_memory_spaces) + default_vector_ = std::move(other.default_vector_); host_space_active_ = other.host_space_active_; /* Reinitialize view to point to the correct vector data: */ @@ -623,7 +627,7 @@ namespace ryujin using DefaultSpace = dealii::MemorySpace::Default; static_assert(std::is_same_v || std::is_same_v, - "Unexpected Kokkos memory space"); + "Unexpected memory space"); return host_space_active_ == std::is_same_v; } @@ -644,13 +648,15 @@ namespace ryujin return; if constexpr (std::is_same_v) { - host_vector_.import_elements(default_vector_, - dealii::VectorOperation::insert); + if constexpr (have_separate_memory_spaces) + host_vector_.import_elements(default_vector_, + dealii::VectorOperation::insert); host_space_active_ = true; } else if constexpr (std::is_same_v) { - default_vector_.import_elements(host_vector_, - dealii::VectorOperation::insert); + if constexpr (have_separate_memory_spaces) + default_vector_.import_elements(host_vector_, + dealii::VectorOperation::insert); host_space_active_ = false; } } @@ -663,15 +669,15 @@ namespace ryujin { using HostSpace = dealii::MemorySpace::Host; using DefaultSpace = dealii::MemorySpace::Default; - static_assert(std::is_same_v || std::is_same_v, "Unexpected memory space"); - if constexpr (std::is_same_v) { - host_vector_.zero_out_ghost_values(); - } else { + if constexpr (have_separate_memory_spaces && + !std::is_same_v) { default_vector_.zero_out_ghost_values(); + } else { + host_vector_.zero_out_ghost_values(); } } @@ -683,15 +689,15 @@ namespace ryujin { using HostSpace = dealii::MemorySpace::Host; using DefaultSpace = dealii::MemorySpace::Default; - static_assert(std::is_same_v || std::is_same_v, "Unexpected memory space"); - if constexpr (std::is_same_v) { - host_vector_.update_ghost_values(); - } else { + if constexpr (have_separate_memory_spaces && + !std::is_same_v) { default_vector_.update_ghost_values(); + } else { + host_vector_.update_ghost_values(); } } @@ -703,15 +709,15 @@ namespace ryujin { using HostSpace = dealii::MemorySpace::Host; using DefaultSpace = dealii::MemorySpace::Default; - static_assert(std::is_same_v || std::is_same_v, "Unexpected memory space"); - if constexpr (std::is_same_v) { - host_vector_.compress(operation); - } else { + if constexpr (have_separate_memory_spaces && + !std::is_same_v) { default_vector_.compress(operation); + } else { + host_vector_.compress(operation); } } @@ -758,15 +764,19 @@ namespace ryujin { using HostSpace = dealii::MemorySpace::Host; using DefaultSpace = dealii::MemorySpace::Default; - static_assert(std::is_same_v || std::is_same_v, "Unexpected memory space"); multi_component_vector_ = &multi_component_vector; - if constexpr (std::is_same_v) { - auto &vector = multi_component_vector_->host_vector_; + constexpr bool have_separate_memory_spaces = + ::ryujin::Vectors::MultiComponentVector:: + have_separate_memory_spaces; + + if constexpr (have_separate_memory_spaces && + !std::is_same_v) { + auto &vector = multi_component_vector_->default_vector_; const auto &partitioner = vector.get_partitioner(); data_ = vector.begin(); @@ -774,8 +784,7 @@ namespace ryujin n_locally_relevant_ = n_locally_owned_ + partitioner->n_ghost_indices(); } else { - - auto &vector = multi_component_vector_->default_vector_; + auto &vector = multi_component_vector_->host_vector_; const auto &partitioner = vector.get_partitioner(); data_ = vector.begin(); @@ -1216,7 +1225,6 @@ namespace ryujin { using HostSpace = dealii::MemorySpace::Host; using DefaultSpace = dealii::MemorySpace::Default; - static_assert(std::is_same_v || std::is_same_v, "Unexpected memory space"); @@ -1242,7 +1250,6 @@ namespace ryujin { using HostSpace = dealii::MemorySpace::Host; using DefaultSpace = dealii::MemorySpace::Default; - static_assert(std::is_same_v || std::is_same_v, "Unexpected memory space"); @@ -1268,7 +1275,6 @@ namespace ryujin { using HostSpace = dealii::MemorySpace::Host; using DefaultSpace = dealii::MemorySpace::Default; - static_assert(std::is_same_v || std::is_same_v, "Unexpected memory space"); From b805418ea203f540b9be98b146cabb2ffb12f4a9 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Sat, 21 Mar 2026 18:35:50 -0500 Subject: [PATCH 10/10] MultiComponentVector: add a test --- tests/common/multicomponent_vector_03.cc | 77 +++++++++++++++++++ .../multicomponent_vector_03.threads=2.output | 31 ++++++++ 2 files changed, 108 insertions(+) create mode 100644 tests/common/multicomponent_vector_03.cc create mode 100644 tests/common/multicomponent_vector_03.threads=2.output diff --git a/tests/common/multicomponent_vector_03.cc b/tests/common/multicomponent_vector_03.cc new file mode 100644 index 000000000..b34642b61 --- /dev/null +++ b/tests/common/multicomponent_vector_03.cc @@ -0,0 +1,77 @@ +#include + +int main(int argc, char *argv[]) +{ + dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); + + /* Set up locally owned and relevant index sets. */ + + dealii::IndexSet locally_owned(12); + dealii::IndexSet locally_relevant(12); + locally_owned.add_range(0, 12); + locally_relevant.add_range(0, 12); + + const auto scalar_partitioner = + std::make_shared( + locally_owned, locally_relevant, MPI_COMM_WORLD); + + const auto vector_partitioner = + ryujin::Vectors::create_vector_partitioner(scalar_partitioner, 4); + + ryujin::Vectors::MultiComponentVector state_vector; + + state_vector.reinit_with_vector_partitioner(vector_partitioner); + + using HostSpace = dealii::MemorySpace::Host; + using DefaultSpace = dealii::MemorySpace::Default; + + const auto print_status = [&]() { + std::cout << "HostSpace active == " + << state_vector.is_active_memory_space() << std::endl + << "DefaultSpace active == " + << state_vector.is_active_memory_space() + << std::endl; + }; + + print_status(); + + for (unsigned int i = 0; i < 12; ++i) { + dealii::Tensor<1, 4, double> tensor{{static_cast(10 * i), + static_cast(10 * i + 1), + static_cast(10 * i + 2), + static_cast(10 * i + 3)}}; + state_vector.write_tensor(tensor, i); + } + + std::cout << "\nState vector:\n"; + for (unsigned int i = 0; i < 8; i += 1) { + std::cout << state_vector.read_tensor(i) << "\n"; + } + + std::cout << "After move to DefaultSpace:" << std::endl; + state_vector.move_to_memory_space(); + print_status(); + std::cout << "After repeated move to DefaultSpace:" << std::endl; + state_vector.move_to_memory_space(); + print_status(); + + const auto &view = state_vector.template get_view(); + using ExecutionSpace = DefaultSpace::kokkos_space::execution_space; + const auto exec = ExecutionSpace{}; + Kokkos::parallel_for("test", + Kokkos::RangePolicy(exec, 0, 3), + [=](std::size_t i) { + auto a = view.read_tensor(i); + a *= 10.; + view.write_tensor(a, i); + }); + + std::cout << "After move to HostSpace:" << std::endl; + state_vector.move_to_memory_space(); + print_status(); + + std::cout << "\nState vector:\n"; + for (unsigned int i = 0; i < 8; i += 1) { + std::cout << state_vector.read_tensor(i) << "\n"; + } +} diff --git a/tests/common/multicomponent_vector_03.threads=2.output b/tests/common/multicomponent_vector_03.threads=2.output new file mode 100644 index 000000000..f6de43dc2 --- /dev/null +++ b/tests/common/multicomponent_vector_03.threads=2.output @@ -0,0 +1,31 @@ +HostSpace active == 1 +DefaultSpace active == 0 + +State vector: +0 1 2 3 +10 11 12 13 +20 21 22 23 +30 31 32 33 +40 41 42 43 +50 51 52 53 +60 61 62 63 +70 71 72 73 +After move to DefaultSpace: +HostSpace active == 0 +DefaultSpace active == 1 +After repeated move to DefaultSpace: +HostSpace active == 0 +DefaultSpace active == 1 +After move to HostSpace: +HostSpace active == 1 +DefaultSpace active == 0 + +State vector: +0 10 20 30 +100 110 120 130 +200 210 220 230 +30 31 32 33 +40 41 42 43 +50 51 52 53 +60 61 62 63 +70 71 72 73