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 : #include <type_traits>
25 :
26 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
27 : #include "DataStructures/DataBox/Subitems.hpp"
28 : #include "DataStructures/DataBox/Tag.hpp"
29 : #include "DataStructures/DataBox/TagName.hpp"
30 : #include "DataStructures/DataBox/TagTraits.hpp"
31 : #include "DataStructures/DataVector.hpp"
32 : #include "DataStructures/MathWrapper.hpp"
33 : #include "DataStructures/SpinWeighted.hpp"
34 : #include "DataStructures/Tensor/IndexType.hpp"
35 : #include "DataStructures/Tensor/Tensor.hpp"
36 : #include "Utilities/EqualWithinRoundoff.hpp"
37 : #include "Utilities/ErrorHandling/Assert.hpp"
38 : #include "Utilities/ForceInline.hpp"
39 : #include "Utilities/Gsl.hpp"
40 : #include "Utilities/Literals.hpp"
41 : #include "Utilities/MakeSignalingNan.hpp"
42 : #include "Utilities/MemoryHelpers.hpp"
43 : #include "Utilities/PrettyType.hpp"
44 : #include "Utilities/Requires.hpp"
45 : #include "Utilities/SetNumberOfGridPoints.hpp"
46 : #include "Utilities/TMPL.hpp"
47 : #include "Utilities/TaggedTuple.hpp"
48 : #include "Utilities/TypeTraits.hpp"
49 : #include "Utilities/TypeTraits/IsA.hpp"
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 <
444 : typename... WrappedTags,
445 : Requires<
446 : sizeof...(WrappedTags) == sizeof...(Tags) and
447 : tmpl2::flat_all_v<std::is_same_v<db::remove_all_prefixes<WrappedTags>,
448 : db::remove_all_prefixes<Tags>>...>
449 : #ifdef __CUDACC__
450 : and tmpl2::flat_all_v<std::is_same_v<typename WrappedTags::type,
451 : typename Tags::type>...>
452 : #endif
453 : > = nullptr>
454 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator+(
455 : const Variables<tmpl::list<WrappedTags...>>& lhs, const Variables& rhs) {
456 : #ifndef __CUDACC__
457 : // nvcc incorrectly hits this static_assert and variants of it, but using
458 : // the additional requires part above is fine.
459 : static_assert(
460 : (std::is_same_v<typename Tags::type, typename WrappedTags::type> and
461 : ...),
462 : "Tensor types do not match!");
463 : #endif
464 : return lhs.get_variable_data() + rhs.variable_data_;
465 : }
466 : template <typename VT, bool VF>
467 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator+(
468 : const blaze::DenseVector<VT, VF>& lhs, const Variables& rhs) {
469 : return *lhs + rhs.variable_data_;
470 : }
471 : template <typename VT, bool VF>
472 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator+(
473 : const Variables& lhs, const blaze::DenseVector<VT, VF>& rhs) {
474 : return lhs.variable_data_ + *rhs;
475 : }
476 :
477 : template <
478 : typename... WrappedTags,
479 : Requires<
480 : sizeof...(WrappedTags) == sizeof...(Tags) and
481 : tmpl2::flat_all_v<std::is_same_v<db::remove_all_prefixes<WrappedTags>,
482 : db::remove_all_prefixes<Tags>>...>
483 : #ifdef __CUDACC__
484 : and tmpl2::flat_all_v<std::is_same_v<typename WrappedTags::type,
485 : typename Tags::type>...>
486 : #endif
487 : > = nullptr>
488 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator-(
489 : const Variables<tmpl::list<WrappedTags...>>& lhs, const Variables& rhs) {
490 : #ifndef __CUDACC__
491 : // nvcc incorrectly hits this static_assert and variants of it, but using
492 : // the additional requires part above is fine.
493 : static_assert(
494 : (std::is_same_v<typename Tags::type, typename WrappedTags::type> and
495 : ...),
496 : "Tensor types do not match!");
497 : #endif
498 : return lhs.get_variable_data() - rhs.variable_data_;
499 : }
500 : template <typename VT, bool VF>
501 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator-(
502 : const blaze::DenseVector<VT, VF>& lhs, const Variables& rhs) {
503 : return *lhs - rhs.variable_data_;
504 : }
505 : template <typename VT, bool VF>
506 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator-(
507 : const Variables& lhs, const blaze::DenseVector<VT, VF>& rhs) {
508 : return lhs.variable_data_ - *rhs;
509 : }
510 :
511 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator*(const Variables& lhs,
512 : const value_type& rhs) {
513 : return lhs.variable_data_ * rhs;
514 : }
515 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator*(const value_type& lhs,
516 : const Variables& rhs) {
517 : return lhs * rhs.variable_data_;
518 : }
519 :
520 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator/(const Variables& lhs,
521 : const value_type& rhs) {
522 : return lhs.variable_data_ / rhs;
523 : }
524 :
525 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator-(const Variables& lhs) {
526 : return -lhs.variable_data_;
527 : }
528 0 : friend SPECTRE_ALWAYS_INLINE decltype(auto) operator+(const Variables& lhs) {
529 : return lhs.variable_data_;
530 : }
531 :
532 : private:
533 : /// @{
534 : /*!
535 : * \brief Subscript operator
536 : *
537 : * The subscript operator is private since it should not be used directly.
538 : * Mathematical operations should be done using the math operators provided.
539 : * Since the internal ordering of variables is implementation defined there
540 : * is no safe way to perform any operation that is not a linear combination of
541 : * Variables. Retrieving a Tensor must be done via the `get()` function.
542 : *
543 : * \requires `i >= 0 and i < size()`
544 : */
545 1 : SPECTRE_ALWAYS_INLINE value_type& operator[](const size_type i) {
546 : return variable_data_[i];
547 : }
548 1 : SPECTRE_ALWAYS_INLINE const value_type& operator[](const size_type i) const {
549 : return variable_data_[i];
550 : }
551 : /// @}
552 :
553 0 : void add_reference_variable_data();
554 :
555 0 : friend bool operator==(const Variables& lhs, const Variables& rhs) {
556 : return blaze::equal<blaze::strict>(lhs.variable_data_, rhs.variable_data_);
557 : }
558 :
559 : template <typename VT, bool TF>
560 0 : friend bool operator==(const Variables& lhs,
561 : const blaze::DenseVector<VT, TF>& rhs) {
562 : return blaze::equal<blaze::strict>(lhs.variable_data_, *rhs);
563 : }
564 :
565 : template <typename VT, bool TF>
566 0 : friend bool operator==(const blaze::DenseVector<VT, TF>& lhs,
567 : const Variables& rhs) {
568 : return blaze::equal<blaze::strict>(*lhs, rhs.variable_data_);
569 : }
570 :
571 : template <class FriendTags>
572 0 : friend class Variables;
573 :
574 : std::array<value_type, number_of_independent_components>
575 0 : variable_data_impl_static_;
576 0 : std::unique_ptr<value_type[]> variable_data_impl_dynamic_{};
577 0 : bool owning_{true};
578 0 : size_t size_ = 0;
579 0 : size_t number_of_grid_points_ = 0;
580 :
581 0 : pointer_type variable_data_;
582 0 : tuples::TaggedTuple<Tags...> reference_variable_data_;
583 : };
584 :
585 : // The above Variables implementation doesn't work for an empty parameter pack,
586 : // so specialize here.
587 : template <>
588 0 : class Variables<tmpl::list<>> {
589 : public:
590 0 : using tags_list = tmpl::list<>;
591 0 : static constexpr size_t number_of_independent_components = 0;
592 0 : Variables() = default;
593 0 : explicit Variables(const size_t /*number_of_grid_points*/) {}
594 : template <typename T>
595 0 : Variables(const T* /*pointer*/, const size_t /*size*/) {}
596 0 : static constexpr size_t size() { return 0; }
597 0 : static constexpr size_t number_of_grid_points() { return 0; }
598 0 : void assign_subset(const Variables<tmpl::list<>>& /*unused*/) {}
599 0 : void assign_subset(const tuples::TaggedTuple<>& /*unused*/) {}
600 : // Initialization for empty variables should ignore any input.
601 0 : void initialize(size_t /*number_of_grid_points*/) {}
602 : };
603 :
604 : // gcc8 screams when the empty Variables has pup as a member function, so we
605 : // declare pup as a free function here.
606 : // clang-tidy: runtime-references
607 0 : SPECTRE_ALWAYS_INLINE void pup(
608 : PUP::er& /*p*/, // NOLINT
609 : Variables<tmpl::list<>>& /* unused */) { // NOLINT
610 : }
611 0 : SPECTRE_ALWAYS_INLINE void operator|(
612 : PUP::er& /*p*/, Variables<tmpl::list<>>& /* unused */) { // NOLINT
613 : }
614 :
615 0 : SPECTRE_ALWAYS_INLINE bool operator==(const Variables<tmpl::list<>>& /*lhs*/,
616 : const Variables<tmpl::list<>>& /*rhs*/) {
617 : return true;
618 : }
619 0 : SPECTRE_ALWAYS_INLINE bool operator!=(const Variables<tmpl::list<>>& /*lhs*/,
620 : const Variables<tmpl::list<>>& /*rhs*/) {
621 : return false;
622 : }
623 :
624 0 : inline std::ostream& operator<<(std::ostream& os,
625 : const Variables<tmpl::list<>>& /*d*/) {
626 : return os << "{}";
627 : }
628 :
629 : template <typename... Tags>
630 : Variables<tmpl::list<Tags...>>::Variables() {
631 : // This makes an assertion trigger if one tries to assign to
632 : // components of a default-constructed Variables.
633 : const auto set_refs = [](auto& var) {
634 : for (auto& dv : var) {
635 : dv.set_data_ref(nullptr, 0);
636 : }
637 : return 0;
638 : };
639 : (void)set_refs;
640 : expand_pack(set_refs(tuples::get<Tags>(reference_variable_data_))...);
641 : }
642 :
643 : template <typename... Tags>
644 : Variables<tmpl::list<Tags...>>::Variables(const size_t number_of_grid_points) {
645 : initialize(number_of_grid_points);
646 : }
647 :
648 : template <typename... Tags>
649 : Variables<tmpl::list<Tags...>>::Variables(const size_t number_of_grid_points,
650 : const value_type value) {
651 : initialize(number_of_grid_points, value);
652 : }
653 :
654 : template <typename... Tags>
655 : void Variables<tmpl::list<Tags...>>::initialize(
656 : const size_t number_of_grid_points) {
657 : if (number_of_grid_points_ == 0) {
658 : variable_data_impl_dynamic_.reset();
659 : size_ = 0;
660 : number_of_grid_points_ = 0;
661 : }
662 : if (number_of_grid_points_ == number_of_grid_points) {
663 : return;
664 : }
665 : if (UNLIKELY(not is_owning())) {
666 : ERROR("Variables::initialize cannot be called on a non-owning Variables. "
667 : "This likely happened because of an attempted resize. The current "
668 : "number of grid points is " << number_of_grid_points_ << " and the "
669 : "requested number is " << number_of_grid_points << ".");
670 : }
671 : number_of_grid_points_ = number_of_grid_points;
672 : size_ = number_of_grid_points * number_of_independent_components;
673 : if (size_ > 0) {
674 : if (number_of_grid_points_ == 1) {
675 : variable_data_impl_dynamic_.reset();
676 : } else {
677 : variable_data_impl_dynamic_ =
678 : cpp20::make_unique_for_overwrite<value_type[]>(size_);
679 : }
680 : add_reference_variable_data();
681 : #if defined(SPECTRE_DEBUG) || defined(SPECTRE_NAN_INIT)
682 : std::fill(variable_data_.data(), variable_data_.data() + size_,
683 : make_signaling_NaN<value_type>());
684 : #endif // SPECTRE_DEBUG
685 : }
686 : }
687 :
688 : template <typename... Tags>
689 : void Variables<tmpl::list<Tags...>>::initialize(
690 : const size_t number_of_grid_points, const value_type value) {
691 : initialize(number_of_grid_points);
692 : std::fill(variable_data_.data(), variable_data_.data() + size_, value);
693 : }
694 :
695 : /// \cond HIDDEN_SYMBOLS
696 : template <typename... Tags>
697 : Variables<tmpl::list<Tags...>>::Variables(
698 : const Variables<tmpl::list<Tags...>>& rhs) {
699 : initialize(rhs.number_of_grid_points());
700 : variable_data_ =
701 : static_cast<const blaze::Vector<pointer_type, transpose_flag>&>(
702 : rhs.variable_data_);
703 : }
704 :
705 : template <typename... Tags>
706 : Variables<tmpl::list<Tags...>>& Variables<tmpl::list<Tags...>>::operator=(
707 : const Variables<tmpl::list<Tags...>>& rhs) {
708 : if (&rhs == this) {
709 : return *this;
710 : }
711 : initialize(rhs.number_of_grid_points());
712 : variable_data_ =
713 : static_cast<const blaze::Vector<pointer_type, transpose_flag>&>(
714 : rhs.variable_data_);
715 : return *this;
716 : }
717 :
718 : template <typename... Tags>
719 : Variables<tmpl::list<Tags...>>::Variables(Variables<tmpl::list<Tags...>>&& rhs)
720 : : variable_data_impl_dynamic_(std::move(rhs.variable_data_impl_dynamic_)),
721 : owning_(rhs.owning_),
722 : size_(rhs.size()),
723 : number_of_grid_points_(rhs.number_of_grid_points()),
724 : variable_data_(std::move(rhs.variable_data_)) {
725 : if (number_of_grid_points_ == 1) {
726 : variable_data_impl_static_ = std::move(rhs.variable_data_impl_static_);
727 : }
728 : rhs.variable_data_impl_dynamic_.reset();
729 : rhs.owning_ = true;
730 : rhs.size_ = 0;
731 : rhs.number_of_grid_points_ = 0;
732 : add_reference_variable_data();
733 : }
734 :
735 : template <typename... Tags>
736 : Variables<tmpl::list<Tags...>>& Variables<tmpl::list<Tags...>>::operator=(
737 : Variables<tmpl::list<Tags...>>&& rhs) {
738 : if (this == &rhs) {
739 : return *this;
740 : }
741 : owning_ = rhs.owning_;
742 : size_ = rhs.size_;
743 : number_of_grid_points_ = std::move(rhs.number_of_grid_points_);
744 : variable_data_ = std::move(rhs.variable_data_);
745 : variable_data_impl_dynamic_ = std::move(rhs.variable_data_impl_dynamic_);
746 : if (number_of_grid_points_ == 1) {
747 : variable_data_impl_static_ = std::move(rhs.variable_data_impl_static_);
748 : }
749 :
750 : rhs.variable_data_impl_dynamic_.reset();
751 : rhs.owning_ = true;
752 : rhs.size_ = 0;
753 : rhs.number_of_grid_points_ = 0;
754 : add_reference_variable_data();
755 : return *this;
756 : }
757 :
758 : template <typename... Tags>
759 : template <typename... WrappedTags,
760 : Requires<tmpl2::flat_all<
761 : std::is_same<db::remove_all_prefixes<WrappedTags>,
762 : db::remove_all_prefixes<Tags>>::value...>::value>>
763 : Variables<tmpl::list<Tags...>>::Variables(
764 : const Variables<tmpl::list<WrappedTags...>>& rhs) {
765 : static_assert(
766 : (std::is_same_v<typename Tags::type, typename WrappedTags::type> and ...),
767 : "Tensor types do not match!");
768 : initialize(rhs.number_of_grid_points());
769 : variable_data_ =
770 : static_cast<const blaze::Vector<pointer_type, transpose_flag>&>(
771 : rhs.variable_data_);
772 : }
773 :
774 : template <typename... Tags>
775 : template <typename... WrappedTags,
776 : Requires<tmpl2::flat_all<
777 : std::is_same<db::remove_all_prefixes<WrappedTags>,
778 : db::remove_all_prefixes<Tags>>::value...>::value>>
779 : Variables<tmpl::list<Tags...>>& Variables<tmpl::list<Tags...>>::operator=(
780 : const Variables<tmpl::list<WrappedTags...>>& rhs) {
781 : static_assert(
782 : (std::is_same_v<typename Tags::type, typename WrappedTags::type> and ...),
783 : "Tensor types do not match!");
784 : initialize(rhs.number_of_grid_points());
785 : variable_data_ =
786 : static_cast<const blaze::Vector<pointer_type, transpose_flag>&>(
787 : rhs.variable_data_);
788 : return *this;
789 : }
790 :
791 : template <typename... Tags>
792 : template <typename... WrappedTags,
793 : Requires<tmpl2::flat_all<
794 : std::is_same<db::remove_all_prefixes<WrappedTags>,
795 : db::remove_all_prefixes<Tags>>::value...>::value>>
796 : Variables<tmpl::list<Tags...>>::Variables(
797 : Variables<tmpl::list<WrappedTags...>>&& rhs)
798 : : variable_data_impl_dynamic_(std::move(rhs.variable_data_impl_dynamic_)),
799 : owning_(rhs.owning_),
800 : size_(rhs.size()),
801 : number_of_grid_points_(rhs.number_of_grid_points()),
802 : variable_data_(std::move(rhs.variable_data_)) {
803 : static_assert(
804 : (std::is_same_v<typename Tags::type, typename WrappedTags::type> and ...),
805 : "Tensor types do not match!");
806 : if (number_of_grid_points_ == 1) {
807 : variable_data_impl_static_ = std::move(rhs.variable_data_impl_static_);
808 : }
809 : rhs.variable_data_impl_dynamic_.reset();
810 : rhs.size_ = 0;
811 : rhs.owning_ = true;
812 : rhs.number_of_grid_points_ = 0;
813 : add_reference_variable_data();
814 : }
815 :
816 : template <typename... Tags>
817 : template <typename... WrappedTags,
818 : Requires<tmpl2::flat_all<
819 : std::is_same<db::remove_all_prefixes<WrappedTags>,
820 : db::remove_all_prefixes<Tags>>::value...>::value>>
821 : Variables<tmpl::list<Tags...>>& Variables<tmpl::list<Tags...>>::operator=(
822 : Variables<tmpl::list<WrappedTags...>>&& rhs) {
823 : static_assert(
824 : (std::is_same_v<typename Tags::type, typename WrappedTags::type> and ...),
825 : "Tensor types do not match!");
826 : variable_data_ = std::move(rhs.variable_data_);
827 : owning_ = rhs.owning_;
828 : size_ = rhs.size_;
829 : number_of_grid_points_ = std::move(rhs.number_of_grid_points_);
830 : variable_data_impl_dynamic_ = std::move(rhs.variable_data_impl_dynamic_);
831 : if (number_of_grid_points_ == 1) {
832 : variable_data_impl_static_ = std::move(rhs.variable_data_impl_static_);
833 : }
834 :
835 : rhs.variable_data_impl_dynamic_.reset();
836 : rhs.size_ = 0;
837 : rhs.owning_ = true;
838 : rhs.number_of_grid_points_ = 0;
839 : add_reference_variable_data();
840 : return *this;
841 : }
842 :
843 : template <typename... Tags>
844 : Variables<tmpl::list<Tags...>>::Variables(const pointer start,
845 : const size_t size) {
846 : set_data_ref(start, size);
847 : }
848 :
849 : template <typename... Tags>
850 : void Variables<tmpl::list<Tags...>>::pup(PUP::er& p) {
851 : ASSERT(
852 : owning_,
853 : "Cannot pup a non-owning Variables! It may be reasonable to pack a "
854 : "non-owning Variables, but not to unpack one. This should be discussed "
855 : "in an issue with the core devs if the feature seems necessary.");
856 : size_t number_of_grid_points = number_of_grid_points_;
857 : p | number_of_grid_points;
858 : if (p.isUnpacking()) {
859 : initialize(number_of_grid_points);
860 : }
861 : PUParray(p, variable_data_.data(), size_);
862 : }
863 : /// \endcond
864 :
865 : /// @{
866 : /*!
867 : * \ingroup DataStructuresGroup
868 : * \brief Return Tag::type pointing into the contiguous array
869 : *
870 : * \tparam Tag the variable to return
871 : */
872 : template <typename Tag, typename TagList>
873 1 : constexpr typename Tag::type& get(Variables<TagList>& v) {
874 : static_assert(tmpl::list_contains_v<TagList, Tag>,
875 : "Could not retrieve Tag from Variables. See the first "
876 : "template parameter of the instantiation for what Tag is "
877 : "being retrieved and the second template parameter for "
878 : "what Tags are available.");
879 : return tuples::get<Tag>(v.reference_variable_data_);
880 : }
881 : template <typename Tag, typename TagList>
882 1 : constexpr const typename Tag::type& get(const Variables<TagList>& v) {
883 : static_assert(tmpl::list_contains_v<TagList, Tag>,
884 : "Could not retrieve Tag from Variables. See the first "
885 : "template parameter of the instantiation for what Tag is "
886 : "being retrieved and the second template parameter for "
887 : "what Tags are available.");
888 : return tuples::get<Tag>(v.reference_variable_data_);
889 : }
890 : /// @}
891 :
892 : template <typename... Tags>
893 : template <typename VT, bool VF>
894 : Variables<tmpl::list<Tags...>>::Variables(
895 : const blaze::DenseVector<VT, VF>& expression) {
896 : ASSERT((*expression).size() % number_of_independent_components == 0,
897 : "Invalid size " << (*expression).size() << " for a Variables with "
898 : << number_of_independent_components << " components.");
899 : initialize((*expression).size() / number_of_independent_components);
900 : variable_data_ = expression;
901 : }
902 :
903 : /// \cond
904 : template <typename... Tags>
905 : template <typename VT, bool VF>
906 : Variables<tmpl::list<Tags...>>& Variables<tmpl::list<Tags...>>::operator=(
907 : const blaze::DenseVector<VT, VF>& expression) {
908 : ASSERT((*expression).size() % number_of_independent_components == 0,
909 : "Invalid size " << (*expression).size() << " for a Variables with "
910 : << number_of_independent_components << " components.");
911 : initialize((*expression).size() / number_of_independent_components);
912 : variable_data_ = expression;
913 : return *this;
914 : }
915 : /// \endcond
916 :
917 : /// \cond HIDDEN_SYMBOLS
918 : template <typename... Tags>
919 : void Variables<tmpl::list<Tags...>>::add_reference_variable_data() {
920 : if (size_ == 0) {
921 : return;
922 : }
923 : if (is_owning()) {
924 : if (number_of_grid_points_ == 1) {
925 : variable_data_.reset(variable_data_impl_static_.data(), size_);
926 : } else {
927 : variable_data_.reset(variable_data_impl_dynamic_.get(), size_);
928 : }
929 : }
930 : ASSERT(variable_data_.size() == size_ and
931 : size_ == number_of_grid_points_ * number_of_independent_components,
932 : "Size mismatch: variable_data_.size() = "
933 : << variable_data_.size() << " size_ = " << size_ << " should be: "
934 : << number_of_grid_points_ * number_of_independent_components
935 : << "\nThis is an internal inconsistency bug in Variables. Please "
936 : "file an issue.");
937 : size_t variable_offset = 0;
938 : tmpl::for_each<tags_list>([this, &variable_offset](auto tag_v) {
939 : using Tag = tmpl::type_from<decltype(tag_v)>;
940 : auto& var = tuples::get<Tag>(reference_variable_data_);
941 : for (size_t i = 0; i < Tag::type::size(); ++i) {
942 : var[i].set_data_ref(
943 : &variable_data_[variable_offset++ * number_of_grid_points_],
944 : number_of_grid_points_);
945 : }
946 : });
947 : }
948 : /// \endcond
949 :
950 : template <typename... Tags>
951 0 : Variables<tmpl::list<Tags...>>& operator*=(
952 : Variables<tmpl::list<Tags...>>& lhs,
953 : const typename Variables<tmpl::list<Tags...>>::vector_type& rhs) {
954 : using value_type = typename Variables<tmpl::list<Tags...>>::value_type;
955 : ASSERT(lhs.number_of_grid_points() == rhs.size(),
956 : "Size mismatch in multiplication: " << lhs.number_of_grid_points()
957 : << " and " << rhs.size());
958 : value_type* const lhs_data = lhs.data();
959 : const value_type* const rhs_data = rhs.data();
960 : for (size_t c = 0; c < lhs.number_of_independent_components; ++c) {
961 : for (size_t s = 0; s < lhs.number_of_grid_points(); ++s) {
962 : // clang-tidy: do not use pointer arithmetic
963 : lhs_data[c * lhs.number_of_grid_points() + s] *= rhs_data[s]; // NOLINT
964 : }
965 : }
966 : return lhs;
967 : }
968 :
969 : template <typename... Tags>
970 0 : Variables<tmpl::list<Tags...>> operator*(
971 : const Variables<tmpl::list<Tags...>>& lhs,
972 : const typename Variables<tmpl::list<Tags...>>::vector_type& rhs) {
973 : auto result = lhs;
974 : result *= rhs;
975 : return result;
976 : }
977 :
978 : template <typename... Tags>
979 0 : Variables<tmpl::list<Tags...>> operator*(
980 : const typename Variables<tmpl::list<Tags...>>::vector_type& lhs,
981 : const Variables<tmpl::list<Tags...>>& rhs) {
982 : auto result = rhs;
983 : result *= lhs;
984 : return result;
985 : }
986 :
987 : template <typename... Tags>
988 0 : Variables<tmpl::list<Tags...>>& operator/=(
989 : Variables<tmpl::list<Tags...>>& lhs,
990 : const typename Variables<tmpl::list<Tags...>>::vector_type& rhs) {
991 : ASSERT(lhs.number_of_grid_points() == rhs.size(),
992 : "Size mismatch in multiplication: " << lhs.number_of_grid_points()
993 : << " and " << rhs.size());
994 : using value_type = typename Variables<tmpl::list<Tags...>>::value_type;
995 : value_type* const lhs_data = lhs.data();
996 : const value_type* const rhs_data = rhs.data();
997 : for (size_t c = 0; c < lhs.number_of_independent_components; ++c) {
998 : for (size_t s = 0; s < lhs.number_of_grid_points(); ++s) {
999 : // clang-tidy: do not use pointer arithmetic
1000 : lhs_data[c * lhs.number_of_grid_points() + s] /= rhs_data[s]; // NOLINT
1001 : }
1002 : }
1003 : return lhs;
1004 : }
1005 :
1006 : template <typename... Tags>
1007 0 : Variables<tmpl::list<Tags...>> operator/(
1008 : const Variables<tmpl::list<Tags...>>& lhs,
1009 : const typename Variables<tmpl::list<Tags...>>::vector_type& rhs) {
1010 : auto result = lhs;
1011 : result /= rhs;
1012 : return result;
1013 : }
1014 :
1015 : template <typename... Tags, Requires<sizeof...(Tags) != 0> = nullptr>
1016 0 : std::ostream& operator<<(std::ostream& os,
1017 : const Variables<tmpl::list<Tags...>>& d) {
1018 : size_t count = 0;
1019 : const auto helper = [&os, &d, &count](auto tag_v) {
1020 : count++;
1021 : using Tag = typename decltype(tag_v)::type;
1022 : os << db::tag_name<Tag>() << ":\n";
1023 : os << get<Tag>(d);
1024 : if (count < sizeof...(Tags)) {
1025 : os << "\n\n";
1026 : }
1027 : };
1028 : EXPAND_PACK_LEFT_TO_RIGHT(helper(tmpl::type_<Tags>{}));
1029 : return os;
1030 : }
1031 :
1032 : template <typename TagsList>
1033 0 : bool operator!=(const Variables<TagsList>& lhs,
1034 : const Variables<TagsList>& rhs) {
1035 : return not(lhs == rhs);
1036 : }
1037 :
1038 : template <typename... TagsLhs, typename... TagsRhs,
1039 : Requires<not std::is_same<tmpl::list<TagsLhs...>,
1040 : tmpl::list<TagsRhs...>>::value> = nullptr>
1041 0 : void swap(Variables<tmpl::list<TagsLhs...>>& lhs,
1042 : Variables<tmpl::list<TagsRhs...>>& rhs) {
1043 : Variables<tmpl::list<TagsLhs...>> temp{std::move(lhs)};
1044 : lhs = std::move(rhs);
1045 : rhs = std::move(temp);
1046 : }
1047 :
1048 : /// \ingroup DataStructuresGroup
1049 : /// Construct a variables from the `Tensor`s in a `TaggedTuple`.
1050 : template <typename... Tags>
1051 1 : Variables<tmpl::list<Tags...>> variables_from_tagged_tuple(
1052 : const tuples::TaggedTuple<Tags...>& tuple) {
1053 : auto result = make_with_value<Variables<tmpl::list<Tags...>>>(
1054 : get<tmpl::front<tmpl::list<Tags...>>>(tuple), 0.0);
1055 : result.assign_subset(tuple);
1056 : return result;
1057 : }
1058 :
1059 : namespace MakeWithValueImpls {
1060 : template <typename TagList>
1061 0 : struct MakeWithSize<Variables<TagList>> {
1062 0 : static SPECTRE_ALWAYS_INLINE Variables<TagList> apply(
1063 : const size_t size, const typename Variables<TagList>::value_type value) {
1064 : return Variables<TagList>(size, value);
1065 : }
1066 : };
1067 :
1068 : template <typename TagList>
1069 0 : struct NumberOfPoints<Variables<TagList>> {
1070 0 : static SPECTRE_ALWAYS_INLINE size_t apply(const Variables<TagList>& input) {
1071 : return input.number_of_grid_points();
1072 : }
1073 : };
1074 : } // namespace MakeWithValueImpls
1075 :
1076 : template <typename TagList>
1077 0 : struct SetNumberOfGridPointsImpls::SetNumberOfGridPointsImpl<
1078 : Variables<TagList>> {
1079 0 : static constexpr bool is_trivial = false;
1080 0 : static SPECTRE_ALWAYS_INLINE void apply(
1081 : const gsl::not_null<Variables<TagList>*> result, const size_t size) {
1082 : result->initialize(size);
1083 : }
1084 : };
1085 :
1086 1 : namespace EqualWithinRoundoffImpls {
1087 : // It would be nice to use `blaze::equal` in these implementations, but it
1088 : // doesn't currently allow to specify a tolerance. See upstream issue:
1089 : // https://bitbucket.org/blaze-lib/blaze/issues/417/adjust-relaxed-equal-accuracy
1090 : template <typename TagList, typename Floating>
1091 : struct EqualWithinRoundoffImpl<Variables<TagList>, Floating,
1092 : Requires<std::is_floating_point_v<Floating>>> {
1093 : static bool apply(const Variables<TagList>& lhs, const Floating& rhs,
1094 : const double eps, const double scale) {
1095 : for (size_t i = 0; i < lhs.size(); ++i) {
1096 : if (not equal_within_roundoff(lhs.data()[i], rhs, eps, scale)) {
1097 : return false;
1098 : }
1099 : }
1100 : return true;
1101 : }
1102 : };
1103 : template <typename TagList, typename Floating>
1104 : struct EqualWithinRoundoffImpl<Floating, Variables<TagList>,
1105 : Requires<std::is_floating_point_v<Floating>>> {
1106 : static SPECTRE_ALWAYS_INLINE bool apply(const Floating& lhs,
1107 : const Variables<TagList>& rhs,
1108 : const double eps,
1109 : const double scale) {
1110 : return equal_within_roundoff(rhs, lhs, eps, scale);
1111 : }
1112 : };
1113 : template <typename LhsTagList, typename RhsTagList>
1114 0 : struct EqualWithinRoundoffImpl<Variables<LhsTagList>, Variables<RhsTagList>> {
1115 0 : static bool apply(const Variables<LhsTagList>& lhs,
1116 : const Variables<RhsTagList>& rhs, const double eps,
1117 : const double scale) {
1118 : ASSERT(lhs.size() == rhs.size(),
1119 : "Can only compare two Variables of the same size, but lhs has size "
1120 : << lhs.size() << " and rhs has size " << rhs.size() << ".");
1121 : for (size_t i = 0; i < lhs.size(); ++i) {
1122 : if (not equal_within_roundoff(lhs.data()[i], rhs.data()[i], eps, scale)) {
1123 : return false;
1124 : }
1125 : }
1126 : return true;
1127 : }
1128 : };
1129 : } // namespace EqualWithinRoundoffImpls
1130 :
1131 : namespace db {
1132 : // Enable subitems for ::Tags::Variables and derived tags (e.g. compute tags).
1133 : // Other tags that hold a `Variables` don't expose the constituent tensors as
1134 : // subitems by default.
1135 : namespace Variables_detail {
1136 : // Check if the argument is a `::Tags::Variables`, or derived from it. Can't use
1137 : // `tt:is_a_v` because we also want to match derived classes. Can't use
1138 : // `std::is_base_of` because `::Tags::Variables` is a template.
1139 : template <typename TagsList>
1140 : constexpr std::true_type is_a_variables_tag(::Tags::Variables<TagsList>&&) {
1141 : return {};
1142 : }
1143 : constexpr std::false_type is_a_variables_tag(...) { return {}; }
1144 : template <typename Tag>
1145 : static constexpr bool is_a_variables_tag_v =
1146 : decltype(is_a_variables_tag(std::declval<Tag>()))::value;
1147 : } // namespace Variables_detail
1148 : template <typename Tag>
1149 : struct Subitems<Tag, Requires<Variables_detail::is_a_variables_tag_v<Tag>>> {
1150 : using type = typename Tag::type::tags_list;
1151 :
1152 : template <typename Subtag, typename LocalTag = Tag>
1153 : static void create_item(
1154 : const gsl::not_null<typename LocalTag::type*> parent_value,
1155 : const gsl::not_null<typename Subtag::type*> sub_value) {
1156 : auto& vars = get<Subtag>(*parent_value);
1157 : // Only update the Tensor if the Variables has changed its allocation
1158 : if constexpr (not is_any_spin_weighted_v<typename Subtag::type::type>) {
1159 : if (vars.begin()->data() != sub_value->begin()->data()) {
1160 : for (auto vars_it = vars.begin(), sub_var_it = sub_value->begin();
1161 : vars_it != vars.end(); ++vars_it, ++sub_var_it) {
1162 : sub_var_it->set_data_ref(make_not_null(&*vars_it));
1163 : }
1164 : }
1165 : } else {
1166 : if (vars.begin()->data().data() != sub_value->begin()->data().data()) {
1167 : for (auto vars_it = vars.begin(), sub_var_it = sub_value->begin();
1168 : vars_it != vars.end(); ++vars_it, ++sub_var_it) {
1169 : sub_var_it->set_data_ref(make_not_null(&*vars_it));
1170 : }
1171 : }
1172 : }
1173 : }
1174 :
1175 : template <typename Subtag>
1176 : static const typename Subtag::type& create_compute_item(
1177 : const typename Tag::type& parent_value) {
1178 : return get<Subtag>(parent_value);
1179 : }
1180 : };
1181 : } // namespace db
1182 :
1183 : namespace Variables_detail {
1184 : template <typename Tags>
1185 : using MathWrapperVectorType = tmpl::conditional_t<
1186 : std::is_same_v<typename Variables<Tags>::value_type, double>, DataVector,
1187 : ComplexDataVector>;
1188 : } // namespace Variables_detail
1189 :
1190 : template <typename Tags>
1191 0 : auto make_math_wrapper(const gsl::not_null<Variables<Tags>*> data) {
1192 : using Vector = Variables_detail::MathWrapperVectorType<Tags>;
1193 : Vector referencing(data->data(), data->size());
1194 : return make_math_wrapper(&referencing);
1195 : }
1196 :
1197 : template <typename Tags>
1198 0 : auto make_math_wrapper(const Variables<Tags>& data) {
1199 : using Vector = Variables_detail::MathWrapperVectorType<Tags>;
1200 : const Vector referencing(
1201 : const_cast<typename Vector::value_type*>(data.data()), data.size());
1202 : return make_math_wrapper(referencing);
1203 : }
1204 :
1205 : /// \cond
1206 : template <size_t I, class... Tags>
1207 : inline constexpr typename tmpl::at_c<tmpl::list<Tags...>, I>::type&& get(
1208 : Variables<tmpl::list<Tags...>>&& t) {
1209 : return get<tmpl::at_c<tmpl::list<Tags...>, I>>(t);
1210 : }
1211 :
1212 : template <size_t I, class... Tags>
1213 : inline constexpr const typename tmpl::at_c<tmpl::list<Tags...>, I>::type& get(
1214 : const Variables<tmpl::list<Tags...>>& t) {
1215 : return get<tmpl::at_c<tmpl::list<Tags...>, I>>(t);
1216 : }
1217 :
1218 : template <size_t I, class... Tags>
1219 : inline constexpr typename tmpl::at_c<tmpl::list<Tags...>, I>::type& get(
1220 : Variables<tmpl::list<Tags...>>& t) {
1221 : return get<tmpl::at_c<tmpl::list<Tags...>, I>>(t);
1222 : }
1223 : /// \endcond
1224 :
1225 : namespace std {
1226 : template <typename... Tags>
1227 : struct tuple_size<Variables<tmpl::list<Tags...>>>
1228 : : std::integral_constant<int, sizeof...(Tags)> {};
1229 : template <size_t I, typename... Tags>
1230 : struct tuple_element<I, Variables<tmpl::list<Tags...>>> {
1231 : using type = typename tmpl::at_c<tmpl::list<Tags...>, I>::type;
1232 : };
1233 : } // namespace std
1234 :
1235 : template <typename TagsList>
1236 0 : bool contains_allocations(const Variables<TagsList>& value) {
1237 : return value.number_of_grid_points() > 1;
1238 : }
1239 :
1240 0 : inline bool contains_allocations(const Variables<tmpl::list<>>& /*value*/) {
1241 : return false;
1242 : }
|
