ryujin 2.1.1 revision 71cdc42292164f8095c0bd62c4ea75ebcfac58fa
Loading...
Searching...
No Matches
multicomponent_vector.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2020 - 2026 by the ryujin authors
4//
5
6#pragma once
7
8#include <compile_time_options.h>
9
10#include "gpu.h"
11#include "loop.h"
12
13#include <deal.II/base/mpi.h>
14#include <deal.II/base/partitioner.h>
15#include <deal.II/base/vectorization.h>
16#include <deal.II/lac/la_parallel_vector.h>
17
18namespace ryujin
19{
20 namespace Vectors
21 {
41 std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
42 create_vector_partitioner(
43 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
44 &scalar_partitioner,
45 const unsigned int n_comp);
46
47
48 template <typename Number,
49 int n_comp,
50 int simd_length = dealii::VectorizedArray<Number>::size(),
51 typename MemorySpace = dealii::MemorySpace::Host,
52 bool writable = true>
53 class MultiComponentVectorView;
54
55
63 template <typename Number,
64 int n_comp,
65 int simd_length = dealii::VectorizedArray<Number>::size()>
66 class MultiComponentVector
67 : public MirroredStorage<
68 MultiComponentVector<Number, n_comp, simd_length>>
69 {
70 public:
75
76 MultiComponentVector() = default;
77
78 MultiComponentVector(const MultiComponentVector &other);
79
80 MultiComponentVector(MultiComponentVector &&other) noexcept;
81
95 void reinit_with_vector_partitioner(
96 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
97 &vector_partitioner,
98 const TransferPolicy transfer_policy =
99 TransferPolicy::explicit_transfers);
100
115 void reinit_with_scalar_partitioner(
116 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
117 &scalar_partitioner,
118 const TransferPolicy transfer_policy =
119 TransferPolicy::explicit_transfers);
120
121 MultiComponentVector &operator=(const MultiComponentVector &other);
122
123 MultiComponentVector &operator=(MultiComponentVector &&other) noexcept;
124
126
130
139 template <typename MemorySpace = dealii::MemorySpace::Host>
140 MultiComponentVectorView<Number, n_comp, simd_length, MemorySpace, true>
141 view();
142
150 template <typename MemorySpace = dealii::MemorySpace::Host>
151 MultiComponentVectorView<Number, n_comp, simd_length, MemorySpace, false>
152 view() const;
153
164 template <typename MemorySpace = dealii::MemorySpace::Host>
165 dealii::LinearAlgebra::distributed::Vector<Number, MemorySpace> &
166 deal_ii_vector()
167 requires(n_comp == 1);
168
176 template <typename MemorySpace = dealii::MemorySpace::Host>
177 const dealii::LinearAlgebra::distributed::Vector<Number, MemorySpace> &
178 deal_ii_vector() const
179 requires(n_comp == 1);
180
181 /*
182 * The is_resident(), copy_to_memory_space(), move_to_memory_space(),
183 * transfer_policy(), and set_transfer_policy() methods are inherited
184 * from the MirroredStorage base class.
185 */
186
188
192
196 template <typename MemorySpace>
197 void zero_out_ghost_values_on_memory_space();
198
203 template <typename MemorySpace>
204 void update_ghost_values_on_memory_space();
205
211 template <typename MemorySpace>
212 void compress_on_memory_space(dealii::VectorOperation::values operation);
213
214 private:
216
220
221 /*
222 * Note: The storage is marked mutable so that the (logically const)
223 * copy_to_memory_space() operation can populate a mirror from within
224 * a const view() under the implicit_transfers policy.
225 *
226 * We avoid setting up the default_vector_ if default and host happen
227 * to be the same memory space (see have_separate_memory_spaces).
228 */
229
230 mutable dealii::LinearAlgebra::distributed::
231 Vector<Number, dealii::MemorySpace::Host>
232 host_vector_;
233
234 mutable dealii::LinearAlgebra::distributed::
235 Vector<Number, dealii::MemorySpace::Default>
236 default_vector_;
237
238 /*
239 * Storage primitives used by the MirroredStorage base class:
240 */
241
242 template <typename MemorySpace>
243 void allocate_storage() const;
244
245 template <typename To, typename From>
246 void deep_copy_storage() const;
247
248 template <typename MemorySpace>
249 void deallocate_storage();
250
251
252 friend class MirroredStorage<
253 MultiComponentVector<Number, n_comp, simd_length>>;
254
255 template <typename, int, int, typename, bool>
256 friend class MultiComponentVectorView;
257
259 };
260
261
276 template <typename Number,
277 int n_comp,
278 int simd_length,
279 typename MemorySpace,
280 bool writable>
281 class MultiComponentVectorView
282 {
283 public:
288
289 MultiComponentVectorView() = default;
290
291 MultiComponentVectorView(MultiComponentVector<Number, n_comp, simd_length>
292 &multi_component_vector)
293 requires(writable);
294
295 MultiComponentVectorView(
296 const MultiComponentVector<Number, n_comp, simd_length>
297 &multi_component_vector)
298 requires(!writable);
299
308 template <bool other_writable>
309 DEAL_II_HOST_DEVICE MultiComponentVectorView(
310 const MultiComponentVectorView<Number,
311 n_comp,
312 simd_length,
313 MemorySpace,
314 other_writable> &other)
315 requires(!writable && other_writable);
316
317 template <typename MultiComponentVector>
318 void reinit(MultiComponentVector &multi_component_vector)
319 requires(writable != std::is_const_v<MultiComponentVector>);
320
322
326
332 using ScalarVector =
333 dealii::LinearAlgebra::distributed::Vector<Number, MemorySpace>;
334
336
340
363 template <typename Functor = std::identity>
364 void extract_component(ScalarVector &scalar_vector,
365 unsigned int component,
366 const Functor &functor = std::identity{}) const;
367
385 template <typename Functor = std::identity>
386 void insert_component(const ScalarVector &scalar_vector,
387 unsigned int component,
388 const Functor &functor = std::identity{}) const
389 requires writable;
390
395 template <typename Functor = std::identity>
396 void insert_component(const dealii::Vector<Number> &scalar_vector,
397 unsigned int component,
398 const Functor &functor = std::identity{}) const
399 requires writable;
400
405 void sadd(const Number s,
406 const Number a,
407 const MultiComponentVectorView<Number,
408 n_comp,
409 simd_length,
410 MemorySpace,
411 /*writable=*/false> &v) const
412 requires writable;
413
415
420
431 template <typename Number2 = Number>
432 DEAL_II_HOST_DEVICE Number2 read_entry(const unsigned int i) const;
433
444 template <typename Number2 = Number>
445 DEAL_II_HOST_DEVICE Number2 read_entry(const unsigned int *js) const;
446
455 template <typename Number2 = Number,
456 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
457 DEAL_II_HOST_DEVICE Tensor read_tensor(const unsigned int i) const;
458
467 template <typename Number2 = Number,
468 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
469 DEAL_II_HOST_DEVICE Tensor read_tensor(const unsigned int *js) const;
470
482 template <typename Number2 = Number>
483 DEAL_II_HOST_DEVICE void write_entry(const Number2 &entry,
484 const unsigned int i) const
485 requires writable;
486
501 template <typename Number2 = Number,
502 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
503 DEAL_II_HOST_DEVICE void write_tensor(const Tensor &tensor,
504 const unsigned int i) const
505 requires writable;
506
521 template <typename Number2 = Number>
522 DEAL_II_HOST_DEVICE void add_entry(const Number2 &entry,
523 const unsigned int i) const
524 requires writable;
525
540 template <typename Number2 = Number,
541 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
542 DEAL_II_HOST_DEVICE void add_tensor(const Tensor &tensor,
543 const unsigned int i) const
544 requires writable;
545
547
551
555 void zero_out_ghost_values() const
556 requires(writable);
557
562 void update_ghost_values() const
563 requires(writable);
564
570 void compress(dealii::VectorOperation::values operation) const
571 requires(writable);
572
573 private:
575
579
580 using MCV = MultiComponentVector<Number, n_comp, simd_length>;
581 std::conditional_t<writable, MCV *, const MCV *> multi_component_vector_;
582
583 Number *data_;
584 unsigned int n_locally_owned_;
585 unsigned int n_locally_relevant_;
586
587 template <typename, int, int, typename, bool>
588 friend class MultiComponentVectorView;
589
591 };
592
593
594#ifndef DOXYGEN
595 /*
596 * -------------------------------------------------------------------------
597 * Inline function definitions
598 * -------------------------------------------------------------------------
599 */
600
601
602 template <typename Number, int n_comp, int simd_length>
603 MultiComponentVector<Number, n_comp, simd_length>::MultiComponentVector(
604 const MultiComponentVector &other)
605 {
606 *this = other;
607 }
608
609
610 template <typename Number, int n_comp, int simd_length>
611 MultiComponentVector<Number, n_comp, simd_length>::MultiComponentVector(
612 MultiComponentVector &&other) noexcept
613 {
614 *this = other;
615 }
616
617
618 template <typename Number, int n_comp, int simd_length>
619 void MultiComponentVector<Number, n_comp, simd_length>::
620 reinit_with_vector_partitioner(
621 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
622 &vector_partitioner,
623 const TransferPolicy transfer_policy)
624 {
625 /*
626 * Drop the transfer policy for the duration of the reinit so that a
627 * pinned memory space of a previous policy does not interfere:
628 */
629 this->set_transfer_policy(TransferPolicy::explicit_transfers);
630
631 /* Special case of a zero component vector */
632 if (n_comp == 0) {
633 /* A zero component vector is trivially resident everywhere: */
634 this->reset_residency(/*host*/ true, /*default*/ true);
635 this->set_transfer_policy(transfer_policy);
636 return;
637 }
638
639 host_vector_.reinit(vector_partitioner);
640
641 /*
642 * The vector is resident on the host memory space only. Device
643 * storage is allocated lazily on the first copy_to_memory_space() /
644 * move_to_memory_space(); drop possibly stale device storage from a
645 * previous reinit:
646 */
647 if constexpr (have_separate_memory_spaces)
648 default_vector_.reinit(0);
649
650 this->reset_residency(/*host*/ true,
651 /*default*/ !have_separate_memory_spaces);
652
653 this->set_transfer_policy(transfer_policy);
654 }
655
656
657 template <typename Number, int n_comp, int simd_length>
658 void MultiComponentVector<Number, n_comp, simd_length>::
659 reinit_with_scalar_partitioner(
660 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
661 &scalar_partitioner,
662 const TransferPolicy transfer_policy)
663 {
664 /*
665 * Drop the transfer policy for the duration of the reinit so that a
666 * pinned memory space of a previous policy does not interfere:
667 */
668 this->set_transfer_policy(TransferPolicy::explicit_transfers);
669
670 /* Special case of a zero component vector: */
671 if (n_comp == 0) {
672 /* A zero component vector is trivially resident everywhere: */
673 this->reset_residency(/*host*/ true, /*default*/ true);
674 this->set_transfer_policy(transfer_policy);
675 return;
676 }
677
678 /* Special case of a scalar vector: */
679 if (n_comp == 1) {
680 host_vector_.reinit(scalar_partitioner);
681 } else {
682 auto vector_partitioner =
683 create_vector_partitioner(scalar_partitioner, n_comp);
684 host_vector_.reinit(vector_partitioner);
685 }
686
687 /*
688 * The vector is resident on the host memory space only. Device
689 * storage is allocated lazily on the first copy_to_memory_space() /
690 * move_to_memory_space(); drop possibly stale device storage from a
691 * previous reinit:
692 */
693 if constexpr (have_separate_memory_spaces)
694 default_vector_.reinit(0);
695
696 /*
697 * If both memory spaces coincide the single host allocation is
698 * trivially resident on both of them:
699 */
700 this->reset_residency(/*host*/ true,
701 /*default*/ !have_separate_memory_spaces);
702
703 this->set_transfer_policy(transfer_policy);
704 }
705
706
707 template <typename Number, int n_comp, int simd_length>
708 auto MultiComponentVector<Number, n_comp, simd_length>::operator=(
709 const MultiComponentVector &other) -> MultiComponentVector &
710 {
711 /* Copy residency state and transfer policy: */
712 static_cast<MirroredStorage<MultiComponentVector> &>(*this) = other;
713
714 host_vector_ = other.host_vector_;
715 if constexpr (have_separate_memory_spaces)
716 default_vector_ = other.default_vector_;
717
718 return *this;
719 }
720
721
722 template <typename Number, int n_comp, int simd_length>
723 auto MultiComponentVector<Number, n_comp, simd_length>::operator=(
724 MultiComponentVector &&other) noexcept -> MultiComponentVector &
725 {
726 /* Copy residency state and transfer policy: */
727 static_cast<MirroredStorage<MultiComponentVector> &>(*this) = other;
728
729 host_vector_ = std::move(other.host_vector_);
730 if constexpr (have_separate_memory_spaces)
731 default_vector_ = std::move(other.default_vector_);
732
733 return *this;
734 }
735
736
737 template <typename Number, int n_comp, int simd_length>
738 template <typename MemorySpace>
739 MultiComponentVectorView<Number, n_comp, simd_length, MemorySpace, true>
740 MultiComponentVector<Number, n_comp, simd_length>::view()
741 {
742 this->template prepare_write_access<MemorySpace>();
743
744 return MultiComponentVectorView<Number,
745 n_comp,
746 simd_length,
747 MemorySpace,
748 true>(*this);
749 }
750
751
752 template <typename Number, int n_comp, int simd_length>
753 template <typename MemorySpace>
754 MultiComponentVectorView<Number, n_comp, simd_length, MemorySpace, false>
755 MultiComponentVector<Number, n_comp, simd_length>::view() const
756 {
757 this->template prepare_read_access<MemorySpace>();
758
759 return MultiComponentVectorView<Number,
760 n_comp,
761 simd_length,
762 MemorySpace,
763 false>(*this);
764 }
765
766
767 template <typename Number, int n_comp, int simd_length>
768 template <typename MemorySpace>
769 dealii::LinearAlgebra::distributed::Vector<Number, MemorySpace> &
770 MultiComponentVector<Number, n_comp, simd_length>::deal_ii_vector()
771 requires(n_comp == 1)
772 {
773 this->template prepare_write_access<MemorySpace>();
774
775 if constexpr (std::is_same_v<MemorySpace, dealii::MemorySpace::Host>)
776 return host_vector_;
777 else
778 return default_vector_;
779 }
780
781
782 template <typename Number, int n_comp, int simd_length>
783 template <typename MemorySpace>
784 const dealii::LinearAlgebra::distributed::Vector<Number, MemorySpace> &
785 MultiComponentVector<Number, n_comp, simd_length>::deal_ii_vector() const
786 requires(n_comp == 1)
787 {
788 this->template prepare_read_access<MemorySpace>();
789
790 if constexpr (std::is_same_v<MemorySpace, dealii::MemorySpace::Host>)
791 return host_vector_;
792 else
793 return default_vector_;
794 }
795
796
797 template <typename Number, int n_comp, int simd_length>
798 template <typename MemorySpace>
799 void
800 MultiComponentVector<Number, n_comp, simd_length>::allocate_storage() const
801 {
802 using HostSpace = dealii::MemorySpace::Host;
803
804 /*
805 * The partitioner is recovered from the other (still resident)
806 * vector:
807 */
808
809 if constexpr (std::is_same_v<MemorySpace, HostSpace>) {
810 Assert(default_vector_.size() != 0, dealii::ExcNotInitialized());
811 host_vector_.reinit(default_vector_.get_partitioner());
812 } else {
813 Assert(host_vector_.size() != 0, dealii::ExcNotInitialized());
814 default_vector_.reinit(host_vector_.get_partitioner());
815 }
816 }
817
818
819 template <typename Number, int n_comp, int simd_length>
820 template <typename To, typename From>
821 void
822 MultiComponentVector<Number, n_comp, simd_length>::deep_copy_storage() const
823 {
824 using HostSpace = dealii::MemorySpace::Host;
825
826 if constexpr (std::is_same_v<To, HostSpace>) {
827 host_vector_.import_elements(default_vector_,
828 dealii::VectorOperation::insert);
829 } else {
830 default_vector_.import_elements(host_vector_,
831 dealii::VectorOperation::insert);
832 }
833 }
834
835
836 template <typename Number, int n_comp, int simd_length>
837 template <typename MemorySpace>
838 void MultiComponentVector<Number, n_comp, simd_length>::deallocate_storage()
839 {
840 using HostSpace = dealii::MemorySpace::Host;
841
842 if constexpr (std::is_same_v<MemorySpace, HostSpace>) {
843 host_vector_.reinit(0);
844
845 } else {
846 default_vector_.reinit(0);
847 }
848 }
849
850
851 template <typename Number, int n_components, int simd_length>
852 template <typename MemorySpace>
853 void MultiComponentVector<Number, n_components, simd_length>::
854 zero_out_ghost_values_on_memory_space()
855 {
856 using HostSpace = dealii::MemorySpace::Host;
857 using DefaultSpace = dealii::MemorySpace::Default;
858 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
859 std::is_same_v<MemorySpace, DefaultSpace>,
860 "Unexpected memory space");
861
862 Assert(this->template is_resident<MemorySpace>(),
863 dealii::ExcMessage("The chosen memory space is not resident."));
864
865 if constexpr (have_separate_memory_spaces &&
866 !std::is_same_v<MemorySpace, HostSpace>) {
867 default_vector_.zero_out_ghost_values();
868 } else {
869 host_vector_.zero_out_ghost_values();
870 }
871 }
872
873
874 template <typename Number, int n_components, int simd_length>
875 template <typename MemorySpace>
876 void MultiComponentVector<Number, n_components, simd_length>::
877 update_ghost_values_on_memory_space()
878 {
879 using HostSpace = dealii::MemorySpace::Host;
880 using DefaultSpace = dealii::MemorySpace::Default;
881 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
882 std::is_same_v<MemorySpace, DefaultSpace>,
883 "Unexpected memory space");
884
885 Assert(this->template is_resident<MemorySpace>(),
886 dealii::ExcMessage("The chosen memory space is not resident."));
887
888 if constexpr (have_separate_memory_spaces &&
889 !std::is_same_v<MemorySpace, HostSpace>) {
890 default_vector_.update_ghost_values();
891 } else {
892 host_vector_.update_ghost_values();
893 }
894 }
895
896
897 template <typename Number, int n_components, int simd_length>
898 template <typename MemorySpace>
899 void MultiComponentVector<Number, n_components, simd_length>::
900 compress_on_memory_space(dealii::VectorOperation::values operation)
901 {
902 using HostSpace = dealii::MemorySpace::Host;
903 using DefaultSpace = dealii::MemorySpace::Default;
904 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
905 std::is_same_v<MemorySpace, DefaultSpace>,
906 "Unexpected memory space");
907
908 Assert(this->template is_resident<MemorySpace>(),
909 dealii::ExcMessage("The chosen memory space is not resident."));
910
911 if constexpr (have_separate_memory_spaces &&
912 !std::is_same_v<MemorySpace, HostSpace>) {
913 default_vector_.compress(operation);
914 } else {
915 host_vector_.compress(operation);
916 }
917 }
918
919
920 template <typename Number,
921 int n_comp,
922 int simd_l,
923 typename MemorySpace,
924 bool writable>
925 MultiComponentVectorView<Number, n_comp, simd_l, MemorySpace, writable>::
926 MultiComponentVectorView(MultiComponentVector<Number, n_comp, simd_l>
927 &multi_component_vector)
928 requires(writable)
929 {
930 reinit(multi_component_vector);
931 }
932
933
934 template <typename Number,
935 int n_comp,
936 int simd_l,
937 typename MemorySpace,
938 bool writable>
939 MultiComponentVectorView<Number, n_comp, simd_l, MemorySpace, writable>::
940 MultiComponentVectorView(
941 const MultiComponentVector<Number, n_comp, simd_l>
942 &multi_component_vector)
943 requires(!writable)
944 {
945 reinit(multi_component_vector);
946 }
947
948
949 template <typename Number,
950 int n_comp,
951 int simd_l,
952 typename MemorySpace,
953 bool writable>
954 template <bool other_writable>
955 DEAL_II_HOST_DEVICE
956 MultiComponentVectorView<Number, n_comp, simd_l, MemorySpace, writable>::
957 MultiComponentVectorView(
958 const MultiComponentVectorView<Number,
959 n_comp,
960 simd_l,
961 MemorySpace,
962 other_writable> &other)
963 requires(!writable && other_writable)
964 : multi_component_vector_(other.multi_component_vector_)
965 , data_(other.data_)
966 , n_locally_owned_(other.n_locally_owned_)
967 , n_locally_relevant_(other.n_locally_relevant_)
968 {
969 }
970
971
972 template <typename Number,
973 int n_comp,
974 int simd_l,
975 typename MemorySpace,
976 bool writable>
977 template <typename MultiComponentVector>
978 void
979 MultiComponentVectorView<Number, n_comp, simd_l, MemorySpace, writable>::
980 reinit(MultiComponentVector &multi_component_vector)
981 requires(writable != std::is_const_v<MultiComponentVector>)
982 {
983 using HostSpace = dealii::MemorySpace::Host;
984 using DefaultSpace = dealii::MemorySpace::Default;
985 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
986 std::is_same_v<MemorySpace, DefaultSpace>,
987 "Unexpected memory space");
988
989 multi_component_vector_ = &multi_component_vector;
990
991 if constexpr (have_separate_memory_spaces &&
992 !std::is_same_v<MemorySpace, HostSpace>) {
993 auto &vector = multi_component_vector_->default_vector_;
994 const auto &partitioner = vector.get_partitioner();
995
996 data_ = vector.begin();
997 n_locally_owned_ = partitioner->locally_owned_size();
998 n_locally_relevant_ = n_locally_owned_ + partitioner->n_ghost_indices();
999
1000 } else {
1001 auto &vector = multi_component_vector_->host_vector_;
1002 const auto &partitioner = vector.get_partitioner();
1003
1004 data_ = vector.begin();
1005 n_locally_owned_ = partitioner->locally_owned_size();
1006 n_locally_relevant_ = n_locally_owned_ + partitioner->n_ghost_indices();
1007 }
1008 }
1009
1010
1011 template <typename Number,
1012 int n_comp,
1013 int simd_length,
1014 typename MemorySpace,
1015 bool writable>
1016 template <typename Functor>
1017 void MultiComponentVectorView<
1018 Number,
1019 n_comp,
1020 simd_length,
1021 MemorySpace,
1022 writable>::extract_component(ScalarVector &scalar_vector,
1023 unsigned int component,
1024 const Functor &functor) const
1025 {
1026 Assert(n_comp > 0,
1027 dealii::ExcMessage(
1028 "Cannot extract from a vector with zero components."));
1029 AssertIndexRange(component, n_comp);
1030
1031 const auto local_size = static_cast<unsigned int>(
1032 scalar_vector.get_partitioner()->locally_owned_size());
1033
1034 Assert(n_comp * local_size == n_locally_owned_,
1035 dealii::ExcMessage("Called with a scalar_vector argument that has "
1036 "incompatible local range."));
1037
1038 const auto *data = data_;
1039 auto *destination = scalar_vector.begin();
1040
1041 const auto body = [=](auto /*sentinel*/, unsigned int i) {
1042 destination[i] = functor(data[i * n_comp + component]);
1043 };
1044
1045 loop<MemorySpace, Number>(
1046 "extract_component", body, 0, /*no vectorization*/ 0, local_size);
1047
1048 scalar_vector.update_ghost_values();
1049 }
1050
1051
1052 template <typename Number,
1053 int n_comp,
1054 int simd_length,
1055 typename MemorySpace,
1056 bool writable>
1057 template <typename Functor>
1058 void MultiComponentVectorView<
1059 Number,
1060 n_comp,
1061 simd_length,
1062 MemorySpace,
1063 writable>::insert_component(const ScalarVector &scalar_vector,
1064 unsigned int component,
1065 const Functor &functor) const
1066 requires writable
1067 {
1068 using HostSpace = dealii::MemorySpace::Host;
1069 AssertThrow((std::is_same_v<MemorySpace, HostSpace>),
1070 dealii::ExcNotImplemented());
1071
1072 Assert(n_comp > 0,
1073 dealii::ExcMessage(
1074 "Cannot insert into a vector with zero components."));
1075 AssertIndexRange(component, n_comp);
1076
1077 const auto local_size =
1078 scalar_vector.get_partitioner()->locally_owned_size();
1079
1080 Assert(n_comp * local_size == n_locally_owned_,
1081 dealii::ExcMessage("Called with a scalar_vector argument that has "
1082 "incompatible local range."));
1083
1084 for (unsigned int i = 0; i < local_size; ++i)
1085 data_[i * n_comp + component] = functor(scalar_vector.local_element(i));
1086 }
1087
1088
1089 template <typename Number,
1090 int n_comp,
1091 int simd_length,
1092 typename MemorySpace,
1093 bool writable>
1094 template <typename Functor>
1095 void MultiComponentVectorView<
1096 Number,
1097 n_comp,
1098 simd_length,
1099 MemorySpace,
1100 writable>::insert_component(const dealii::Vector<Number> &scalar_vector,
1101 unsigned int component,
1102 const Functor &functor) const
1103 requires writable
1104 {
1105 using HostSpace = dealii::MemorySpace::Host;
1106 AssertThrow((std::is_same_v<MemorySpace, HostSpace>),
1107 dealii::ExcInternalError());
1108
1109 Assert(n_comp > 0,
1110 dealii::ExcMessage(
1111 "Cannot insert into a vector with zero components."));
1112 AssertIndexRange(component, n_comp);
1113
1114 const auto local_size = scalar_vector.size();
1115
1116 Assert(n_comp * local_size >= n_locally_owned_,
1117 dealii::ExcMessage("Called with a scalar_vector argument that has "
1118 "incompatible local range."));
1119
1120 for (unsigned int i = 0; i < local_size; ++i)
1121 data_[i * n_comp + component] = functor(scalar_vector[i]);
1122 }
1123
1124
1125 template <typename Num, int n_comp, int simd_l, typename MS, bool writable>
1126 void MultiComponentVectorView<Num, n_comp, simd_l, MS, writable>::sadd(
1127 const Num s,
1128 const Num a,
1129 const MultiComponentVectorView<Num,
1130 n_comp,
1131 simd_l,
1132 MS,
1133 /*writable=*/false> &v) const
1134 requires writable
1135 {
1136 using HS = dealii::MemorySpace::Host;
1137 using DS = dealii::MemorySpace::Default;
1138 static_assert(std::is_same_v<MS, HS> || std::is_same_v<MS, DS>,
1139 "Unexpected memory space");
1140
1141 /*
1142 * Note: If the host and default memory spaces coincide all views
1143 * reference the host storage.
1144 */
1145 if constexpr (have_separate_memory_spaces && !std::is_same_v<MS, HS>) {
1146 const auto &other = v.multi_component_vector_->default_vector_;
1147 multi_component_vector_->default_vector_.sadd(s, a, other);
1148 } else {
1149 const auto &other = v.multi_component_vector_->host_vector_;
1150 multi_component_vector_->host_vector_.sadd(s, a, other);
1151 }
1152 }
1153
1154
1155 template <typename Number,
1156 int n_comp,
1157 int simd_length,
1158 typename MemorySpace,
1159 bool writable>
1160 template <typename Number2>
1161 DEAL_II_HOST_DEVICE_ALWAYS_INLINE Number2
1162 MultiComponentVectorView<Number,
1163 n_comp,
1164 simd_length,
1165 MemorySpace,
1166 writable>::read_entry(const unsigned int i) const
1167 {
1168 static_assert(
1169 n_comp == 1,
1170 "Attempted to read a scalar value from a tensor-valued vector entry");
1171
1172 AssertIndexRange(i, n_locally_relevant_);
1173
1174 const auto result = read_tensor<Number2>(i);
1175 return result[0];
1176 }
1177
1178
1179 template <typename Number,
1180 int n_comp,
1181 int simd_length,
1182 typename MemorySpace,
1183 bool writable>
1184 template <typename Number2, typename Tensor>
1185 DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor
1186 MultiComponentVectorView<Number,
1187 n_comp,
1188 simd_length,
1189 MemorySpace,
1190 writable>::read_tensor(const unsigned int i) const
1191 {
1192 static_assert(std::is_same_v<Number2, typename Tensor::value_type>,
1193 "type mismatch");
1194
1195 AssertIndexRange(i, n_locally_relevant_);
1196
1197 Tensor tensor;
1198
1199 /* Special case of a zero component vector */
1200 if constexpr (n_comp == 0)
1201 return tensor;
1202
1203 using VA = dealii::VectorizedArray<Number>;
1204 if constexpr (std::is_same_v<VA, Number2>) {
1205 /* Vectorized fast access. index must be divisible by simd_length */
1206 std::array<unsigned int, VA::size()> indices;
1207 for (unsigned int k = 0; k < VA::size(); ++k)
1208 indices[k] = k * n_comp;
1209
1210 dealii::vectorized_load_and_transpose(
1211 n_comp, data_ + i * n_comp, indices.data(), &tensor[0]);
1212
1213 } else {
1214 /* Non-vectorized sequential access. */
1215 for (unsigned int d = 0; d < n_comp; ++d)
1216 tensor[d] = data_[i * n_comp + d];
1217 }
1218
1219 return tensor;
1220 }
1221
1222
1223 template <typename Number,
1224 int n_comp,
1225 int simd_length,
1226 typename MemorySpace,
1227 bool writable>
1228 template <typename Number2>
1229 DEAL_II_HOST_DEVICE_ALWAYS_INLINE Number2
1230 MultiComponentVectorView<Number,
1231 n_comp,
1232 simd_length,
1233 MemorySpace,
1234 writable>::read_entry(const unsigned int *js) const
1235 {
1236 static_assert(
1237 n_comp == 1,
1238 "Attempted to read a scalar value from a tensor-valued vector entry");
1239
1240 const auto result = read_tensor<Number2>(js);
1241 return result[0];
1242 }
1243
1244
1245 template <typename Number,
1246 int n_comp,
1247 int simd_length,
1248 typename MemorySpace,
1249 bool writable>
1250 template <typename Number2, typename Tensor>
1251 DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor
1252 MultiComponentVectorView<Number,
1253 n_comp,
1254 simd_length,
1255 MemorySpace,
1256 writable>::read_tensor(const unsigned int *js)
1257 const
1258 {
1259 static_assert(std::is_same_v<Number2, typename Tensor::value_type>,
1260 "type mismatch");
1261
1262 Tensor tensor;
1263
1264 /* Special case of a zero component vector */
1265 if constexpr (n_comp == 0)
1266 return tensor;
1267
1268 using VA = dealii::VectorizedArray<Number>;
1269 if constexpr (std::is_same_v<VA, Number2>) {
1270 /* Vectorized fast access. index must be divisible by simd_length */
1271
1272 std::array<unsigned int, VA::size()> indices;
1273 for (unsigned int k = 0; k < VA::size(); ++k) {
1274 AssertIndexRange(js[k], n_locally_relevant_);
1275 indices[k] = js[k] * n_comp;
1276 }
1277
1278 dealii::vectorized_load_and_transpose(
1279 n_comp, data_, indices.data(), &tensor[0]);
1280
1281 } else {
1282 /* Non-vectorized sequential access. */
1283
1284 AssertIndexRange(*js, n_locally_relevant_);
1285
1286 for (unsigned int d = 0; d < n_comp; ++d)
1287 tensor[d] = data_[js[0] * n_comp + d];
1288 }
1289
1290 return tensor;
1291 }
1292
1293
1294 template <typename Number,
1295 int n_comp,
1296 int simd_length,
1297 typename MemorySpace,
1298 bool writable>
1299 template <typename Number2>
1300 DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1301 MultiComponentVectorView<Number,
1302 n_comp,
1303 simd_length,
1304 MemorySpace,
1305 writable>::write_entry(const Number2 &entry,
1306 const unsigned int i) const
1307 requires writable
1308 {
1309 static_assert(n_comp == 1,
1310 "Attempted to write a scalar value into a tensor-valued "
1311 "vector entry");
1312
1313 AssertIndexRange(i, n_locally_relevant_);
1314
1315 dealii::Tensor<1, n_comp, Number2> tensor;
1316 tensor[0] = entry;
1317
1318 write_tensor<Number2>(tensor, i);
1319 }
1320
1321
1322 template <typename Number,
1323 int n_comp,
1324 int simd_length,
1325 typename MemorySpace,
1326 bool writable>
1327 template <typename Number2, typename Tensor>
1328 DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1329 MultiComponentVectorView<Number,
1330 n_comp,
1331 simd_length,
1332 MemorySpace,
1333 writable>::write_tensor(const Tensor &tensor,
1334 const unsigned int i) const
1335 requires writable
1336 {
1337 static_assert(std::is_same_v<Number2, typename Tensor::value_type>,
1338 "type mismatch");
1339
1340 AssertIndexRange(i, n_locally_relevant_);
1341
1342 /* Special case of a zero component vector */
1343 if constexpr (n_comp == 0)
1344 return;
1345
1346 using VA = dealii::VectorizedArray<Number>;
1347 if constexpr (std::is_same_v<VA, Number2>) {
1348 /* Vectorized fast access. index must be divisible by simd_length */
1349
1350 std::array<unsigned int, VA::size()> indices;
1351 for (unsigned int k = 0; k < VA::size(); ++k)
1352 indices[k] = k * n_comp;
1353
1354 dealii::vectorized_transpose_and_store(/*add into*/ false,
1355 n_comp,
1356 &tensor[0],
1357 indices.data(),
1358 data_ + i * n_comp);
1359
1360 } else {
1361 /* Non-vectorized sequential access. */
1362
1363 for (unsigned int d = 0; d < n_comp; ++d)
1364 data_[i * n_comp + d] = tensor[d];
1365 }
1366 }
1367
1368
1369 template <typename Number,
1370 int n_comp,
1371 int simd_length,
1372 typename MemorySpace,
1373 bool writable>
1374 template <typename Number2>
1375 DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1376 MultiComponentVectorView<Number,
1377 n_comp,
1378 simd_length,
1379 MemorySpace,
1380 writable>::add_entry(const Number2 &entry,
1381 const unsigned int i) const
1382 requires writable
1383 {
1384 static_assert(n_comp == 1,
1385 "Attempted to write a scalar value into a tensor-valued "
1386 "matrix entry");
1387
1388 AssertIndexRange(i, n_locally_relevant_);
1389
1390 dealii::Tensor<1, n_comp, Number2> tensor;
1391 tensor[0] = entry;
1392
1393 add_tensor<Number2>(tensor, i);
1394 }
1395
1396
1397 template <typename Number,
1398 int n_comp,
1399 int simd_length,
1400 typename MemorySpace,
1401 bool writable>
1402 template <typename Number2, typename Tensor>
1403 DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1404 MultiComponentVectorView<Number,
1405 n_comp,
1406 simd_length,
1407 MemorySpace,
1408 writable>::add_tensor(const Tensor &tensor,
1409 const unsigned int i) const
1410 requires writable
1411 {
1412 static_assert(std::is_same_v<Number2, typename Tensor::value_type>,
1413 "type mismatch");
1414
1415 AssertIndexRange(i, n_locally_relevant_);
1416
1417 /* Special case of a zero component vector */
1418 if constexpr (n_comp == 0)
1419 return;
1420
1421 using VA = dealii::VectorizedArray<Number>;
1422 if constexpr (std::is_same_v<VA, Number2>) {
1423 /* Vectorized fast access. index must be divisible by simd_length */
1424
1425 std::array<unsigned int, VA::size()> indices;
1426 for (unsigned int k = 0; k < VA::size(); ++k)
1427 indices[k] = k * n_comp;
1428
1429 dealii::vectorized_transpose_and_store(/*add into*/ true,
1430 n_comp,
1431 &tensor[0],
1432 indices.data(),
1433 data_ + i * n_comp);
1434
1435 } else {
1436 /* Non-vectorized sequential access. */
1437
1438 for (unsigned int d = 0; d < n_comp; ++d)
1439 data_[i * n_comp + d] += tensor[d];
1440 }
1441 }
1442
1443
1444 template <typename Number,
1445 int n_comp,
1446 int simd_l,
1447 typename MemorySpace,
1448 bool writable>
1449 void
1450 MultiComponentVectorView<Number, n_comp, simd_l, MemorySpace, writable>::
1451 zero_out_ghost_values() const
1452 requires(writable)
1453 {
1454 using HostSpace = dealii::MemorySpace::Host;
1455 using DefaultSpace = dealii::MemorySpace::Default;
1456 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
1457 std::is_same_v<MemorySpace, DefaultSpace>,
1458 "Unexpected memory space");
1459
1460 Assert(multi_component_vector_->template is_resident<MemorySpace>(),
1461 dealii::ExcMessage("The chosen memory space is not resident."));
1462
1463 multi_component_vector_
1464 ->template zero_out_ghost_values_on_memory_space<MemorySpace>();
1465 }
1466
1467
1468 template <typename Number,
1469 int n_comp,
1470 int simd_l,
1471 typename MemorySpace,
1472 bool writable>
1473 void
1474 MultiComponentVectorView<Number, n_comp, simd_l, MemorySpace, writable>::
1475 update_ghost_values() const
1476 requires(writable)
1477 {
1478 using HostSpace = dealii::MemorySpace::Host;
1479 using DefaultSpace = dealii::MemorySpace::Default;
1480 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
1481 std::is_same_v<MemorySpace, DefaultSpace>,
1482 "Unexpected memory space");
1483
1484 Assert(multi_component_vector_->template is_resident<MemorySpace>(),
1485 dealii::ExcMessage("The chosen memory space is not resident."));
1486
1487 multi_component_vector_
1488 ->template update_ghost_values_on_memory_space<MemorySpace>();
1489 }
1490
1491
1492 template <typename Number,
1493 int n_comp,
1494 int simd_l,
1495 typename MemorySpace,
1496 bool writable>
1497 void
1498 MultiComponentVectorView<Number, n_comp, simd_l, MemorySpace, writable>::
1499 compress(dealii::VectorOperation::values operation) const
1500 requires(writable)
1501 {
1502 using HostSpace = dealii::MemorySpace::Host;
1503 using DefaultSpace = dealii::MemorySpace::Default;
1504 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
1505 std::is_same_v<MemorySpace, DefaultSpace>,
1506 "Unexpected memory space");
1507
1508 Assert(multi_component_vector_->template is_resident<MemorySpace>(),
1509 dealii::ExcMessage("The chosen memory space is not resident."));
1510
1511 multi_component_vector_->template compress_on_memory_space<MemorySpace>(
1512 operation);
1513 }
1514
1515#endif
1516 } // namespace Vectors
1517} // namespace ryujin
ryujin::MirroredStorage::set_transfer_policy
void set_transfer_policy(const TransferPolicy transfer_policy)
ryujin::MirroredStorage::reset_residency
void reset_residency(const bool host_resident, const bool default_resident)
ryujin::MirroredStorage::transfer_policy
TransferPolicy transfer_policy() const
ryujin::Vectors::MultiComponentVectorView::write_entry
DEAL_II_HOST_DEVICE void write_entry(const Number2 &entry, const unsigned int i) const
ryujin::Vectors::MultiComponentVectorView::extract_component
void extract_component(ScalarVector &scalar_vector, unsigned int component, const Functor &functor=std::identity{}) const
ryujin::Vectors::MultiComponentVectorView::MultiComponentVectorView
DEAL_II_HOST_DEVICE MultiComponentVectorView(const MultiComponentVectorView< Number, n_comp, simd_length, MemorySpace, other_writable > &other)
ryujin::Vectors::MultiComponentVectorView::MultiComponentVectorView
MultiComponentVectorView(const MultiComponentVector< Number, n_comp, simd_length > &multi_component_vector)
ryujin::Vectors::MultiComponentVectorView::sadd
void sadd(const Number s, const Number a, const MultiComponentVectorView< Number, n_comp, simd_length, MemorySpace, false > &v) const
ryujin::Vectors::MultiComponentVectorView::write_tensor
DEAL_II_HOST_DEVICE void write_tensor(const Tensor &tensor, const unsigned int i) const
ryujin::Vectors::MultiComponentVectorView::add_tensor
DEAL_II_HOST_DEVICE void add_tensor(const Tensor &tensor, const unsigned int i) const
ryujin::Vectors::MultiComponentVectorView::ScalarVector
dealii::LinearAlgebra::distributed::Vector< Number, MemorySpace > ScalarVector
Definition multicomponent_vector.h:333
ryujin::Vectors::MultiComponentVectorView::read_tensor
DEAL_II_HOST_DEVICE Tensor read_tensor(const unsigned int i) const
ryujin::Vectors::MultiComponentVectorView::compress
void compress(dealii::VectorOperation::values operation) const
ryujin::Vectors::MultiComponentVectorView::MultiComponentVectorView
MultiComponentVectorView(MultiComponentVector< Number, n_comp, simd_length > &multi_component_vector)
ryujin::Vectors::MultiComponentVectorView::add_entry
DEAL_II_HOST_DEVICE void add_entry(const Number2 &entry, const unsigned int i) const
ryujin::Vectors::MultiComponentVectorView::reinit
void reinit(MultiComponentVector &multi_component_vector)
ryujin::Vectors::MultiComponentVectorView::read_entry
DEAL_II_HOST_DEVICE Number2 read_entry(const unsigned int i) const
ryujin::Vectors::MultiComponentVectorView::read_tensor
DEAL_II_HOST_DEVICE Tensor read_tensor(const unsigned int *js) const
ryujin::Vectors::MultiComponentVectorView::insert_component
void insert_component(const dealii::Vector< Number > &scalar_vector, unsigned int component, const Functor &functor=std::identity{}) const
ryujin::Vectors::MultiComponentVectorView::insert_component
void insert_component(const ScalarVector &scalar_vector, unsigned int component, const Functor &functor=std::identity{}) const
ryujin::Vectors::MultiComponentVectorView::read_entry
DEAL_II_HOST_DEVICE Number2 read_entry(const unsigned int *js) const
ryujin::Vectors::MultiComponentVector::deal_ii_vector
dealii::LinearAlgebra::distributed::Vector< Number, MemorySpace > & deal_ii_vector()
ryujin::Vectors::MultiComponentVector::deal_ii_vector
const dealii::LinearAlgebra::distributed::Vector< Number, MemorySpace > & deal_ii_vector() const
ryujin::Vectors::MultiComponentVector::reinit_with_vector_partitioner
void reinit_with_vector_partitioner(const std::shared_ptr< const dealii::Utilities::MPI::Partitioner > &vector_partitioner, const TransferPolicy transfer_policy=TransferPolicy::explicit_transfers)
ryujin::Vectors::MultiComponentVector::view
MultiComponentVectorView< Number, n_comp, simd_length, MemorySpace, false > view() const
ryujin::Vectors::MultiComponentVector::operator=
MultiComponentVector & operator=(const MultiComponentVector &other)
ryujin::Vectors::MultiComponentVector::compress_on_memory_space
void compress_on_memory_space(dealii::VectorOperation::values operation)
ryujin::Vectors::MultiComponentVector::MultiComponentVector
MultiComponentVector(const MultiComponentVector &other)
ryujin::Vectors::MultiComponentVector::MultiComponentVector
MultiComponentVector(MultiComponentVector &&other) noexcept
ryujin::Vectors::MultiComponentVector::reinit_with_scalar_partitioner
void reinit_with_scalar_partitioner(const std::shared_ptr< const dealii::Utilities::MPI::Partitioner > &scalar_partitioner, const TransferPolicy transfer_policy=TransferPolicy::explicit_transfers)
ryujin::Vectors::MultiComponentVector::view
MultiComponentVectorView< Number, n_comp, simd_length, MemorySpace, true > view()
ryujin::Vectors::MultiComponentVector::operator=
MultiComponentVector & operator=(MultiComponentVector &&other) noexcept
ryujin::TransferPolicy
TransferPolicy
Definition gpu.h:88
ryujin::have_separate_memory_spaces
constexpr bool have_separate_memory_spaces
Definition gpu.h:29
ryujin::TransferPolicy::explicit_transfers
@ explicit_transfers
ryujin::Vectors::create_vector_partitioner
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_comp)
ryujin::read_entry
DEAL_II_ALWAYS_INLINE T read_entry(const V &vector, unsigned int i)
Definition simd.h:400
ryujin::write_entry
DEAL_II_ALWAYS_INLINE void write_entry(V &vector, const T &values, unsigned int i)
Definition simd.h:513
ryujin::Vectors::ScalarVector
MultiComponentVector< Number, 1 > ScalarVector
Definition state_vector.h:43
Generated by doxygen 1.9.8