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