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