Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines class Variables
6 :
7 : #pragma once
8 :
9 : #include <algorithm>
10 : #include <array>
11 : #include <blaze/math/AlignmentFlag.h>
12 : #include <blaze/math/CustomVector.h>
13 : #include <blaze/math/DenseVector.h>
14 : #include <blaze/math/GroupTag.h>
15 : #include <blaze/math/PaddingFlag.h>
16 : #include <blaze/math/TransposeFlag.h>
17 : #include <blaze/math/Vector.h>
18 : #include <limits>
19 : #include <memory>
20 : #include <ostream>
21 : #include <pup.h>
22 : #include <string>
23 : #include <tuple>
24 :
25 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
26 : #include "DataStructures/DataBox/Subitems.hpp"
27 : #include "DataStructures/DataBox/Tag.hpp"
28 : #include "DataStructures/DataBox/TagName.hpp"
29 : #include "DataStructures/DataBox/TagTraits.hpp"
30 : #include "DataStructures/DataVector.hpp"
31 : #include "DataStructures/MathWrapper.hpp"
32 : #include "DataStructures/Tensor/IndexType.hpp"
33 : #include "DataStructures/Tensor/Tensor.hpp"
34 : #include "Utilities/EqualWithinRoundoff.hpp"
35 : #include "Utilities/ErrorHandling/Assert.hpp"
36 : #include "Utilities/ForceInline.hpp"
37 : #include "Utilities/Gsl.hpp"
38 : #include "Utilities/Literals.hpp"
39 : #include "Utilities/MakeSignalingNan.hpp"
40 : #include "Utilities/MemoryHelpers.hpp"
41 : #include "Utilities/PrettyType.hpp"
42 : #include "Utilities/Requires.hpp"
43 : #include "Utilities/TMPL.hpp"
44 : #include "Utilities/TaggedTuple.hpp"
45 : #include "Utilities/TypeTraits.hpp"
46 : #include "Utilities/TypeTraits/IsA.hpp"
47 :
48 : // IWYU pragma: no_forward_declare MakeWithValueImpl
49 : // IWYU pragma: no_forward_declare Variables
50 :
51 : /// \cond
52 : template <typename X, typename Symm, typename IndexList>
53 : class Tensor;
54 : template <typename TagsList>
55 : class Variables;
56 : template <typename T, typename VectorType, size_t StaticSize>
57 : class VectorImpl;
58 : namespace Tags {
59 : template <typename TagsList>
60 : class Variables;
61 : } // namespace Tags
62 : /// \endcond
63 :
64 : /*!
65 : * \ingroup DataStructuresGroup
66 : * \brief A Variables holds a contiguous memory block with Tensors pointing
67 : * into it.
68 : *
69 : * The `Tags` are `struct`s that must have a public type alias `type` whose
70 : * value must be a `Tensor<DataVector, ...>`, a `static' method `name()` that
71 : * returns a `std::string` of the tag name, and must derive off of
72 : * `db::SimpleTag`. In general, they should be DataBoxTags that are not compute
73 : * items. For example,
74 : *
75 : * \snippet Helpers/DataStructures/TestTags.hpp simple_variables_tag
76 : *
77 : * Prefix tags can also be stored and their format is:
78 : *
79 : * \snippet Helpers/DataStructures/TestTags.hpp prefix_variables_tag
80 : *
81 : * #### Design Decisions
82 : *
83 : * The `Variables` class is designed to hold several different `Tensor`s
84 : * performing one memory allocation for all the `Tensor`s. The advantage is that
85 : * memory allocations are quite expensive, especially in a parallel environment.
86 : *
87 : * In Debug mode, or if the macro `SPECTRE_NAN_INIT` is defined, the contents
88 : * are initialized with `NaN`s.
89 : *
90 : * `Variables` stores the data it owns in a `std::unique_ptr<double[]>`
91 : * instead of a `std::vector` because `std::vector` value-initializes its
92 : * contents, which is very slow.
93 : */
94 : template <typename... Tags>
95 1 : class Variables<tmpl::list<Tags...>> {
96 : public:
97 0 : using size_type = size_t;
98 0 : using difference_type = std::ptrdiff_t;
99 0 : static constexpr auto transpose_flag = blaze::defaultTransposeFlag;
100 :
101 : /// A typelist of the Tags whose variables are held
102 1 : using tags_list = tmpl::list<Tags...>;
103 : static_assert(sizeof...(Tags) > 0,
104 : "You must provide at least one tag to the Variables "
105 : "for type inference");
106 : static_assert((db::is_simple_tag_v<Tags> and ...));
107 : static_assert(tmpl2::flat_all_v<tt::is_a_v<Tensor, typename Tags::type>...>);
108 :
109 : private:
110 0 : using first_tensors_type = typename tmpl::front<tags_list>::type::type;
111 :
112 : public:
113 : static_assert(tmpl2::flat_all_v<std::is_same_v<typename Tags::type::type,
114 : first_tensors_type>...> or
115 : tmpl2::flat_all_v<is_spin_weighted_of_same_type_v<
116 : typename Tags::type::type, first_tensors_type>...>,
117 : "All tensors stored in a single Variables must "
118 : "have the same internal storage type.");
119 :
120 0 : using vector_type =
121 : tmpl::conditional_t<is_any_spin_weighted_v<first_tensors_type>,
122 : typename first_tensors_type::value_type,
123 : first_tensors_type>;
124 0 : using value_type = typename vector_type::value_type;
125 0 : using pointer = value_type*;
126 0 : using const_pointer = const value_type*;
127 0 : using allocator_type = std::allocator<value_type>;
128 0 : using pointer_type =
129 : blaze::CustomVector<value_type, blaze::AlignmentFlag::unaligned,
130 : blaze::PaddingFlag::unpadded, transpose_flag,
131 : blaze::GroupTag<0>, vector_type>;
132 :
133 : static_assert(
134 : std::is_fundamental_v<value_type> or tt::is_a_v<std::complex, value_type>,
135 : "`value_type` of the Variables (so the storage type of the vector type "
136 : "within the tensors in the Variables) must be either a fundamental type "
137 : "or a std::complex. If this constraint is relaxed, the value_type "
138 : "should be handled differently in the Variables, including pass by "
139 : "reference.");
140 :
141 : /// The number of variables of the Variables object is holding. E.g.
142 : /// \f$\psi_{ab}\f$ would be counted as one variable.
143 1 : static constexpr auto number_of_variables = sizeof...(Tags);
144 :
145 : /// The total number of independent components of all the variables. E.g.
146 : /// a rank-2 symmetric spacetime Tensor \f$\psi_{ab}\f$ in 3 spatial
147 : /// dimensions would have 10 independent components.
148 1 : static constexpr size_t number_of_independent_components =
149 : (... + Tags::type::size());
150 :
151 : /// Default construct an empty Variables class, Charm++ needs this
152 1 : Variables();
153 :
154 0 : explicit Variables(size_t number_of_grid_points);
155 :
156 0 : Variables(size_t number_of_grid_points, value_type value);
157 :
158 : /// Construct a non-owning Variables that points to `start`. `size` is the
159 : /// size of the allocation, which must be
160 : /// `number_of_grid_points * Variables::number_of_independent_components`
161 1 : Variables(pointer start, size_t size);
162 :
163 0 : Variables(Variables&& rhs);
164 0 : Variables& operator=(Variables&& rhs);
165 :
166 0 : Variables(const Variables& rhs);
167 0 : Variables& operator=(const Variables& rhs);
168 :
169 : /// @{
170 : /// Copy and move semantics for wrapped variables
171 : template <typename... WrappedTags,
172 : Requires<tmpl2::flat_all<std::is_same<
173 : db::remove_all_prefixes<WrappedTags>,
174 : db::remove_all_prefixes<Tags>>::value...>::value> = nullptr>
175 1 : explicit Variables(Variables<tmpl::list<WrappedTags...>>&& rhs);
176 : template <typename... WrappedTags,
177 : Requires<tmpl2::flat_all<std::is_same<
178 : db::remove_all_prefixes<WrappedTags>,
179 : db::remove_all_prefixes<Tags>>::value...>::value> = nullptr>
180 1 : Variables& operator=(Variables<tmpl::list<WrappedTags...>>&& rhs);
181 :
182 : template <typename... WrappedTags,
183 : Requires<tmpl2::flat_all<std::is_same<
184 : db::remove_all_prefixes<WrappedTags>,
185 : db::remove_all_prefixes<Tags>>::value...>::value> = nullptr>
186 1 : explicit Variables(const Variables<tmpl::list<WrappedTags...>>& rhs);
187 : template <typename... WrappedTags,
188 : Requires<tmpl2::flat_all<std::is_same<
189 : db::remove_all_prefixes<WrappedTags>,
190 : db::remove_all_prefixes<Tags>>::value...>::value> = nullptr>
191 1 : Variables& operator=(const Variables<tmpl::list<WrappedTags...>>& rhs);
192 : /// @}
193 :
194 : /// \cond HIDDEN_SYMBOLS
195 : ~Variables() = default;
196 : /// \endcond
197 :
198 : /// @{
199 : /// Initialize a Variables to the state it would have after calling
200 : /// the constructor with the same arguments.
201 : // this should be updated if we ever use a variables which has a `value_type`
202 : // larger than ~2 doubles in size.
203 1 : void initialize(size_t number_of_grid_points);
204 1 : void initialize(size_t number_of_grid_points, value_type value);
205 : /// @}
206 :
207 : /// @{
208 : /// Set the VectorImpl to be a reference to another VectorImpl object
209 1 : void set_data_ref(const gsl::not_null<Variables*> rhs) {
210 : set_data_ref(rhs->data(), rhs->size());
211 : }
212 :
213 1 : void set_data_ref(pointer const start, const size_t size) {
214 : variable_data_impl_dynamic_.reset();
215 : if (start == nullptr) {
216 : variable_data_ = pointer_type{};
217 : size_ = 0;
218 : number_of_grid_points_ = 0;
219 : } else {
220 : size_ = size;
221 : variable_data_.reset(start, size_);
222 : ASSERT(
223 : size_ % number_of_independent_components == 0,
224 : "The size (" << size_
225 : << ") must be a multiple of the number of independent "
226 : "components ("
227 : << number_of_independent_components
228 : << ") since we calculate the number of grid points from "
229 : "the size and number of independent components.");
230 : number_of_grid_points_ = size_ / number_of_independent_components;
231 : }
232 : owning_ = false;
233 : add_reference_variable_data();
234 : }
235 : /// @}
236 :
237 0 : constexpr SPECTRE_ALWAYS_INLINE size_t number_of_grid_points() const {
238 : return number_of_grid_points_;
239 : }
240 :
241 : /// Number of grid points * number of independent components
242 1 : constexpr SPECTRE_ALWAYS_INLINE size_type size() const { return size_; }
243 :
244 : /// @{
245 : /// Access pointer to underlying data
246 1 : pointer data() { return variable_data_.data(); }
247 1 : const_pointer data() const { return variable_data_.data(); }
248 : /// @}
249 :
250 : /// \cond HIDDEN_SYMBOLS
251 : /// Needed because of limitations and inconsistency between compiler
252 : /// implementations of friend function templates with auto return type of
253 : /// class templates
254 : const auto& get_variable_data() const { return variable_data_; }
255 : /// \endcond
256 :
257 : /// Returns true if the class owns the data
258 1 : bool is_owning() const { return owning_; }
259 :
260 : // clang-tidy: redundant-declaration
261 : template <typename Tag, typename TagList>
262 1 : friend constexpr typename Tag::type& get( // NOLINT
263 : Variables<TagList>& v);
264 : template <typename Tag, typename TagList>
265 0 : friend constexpr const typename Tag::type& get( // NOLINT
266 : const Variables<TagList>& v);
267 :
268 : /// Serialization for Charm++.
269 : // NOLINTNEXTLINE(google-runtime-references)
270 1 : void pup(PUP::er& p);
271 :
272 : /// @{
273 : /// \brief Assign a subset of the `Tensor`s from another Variables or a
274 : /// tuples::TaggedTuple. Any tags that aren't in both containers are
275 : /// ignored.
276 : ///
277 : /// \note There is no need for an rvalue overload because we need to copy into
278 : /// the contiguous array anyway
279 : template <typename... SubsetOfTags>
280 1 : void assign_subset(const Variables<tmpl::list<SubsetOfTags...>>& vars) {
281 : EXPAND_PACK_LEFT_TO_RIGHT([this, &vars]() {
282 : if constexpr (tmpl::list_contains_v<tmpl::list<Tags...>, SubsetOfTags>) {
283 : get<SubsetOfTags>(*this) = get<SubsetOfTags>(vars);
284 : } else {
285 : (void)this;
286 : (void)vars;
287 : }
288 : }());
289 : }
290 :
291 : template <typename... SubsetOfTags>
292 1 : void assign_subset(const tuples::TaggedTuple<SubsetOfTags...>& vars) {
293 : EXPAND_PACK_LEFT_TO_RIGHT([this, &vars]() {
294 : if constexpr (tmpl::list_contains_v<tmpl::list<Tags...>, SubsetOfTags>) {
295 : get<SubsetOfTags>(*this) = get<SubsetOfTags>(vars);
296 : } else {
297 : (void)this;
298 : (void)vars;
299 : }
300 : }());
301 : }
302 : /// @}
303 :
304 : /// Create a Variables from a subset of the `Tensor`s in this
305 : /// Variables
306 : template <typename SubsetOfTags>
307 1 : Variables<SubsetOfTags> extract_subset() const {
308 : Variables<SubsetOfTags> sub_vars(number_of_grid_points());
309 : tmpl::for_each<SubsetOfTags>([&](const auto tag_v) {
310 : using tag = tmpl::type_from<decltype(tag_v)>;
311 : get<tag>(sub_vars) = get<tag>(*this);
312 : });
313 : return sub_vars;
314 : }
315 :
316 : /// Create a non-owning Variables referencing a subset of the
317 : /// `Tensor`s in this Variables. The referenced tensors must be
318 : /// consecutive in this Variables's tags list.
319 : ///
320 : /// \warning As with other appearances of non-owning variables, this
321 : /// method can cast away constness.
322 : template <typename SubsetOfTags>
323 1 : Variables<SubsetOfTags> reference_subset() const {
324 : if constexpr (std::is_same_v<SubsetOfTags, tmpl::list<>>) {
325 : return {};
326 : } else {
327 : using subset_from_tags_list = tmpl::reverse_find<
328 : tmpl::find<
329 : tags_list,
330 : std::is_same<tmpl::_1, tmpl::pin<tmpl::front<SubsetOfTags>>>>,
331 : std::is_same<tmpl::_1, tmpl::pin<tmpl::back<SubsetOfTags>>>>;
332 : static_assert(std::is_same_v<subset_from_tags_list, SubsetOfTags>,
333 : "Tags passed to reference_subset must be consecutive tags "
334 : "in the original tags_list.");
335 :
336 : using preceeding_tags = tmpl::pop_back<tmpl::reverse_find<
337 : tags_list,
338 : std::is_same<tmpl::_1, tmpl::pin<tmpl::front<SubsetOfTags>>>>>;
339 :
340 : static constexpr auto count_components = [](auto... tags) {
341 : return (0_st + ... + tmpl::type_from<decltype(tags)>::type::size());
342 : };
343 : static constexpr size_t number_of_preceeding_components =
344 : tmpl::as_pack<preceeding_tags>(count_components);
345 : static constexpr size_t number_of_components_in_subset =
346 : tmpl::as_pack<SubsetOfTags>(count_components);
347 :
348 : return {const_cast<value_type*>(data()) +
349 : number_of_grid_points() * number_of_preceeding_components,
350 : number_of_grid_points() * number_of_components_in_subset};
351 : }
352 : }
353 :
354 : /// Create a non-owning version of this Variables with different
355 : /// prefixes on the tensors. Both sets of prefixes must have the
356 : /// same tensor types.
357 : ///
358 : /// \warning As with other appearances of non-owning variables, this
359 : /// method can cast away constness.
360 : template <typename WrappedVariables>
361 1 : WrappedVariables reference_with_different_prefixes() const {
362 : static_assert(
363 : tmpl::all<tmpl::transform<
364 : typename WrappedVariables::tags_list, tmpl::list<Tags...>,
365 : std::is_same<tmpl::bind<db::remove_all_prefixes, tmpl::_1>,
366 : tmpl::bind<db::remove_all_prefixes, tmpl::_2>>>>::
367 : value,
368 : "Unprefixed tags do not match!");
369 : static_assert(
370 : tmpl::all<tmpl::transform<
371 : typename WrappedVariables::tags_list, tmpl::list<Tags...>,
372 : std::is_same<tmpl::bind<tmpl::type_from, tmpl::_1>,
373 : tmpl::bind<tmpl::type_from, tmpl::_2>>>>::value,
374 : "Tensor types do not match!");
375 : return {const_cast<value_type*>(data()), size()};
376 : }
377 :
378 : /// Converting constructor for an expression to a Variables class
379 : // clang-tidy: mark as explicit (we want conversion to Variables)
380 : template <typename VT, bool VF>
381 1 : Variables(const blaze::DenseVector<VT, VF>& expression); // NOLINT
382 :
383 : template <typename VT, bool VF>
384 0 : Variables& operator=(const blaze::DenseVector<VT, VF>& expression);
385 :
386 : // Prevent the previous two declarations from implicitly converting
387 : // DataVectors and similar to Variables.
388 : template <typename T, typename VectorType, size_t StaticSize>
389 0 : Variables(const VectorImpl<T, VectorType, StaticSize>&) = delete;
390 : template <typename T, typename VectorType, size_t StaticSize>
391 0 : Variables& operator=(const VectorImpl<T, VectorType, StaticSize>&) = delete;
392 :
393 : template <typename... WrappedTags,
394 : Requires<tmpl2::flat_all<std::is_same_v<
395 : db::remove_all_prefixes<WrappedTags>,
396 : db::remove_all_prefixes<Tags>>...>::value> = nullptr>
397 0 : SPECTRE_ALWAYS_INLINE Variables& operator+=(
398 : const Variables<tmpl::list<WrappedTags...>>& rhs) {
399 : static_assert(
400 : (std::is_same_v<typename Tags::type, typename WrappedTags::type> and
401 : ...),
402 : "Tensor types do not match!");
403 : variable_data_ += rhs.variable_data_;
404 : return *this;
405 : }
406 : template <typename VT, bool VF>
407 0 : SPECTRE_ALWAYS_INLINE Variables& operator+=(
408 : const blaze::Vector<VT, VF>& rhs) {
409 : variable_data_ += rhs;
410 : return *this;
411 : }
412 :
413 : template <typename... WrappedTags,
414 : Requires<tmpl2::flat_all<std::is_same_v<
415 : db::remove_all_prefixes<WrappedTags>,
416 : db::remove_all_prefixes<Tags>>...>::value> = nullptr>
417 0 : SPECTRE_ALWAYS_INLINE Variables& operator-=(
418 : const Variables<tmpl::list<WrappedTags...>>& rhs) {
419 : static_assert(
420 : (std::is_same_v<typename Tags::type, typename WrappedTags::type> and
421 : ...),
422 : "Tensor types do not match!");
423 : variable_data_ -= rhs.variable_data_;
424 : return *this;
425 : }
426 : template <typename VT, bool VF>
427 0 : SPECTRE_ALWAYS_INLINE Variables& operator-=(
428 : const blaze::Vector<VT, VF>& rhs) {
429 : variable_data_ -= rhs;
430 : return *this;
431 : }
432 :
433 0 : SPECTRE_ALWAYS_INLINE Variables& operator*=(const value_type& rhs) {
434 : variable_data_ *= rhs;
435 : return *this;
436 : }
437 :
438 0 : SPECTRE_ALWAYS_INLINE Variables& operator/=(const value_type& rhs) {
439 : variable_data_ /= rhs;
440 : return *this;
441 : }
442 :
443 : template <typename... WrappedTags,
444 : Requires<tmpl2::flat_all<std::is_same_v<
445 : db::remove_all_prefixes<WrappedTags>,
446 : db::remove_all_prefixes<Tags>>...>::value> = nullptr>
447 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator+(
448 : const Variables<tmpl::list<WrappedTags...>>& lhs, const Variables& rhs) {
449 : static_assert(
450 : (std::is_same_v<typename Tags::type, typename WrappedTags::type> and
451 : ...),
452 : "Tensor types do not match!");
453 : return lhs.get_variable_data() + rhs.variable_data_;
454 : }
455 : template <typename VT, bool VF>
456 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator+(
457 : const blaze::DenseVector<VT, VF>& lhs, const Variables& rhs) {
458 : return *lhs + rhs.variable_data_;
459 : }
460 : template <typename VT, bool VF>
461 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator+(
462 : const Variables& lhs, const blaze::DenseVector<VT, VF>& rhs) {
463 : return lhs.variable_data_ + *rhs;
464 : }
465 :
466 : template <typename... WrappedTags,
467 : Requires<tmpl2::flat_all<std::is_same_v<
468 : db::remove_all_prefixes<WrappedTags>,
469 : db::remove_all_prefixes<Tags>>...>::value> = nullptr>
470 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator-(
471 : const Variables<tmpl::list<WrappedTags...>>& lhs, const Variables& rhs) {
472 : static_assert(
473 : (std::is_same_v<typename Tags::type, typename WrappedTags::type> and
474 : ...),
475 : "Tensor types do not match!");
476 : return lhs.get_variable_data() - rhs.variable_data_;
477 : }
478 : template <typename VT, bool VF>
479 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator-(
480 : const blaze::DenseVector<VT, VF>& lhs, const Variables& rhs) {
481 : return *lhs - rhs.variable_data_;
482 : }
483 : template <typename VT, bool VF>
484 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator-(
485 : const Variables& lhs, const blaze::DenseVector<VT, VF>& rhs) {
486 : return lhs.variable_data_ - *rhs;
487 : }
488 :
489 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator*(const Variables& lhs,
490 : const value_type& rhs) {
491 : return lhs.variable_data_ * rhs;
492 : }
493 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator*(const value_type& lhs,
494 : const Variables& rhs) {
495 : return lhs * rhs.variable_data_;
496 : }
497 :
498 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator/(const Variables& lhs,
499 : const value_type& rhs) {
500 : return lhs.variable_data_ / rhs;
501 : }
502 :
503 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator-(const Variables& lhs) {
504 : return -lhs.variable_data_;
505 : }
506 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator+(const Variables& lhs) {
507 : return lhs.variable_data_;
508 : }
509 :
510 : private:
511 : /// @{
512 : /*!
513 : * \brief Subscript operator
514 : *
515 : * The subscript operator is private since it should not be used directly.
516 : * Mathematical operations should be done using the math operators provided.
517 : * Since the internal ordering of variables is implementation defined there
518 : * is no safe way to perform any operation that is not a linear combination of
519 : * Variables. Retrieving a Tensor must be done via the `get()` function.
520 : *
521 : * \requires `i >= 0 and i < size()`
522 : */
523 1 : SPECTRE_ALWAYS_INLINE value_type& operator[](const size_type i) {
524 : return variable_data_[i];
525 : }
526 1 : SPECTRE_ALWAYS_INLINE const value_type& operator[](const size_type i) const {
527 : return variable_data_[i];
528 : }
529 : /// @}
530 :
531 0 : void add_reference_variable_data();
532 :
533 0 : friend bool operator==(const Variables& lhs, const Variables& rhs) {
534 : return blaze::equal<blaze::strict>(lhs.variable_data_, rhs.variable_data_);
535 : }
536 :
537 : template <typename VT, bool TF>
538 0 : friend bool operator==(const Variables& lhs,
539 : const blaze::DenseVector<VT, TF>& rhs) {
540 : return blaze::equal<blaze::strict>(lhs.variable_data_, *rhs);
541 : }
542 :
543 : template <typename VT, bool TF>
544 0 : friend bool operator==(const blaze::DenseVector<VT, TF>& lhs,
545 : const Variables& rhs) {
546 : return blaze::equal<blaze::strict>(*lhs, rhs.variable_data_);
547 : }
548 :
549 : template <class FriendTags>
550 0 : friend class Variables;
551 :
552 : std::array<value_type, number_of_independent_components>
553 0 : variable_data_impl_static_;
554 0 : std::unique_ptr<value_type[]> variable_data_impl_dynamic_{};
555 0 : bool owning_{true};
556 0 : size_t size_ = 0;
557 0 : size_t number_of_grid_points_ = 0;
558 :
559 0 : pointer_type variable_data_;
560 0 : tuples::TaggedTuple<Tags...> reference_variable_data_;
561 : };
562 :
563 : // The above Variables implementation doesn't work for an empty parameter pack,
564 : // so specialize here.
565 : template <>
566 0 : class Variables<tmpl::list<>> {
567 : public:
568 0 : using tags_list = tmpl::list<>;
569 0 : static constexpr size_t number_of_independent_components = 0;
570 0 : Variables() = default;
571 0 : explicit Variables(const size_t /*number_of_grid_points*/) {}
572 : template <typename T>
573 0 : Variables(const T* /*pointer*/, const size_t /*size*/) {}
574 0 : static constexpr size_t size() { return 0; }
575 : };
576 :
577 : // gcc8 screams when the empty Variables has pup as a member function, so we
578 : // declare pup as a free function here.
579 : // clang-tidy: runtime-references
580 0 : SPECTRE_ALWAYS_INLINE void pup(
581 : PUP::er& /*p*/, // NOLINT
582 : Variables<tmpl::list<>>& /* unused */) { // NOLINT
583 : }
584 0 : SPECTRE_ALWAYS_INLINE void operator|(
585 : PUP::er& /*p*/, Variables<tmpl::list<>>& /* unused */) { // NOLINT
586 : }
587 :
588 0 : SPECTRE_ALWAYS_INLINE bool operator==(const Variables<tmpl::list<>>& /*lhs*/,
589 : const Variables<tmpl::list<>>& /*rhs*/) {
590 : return true;
591 : }
592 0 : SPECTRE_ALWAYS_INLINE bool operator!=(const Variables<tmpl::list<>>& /*lhs*/,
593 : const Variables<tmpl::list<>>& /*rhs*/) {
594 : return false;
595 : }
596 :
597 0 : inline std::ostream& operator<<(std::ostream& os,
598 : const Variables<tmpl::list<>>& /*d*/) {
599 : return os << "{}";
600 : }
601 :
602 : template <typename... Tags>
603 : Variables<tmpl::list<Tags...>>::Variables() {
604 : // This makes an assertion trigger if one tries to assign to
605 : // components of a default-constructed Variables.
606 : const auto set_refs = [](auto& var) {
607 : for (auto& dv : var) {
608 : dv.set_data_ref(nullptr, 0);
609 : }
610 : return 0;
611 : };
612 : (void)set_refs;
613 : expand_pack(set_refs(tuples::get<Tags>(reference_variable_data_))...);
614 : }
615 :
616 : template <typename... Tags>
617 : Variables<tmpl::list<Tags...>>::Variables(const size_t number_of_grid_points) {
618 : initialize(number_of_grid_points);
619 : }
620 :
621 : template <typename... Tags>
622 : Variables<tmpl::list<Tags...>>::Variables(const size_t number_of_grid_points,
623 : const value_type value) {
624 : initialize(number_of_grid_points, value);
625 : }
626 :
627 : template <typename... Tags>
628 : void Variables<tmpl::list<Tags...>>::initialize(
629 : const size_t number_of_grid_points) {
630 : if (number_of_grid_points_ == 0) {
631 : variable_data_impl_dynamic_.reset();
632 : size_ = 0;
633 : number_of_grid_points_ = 0;
634 : }
635 : if (number_of_grid_points_ == number_of_grid_points) {
636 : return;
637 : }
638 : if (UNLIKELY(not is_owning())) {
639 : ERROR("Variables::initialize cannot be called on a non-owning Variables. "
640 : "This likely happened because of an attempted resize. The current "
641 : "number of grid points is " << number_of_grid_points_ << " and the "
642 : "requested number is " << number_of_grid_points << ".");
643 : }
644 : number_of_grid_points_ = number_of_grid_points;
645 : size_ = number_of_grid_points * number_of_independent_components;
646 : if (size_ > 0) {
647 : if (number_of_grid_points_ == 1) {
648 : variable_data_impl_dynamic_.reset();
649 : } else {
650 : variable_data_impl_dynamic_ =
651 : cpp20::make_unique_for_overwrite<value_type[]>(size_);
652 : }
653 : add_reference_variable_data();
654 : #if defined(SPECTRE_DEBUG) || defined(SPECTRE_NAN_INIT)
655 : std::fill(variable_data_.data(), variable_data_.data() + size_,
656 : make_signaling_NaN<value_type>());
657 : #endif // SPECTRE_DEBUG
658 : }
659 : }
660 :
661 : template <typename... Tags>
662 : void Variables<tmpl::list<Tags...>>::initialize(
663 : const size_t number_of_grid_points, const value_type value) {
664 : initialize(number_of_grid_points);
665 : std::fill(variable_data_.data(), variable_data_.data() + size_, value);
666 : }
667 :
668 : /// \cond HIDDEN_SYMBOLS
669 : template <typename... Tags>
670 : Variables<tmpl::list<Tags...>>::Variables(
671 : const Variables<tmpl::list<Tags...>>& rhs) {
672 : initialize(rhs.number_of_grid_points());
673 : variable_data_ =
674 : static_cast<const blaze::Vector<pointer_type, transpose_flag>&>(
675 : rhs.variable_data_);
676 : }
677 :
678 : template <typename... Tags>
679 : Variables<tmpl::list<Tags...>>& Variables<tmpl::list<Tags...>>::operator=(
680 : const Variables<tmpl::list<Tags...>>& rhs) {
681 : if (&rhs == this) {
682 : return *this;
683 : }
684 : initialize(rhs.number_of_grid_points());
685 : variable_data_ =
686 : static_cast<const blaze::Vector<pointer_type, transpose_flag>&>(
687 : rhs.variable_data_);
688 : return *this;
689 : }
690 :
691 : template <typename... Tags>
692 : Variables<tmpl::list<Tags...>>::Variables(Variables<tmpl::list<Tags...>>&& rhs)
693 : : variable_data_impl_dynamic_(std::move(rhs.variable_data_impl_dynamic_)),
694 : owning_(rhs.owning_),
695 : size_(rhs.size()),
696 : number_of_grid_points_(rhs.number_of_grid_points()),
697 : variable_data_(std::move(rhs.variable_data_)) {
698 : if (number_of_grid_points_ == 1) {
699 : #if defined(__GNUC__) and not defined(__clang__)
700 : #pragma GCC diagnostic push
701 : #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
702 : #endif // defined(__GNUC__) and not defined(__clang__)
703 : variable_data_impl_static_ = std::move(rhs.variable_data_impl_static_);
704 : #if defined(__GNUC__) and not defined(__clang__)
705 : #pragma GCC diagnostic pop
706 : #endif // defined(__GNUC__) and not defined(__clang__)
707 : }
708 : rhs.variable_data_impl_dynamic_.reset();
709 : rhs.owning_ = true;
710 : rhs.size_ = 0;
711 : rhs.number_of_grid_points_ = 0;
712 : add_reference_variable_data();
713 : }
714 :
715 : template <typename... Tags>
716 : Variables<tmpl::list<Tags...>>& Variables<tmpl::list<Tags...>>::operator=(
717 : Variables<tmpl::list<Tags...>>&& rhs) {
718 : if (this == &rhs) {
719 : return *this;
720 : }
721 : owning_ = rhs.owning_;
722 : size_ = rhs.size_;
723 : number_of_grid_points_ = std::move(rhs.number_of_grid_points_);
724 : variable_data_ = std::move(rhs.variable_data_);
725 : variable_data_impl_dynamic_ = std::move(rhs.variable_data_impl_dynamic_);
726 : if (number_of_grid_points_ == 1) {
727 : #if defined(__GNUC__) and not defined(__clang__)
728 : #pragma GCC diagnostic push
729 : #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
730 : #endif // defined(__GNUC__) and not defined(__clang__)
731 : variable_data_impl_static_ = std::move(rhs.variable_data_impl_static_);
732 : #if defined(__GNUC__) and not defined(__clang__)
733 : #pragma GCC diagnostic pop
734 : #endif // defined(__GNUC__) and not defined(__clang__)
735 : }
736 :
737 : rhs.variable_data_impl_dynamic_.reset();
738 : rhs.owning_ = true;
739 : rhs.size_ = 0;
740 : rhs.number_of_grid_points_ = 0;
741 : add_reference_variable_data();
742 : return *this;
743 : }
744 :
745 : template <typename... Tags>
746 : template <typename... WrappedTags,
747 : Requires<tmpl2::flat_all<
748 : std::is_same<db::remove_all_prefixes<WrappedTags>,
749 : db::remove_all_prefixes<Tags>>::value...>::value>>
750 : Variables<tmpl::list<Tags...>>::Variables(
751 : const Variables<tmpl::list<WrappedTags...>>& rhs) {
752 : static_assert(
753 : (std::is_same_v<typename Tags::type, typename WrappedTags::type> and ...),
754 : "Tensor types do not match!");
755 : initialize(rhs.number_of_grid_points());
756 : variable_data_ =
757 : static_cast<const blaze::Vector<pointer_type, transpose_flag>&>(
758 : rhs.variable_data_);
759 : }
760 :
761 : template <typename... Tags>
762 : template <typename... WrappedTags,
763 : Requires<tmpl2::flat_all<
764 : std::is_same<db::remove_all_prefixes<WrappedTags>,
765 : db::remove_all_prefixes<Tags>>::value...>::value>>
766 : Variables<tmpl::list<Tags...>>& Variables<tmpl::list<Tags...>>::operator=(
767 : const Variables<tmpl::list<WrappedTags...>>& rhs) {
768 : static_assert(
769 : (std::is_same_v<typename Tags::type, typename WrappedTags::type> and ...),
770 : "Tensor types do not match!");
771 : initialize(rhs.number_of_grid_points());
772 : variable_data_ =
773 : static_cast<const blaze::Vector<pointer_type, transpose_flag>&>(
774 : rhs.variable_data_);
775 : return *this;
776 : }
777 :
778 : template <typename... Tags>
779 : template <typename... WrappedTags,
780 : Requires<tmpl2::flat_all<
781 : std::is_same<db::remove_all_prefixes<WrappedTags>,
782 : db::remove_all_prefixes<Tags>>::value...>::value>>
783 : Variables<tmpl::list<Tags...>>::Variables(
784 : Variables<tmpl::list<WrappedTags...>>&& rhs)
785 : : variable_data_impl_dynamic_(std::move(rhs.variable_data_impl_dynamic_)),
786 : owning_(rhs.owning_),
787 : size_(rhs.size()),
788 : number_of_grid_points_(rhs.number_of_grid_points()),
789 : variable_data_(std::move(rhs.variable_data_)) {
790 : static_assert(
791 : (std::is_same_v<typename Tags::type, typename WrappedTags::type> and ...),
792 : "Tensor types do not match!");
793 : if (number_of_grid_points_ == 1) {
794 : #if defined(__GNUC__) and not defined(__clang__)
795 : #pragma GCC diagnostic push
796 : #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
797 : #endif // defined(__GNUC__) and not defined(__clang__)
798 : variable_data_impl_static_ = std::move(rhs.variable_data_impl_static_);
799 : #if defined(__GNUC__) and not defined(__clang__)
800 : #pragma GCC diagnostic pop
801 : #endif // defined(__GNUC__) and not defined(__clang__)
802 : }
803 : rhs.variable_data_impl_dynamic_.reset();
804 : rhs.size_ = 0;
805 : rhs.owning_ = true;
806 : rhs.number_of_grid_points_ = 0;
807 : add_reference_variable_data();
808 : }
809 :
810 : template <typename... Tags>
811 : template <typename... WrappedTags,
812 : Requires<tmpl2::flat_all<
813 : std::is_same<db::remove_all_prefixes<WrappedTags>,
814 : db::remove_all_prefixes<Tags>>::value...>::value>>
815 : Variables<tmpl::list<Tags...>>& Variables<tmpl::list<Tags...>>::operator=(
816 : Variables<tmpl::list<WrappedTags...>>&& rhs) {
817 : static_assert(
818 : (std::is_same_v<typename Tags::type, typename WrappedTags::type> and ...),
819 : "Tensor types do not match!");
820 : variable_data_ = std::move(rhs.variable_data_);
821 : owning_ = rhs.owning_;
822 : size_ = rhs.size_;
823 : number_of_grid_points_ = std::move(rhs.number_of_grid_points_);
824 : variable_data_impl_dynamic_ = std::move(rhs.variable_data_impl_dynamic_);
825 : if (number_of_grid_points_ == 1) {
826 : #if defined(__GNUC__) and not defined(__clang__)
827 : #pragma GCC diagnostic push
828 : #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
829 : #endif // defined(__GNUC__) and not defined(__clang__)
830 : variable_data_impl_static_ = std::move(rhs.variable_data_impl_static_);
831 : #if defined(__GNUC__) and not defined(__clang__)
832 : #pragma GCC diagnostic pop
833 : #endif // defined(__GNUC__) and not defined(__clang__)
834 : }
835 :
836 : rhs.variable_data_impl_dynamic_.reset();
837 : rhs.size_ = 0;
838 : rhs.owning_ = true;
839 : rhs.number_of_grid_points_ = 0;
840 : add_reference_variable_data();
841 : return *this;
842 : }
843 :
844 : template <typename... Tags>
845 : Variables<tmpl::list<Tags...>>::Variables(const pointer start,
846 : const size_t size) {
847 : set_data_ref(start, size);
848 : }
849 :
850 : template <typename... Tags>
851 : void Variables<tmpl::list<Tags...>>::pup(PUP::er& p) {
852 : ASSERT(
853 : owning_,
854 : "Cannot pup a non-owning Variables! It may be reasonable to pack a "
855 : "non-owning Variables, but not to unpack one. This should be discussed "
856 : "in an issue with the core devs if the feature seems necessary.");
857 : size_t number_of_grid_points = number_of_grid_points_;
858 : p | number_of_grid_points;
859 : if (p.isUnpacking()) {
860 : initialize(number_of_grid_points);
861 : }
862 : PUParray(p, variable_data_.data(), size_);
863 : }
864 : /// \endcond
865 :
866 : /// @{
867 : /*!
868 : * \ingroup DataStructuresGroup
869 : * \brief Return Tag::type pointing into the contiguous array
870 : *
871 : * \tparam Tag the variable to return
872 : */
873 : template <typename Tag, typename TagList>
874 1 : constexpr typename Tag::type& get(Variables<TagList>& v) {
875 : static_assert(tmpl::list_contains_v<TagList, Tag>,
876 : "Could not retrieve Tag from Variables. See the first "
877 : "template parameter of the instantiation for what Tag is "
878 : "being retrieved and the second template parameter for "
879 : "what Tags are available.");
880 : return tuples::get<Tag>(v.reference_variable_data_);
881 : }
882 : template <typename Tag, typename TagList>
883 1 : constexpr const typename Tag::type& get(const Variables<TagList>& v) {
884 : static_assert(tmpl::list_contains_v<TagList, Tag>,
885 : "Could not retrieve Tag from Variables. See the first "
886 : "template parameter of the instantiation for what Tag is "
887 : "being retrieved and the second template parameter for "
888 : "what Tags are available.");
889 : return tuples::get<Tag>(v.reference_variable_data_);
890 : }
891 : /// @}
892 :
893 : template <typename... Tags>
894 : template <typename VT, bool VF>
895 : Variables<tmpl::list<Tags...>>::Variables(
896 : const blaze::DenseVector<VT, VF>& expression) {
897 : ASSERT((*expression).size() % number_of_independent_components == 0,
898 : "Invalid size " << (*expression).size() << " for a Variables with "
899 : << number_of_independent_components << " components.");
900 : initialize((*expression).size() / number_of_independent_components);
901 : variable_data_ = expression;
902 : }
903 :
904 : /// \cond
905 : template <typename... Tags>
906 : template <typename VT, bool VF>
907 : Variables<tmpl::list<Tags...>>& Variables<tmpl::list<Tags...>>::operator=(
908 : const blaze::DenseVector<VT, VF>& expression) {
909 : ASSERT((*expression).size() % number_of_independent_components == 0,
910 : "Invalid size " << (*expression).size() << " for a Variables with "
911 : << number_of_independent_components << " components.");
912 : initialize((*expression).size() / number_of_independent_components);
913 : variable_data_ = expression;
914 : return *this;
915 : }
916 : /// \endcond
917 :
918 : /// \cond HIDDEN_SYMBOLS
919 : template <typename... Tags>
920 : void Variables<tmpl::list<Tags...>>::add_reference_variable_data() {
921 : if (size_ == 0) {
922 : return;
923 : }
924 : if (is_owning()) {
925 : if (number_of_grid_points_ == 1) {
926 : variable_data_.reset(variable_data_impl_static_.data(), size_);
927 : } else {
928 : variable_data_.reset(variable_data_impl_dynamic_.get(), size_);
929 : }
930 : }
931 : ASSERT(variable_data_.size() == size_ and
932 : size_ == number_of_grid_points_ * number_of_independent_components,
933 : "Size mismatch: variable_data_.size() = "
934 : << variable_data_.size() << " size_ = " << size_ << " should be: "
935 : << number_of_grid_points_ * number_of_independent_components
936 : << "\nThis is an internal inconsistency bug in Variables. Please "
937 : "file an issue.");
938 : size_t variable_offset = 0;
939 : tmpl::for_each<tags_list>([this, &variable_offset](auto tag_v) {
940 : using Tag = tmpl::type_from<decltype(tag_v)>;
941 : auto& var = tuples::get<Tag>(reference_variable_data_);
942 : for (size_t i = 0; i < Tag::type::size(); ++i) {
943 : var[i].set_data_ref(
944 : &variable_data_[variable_offset++ * number_of_grid_points_],
945 : number_of_grid_points_);
946 : }
947 : });
948 : }
949 : /// \endcond
950 :
951 : template <typename... Tags>
952 0 : Variables<tmpl::list<Tags...>>& operator*=(
953 : Variables<tmpl::list<Tags...>>& lhs,
954 : const typename Variables<tmpl::list<Tags...>>::vector_type& rhs) {
955 : using value_type = typename Variables<tmpl::list<Tags...>>::value_type;
956 : ASSERT(lhs.number_of_grid_points() == rhs.size(),
957 : "Size mismatch in multiplication: " << lhs.number_of_grid_points()
958 : << " and " << rhs.size());
959 : value_type* const lhs_data = lhs.data();
960 : const value_type* const rhs_data = rhs.data();
961 : for (size_t c = 0; c < lhs.number_of_independent_components; ++c) {
962 : for (size_t s = 0; s < lhs.number_of_grid_points(); ++s) {
963 : // clang-tidy: do not use pointer arithmetic
964 : lhs_data[c * lhs.number_of_grid_points() + s] *= rhs_data[s]; // NOLINT
965 : }
966 : }
967 : return lhs;
968 : }
969 :
970 : template <typename... Tags>
971 0 : Variables<tmpl::list<Tags...>> operator*(
972 : const Variables<tmpl::list<Tags...>>& lhs,
973 : const typename Variables<tmpl::list<Tags...>>::vector_type& rhs) {
974 : auto result = lhs;
975 : result *= rhs;
976 : return result;
977 : }
978 :
979 : template <typename... Tags>
980 0 : Variables<tmpl::list<Tags...>> operator*(
981 : const typename Variables<tmpl::list<Tags...>>::vector_type& lhs,
982 : const Variables<tmpl::list<Tags...>>& rhs) {
983 : auto result = rhs;
984 : result *= lhs;
985 : return result;
986 : }
987 :
988 : template <typename... Tags>
989 0 : Variables<tmpl::list<Tags...>>& operator/=(
990 : Variables<tmpl::list<Tags...>>& lhs,
991 : const typename Variables<tmpl::list<Tags...>>::vector_type& rhs) {
992 : ASSERT(lhs.number_of_grid_points() == rhs.size(),
993 : "Size mismatch in multiplication: " << lhs.number_of_grid_points()
994 : << " and " << rhs.size());
995 : using value_type = typename Variables<tmpl::list<Tags...>>::value_type;
996 : value_type* const lhs_data = lhs.data();
997 : const value_type* const rhs_data = rhs.data();
998 : for (size_t c = 0; c < lhs.number_of_independent_components; ++c) {
999 : for (size_t s = 0; s < lhs.number_of_grid_points(); ++s) {
1000 : // clang-tidy: do not use pointer arithmetic
1001 : lhs_data[c * lhs.number_of_grid_points() + s] /= rhs_data[s]; // NOLINT
1002 : }
1003 : }
1004 : return lhs;
1005 : }
1006 :
1007 : template <typename... Tags>
1008 0 : Variables<tmpl::list<Tags...>> operator/(
1009 : const Variables<tmpl::list<Tags...>>& lhs,
1010 : const typename Variables<tmpl::list<Tags...>>::vector_type& rhs) {
1011 : auto result = lhs;
1012 : result /= rhs;
1013 : return result;
1014 : }
1015 :
1016 : template <typename... Tags, Requires<sizeof...(Tags) != 0> = nullptr>
1017 0 : std::ostream& operator<<(std::ostream& os,
1018 : const Variables<tmpl::list<Tags...>>& d) {
1019 : size_t count = 0;
1020 : const auto helper = [&os, &d, &count](auto tag_v) {
1021 : count++;
1022 : using Tag = typename decltype(tag_v)::type;
1023 : os << db::tag_name<Tag>() << ":\n";
1024 : os << get<Tag>(d);
1025 : if (count < sizeof...(Tags)) {
1026 : os << "\n\n";
1027 : }
1028 : };
1029 : EXPAND_PACK_LEFT_TO_RIGHT(helper(tmpl::type_<Tags>{}));
1030 : return os;
1031 : }
1032 :
1033 : template <typename TagsList>
1034 0 : bool operator!=(const Variables<TagsList>& lhs,
1035 : const Variables<TagsList>& rhs) {
1036 : return not(lhs == rhs);
1037 : }
1038 :
1039 : template <typename... TagsLhs, typename... TagsRhs,
1040 : Requires<not std::is_same<tmpl::list<TagsLhs...>,
1041 : tmpl::list<TagsRhs...>>::value> = nullptr>
1042 0 : void swap(Variables<tmpl::list<TagsLhs...>>& lhs,
1043 : Variables<tmpl::list<TagsRhs...>>& rhs) {
1044 : Variables<tmpl::list<TagsLhs...>> temp{std::move(lhs)};
1045 : lhs = std::move(rhs);
1046 : rhs = std::move(temp);
1047 : }
1048 :
1049 : /// \ingroup DataStructuresGroup
1050 : /// Construct a variables from the `Tensor`s in a `TaggedTuple`.
1051 : template <typename... Tags>
1052 1 : Variables<tmpl::list<Tags...>> variables_from_tagged_tuple(
1053 : const tuples::TaggedTuple<Tags...>& tuple) {
1054 : auto result = make_with_value<Variables<tmpl::list<Tags...>>>(
1055 : get<tmpl::front<tmpl::list<Tags...>>>(tuple), 0.0);
1056 : result.assign_subset(tuple);
1057 : return result;
1058 : }
1059 :
1060 : namespace MakeWithValueImpls {
1061 : template <typename TagList>
1062 0 : struct MakeWithSize<Variables<TagList>> {
1063 0 : static SPECTRE_ALWAYS_INLINE Variables<TagList> apply(
1064 : const size_t size, const typename Variables<TagList>::value_type value) {
1065 : return Variables<TagList>(size, value);
1066 : }
1067 : };
1068 :
1069 : template <typename TagList>
1070 0 : struct NumberOfPoints<Variables<TagList>> {
1071 0 : static SPECTRE_ALWAYS_INLINE size_t apply(const Variables<TagList>& input) {
1072 : return input.number_of_grid_points();
1073 : }
1074 : };
1075 : } // namespace MakeWithValueImpls
1076 :
1077 1 : namespace EqualWithinRoundoffImpls {
1078 : // It would be nice to use `blaze::equal` in these implementations, but it
1079 : // doesn't currently allow to specify a tolerance. See upstream issue:
1080 : // https://bitbucket.org/blaze-lib/blaze/issues/417/adjust-relaxed-equal-accuracy
1081 : template <typename TagList, typename Floating>
1082 : struct EqualWithinRoundoffImpl<Variables<TagList>, Floating,
1083 : Requires<std::is_floating_point_v<Floating>>> {
1084 : static bool apply(const Variables<TagList>& lhs, const Floating& rhs,
1085 : const double eps, const double scale) {
1086 : for (size_t i = 0; i < lhs.size(); ++i) {
1087 : if (not equal_within_roundoff(lhs.data()[i], rhs, eps, scale)) {
1088 : return false;
1089 : }
1090 : }
1091 : return true;
1092 : }
1093 : };
1094 : template <typename TagList, typename Floating>
1095 : struct EqualWithinRoundoffImpl<Floating, Variables<TagList>,
1096 : Requires<std::is_floating_point_v<Floating>>> {
1097 : static SPECTRE_ALWAYS_INLINE bool apply(const Floating& lhs,
1098 : const Variables<TagList>& rhs,
1099 : const double eps,
1100 : const double scale) {
1101 : return equal_within_roundoff(rhs, lhs, eps, scale);
1102 : }
1103 : };
1104 : template <typename LhsTagList, typename RhsTagList>
1105 0 : struct EqualWithinRoundoffImpl<Variables<LhsTagList>, Variables<RhsTagList>> {
1106 0 : static bool apply(const Variables<LhsTagList>& lhs,
1107 : const Variables<RhsTagList>& rhs, const double eps,
1108 : const double scale) {
1109 : ASSERT(lhs.size() == rhs.size(),
1110 : "Can only compare two Variables of the same size, but lhs has size "
1111 : << lhs.size() << " and rhs has size " << rhs.size() << ".");
1112 : for (size_t i = 0; i < lhs.size(); ++i) {
1113 : if (not equal_within_roundoff(lhs.data()[i], rhs.data()[i], eps, scale)) {
1114 : return false;
1115 : }
1116 : }
1117 : return true;
1118 : }
1119 : };
1120 : } // namespace EqualWithinRoundoffImpls
1121 :
1122 : namespace db {
1123 : // Enable subitems for ::Tags::Variables and derived tags (e.g. compute tags).
1124 : // Other tags that hold a `Variables` don't expose the constituent tensors as
1125 : // subitems by default.
1126 : namespace Variables_detail {
1127 : // Check if the argument is a `::Tags::Variables`, or derived from it. Can't use
1128 : // `tt:is_a_v` because we also want to match derived classes. Can't use
1129 : // `std::is_base_of` because `::Tags::Variables` is a template.
1130 : template <typename TagsList>
1131 : constexpr std::true_type is_a_variables_tag(::Tags::Variables<TagsList>&&) {
1132 : return {};
1133 : }
1134 : constexpr std::false_type is_a_variables_tag(...) { return {}; }
1135 : template <typename Tag>
1136 : static constexpr bool is_a_variables_tag_v =
1137 : decltype(is_a_variables_tag(std::declval<Tag>()))::value;
1138 : } // namespace Variables_detail
1139 : template <typename Tag>
1140 : struct Subitems<Tag, Requires<Variables_detail::is_a_variables_tag_v<Tag>>> {
1141 : using type = typename Tag::type::tags_list;
1142 :
1143 : template <typename Subtag, typename LocalTag = Tag>
1144 : static void create_item(
1145 : const gsl::not_null<typename LocalTag::type*> parent_value,
1146 : const gsl::not_null<typename Subtag::type*> sub_value) {
1147 : auto& vars = get<Subtag>(*parent_value);
1148 : // Only update the Tensor if the Variables has changed its allocation
1149 : if constexpr (not is_any_spin_weighted_v<typename Subtag::type::type>) {
1150 : if (vars.begin()->data() != sub_value->begin()->data()) {
1151 : for (auto vars_it = vars.begin(), sub_var_it = sub_value->begin();
1152 : vars_it != vars.end(); ++vars_it, ++sub_var_it) {
1153 : sub_var_it->set_data_ref(make_not_null(&*vars_it));
1154 : }
1155 : }
1156 : } else {
1157 : if (vars.begin()->data().data() != sub_value->begin()->data().data()) {
1158 : for (auto vars_it = vars.begin(), sub_var_it = sub_value->begin();
1159 : vars_it != vars.end(); ++vars_it, ++sub_var_it) {
1160 : sub_var_it->set_data_ref(make_not_null(&*vars_it));
1161 : }
1162 : }
1163 : }
1164 : }
1165 :
1166 : template <typename Subtag>
1167 : static const typename Subtag::type& create_compute_item(
1168 : const typename Tag::type& parent_value) {
1169 : return get<Subtag>(parent_value);
1170 : }
1171 : };
1172 : } // namespace db
1173 :
1174 : namespace Variables_detail {
1175 : template <typename Tags>
1176 : using MathWrapperVectorType = tmpl::conditional_t<
1177 : std::is_same_v<typename Variables<Tags>::value_type, double>, DataVector,
1178 : ComplexDataVector>;
1179 : } // namespace Variables_detail
1180 :
1181 : template <typename Tags>
1182 0 : auto make_math_wrapper(const gsl::not_null<Variables<Tags>*> data) {
1183 : using Vector = Variables_detail::MathWrapperVectorType<Tags>;
1184 : Vector referencing(data->data(), data->size());
1185 : return make_math_wrapper(&referencing);
1186 : }
1187 :
1188 : template <typename Tags>
1189 0 : auto make_math_wrapper(const Variables<Tags>& data) {
1190 : using Vector = Variables_detail::MathWrapperVectorType<Tags>;
1191 : const Vector referencing(
1192 : const_cast<typename Vector::value_type*>(data.data()), data.size());
1193 : return make_math_wrapper(referencing);
1194 : }
1195 :
1196 : /// \cond
1197 : template <size_t I, class... Tags>
1198 : inline constexpr typename tmpl::at_c<tmpl::list<Tags...>, I>::type&& get(
1199 : Variables<tmpl::list<Tags...>>&& t) {
1200 : return get<tmpl::at_c<tmpl::list<Tags...>, I>>(t);
1201 : }
1202 :
1203 : template <size_t I, class... Tags>
1204 : inline constexpr const typename tmpl::at_c<tmpl::list<Tags...>, I>::type& get(
1205 : const Variables<tmpl::list<Tags...>>& t) {
1206 : return get<tmpl::at_c<tmpl::list<Tags...>, I>>(t);
1207 : }
1208 :
1209 : template <size_t I, class... Tags>
1210 : inline constexpr typename tmpl::at_c<tmpl::list<Tags...>, I>::type& get(
1211 : Variables<tmpl::list<Tags...>>& t) {
1212 : return get<tmpl::at_c<tmpl::list<Tags...>, I>>(t);
1213 : }
1214 : /// \endcond
1215 :
1216 : namespace std {
1217 : template <typename... Tags>
1218 : struct tuple_size<Variables<tmpl::list<Tags...>>>
1219 : : std::integral_constant<int, sizeof...(Tags)> {};
1220 : template <size_t I, typename... Tags>
1221 : struct tuple_element<I, Variables<tmpl::list<Tags...>>> {
1222 : using type = typename tmpl::at_c<tmpl::list<Tags...>, I>::type;
1223 : };
1224 : } // namespace std
1225 :
1226 : template <typename TagsList>
1227 0 : bool contains_allocations(const Variables<TagsList>& value) {
1228 : return value.number_of_grid_points() > 1;
1229 : }
1230 :
1231 0 : inline bool contains_allocations(const Variables<tmpl::list<>>& /*value*/) {
1232 : return false;
1233 : }
|
