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