CceComputationTestHelpers.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 
8 #include "DataStructures/ComplexDataVector.hpp"
9 #include "DataStructures/ComplexModalVector.hpp"
12 #include "DataStructures/DataBox/Tag.hpp"
13 #include "DataStructures/Tags.hpp"
15 #include "Evolution/Systems/Cce/IntegrandInputSteps.hpp"
16 #include "Evolution/Systems/Cce/Tags.hpp"
20 #include "NumericalAlgorithms/Spectral/SwshCollocation.hpp"
21 #include "NumericalAlgorithms/Spectral/SwshDerivatives.hpp"
22 #include "NumericalAlgorithms/Spectral/SwshFiltering.hpp"
23 #include "Utilities/Gsl.hpp"
24 #include "Utilities/TMPL.hpp"
25 #include "Utilities/VectorAlgebra.hpp"
26 
27 namespace Cce {
28 namespace TestHelpers {
29 
30 // For representing a primitive series of powers in inverse r for diagnostic
31 // computations
32 template <typename Tag>
34  using type = Scalar<ComplexModalVector>;
35  using tag = Tag;
36 };
37 
38 // For representing the angular function in a separable quantity.
39 template <typename Tag>
41  using type = typename Tag::type;
42  using tag = Tag;
43 };
44 
45 // shortcut method for evaluating a frequently used radial quantity in CCE
46 void volume_one_minus_y(
48  size_t l_max) noexcept;
49 
50 // explicit power method avoids behavior of Blaze to occasionally FPE on
51 // powers of complex that operate fine when repeated multiplication is used
52 // instead.
53 ComplexDataVector power(const ComplexDataVector& value,
54  size_t exponent) noexcept;
55 
56 // A utility for copying a set of tags from one DataBox to another. Useful in
57 // the CCE tests, where often an expected input needs to be copied to the box
58 // which is to be tested.
59 template <typename... Tags>
61  template <typename FromDataBox, typename ToDataBox>
62  static void apply(const gsl::not_null<ToDataBox*> to_data_box,
63  const FromDataBox& from_data_box) noexcept {
64  db::mutate<Tags...>(
65  to_data_box,
66  [](const gsl::not_null<db::item_type<Tags>*>... to_value,
67  const typename Tags::type&... from_value) noexcept {
68  const auto assign = [](auto to, const auto& from) noexcept {
69  *to = from;
70  return 0;
71  };
72  expand_pack(assign(to_value, from_value)...);
73  },
74  db::get<Tags>(from_data_box)...);
75  }
76 };
77 
78 // Given the angular and radial dependence separately, assembles the volume data
79 // by computing the tensor product.
80 void generate_volume_data_from_separated_values(
82  gsl::not_null<ComplexDataVector*> one_divided_by_r,
83  const ComplexDataVector& angular_collocation,
84  const ComplexModalVector& radial_coefficients, size_t l_max,
85  size_t number_of_radial_grid_points) noexcept;
86 
87 // A utility for separately verifying the values in several computation
88 // routines in the CCE quantity derivations. These are an independent
89 // computation of similar quantities needed in the CCE systems under the
90 // assumption that the values are separable, i.e. can be written as F(theta,
91 // phi, r) = f(theta, phi) * g(r), and performing semi-analytic manipulations
92 // using simplifications from separability and an explicit decomposition of g in
93 // inverse powers of r.
94 template <typename Tag>
96 
97 template <typename Tag>
98 struct CalculateSeparatedTag<Tags::Dy<Tag>> {
99  // d x/dy = d x /dr * (2 R / (1-y)^2) = dx/dr * r^2 / (2 R)
100  // So, for x = r^(-n),
101  // d (r^(-n))/dy = -n * r^(-n - 1) * r^2 / (2 R) = -n / (2 R) r^(-n + 1)
102  // except when n=0, then the result is zero. So, the entire process actually
103  // moves the polynomials down in order
104  template <typename AngularTagList, typename RadialCoefficientTagList>
105  void operator()(
106  const gsl::not_null<Variables<AngularTagList>*> angular_collocation,
107  const gsl::not_null<Variables<RadialCoefficientTagList>*>
108  radial_coefficients,
109  const ComplexDataVector& /*one_divided_by_r*/,
110  const ComplexDataVector& boundary_r, const size_t /*l_max*/) noexcept {
111  get(get<AngularCollocationsFor<Tags::Dy<Tag>>>(*angular_collocation))
112  .data() =
113  get(get<AngularCollocationsFor<Tag>>(*angular_collocation)).data() /
114  (2.0 * boundary_r);
115 
116  ComplexModalVector& dy_radial_values = get(
117  get<RadialPolyCoefficientsFor<Tags::Dy<Tag>>>(*radial_coefficients));
118  const ComplexModalVector& radial_values =
119  get(get<RadialPolyCoefficientsFor<Tag>>(*radial_coefficients));
120 
121  for (size_t radial_power = 1;
122  radial_power < radial_coefficients->number_of_grid_points();
123  ++radial_power) {
124  dy_radial_values[radial_power - 1] =
125  -static_cast<double>(radial_power) * radial_values[radial_power];
126  }
127  dy_radial_values[radial_coefficients->number_of_grid_points() - 1] = 0.0;
128  }
129 };
130 
131 template <typename Tag, typename DerivKind>
132 struct CalculateSeparatedTag<Spectral::Swsh::Tags::Derivative<Tag, DerivKind>> {
133  template <typename AngularTagList, typename RadialCoefficientTagList>
134  void operator()(
135  const gsl::not_null<Variables<AngularTagList>*> angular_collocation,
136  const gsl::not_null<Variables<RadialCoefficientTagList>*>
137  radial_coefficients,
138  const ComplexDataVector& /*one_divided_by_r*/,
139  const ComplexDataVector& /*boundary_r*/, const size_t l_max) noexcept {
142  *angular_collocation)) =
143  Spectral::Swsh::angular_derivative<DerivKind>(
144  l_max, 1,
145  get(get<AngularCollocationsFor<Tag>>(*angular_collocation)));
146  // The spin-weighted derivatives are evaluated at constant r, so the radial
147  // coefficients are unaltered.
150  *radial_coefficients)) =
151  get(get<RadialPolyCoefficientsFor<Tag>>(*radial_coefficients));
152  }
153 };
154 
155 template <typename LhsTag, typename RhsTag>
156 struct CalculateSeparatedTag<::Tags::Multiplies<LhsTag, RhsTag>> {
157  template <typename AngularTagList, typename RadialCoefficientTagList>
158  void operator()(
159  const gsl::not_null<Variables<AngularTagList>*> angular_collocation,
160  const gsl::not_null<Variables<RadialCoefficientTagList>*>
161  radial_coefficients,
162  const ComplexDataVector& /*one_divided_by_r*/,
163  const ComplexDataVector& /*boundary_r*/,
164  const size_t /*l_max*/) noexcept {
166  *angular_collocation)) =
167  get(get<AngularCollocationsFor<LhsTag>>(*angular_collocation)) *
168  get(get<AngularCollocationsFor<RhsTag>>(*angular_collocation));
169 
170  ComplexModalVector& multiplied_radial_values =
172  *radial_coefficients));
173  const ComplexModalVector& lhs_radial_values =
174  get(get<RadialPolyCoefficientsFor<LhsTag>>(*radial_coefficients));
175  const ComplexModalVector& rhs_radial_values =
176  get(get<RadialPolyCoefficientsFor<RhsTag>>(*radial_coefficients));
177  multiplied_radial_values = 0.0;
178  for (size_t lhs_radial_power = 0;
179  lhs_radial_power < radial_coefficients->number_of_grid_points();
180  ++lhs_radial_power) {
181  for (size_t rhs_radial_power = 0;
182  rhs_radial_power <
183  radial_coefficients->number_of_grid_points() - lhs_radial_power;
184  ++rhs_radial_power) {
185  multiplied_radial_values[lhs_radial_power + rhs_radial_power] +=
186  lhs_radial_values[lhs_radial_power] *
187  rhs_radial_values[rhs_radial_power];
188  }
189  }
190  }
191 };
192 
193 template <>
194 struct CalculateSeparatedTag<Tags::BondiJbar> {
195  template <typename AngularTagList, typename RadialCoefficientTagList>
196  void operator()(
197  const gsl::not_null<Variables<AngularTagList>*> angular_collocation,
198  const gsl::not_null<Variables<RadialCoefficientTagList>*>
199  radial_coefficients,
200  const ComplexDataVector& /*one_divided_by_r*/,
201  const ComplexDataVector& /*boundary_r*/,
202  const size_t /*l_max*/) noexcept {
203  get(get<AngularCollocationsFor<Tags::BondiJbar>>(*angular_collocation)) =
204  conj(get(
205  get<AngularCollocationsFor<Tags::BondiJ>>(*angular_collocation)));
206  get(get<RadialPolyCoefficientsFor<Tags::BondiJbar>>(*radial_coefficients)) =
208  *radial_coefficients)));
209  }
210 };
211 
212 template <>
213 struct CalculateSeparatedTag<Tags::BondiUbar> {
214  template <typename AngularTagList, typename RadialCoefficientTagList>
215  void operator()(
216  const gsl::not_null<Variables<AngularTagList>*> angular_collocation,
217  const gsl::not_null<Variables<RadialCoefficientTagList>*>
218  radial_coefficients,
219  const ComplexDataVector& /*one_divided_by_r*/,
220  const ComplexDataVector& /*boundary_r*/,
221  const size_t /*l_max*/) noexcept {
222  get(get<AngularCollocationsFor<Tags::BondiUbar>>(*angular_collocation)) =
223  conj(get(
224  get<AngularCollocationsFor<Tags::BondiU>>(*angular_collocation)));
225  get(get<RadialPolyCoefficientsFor<Tags::BondiUbar>>(*radial_coefficients)) =
227  *radial_coefficients)));
228  }
229 };
230 
231 template <>
232 struct CalculateSeparatedTag<Tags::BondiQbar> {
233  template <typename AngularTagList, typename RadialCoefficientTagList>
234  void operator()(
235  const gsl::not_null<Variables<AngularTagList>*> angular_collocation,
236  const gsl::not_null<Variables<RadialCoefficientTagList>*>
237  radial_coefficients,
238  const ComplexDataVector& /*one_divided_by_r*/,
239  const ComplexDataVector& /*boundary_r*/,
240  const size_t /*l_max*/) noexcept {
241  get(get<AngularCollocationsFor<Tags::BondiQbar>>(*angular_collocation)) =
242  conj(get(
243  get<AngularCollocationsFor<Tags::BondiQ>>(*angular_collocation)));
244  get(get<RadialPolyCoefficientsFor<Tags::BondiQbar>>(*radial_coefficients)) =
246  *radial_coefficients)));
247  }
248 };
249 
250 // A generation function used in the tests of the CCE computations for the
251 // inputs to the integrand computation (`Test_PreSwshDerivatives.cpp`,
252 // and `Test_ComputeSwshDerivatives.cpp`). This function first generates the set
253 // of inputs specified by `InputTagList`, then emulates the cascaded
254 // computations of those utilities, using the tag lists as though the integrands
255 // that we wanted to compute were of the tags in `TargetTagList`. Instead of
256 // computing the values using the utilities in the main code base, though, the
257 // generation creates the expected values using the separable computations
258 // above. While this is also a fairly complicated computation (though simpler
259 // due to the assumption of separability), it is a completely independent method
260 // so can be regarded as a robust way of verifying the utilities in the main
261 // code base.
262 template <typename InputTagList, typename TargetTagList,
263  typename PreSwshDerivativesTagList, typename SwshDerivativesTagList,
264  typename AngularCollocationTagList,
265  typename PreSwshDerivativesRadialModeTagList, typename Generator>
266 void generate_separable_expected(
267  const gsl::not_null<Variables<PreSwshDerivativesTagList>*>
268  pre_swsh_derivatives,
269  const gsl::not_null<Variables<SwshDerivativesTagList>*> swsh_derivatives,
270  const gsl::not_null<Variables<AngularCollocationTagList>*>
271  angular_collocations,
272  const gsl::not_null<Variables<PreSwshDerivativesRadialModeTagList>*>
273  radial_modes,
274  const gsl::not_null<Generator*> generator,
275  const SpinWeighted<ComplexDataVector, 0>& boundary_r, const size_t l_max,
276  const size_t number_of_radial_points) noexcept {
280  Spectral::collocation_points<Spectral::Basis::Legendre,
281  Spectral::Quadrature::GaussLobatto>(
282  number_of_radial_points));
283 
284  // generate the separable variables
285  UniformCustomDistribution<double> dist(0.1, 1.0);
286 
287  ComplexDataVector one_divided_by_r =
288  (1.0 - y) / (2.0 * create_vector_of_n_copies(boundary_r.data(),
289  number_of_radial_points));
290 
291  // the generation step for the 'input' tags
292  tmpl::for_each<InputTagList>([
293  &angular_collocations, &radial_modes, &pre_swsh_derivatives, &dist,
294  &generator, &l_max, &number_of_radial_points, &one_divided_by_r
295  ](auto tag_v) noexcept {
296  using tag = typename decltype(tag_v)::type;
297  using radial_polynomial_tag = RadialPolyCoefficientsFor<tag>;
298  using angular_collocation_tag = AngularCollocationsFor<tag>;
299  // generate the angular part randomly
300  get(get<angular_collocation_tag>(*angular_collocations)).data() =
301  make_with_random_values<ComplexDataVector>(
302  generator, make_not_null(&dist),
304  Spectral::Swsh::filter_swsh_boundary_quantity(
306  &get(get<angular_collocation_tag>(*angular_collocations))),
307  l_max, 2);
308  // generate the radial part
309  auto& radial_polynomial = get(get<radial_polynomial_tag>(*radial_modes));
310  fill_with_random_values(make_not_null(&radial_polynomial), generator,
311  make_not_null(&dist));
312  // filtering to make sure the results are reasonable,
313  // filter function: exp(-10.0 * (i / N)**2)
314  for (size_t i = 0; i < radial_polynomial.size(); ++i) {
315  radial_polynomial[i] *= exp(
316  -10.0 * square(static_cast<double>(i) /
317  static_cast<double>(radial_polynomial.size() - 1)));
318  // aggressive Heaviside needed to get acceptable test precision at the low
319  // resolutions needed to keep tests fast.
320  if (i > 3) {
321  radial_polynomial[i] = 0.0;
322  }
323  }
324 
325  generate_volume_data_from_separated_values(
326  make_not_null(&get(get<tag>(*pre_swsh_derivatives)).data()),
327  make_not_null(&one_divided_by_r),
328  get(get<angular_collocation_tag>(*angular_collocations)).data(),
329  radial_polynomial, l_max, number_of_radial_points);
330  });
331 
332  // Compute the expected versions using the separable mathematics from above
333  // utilities
334  const auto calculate_separable_for_pre_swsh_derivative = [
335  &angular_collocations, &radial_modes, &pre_swsh_derivatives, &l_max,
336  &number_of_radial_points, &one_divided_by_r, &boundary_r
337  ](auto pre_swsh_derivative_tag_v) noexcept {
338  using pre_swsh_derivative_tag =
339  typename decltype(pre_swsh_derivative_tag_v)::type;
340  CalculateSeparatedTag<pre_swsh_derivative_tag>{}(
341  angular_collocations, radial_modes, one_divided_by_r, boundary_r.data(),
342  l_max);
343  generate_volume_data_from_separated_values(
345  &get(get<pre_swsh_derivative_tag>(*pre_swsh_derivatives)).data()),
346  make_not_null(&one_divided_by_r),
347  get(get<AngularCollocationsFor<pre_swsh_derivative_tag>>(
348  *angular_collocations))
349  .data(),
350  get(get<RadialPolyCoefficientsFor<pre_swsh_derivative_tag>>(
351  *radial_modes)),
352  l_max, number_of_radial_points);
353  };
354 
355  const auto calculate_separable_for_swsh_derivative = [
356  &angular_collocations, &radial_modes, &swsh_derivatives, &l_max,
357  &number_of_radial_points, &one_divided_by_r, &boundary_r
358  ](auto swsh_derivative_tag_v) noexcept {
359  using swsh_derivative_tag = typename decltype(swsh_derivative_tag_v)::type;
360  CalculateSeparatedTag<swsh_derivative_tag>{}(angular_collocations,
361  radial_modes, one_divided_by_r,
362  boundary_r.data(), l_max);
363  generate_volume_data_from_separated_values(
364  make_not_null(&get(get<swsh_derivative_tag>(*swsh_derivatives)).data()),
365  make_not_null(&one_divided_by_r),
366  get(get<AngularCollocationsFor<swsh_derivative_tag>>(
367  *angular_collocations))
368  .data(),
369  get(get<RadialPolyCoefficientsFor<swsh_derivative_tag>>(*radial_modes)),
370  l_max, number_of_radial_points);
371  };
372 
373  tmpl::for_each<TargetTagList>([
374  &calculate_separable_for_pre_swsh_derivative, &
375  calculate_separable_for_swsh_derivative
376  ](auto target_tag_v) noexcept {
377  using target_tag = typename decltype(target_tag_v)::type;
378  tmpl::for_each<typename detail::TagsToComputeForImpl<
379  target_tag>::pre_swsh_derivative_tags>(
380  calculate_separable_for_pre_swsh_derivative);
381  tmpl::for_each<typename detail::TagsToComputeForImpl<
382  target_tag>::swsh_derivative_tags>(
383  calculate_separable_for_swsh_derivative);
384  tmpl::for_each<typename detail::TagsToComputeForImpl<
385  target_tag>::second_swsh_derivative_tags>(
386  calculate_separable_for_swsh_derivative);
387  });
388 }
389 } // namespace TestHelpers
390 } // namespace Cce
expand_pack
constexpr void expand_pack(Ts &&...) noexcept
Allows zero-cost unordered expansion of a parameter.
Definition: TMPL.hpp:547
DataBoxTag.hpp
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:689
Spectral::collocation_points
const DataVector & collocation_points(const size_t num_points) noexcept
Collocation points.
Definition: Spectral.cpp:280
db::PrefixTag
Marks an item as being a prefix to another tag.
Definition: Tag.hpp:66
square
constexpr decltype(auto) square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:55
TestingFramework.hpp
fill_with_random_values
void fill_with_random_values(const gsl::not_null< T * > data, const gsl::not_null< UniformRandomBitGenerator * > generator, const gsl::not_null< RandomNumberDistribution * > distribution) noexcept
Fill an existing data structure with random values.
Definition: MakeWithRandomValues.hpp:144
MakeWithRandomValues.hpp
Cce::TestHelpers::RadialPolyCoefficientsFor
Definition: CceComputationTestHelpers.hpp:33
Spectral::Swsh::Tags::Derivative
Prefix tag representing the spin-weighted derivative of a spin-weighted scalar.
Definition: SwshTags.hpp:152
db::SimpleTag
Tags for the DataBox inherit from this type.
Definition: Tag.hpp:23
Spectral.hpp
SpinWeighted
Make a spin-weighted type T with spin-weight Spin. Mathematical operators are restricted to addition,...
Definition: SpinWeighted.hpp:24
Cce::TestHelpers::CalculateSeparatedTag
Definition: CceComputationTestHelpers.hpp:95
db::mutate
void mutate(const gsl::not_null< DataBox< TagList > * > box, Invokable &&invokable, Args &&... args) noexcept
Allows changing the state of one or more non-computed elements in the DataBox.
Definition: DataBox.hpp:1020
UniformCustomDistribution
A uniform distribution function object which redirects appropriately to either the std::uniform_int_d...
Definition: MakeWithRandomValues.hpp:117
DataBox.hpp
cstddef
Cce::Tags::Dy
The derivative with respect to the numerical coordinate , where is Bondi radius of the worldtube.
Definition: Tags.hpp:115
Cce::TestHelpers::AngularCollocationsFor
Definition: CceComputationTestHelpers.hpp:40
Cce::TestHelpers::CopyDataBoxTags
Definition: CceComputationTestHelpers.hpp:60
ComplexModalVector
A class for storing complex spectral coefficients on a spectral grid.
Definition: ComplexModalVector.hpp:39
Variables.hpp
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Cce
The set of utilities for performing Cauchy characteristic evolution and Cauchy characteristic matchin...
Definition: BoundaryComputeAndSendToEvolution.hpp:28
Gsl.hpp
make_not_null
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion,...
Definition: Gsl.hpp:880
ComplexDataVector
Stores a collection of complex function values.
Definition: ComplexDataVector.hpp:47
std::data
T data(T... args)
outer_product
void outer_product(const gsl::not_null< ResultVectorType * > result, const LhsVectorType &lhs, const RhsVectorType &rhs) noexcept
Computes the outer product between two vectors.
Definition: VectorAlgebra.hpp:19
TMPL.hpp
Tags::Multiplies
A prefix tag representing the product of two other tags. Note that if non-spin-weighted types are nee...
Definition: Tags.hpp:86
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183
Spectral
Functionality associated with a particular choice of basis functions and quadrature for spectral oper...
Definition: Chebyshev.cpp:19
Spectral::Swsh::number_of_swsh_collocation_points
constexpr size_t number_of_swsh_collocation_points(const size_t l_max) noexcept
Convenience function for determining the number of spin-weighted spherical harmonic collocation value...
Definition: SwshCollocation.hpp:25