ryujin 2.1.1 revision 2062b43aae59b7ee510dec4ae666e29d56c4f0be
multicomponent_vector.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: MIT
3// Copyright (C) 2020 - 2023 by the ryujin authors
4//
5
6#pragma once
7
8#include "simd.h"
9
10#include <deal.II/base/mpi.h>
11#include <deal.II/base/partitioner.h>
12#include <deal.II/base/vectorization.h>
13#include <deal.II/lac/la_parallel_vector.h>
14
15namespace ryujin
16{
36 std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
38 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
39 &scalar_partitioner,
40 const unsigned int n_components);
41
42
53 template <typename Number,
54 int n_comp,
55 int simd_length = dealii::VectorizedArray<Number>::size()>
57 : public dealii::LinearAlgebra::distributed::Vector<Number>
58 {
59 public:
65 using VectorizedArray = dealii::VectorizedArray<Number, simd_length>;
66
72 using scalar_type = dealii::LinearAlgebra::distributed::Vector<Number>;
73
78 using scalar_type::operator=;
79
87 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
88 &scalar_partitioner);
89
104 void extract_component(scalar_type &scalar_vector,
105 unsigned int component) const;
106
120 void insert_component(const scalar_type &scalar_vector,
121 unsigned int component);
122
132 template <typename Number2 = Number,
133 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
134 Tensor get_tensor(const unsigned int i) const;
135
144 template <typename Number2 = Number,
145 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
146 Tensor get_tensor(const unsigned int *js) const;
147
161 template <typename Number2 = Number,
162 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
163 void write_tensor(const Tensor &tensor, const unsigned int i);
164 };
165
166
167#ifndef DOXYGEN
168 /* Template definitions: */
169
170 template <typename Number, int n_comp, int simd_length>
173 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
174 &scalar_partitioner)
175 {
176 /* Special case of a zero component vector */
177 if (n_comp == 0)
178 return;
179
180 auto vector_partitioner =
181 create_vector_partitioner(scalar_partitioner, n_comp);
182
183 dealii::LinearAlgebra::distributed::Vector<Number>::reinit(
184 vector_partitioner);
185 }
186
187
188 template <typename Number, int n_comp, int simd_length>
190 scalar_type &scalar_vector, unsigned int component) const
191 {
192 Assert(n_comp > 0,
193 dealii::ExcMessage(
194 "Cannot extract from a vector with zero components."));
195
196 Assert(n_comp * scalar_vector.get_partitioner()->locally_owned_size() ==
197 this->get_partitioner()->locally_owned_size(),
198 dealii::ExcMessage("Called with a scalar_vector argument that has "
199 "incompatible local range."));
200 const auto local_size =
201 scalar_vector.get_partitioner()->locally_owned_size();
202 for (unsigned int i = 0; i < local_size; ++i)
203 scalar_vector.local_element(i) =
204 this->local_element(i * n_comp + component);
205 scalar_vector.update_ghost_values();
206 }
207
208
209 template <typename Number, int n_comp, int simd_length>
211 const scalar_type &scalar_vector, unsigned int component)
212 {
213 Assert(n_comp > 0,
214 dealii::ExcMessage(
215 "Cannot insert into a vector with zero components."));
216
217 Assert(n_comp * scalar_vector.get_partitioner()->locally_owned_size() ==
218 this->get_partitioner()->locally_owned_size(),
219 dealii::ExcMessage("Called with a scalar_vector argument that has "
220 "incompatible local range."));
221 const auto local_size =
222 scalar_vector.get_partitioner()->locally_owned_size();
223 for (unsigned int i = 0; i < local_size; ++i)
224 this->local_element(i * n_comp + component) =
225 scalar_vector.local_element(i);
226 }
227
228 /* Inline function definitions: */
229
230 template <typename Number, int n_comp, int simd_length>
231 template <typename Number2, typename Tensor>
232 DEAL_II_ALWAYS_INLINE inline Tensor
234 const unsigned int i) const
235 {
236 static_assert(std::is_same<Number2, typename Tensor::value_type>::value,
237 "dummy type mismatch");
238 Tensor tensor;
239
240 /* Special case of a zero component vector */
241 if constexpr (n_comp == 0)
242 return tensor;
243
244 if constexpr (std::is_same<Number, Number2>::value) {
245 /* Non-vectorized sequential access. */
246
247 for (unsigned int d = 0; d < n_comp; ++d)
248 tensor[d] = this->local_element(i * n_comp + d);
249
250 } else if constexpr (std::is_same<VectorizedArray, Number2>::value) {
251
252 /* Vectorized fast access. index must be divisible by simd_length */
253 std::array<unsigned int, VectorizedArray::size()> indices;
254 for (unsigned int k = 0; k < VectorizedArray::size(); ++k)
255 indices[k] = k * n_comp;
256
257 dealii::vectorized_load_and_transpose(
258 n_comp, this->begin() + i * n_comp, indices.data(), &tensor[0]);
259
260 } else {
261 /* not implemented */
262 __builtin_trap();
263 }
264
265 return tensor;
266 }
267
268
269 template <typename Number, int n_comp, int simd_length>
270 template <typename Number2, typename Tensor>
271 DEAL_II_ALWAYS_INLINE inline Tensor
273 const unsigned int *js) const
274 {
275 static_assert(std::is_same<Number2, typename Tensor::value_type>::value,
276 "dummy type mismatch");
277 Tensor tensor;
278
279 /* Special case of a zero component vector */
280 if constexpr (n_comp == 0)
281 return tensor;
282
283 if constexpr (std::is_same<Number, Number2>::value) {
284 /* Non-vectorized sequential access. */
285
286 for (unsigned int d = 0; d < n_comp; ++d)
287 tensor[d] = this->local_element(js[0] * n_comp + d);
288
289 } else if constexpr (std::is_same<VectorizedArray, Number2>::value) {
290 /* Vectorized fast access. index must be divisible by simd_length */
291
292 std::array<unsigned int, VectorizedArray::size()> indices;
293 for (unsigned int k = 0; k < VectorizedArray::size(); ++k)
294 indices[k] = js[k] * n_comp;
295
296 dealii::vectorized_load_and_transpose(
297 n_comp, this->begin(), indices.data(), &tensor[0]);
298
299 } else {
300 /* not implemented */
301 __builtin_trap();
302 }
303
304 return tensor;
305 }
306
307
308 template <typename Number, int n_comp, int simd_length>
309 template <typename Number2, typename Tensor>
310 DEAL_II_ALWAYS_INLINE inline void
312 const Tensor &tensor, const unsigned int i)
313 {
314 static_assert(std::is_same<Number2, typename Tensor::value_type>::value,
315 "dummy type mismatch");
316
317 /* Special case of a zero component vector */
318 if constexpr (n_comp == 0)
319 return;
320
321 if constexpr (std::is_same<Number, Number2>::value) {
322 /* Non-vectorized sequential access. */
323
324 for (unsigned int d = 0; d < n_comp; ++d)
325 this->local_element(i * n_comp + d) = tensor[d];
326
327 } else if constexpr (std::is_same<VectorizedArray, Number2>::value) {
328 /* Vectorized fast access. index must be divisible by simd_length */
329
330 std::array<unsigned int, VectorizedArray::size()> indices;
331 for (unsigned int k = 0; k < VectorizedArray::size(); ++k)
332 indices[k] = k * n_comp;
333
334 dealii::vectorized_transpose_and_store(false,
335 n_comp,
336 &tensor[0],
337 indices.data(),
338 this->begin() + i * n_comp);
339
340 } else {
341 /* not implemented */
342 __builtin_trap();
343 }
344 }
345#endif
346
347} // namespace ryujin
void reinit_with_scalar_partitioner(const std::shared_ptr< const dealii::Utilities::MPI::Partitioner > &scalar_partitioner)
void insert_component(const scalar_type &scalar_vector, unsigned int component)
void extract_component(scalar_type &scalar_vector, unsigned int component) const
Tensor get_tensor(const unsigned int i) const
dealii::LinearAlgebra::distributed::Vector< Number > scalar_type
dealii::VectorizedArray< Number, simd_length > VectorizedArray
void write_tensor(const Tensor &tensor, const unsigned int i)
Tensor get_tensor(const unsigned int *js) const
std::shared_ptr< const dealii::Utilities::MPI::Partitioner > create_vector_partitioner(const std::shared_ptr< const dealii::Utilities::MPI::Partitioner > &scalar_partitioner, const unsigned int n_components)