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