ryujin 2.1.1 revision ee5cbcbf2346c1299c942d0e1f13b46449973c18
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
12#include <deal.II/base/mpi.h>
13#include <deal.II/base/partitioner.h>
14#include <deal.II/base/vectorization.h>
15#include <deal.II/lac/la_parallel_vector.h>
16
17namespace ryujin
18{
19 namespace Vectors
20 {
40 std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
41 create_vector_partitioner(
42 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
43 &scalar_partitioner,
44 const unsigned int n_comp);
45
46
47 template <typename Number,
48 int n_comp,
49 int simd_length = dealii::VectorizedArray<Number>::size(),
50 typename MemorySpace = dealii::MemorySpace::Host,
51 bool writable = true>
52 class MultiComponentVectorView;
53
54
62 template <typename Number,
63 int n_comp,
64 int simd_length = dealii::VectorizedArray<Number>::size()>
65 class MultiComponentVector
66 : public MirroredStorage<
67 MultiComponentVector<Number, n_comp, simd_length>>
68 {
69 public:
74
75 MultiComponentVector() = default;
76
77 MultiComponentVector(const MultiComponentVector &other);
78
79 MultiComponentVector(MultiComponentVector &&other) noexcept;
80
94 void reinit_with_vector_partitioner(
95 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
96 &vector_partitioner,
97 const TransferPolicy transfer_policy =
98 TransferPolicy::explicit_transfers);
99
114 void reinit_with_scalar_partitioner(
115 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
116 &scalar_partitioner,
117 const TransferPolicy transfer_policy =
118 TransferPolicy::explicit_transfers);
119
120 MultiComponentVector &operator=(const MultiComponentVector &other);
121
122 MultiComponentVector &operator=(MultiComponentVector &&other) noexcept;
123
125
129
138 template <typename MemorySpace = dealii::MemorySpace::Host>
139 MultiComponentVectorView<Number, n_comp, simd_length, MemorySpace, true>
140 view();
141
149 template <typename MemorySpace = dealii::MemorySpace::Host>
150 MultiComponentVectorView<Number, n_comp, simd_length, MemorySpace, false>
151 view() const;
152
153 /*
154 * The is_resident(), copy_to_memory_space(), move_to_memory_space(),
155 * transfer_policy(), and set_transfer_policy() methods are inherited
156 * from the MirroredStorage base class.
157 */
158
160
164
168 template <typename MemorySpace>
169 void zero_out_ghost_values_on_memory_space();
170
175 template <typename MemorySpace>
176 void update_ghost_values_on_memory_space();
177
183 template <typename MemorySpace>
184 void compress_on_memory_space(dealii::VectorOperation::values operation);
185
186 private:
188
192
193 /*
194 * Note: The storage is marked mutable so that the (logically const)
195 * copy_to_memory_space() operation can populate a mirror from within
196 * a const view() under the implicit_transfers policy.
197 *
198 * We avoid setting up the default_vector_ if default and host happen
199 * to be the same memory space (see have_separate_memory_spaces).
200 */
201
202 mutable dealii::LinearAlgebra::distributed::
203 Vector<Number, dealii::MemorySpace::Host>
204 host_vector_;
205
206 mutable dealii::LinearAlgebra::distributed::
207 Vector<Number, dealii::MemorySpace::Default>
208 default_vector_;
209
210 /*
211 * Storage primitives used by the MirroredStorage base class:
212 */
213
214 template <typename MemorySpace>
215 void allocate_storage() const;
216
217 template <typename To, typename From>
218 void deep_copy_storage() const;
219
220 template <typename MemorySpace>
221 void deallocate_storage();
222
223
224 friend class MirroredStorage<
225 MultiComponentVector<Number, n_comp, simd_length>>;
226
227 template <typename, int, int, typename, bool>
228 friend class MultiComponentVectorView;
229
231 };
232
233
248 template <typename Number,
249 int n_comp,
250 int simd_length,
251 typename MemorySpace,
252 bool writable>
253 class MultiComponentVectorView
254 {
255 public:
260
261 MultiComponentVectorView() = default;
262
263 MultiComponentVectorView(MultiComponentVector<Number, n_comp, simd_length>
264 &multi_component_vector)
265 requires(writable);
266
267 MultiComponentVectorView(
268 const MultiComponentVector<Number, n_comp, simd_length>
269 &multi_component_vector)
270 requires(!writable);
271
280 template <bool other_writable>
281 DEAL_II_HOST_DEVICE MultiComponentVectorView(
282 const MultiComponentVectorView<Number,
283 n_comp,
284 simd_length,
285 MemorySpace,
286 other_writable> &other)
287 requires(!writable && other_writable);
288
289 template <typename MultiComponentVector>
290 void reinit(MultiComponentVector &multi_component_vector)
291 requires(writable != std::is_const_v<MultiComponentVector>);
292
294
298
304 using ScalarVector =
305 dealii::LinearAlgebra::distributed::Vector<Number, MemorySpace>;
306
308
312
331 template <typename Functor = std::identity>
332 void extract_component(ScalarVector &scalar_vector,
333 unsigned int component,
334 const Functor &functor = std::identity{}) const;
335
353 template <typename Functor = std::identity>
354 void insert_component(const ScalarVector &scalar_vector,
355 unsigned int component,
356 const Functor &functor = std::identity{}) const
357 requires writable;
358
363 template <typename Functor = std::identity>
364 void insert_component(const dealii::Vector<Number> &scalar_vector,
365 unsigned int component,
366 const Functor &functor = std::identity{}) const
367 requires writable;
368
373 void sadd(const Number s,
374 const Number a,
375 const MultiComponentVectorView<Number,
376 n_comp,
377 simd_length,
378 MemorySpace,
379 /*writable=*/false> &v) const
380 requires writable;
381
383
388
399 template <typename Number2 = Number>
400 DEAL_II_HOST_DEVICE Number2 read_entry(const unsigned int i) const;
401
412 template <typename Number2 = Number>
413 DEAL_II_HOST_DEVICE Number2 read_entry(const unsigned int *js) const;
414
423 template <typename Number2 = Number,
424 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
425 DEAL_II_HOST_DEVICE Tensor read_tensor(const unsigned int i) const;
426
435 template <typename Number2 = Number,
436 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
437 DEAL_II_HOST_DEVICE Tensor read_tensor(const unsigned int *js) const;
438
450 template <typename Number2 = Number>
451 DEAL_II_HOST_DEVICE void write_entry(const Number2 &entry,
452 const unsigned int i) const
453 requires writable;
454
469 template <typename Number2 = Number,
470 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
471 DEAL_II_HOST_DEVICE void write_tensor(const Tensor &tensor,
472 const unsigned int i) const
473 requires writable;
474
489 template <typename Number2 = Number>
490 DEAL_II_HOST_DEVICE void add_entry(const Number2 &entry,
491 const unsigned int i) const
492 requires writable;
493
508 template <typename Number2 = Number,
509 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
510 DEAL_II_HOST_DEVICE void add_tensor(const Tensor &tensor,
511 const unsigned int i) const
512 requires writable;
513
515
519
523 void zero_out_ghost_values() const
524 requires(writable);
525
530 void update_ghost_values() const
531 requires(writable);
532
538 void compress(dealii::VectorOperation::values operation) const
539 requires(writable);
540
541 private:
543
547
548 using MCV = MultiComponentVector<Number, n_comp, simd_length>;
549 std::conditional_t<writable, MCV *, const MCV *> multi_component_vector_;
550
551 Number *data_;
552 unsigned int n_locally_owned_;
553 unsigned int n_locally_relevant_;
554
555 template <typename, int, int, typename, bool>
556 friend class MultiComponentVectorView;
557
559 };
560
561
562#ifndef DOXYGEN
563 /*
564 * -------------------------------------------------------------------------
565 * Inline function definitions
566 * -------------------------------------------------------------------------
567 */
568
569
570 template <typename Number, int n_comp, int simd_length>
571 MultiComponentVector<Number, n_comp, simd_length>::MultiComponentVector(
572 const MultiComponentVector &other)
573 {
574 *this = other;
575 }
576
577
578 template <typename Number, int n_comp, int simd_length>
579 MultiComponentVector<Number, n_comp, simd_length>::MultiComponentVector(
580 MultiComponentVector &&other) noexcept
581 {
582 *this = other;
583 }
584
585
586 template <typename Number, int n_comp, int simd_length>
587 void MultiComponentVector<Number, n_comp, simd_length>::
588 reinit_with_vector_partitioner(
589 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
590 &vector_partitioner,
591 const TransferPolicy transfer_policy)
592 {
593 /*
594 * Drop the transfer policy for the duration of the reinit so that a
595 * pinned memory space of a previous policy does not interfere:
596 */
597 this->set_transfer_policy(TransferPolicy::explicit_transfers);
598
599 /* Special case of a zero component vector */
600 if (n_comp == 0) {
601 /* A zero component vector is trivially resident everywhere: */
602 this->reset_residency(/*host*/ true, /*default*/ true);
603 this->set_transfer_policy(transfer_policy);
604 return;
605 }
606
607 host_vector_.reinit(vector_partitioner);
608
609 /*
610 * The vector is resident on the host memory space only. Device
611 * storage is allocated lazily on the first copy_to_memory_space() /
612 * move_to_memory_space(); drop possibly stale device storage from a
613 * previous reinit:
614 */
615 if constexpr (have_separate_memory_spaces)
616 default_vector_.reinit(0);
617
618 this->reset_residency(/*host*/ true,
619 /*default*/ !have_separate_memory_spaces);
620
621 this->set_transfer_policy(transfer_policy);
622 }
623
624
625 template <typename Number, int n_comp, int simd_length>
626 void MultiComponentVector<Number, n_comp, simd_length>::
627 reinit_with_scalar_partitioner(
628 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
629 &scalar_partitioner,
630 const TransferPolicy transfer_policy)
631 {
632 /*
633 * Drop the transfer policy for the duration of the reinit so that a
634 * pinned memory space of a previous policy does not interfere:
635 */
636 this->set_transfer_policy(TransferPolicy::explicit_transfers);
637
638 /* Special case of a zero component vector: */
639 if (n_comp == 0) {
640 /* A zero component vector is trivially resident everywhere: */
641 this->reset_residency(/*host*/ true, /*default*/ true);
642 this->set_transfer_policy(transfer_policy);
643 return;
644 }
645
646 /* Special case of a scalar vector: */
647 if (n_comp == 1)
648 host_vector_.reinit(scalar_partitioner);
649
650 auto vector_partitioner =
651 create_vector_partitioner(scalar_partitioner, n_comp);
652
653 host_vector_.reinit(vector_partitioner);
654
655 /*
656 * The vector is resident on the host memory space only. Device
657 * storage is allocated lazily on the first copy_to_memory_space() /
658 * move_to_memory_space(); drop possibly stale device storage from a
659 * previous reinit:
660 */
661 if constexpr (have_separate_memory_spaces)
662 default_vector_.reinit(0);
663
664 /*
665 * If both memory spaces coincide the single host allocation is
666 * trivially resident on both of them:
667 */
668 this->reset_residency(/*host*/ true,
669 /*default*/ !have_separate_memory_spaces);
670
671 this->set_transfer_policy(transfer_policy);
672 }
673
674
675 template <typename Number, int n_comp, int simd_length>
676 auto MultiComponentVector<Number, n_comp, simd_length>::operator=(
677 const MultiComponentVector &other) -> MultiComponentVector &
678 {
679 /* Copy residency state and transfer policy: */
680 static_cast<MirroredStorage<MultiComponentVector> &>(*this) = other;
681
682 host_vector_ = other.host_vector_;
683 if constexpr (have_separate_memory_spaces)
684 default_vector_ = other.default_vector_;
685
686 return *this;
687 }
688
689
690 template <typename Number, int n_comp, int simd_length>
691 auto MultiComponentVector<Number, n_comp, simd_length>::operator=(
692 MultiComponentVector &&other) noexcept -> MultiComponentVector &
693 {
694 /* Copy residency state and transfer policy: */
695 static_cast<MirroredStorage<MultiComponentVector> &>(*this) = other;
696
697 host_vector_ = std::move(other.host_vector_);
698 if constexpr (have_separate_memory_spaces)
699 default_vector_ = std::move(other.default_vector_);
700
701 return *this;
702 }
703
704
705 template <typename Number, int n_comp, int simd_length>
706 template <typename MemorySpace>
707 MultiComponentVectorView<Number, n_comp, simd_length, MemorySpace, true>
708 MultiComponentVector<Number, n_comp, simd_length>::view()
709 {
710 this->template prepare_write_access<MemorySpace>();
711
712 return MultiComponentVectorView<Number,
713 n_comp,
714 simd_length,
715 MemorySpace,
716 true>(*this);
717 }
718
719
720 template <typename Number, int n_comp, int simd_length>
721 template <typename MemorySpace>
722 MultiComponentVectorView<Number, n_comp, simd_length, MemorySpace, false>
723 MultiComponentVector<Number, n_comp, simd_length>::view() const
724 {
725 this->template prepare_read_access<MemorySpace>();
726
727 return MultiComponentVectorView<Number,
728 n_comp,
729 simd_length,
730 MemorySpace,
731 false>(*this);
732 }
733
734
735 template <typename Number, int n_comp, int simd_length>
736 template <typename MemorySpace>
737 void
738 MultiComponentVector<Number, n_comp, simd_length>::allocate_storage() const
739 {
740 using HostSpace = dealii::MemorySpace::Host;
741
742 /*
743 * The partitioner is recovered from the other (still resident)
744 * vector:
745 */
746
747 if constexpr (std::is_same_v<MemorySpace, HostSpace>) {
748 Assert(default_vector_.size() != 0, dealii::ExcNotInitialized());
749 host_vector_.reinit(default_vector_.get_partitioner());
750 } else {
751 Assert(host_vector_.size() != 0, dealii::ExcNotInitialized());
752 default_vector_.reinit(host_vector_.get_partitioner());
753 }
754 }
755
756
757 template <typename Number, int n_comp, int simd_length>
758 template <typename To, typename From>
759 void
760 MultiComponentVector<Number, n_comp, simd_length>::deep_copy_storage() const
761 {
762 using HostSpace = dealii::MemorySpace::Host;
763
764 if constexpr (std::is_same_v<To, HostSpace>) {
765 host_vector_.import_elements(default_vector_,
766 dealii::VectorOperation::insert);
767 } else {
768 default_vector_.import_elements(host_vector_,
769 dealii::VectorOperation::insert);
770 }
771 }
772
773
774 template <typename Number, int n_comp, int simd_length>
775 template <typename MemorySpace>
776 void MultiComponentVector<Number, n_comp, simd_length>::deallocate_storage()
777 {
778 using HostSpace = dealii::MemorySpace::Host;
779
780 if constexpr (std::is_same_v<MemorySpace, HostSpace>) {
781 host_vector_.reinit(0);
782
783 } else {
784 default_vector_.reinit(0);
785 }
786 }
787
788
789 template <typename Number, int n_components, int simd_length>
790 template <typename MemorySpace>
791 void MultiComponentVector<Number, n_components, simd_length>::
792 zero_out_ghost_values_on_memory_space()
793 {
794 using HostSpace = dealii::MemorySpace::Host;
795 using DefaultSpace = dealii::MemorySpace::Default;
796 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
797 std::is_same_v<MemorySpace, DefaultSpace>,
798 "Unexpected memory space");
799
800 Assert(this->template is_resident<MemorySpace>(),
801 dealii::ExcMessage("The chosen memory space is not resident."));
802
803 if constexpr (have_separate_memory_spaces &&
804 !std::is_same_v<MemorySpace, HostSpace>) {
805 default_vector_.zero_out_ghost_values();
806 } else {
807 host_vector_.zero_out_ghost_values();
808 }
809 }
810
811
812 template <typename Number, int n_components, int simd_length>
813 template <typename MemorySpace>
814 void MultiComponentVector<Number, n_components, simd_length>::
815 update_ghost_values_on_memory_space()
816 {
817 using HostSpace = dealii::MemorySpace::Host;
818 using DefaultSpace = dealii::MemorySpace::Default;
819 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
820 std::is_same_v<MemorySpace, DefaultSpace>,
821 "Unexpected memory space");
822
823 Assert(this->template is_resident<MemorySpace>(),
824 dealii::ExcMessage("The chosen memory space is not resident."));
825
826 if constexpr (have_separate_memory_spaces &&
827 !std::is_same_v<MemorySpace, HostSpace>) {
828 default_vector_.update_ghost_values();
829 } else {
830 host_vector_.update_ghost_values();
831 }
832 }
833
834
835 template <typename Number, int n_components, int simd_length>
836 template <typename MemorySpace>
837 void MultiComponentVector<Number, n_components, simd_length>::
838 compress_on_memory_space(dealii::VectorOperation::values operation)
839 {
840 using HostSpace = dealii::MemorySpace::Host;
841 using DefaultSpace = dealii::MemorySpace::Default;
842 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
843 std::is_same_v<MemorySpace, DefaultSpace>,
844 "Unexpected memory space");
845
846 Assert(this->template is_resident<MemorySpace>(),
847 dealii::ExcMessage("The chosen memory space is not resident."));
848
849 if constexpr (have_separate_memory_spaces &&
850 !std::is_same_v<MemorySpace, HostSpace>) {
851 default_vector_.compress(operation);
852 } else {
853 host_vector_.compress(operation);
854 }
855 }
856
857
858 template <typename Number,
859 int n_comp,
860 int simd_l,
861 typename MemorySpace,
862 bool writable>
863 MultiComponentVectorView<Number, n_comp, simd_l, MemorySpace, writable>::
864 MultiComponentVectorView(MultiComponentVector<Number, n_comp, simd_l>
865 &multi_component_vector)
866 requires(writable)
867 {
868 reinit(multi_component_vector);
869 }
870
871
872 template <typename Number,
873 int n_comp,
874 int simd_l,
875 typename MemorySpace,
876 bool writable>
877 MultiComponentVectorView<Number, n_comp, simd_l, MemorySpace, writable>::
878 MultiComponentVectorView(
879 const MultiComponentVector<Number, n_comp, simd_l>
880 &multi_component_vector)
881 requires(!writable)
882 {
883 reinit(multi_component_vector);
884 }
885
886
887 template <typename Number,
888 int n_comp,
889 int simd_l,
890 typename MemorySpace,
891 bool writable>
892 template <bool other_writable>
893 DEAL_II_HOST_DEVICE
894 MultiComponentVectorView<Number, n_comp, simd_l, MemorySpace, writable>::
895 MultiComponentVectorView(
896 const MultiComponentVectorView<Number,
897 n_comp,
898 simd_l,
899 MemorySpace,
900 other_writable> &other)
901 requires(!writable && other_writable)
902 : multi_component_vector_(other.multi_component_vector_)
903 , data_(other.data_)
904 , n_locally_owned_(other.n_locally_owned_)
905 , n_locally_relevant_(other.n_locally_relevant_)
906 {
907 }
908
909
910 template <typename Number,
911 int n_comp,
912 int simd_l,
913 typename MemorySpace,
914 bool writable>
915 template <typename MultiComponentVector>
916 void
917 MultiComponentVectorView<Number, n_comp, simd_l, MemorySpace, writable>::
918 reinit(MultiComponentVector &multi_component_vector)
919 requires(writable != std::is_const_v<MultiComponentVector>)
920 {
921 using HostSpace = dealii::MemorySpace::Host;
922 using DefaultSpace = dealii::MemorySpace::Default;
923 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
924 std::is_same_v<MemorySpace, DefaultSpace>,
925 "Unexpected memory space");
926
927 multi_component_vector_ = &multi_component_vector;
928
929 if constexpr (have_separate_memory_spaces &&
930 !std::is_same_v<MemorySpace, HostSpace>) {
931 auto &vector = multi_component_vector_->default_vector_;
932 const auto &partitioner = vector.get_partitioner();
933
934 data_ = vector.begin();
935 n_locally_owned_ = partitioner->locally_owned_size();
936 n_locally_relevant_ = n_locally_owned_ + partitioner->n_ghost_indices();
937
938 } else {
939 auto &vector = multi_component_vector_->host_vector_;
940 const auto &partitioner = vector.get_partitioner();
941
942 data_ = vector.begin();
943 n_locally_owned_ = partitioner->locally_owned_size();
944 n_locally_relevant_ = n_locally_owned_ + partitioner->n_ghost_indices();
945 }
946 }
947
948
949 template <typename Number,
950 int n_comp,
951 int simd_length,
952 typename MemorySpace,
953 bool writable>
954 template <typename Functor>
955 void MultiComponentVectorView<
956 Number,
957 n_comp,
958 simd_length,
959 MemorySpace,
960 writable>::extract_component(ScalarVector &scalar_vector,
961 unsigned int component,
962 const Functor &functor) const
963 {
964 using HostSpace = dealii::MemorySpace::Host;
965 AssertThrow((std::is_same_v<MemorySpace, HostSpace>),
966 dealii::ExcNotImplemented());
967
968 Assert(n_comp > 0,
969 dealii::ExcMessage(
970 "Cannot extract from a vector with zero components."));
971 AssertIndexRange(component, n_comp);
972
973 const auto local_size =
974 scalar_vector.get_partitioner()->locally_owned_size();
975
976 Assert(n_comp * local_size == n_locally_owned_,
977 dealii::ExcMessage("Called with a scalar_vector argument that has "
978 "incompatible local range."));
979
980 for (unsigned int i = 0; i < local_size; ++i)
981 scalar_vector.local_element(i) = functor(data_[i * n_comp + component]);
982 scalar_vector.update_ghost_values();
983 }
984
985
986 template <typename Number,
987 int n_comp,
988 int simd_length,
989 typename MemorySpace,
990 bool writable>
991 template <typename Functor>
992 void MultiComponentVectorView<
993 Number,
994 n_comp,
995 simd_length,
996 MemorySpace,
997 writable>::insert_component(const ScalarVector &scalar_vector,
998 unsigned int component,
999 const Functor &functor) const
1000 requires writable
1001 {
1002 using HostSpace = dealii::MemorySpace::Host;
1003 AssertThrow((std::is_same_v<MemorySpace, HostSpace>),
1004 dealii::ExcNotImplemented());
1005
1006 Assert(n_comp > 0,
1007 dealii::ExcMessage(
1008 "Cannot insert into a vector with zero components."));
1009 AssertIndexRange(component, n_comp);
1010
1011 const auto local_size =
1012 scalar_vector.get_partitioner()->locally_owned_size();
1013
1014 Assert(n_comp * local_size == n_locally_owned_,
1015 dealii::ExcMessage("Called with a scalar_vector argument that has "
1016 "incompatible local range."));
1017
1018 for (unsigned int i = 0; i < local_size; ++i)
1019 data_[i * n_comp + component] = functor(scalar_vector.local_element(i));
1020 }
1021
1022
1023 template <typename Number,
1024 int n_comp,
1025 int simd_length,
1026 typename MemorySpace,
1027 bool writable>
1028 template <typename Functor>
1029 void MultiComponentVectorView<
1030 Number,
1031 n_comp,
1032 simd_length,
1033 MemorySpace,
1034 writable>::insert_component(const dealii::Vector<Number> &scalar_vector,
1035 unsigned int component,
1036 const Functor &functor) const
1037 requires writable
1038 {
1039 using HostSpace = dealii::MemorySpace::Host;
1040 AssertThrow((std::is_same_v<MemorySpace, HostSpace>),
1041 dealii::ExcInternalError());
1042
1043 Assert(n_comp > 0,
1044 dealii::ExcMessage(
1045 "Cannot insert into a vector with zero components."));
1046 AssertIndexRange(component, n_comp);
1047
1048 const auto local_size = scalar_vector.size();
1049
1050 Assert(n_comp * local_size >= n_locally_owned_,
1051 dealii::ExcMessage("Called with a scalar_vector argument that has "
1052 "incompatible local range."));
1053
1054 for (unsigned int i = 0; i < local_size; ++i)
1055 data_[i * n_comp + component] = functor(scalar_vector[i]);
1056 }
1057
1058
1059 template <typename Num, int n_comp, int simd_l, typename MS, bool writable>
1060 void MultiComponentVectorView<Num, n_comp, simd_l, MS, writable>::sadd(
1061 const Num s,
1062 const Num a,
1063 const MultiComponentVectorView<Num,
1064 n_comp,
1065 simd_l,
1066 MS,
1067 /*writable=*/false> &v) const
1068 requires writable
1069 {
1070 using HS = dealii::MemorySpace::Host;
1071 using DS = dealii::MemorySpace::Default;
1072 static_assert(std::is_same_v<MS, HS> || std::is_same_v<MS, DS>,
1073 "Unexpected memory space");
1074
1075 /*
1076 * Note: If the host and default memory spaces coincide all views
1077 * reference the host storage.
1078 */
1079 if constexpr (have_separate_memory_spaces && !std::is_same_v<MS, HS>) {
1080 const auto &other = v.multi_component_vector_->default_vector_;
1081 multi_component_vector_->default_vector_.sadd(s, a, other);
1082 } else {
1083 const auto &other = v.multi_component_vector_->host_vector_;
1084 multi_component_vector_->host_vector_.sadd(s, a, other);
1085 }
1086 }
1087
1088
1089 template <typename Number,
1090 int n_comp,
1091 int simd_length,
1092 typename MemorySpace,
1093 bool writable>
1094 template <typename Number2>
1095 DEAL_II_HOST_DEVICE_ALWAYS_INLINE Number2
1096 MultiComponentVectorView<Number,
1097 n_comp,
1098 simd_length,
1099 MemorySpace,
1100 writable>::read_entry(const unsigned int i) const
1101 {
1102 static_assert(
1103 n_comp == 1,
1104 "Attempted to read a scalar value from a tensor-valued vector entry");
1105
1106 AssertIndexRange(i, n_locally_relevant_);
1107
1108 const auto result = read_tensor<Number2>(i);
1109 return result[0];
1110 }
1111
1112
1113 template <typename Number,
1114 int n_comp,
1115 int simd_length,
1116 typename MemorySpace,
1117 bool writable>
1118 template <typename Number2, typename Tensor>
1119 DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor
1120 MultiComponentVectorView<Number,
1121 n_comp,
1122 simd_length,
1123 MemorySpace,
1124 writable>::read_tensor(const unsigned int i) const
1125 {
1126 static_assert(std::is_same_v<Number2, typename Tensor::value_type>,
1127 "type mismatch");
1128
1129 AssertIndexRange(i, n_locally_relevant_);
1130
1131 Tensor tensor;
1132
1133 /* Special case of a zero component vector */
1134 if constexpr (n_comp == 0)
1135 return tensor;
1136
1137 using VA = dealii::VectorizedArray<Number>;
1138 if constexpr (std::is_same_v<VA, Number2>) {
1139 /* Vectorized fast access. index must be divisible by simd_length */
1140 std::array<unsigned int, VA::size()> indices;
1141 for (unsigned int k = 0; k < VA::size(); ++k)
1142 indices[k] = k * n_comp;
1143
1144 dealii::vectorized_load_and_transpose(
1145 n_comp, data_ + i * n_comp, indices.data(), &tensor[0]);
1146
1147 } else {
1148 /* Non-vectorized sequential access. */
1149 for (unsigned int d = 0; d < n_comp; ++d)
1150 tensor[d] = data_[i * n_comp + d];
1151 }
1152
1153 return tensor;
1154 }
1155
1156
1157 template <typename Number,
1158 int n_comp,
1159 int simd_length,
1160 typename MemorySpace,
1161 bool writable>
1162 template <typename Number2>
1163 DEAL_II_HOST_DEVICE_ALWAYS_INLINE Number2
1164 MultiComponentVectorView<Number,
1165 n_comp,
1166 simd_length,
1167 MemorySpace,
1168 writable>::read_entry(const unsigned int *js) const
1169 {
1170 static_assert(
1171 n_comp == 1,
1172 "Attempted to read a scalar value from a tensor-valued vector entry");
1173
1174 const auto result = read_tensor<Number2>(js);
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 *js)
1191 const
1192 {
1193 static_assert(std::is_same_v<Number2, typename Tensor::value_type>,
1194 "type mismatch");
1195
1196 Tensor tensor;
1197
1198 /* Special case of a zero component vector */
1199 if constexpr (n_comp == 0)
1200 return tensor;
1201
1202 using VA = dealii::VectorizedArray<Number>;
1203 if constexpr (std::is_same_v<VA, Number2>) {
1204 /* Vectorized fast access. index must be divisible by simd_length */
1205
1206 std::array<unsigned int, VA::size()> indices;
1207 for (unsigned int k = 0; k < VA::size(); ++k) {
1208 AssertIndexRange(js[k], n_locally_relevant_);
1209 indices[k] = js[k] * n_comp;
1210 }
1211
1212 dealii::vectorized_load_and_transpose(
1213 n_comp, data_, indices.data(), &tensor[0]);
1214
1215 } else {
1216 /* Non-vectorized sequential access. */
1217
1218 AssertIndexRange(*js, n_locally_relevant_);
1219
1220 for (unsigned int d = 0; d < n_comp; ++d)
1221 tensor[d] = data_[js[0] * n_comp + d];
1222 }
1223
1224 return tensor;
1225 }
1226
1227
1228 template <typename Number,
1229 int n_comp,
1230 int simd_length,
1231 typename MemorySpace,
1232 bool writable>
1233 template <typename Number2>
1234 DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1235 MultiComponentVectorView<Number,
1236 n_comp,
1237 simd_length,
1238 MemorySpace,
1239 writable>::write_entry(const Number2 &entry,
1240 const unsigned int i) const
1241 requires writable
1242 {
1243 static_assert(n_comp == 1,
1244 "Attempted to write a scalar value into a tensor-valued "
1245 "vector entry");
1246
1247 AssertIndexRange(i, n_locally_relevant_);
1248
1249 dealii::Tensor<1, n_comp, Number2> tensor;
1250 tensor[0] = entry;
1251
1252 write_tensor<Number2>(tensor, i);
1253 }
1254
1255
1256 template <typename Number,
1257 int n_comp,
1258 int simd_length,
1259 typename MemorySpace,
1260 bool writable>
1261 template <typename Number2, typename Tensor>
1262 DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1263 MultiComponentVectorView<Number,
1264 n_comp,
1265 simd_length,
1266 MemorySpace,
1267 writable>::write_tensor(const Tensor &tensor,
1268 const unsigned int i) const
1269 requires writable
1270 {
1271 static_assert(std::is_same_v<Number2, typename Tensor::value_type>,
1272 "type mismatch");
1273
1274 AssertIndexRange(i, n_locally_relevant_);
1275
1276 /* Special case of a zero component vector */
1277 if constexpr (n_comp == 0)
1278 return;
1279
1280 using VA = dealii::VectorizedArray<Number>;
1281 if constexpr (std::is_same_v<VA, Number2>) {
1282 /* Vectorized fast access. index must be divisible by simd_length */
1283
1284 std::array<unsigned int, VA::size()> indices;
1285 for (unsigned int k = 0; k < VA::size(); ++k)
1286 indices[k] = k * n_comp;
1287
1288 dealii::vectorized_transpose_and_store(/*add into*/ false,
1289 n_comp,
1290 &tensor[0],
1291 indices.data(),
1292 data_ + i * n_comp);
1293
1294 } else {
1295 /* Non-vectorized sequential access. */
1296
1297 for (unsigned int d = 0; d < n_comp; ++d)
1298 data_[i * n_comp + d] = tensor[d];
1299 }
1300 }
1301
1302
1303 template <typename Number,
1304 int n_comp,
1305 int simd_length,
1306 typename MemorySpace,
1307 bool writable>
1308 template <typename Number2>
1309 DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1310 MultiComponentVectorView<Number,
1311 n_comp,
1312 simd_length,
1313 MemorySpace,
1314 writable>::add_entry(const Number2 &entry,
1315 const unsigned int i) const
1316 requires writable
1317 {
1318 static_assert(n_comp == 1,
1319 "Attempted to write a scalar value into a tensor-valued "
1320 "matrix entry");
1321
1322 AssertIndexRange(i, n_locally_relevant_);
1323
1324 dealii::Tensor<1, n_comp, Number2> tensor;
1325 tensor[0] = entry;
1326
1327 add_tensor<Number2>(tensor, i);
1328 }
1329
1330
1331 template <typename Number,
1332 int n_comp,
1333 int simd_length,
1334 typename MemorySpace,
1335 bool writable>
1336 template <typename Number2, typename Tensor>
1337 DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1338 MultiComponentVectorView<Number,
1339 n_comp,
1340 simd_length,
1341 MemorySpace,
1342 writable>::add_tensor(const Tensor &tensor,
1343 const unsigned int i) const
1344 requires writable
1345 {
1346 static_assert(std::is_same_v<Number2, typename Tensor::value_type>,
1347 "type mismatch");
1348
1349 AssertIndexRange(i, n_locally_relevant_);
1350
1351 /* Special case of a zero component vector */
1352 if constexpr (n_comp == 0)
1353 return;
1354
1355 using VA = dealii::VectorizedArray<Number>;
1356 if constexpr (std::is_same_v<VA, Number2>) {
1357 /* Vectorized fast access. index must be divisible by simd_length */
1358
1359 std::array<unsigned int, VA::size()> indices;
1360 for (unsigned int k = 0; k < VA::size(); ++k)
1361 indices[k] = k * n_comp;
1362
1363 dealii::vectorized_transpose_and_store(/*add into*/ true,
1364 n_comp,
1365 &tensor[0],
1366 indices.data(),
1367 data_ + i * n_comp);
1368
1369 } else {
1370 /* Non-vectorized sequential access. */
1371
1372 for (unsigned int d = 0; d < n_comp; ++d)
1373 data_[i * n_comp + d] += tensor[d];
1374 }
1375 }
1376
1377
1378 template <typename Number,
1379 int n_comp,
1380 int simd_l,
1381 typename MemorySpace,
1382 bool writable>
1383 void
1384 MultiComponentVectorView<Number, n_comp, simd_l, MemorySpace, writable>::
1385 zero_out_ghost_values() const
1386 requires(writable)
1387 {
1388 using HostSpace = dealii::MemorySpace::Host;
1389 using DefaultSpace = dealii::MemorySpace::Default;
1390 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
1391 std::is_same_v<MemorySpace, DefaultSpace>,
1392 "Unexpected memory space");
1393
1394 Assert(multi_component_vector_->template is_resident<MemorySpace>(),
1395 dealii::ExcMessage("The chosen memory space is not resident."));
1396
1397 multi_component_vector_
1398 ->template zero_out_ghost_values_on_memory_space<MemorySpace>();
1399 }
1400
1401
1402 template <typename Number,
1403 int n_comp,
1404 int simd_l,
1405 typename MemorySpace,
1406 bool writable>
1407 void
1408 MultiComponentVectorView<Number, n_comp, simd_l, MemorySpace, writable>::
1409 update_ghost_values() const
1410 requires(writable)
1411 {
1412 using HostSpace = dealii::MemorySpace::Host;
1413 using DefaultSpace = dealii::MemorySpace::Default;
1414 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
1415 std::is_same_v<MemorySpace, DefaultSpace>,
1416 "Unexpected memory space");
1417
1418 Assert(multi_component_vector_->template is_resident<MemorySpace>(),
1419 dealii::ExcMessage("The chosen memory space is not resident."));
1420
1421 multi_component_vector_
1422 ->template update_ghost_values_on_memory_space<MemorySpace>();
1423 }
1424
1425
1426 template <typename Number,
1427 int n_comp,
1428 int simd_l,
1429 typename MemorySpace,
1430 bool writable>
1431 void
1432 MultiComponentVectorView<Number, n_comp, simd_l, MemorySpace, writable>::
1433 compress(dealii::VectorOperation::values operation) const
1434 requires(writable)
1435 {
1436 using HostSpace = dealii::MemorySpace::Host;
1437 using DefaultSpace = dealii::MemorySpace::Default;
1438 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
1439 std::is_same_v<MemorySpace, DefaultSpace>,
1440 "Unexpected memory space");
1441
1442 Assert(multi_component_vector_->template is_resident<MemorySpace>(),
1443 dealii::ExcMessage("The chosen memory space is not resident."));
1444
1445 multi_component_vector_->template compress_on_memory_space<MemorySpace>(
1446 operation);
1447 }
1448
1449#endif
1450 } // namespace Vectors
1451} // namespace ryujin
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:305
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::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)
Generated by doxygen 1.9.8