BoundaryTestHelpers.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <boost/iterator/zip_iterator.hpp>
7 #include <cstddef>
8 
9 #include "DataStructures/ComplexDataVector.hpp"
10 #include "DataStructures/DataVector.hpp"
13 #include "Evolution/Systems/Cce/BoundaryData.hpp"
14 #include "Evolution/Systems/Cce/ReadBoundaryDataH5.hpp"
15 #include "Helpers/Evolution/Systems/Cce/WriteToWorldtubeH5.hpp"
16 #include "NumericalAlgorithms/Spectral/SwshCollocation.hpp"
17 #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp"
20 #include "Utilities/FileSystem.hpp"
21 #include "Utilities/Gsl.hpp"
22 
23 namespace Cce {
24 namespace TestHelpers {
25 template <typename... Structure>
26 Tensor<ComplexModalVector, Structure...> tensor_to_goldberg_coefficients(
27  const Tensor<DataVector, Structure...>& nodal_data, size_t l_max) noexcept {
28  Tensor<ComplexModalVector, Structure...> goldberg_modal_data{
29  square(l_max + 1)};
30  SpinWeighted<ComplexDataVector, 0> transform_buffer{
32  for (size_t i = 0; i < nodal_data.size(); ++i) {
33  transform_buffer.data() = std::complex<double>(1.0, 0.0) * nodal_data[i];
34  goldberg_modal_data[i] =
36  Spectral::Swsh::swsh_transform(l_max, 1, transform_buffer), l_max)
37  .data();
38  }
39  return goldberg_modal_data;
40 }
41 
42 template <typename... Structure>
43 Tensor<ComplexModalVector, Structure...> tensor_to_libsharp_coefficients(
44  const Tensor<DataVector, Structure...>& nodal_data,
45  const size_t l_max) // NOLINT(readability-avoid-const-params-in-decls)
46  noexcept {
47  Tensor<ComplexModalVector, Structure...> libsharp_modal_data{
49  SpinWeighted<ComplexDataVector, 0> transform_buffer{
51  for (size_t i = 0; i < nodal_data.size(); ++i) {
52  transform_buffer.data() = std::complex<double>(1.0, 0.0) * nodal_data[i];
53  libsharp_modal_data[i] =
54  Spectral::Swsh::swsh_transform(l_max, 1, transform_buffer).data();
55  }
56  return libsharp_modal_data;
57 }
58 
59 template <typename AnalyticSolution>
60 void create_fake_time_varying_modal_data(
61  const gsl::not_null<tnsr::ii<ComplexModalVector, 3>*>
62  spatial_metric_coefficients,
63  const gsl::not_null<tnsr::ii<ComplexModalVector, 3>*>
64  dt_spatial_metric_coefficients,
65  const gsl::not_null<tnsr::ii<ComplexModalVector, 3>*>
66  dr_spatial_metric_coefficients,
67  const gsl::not_null<tnsr::I<ComplexModalVector, 3>*> shift_coefficients,
68  const gsl::not_null<tnsr::I<ComplexModalVector, 3>*> dt_shift_coefficients,
69  const gsl::not_null<tnsr::I<ComplexModalVector, 3>*> dr_shift_coefficients,
70  const gsl::not_null<Scalar<ComplexModalVector>*> lapse_coefficients,
71  const gsl::not_null<Scalar<ComplexModalVector>*> dt_lapse_coefficients,
72  const gsl::not_null<Scalar<ComplexModalVector>*> dr_lapse_coefficients,
73  const AnalyticSolution& solution, const double extraction_radius,
74  const double amplitude, const double frequency, const double time,
75  const size_t l_max, const bool convert_to_goldberg = true,
76  const bool apply_normalization_bug = false) noexcept {
77  const size_t number_of_angular_points =
79  // create the vector of collocation points that we want to interpolate to
80 
81  tnsr::I<DataVector, 3> collocation_points{number_of_angular_points};
82  const auto& collocation = Spectral::Swsh::cached_collocation_metadata<
83  Spectral::Swsh::ComplexRepresentation::Interleaved>(l_max);
84  for (const auto& collocation_point : collocation) {
85  get<0>(collocation_points)[collocation_point.offset] =
86  extraction_radius * (1.0 + amplitude * sin(frequency * time)) *
87  sin(collocation_point.theta) * cos(collocation_point.phi);
88  get<1>(collocation_points)[collocation_point.offset] =
89  extraction_radius * (1.0 + amplitude * sin(frequency * time)) *
90  sin(collocation_point.theta) * sin(collocation_point.phi);
91  get<2>(collocation_points)[collocation_point.offset] =
92  extraction_radius * (1.0 + amplitude * sin(frequency * time)) *
93  cos(collocation_point.theta);
94  }
95 
96  const auto kerr_schild_variables = solution.variables(
97  collocation_points, 0.0, gr::Solutions::KerrSchild::tags<DataVector>{});
98 
99  const Scalar<DataVector>& lapse =
100  get<gr::Tags::Lapse<DataVector>>(kerr_schild_variables);
101  const Scalar<DataVector>& dt_lapse =
102  get<::Tags::dt<gr::Tags::Lapse<DataVector>>>(kerr_schild_variables);
103  const auto& d_lapse = get<gr::Solutions::KerrSchild::DerivLapse<DataVector>>(
104  kerr_schild_variables);
105 
106  const auto& shift = get<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>(
107  kerr_schild_variables);
108  const auto& dt_shift =
109  get<::Tags::dt<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>>(
110  kerr_schild_variables);
111  const auto& d_shift = get<gr::Solutions::KerrSchild::DerivShift<DataVector>>(
112  kerr_schild_variables);
113 
114  const auto& spatial_metric =
115  get<gr::Tags::SpatialMetric<3, ::Frame::Inertial, DataVector>>(
116  kerr_schild_variables);
117  const auto& dt_spatial_metric = get<
119  kerr_schild_variables);
120  const auto& d_spatial_metric =
121  get<gr::Solutions::KerrSchild::DerivSpatialMetric<DataVector>>(
122  kerr_schild_variables);
123 
124  DataVector normalization_factor{number_of_angular_points, 1.0};
125  if (apply_normalization_bug) {
126  normalization_factor = 0.0;
127  const auto inverse_spatial_metric =
128  determinant_and_inverse(spatial_metric).second;
129  for (size_t i = 0; i < 3; ++i) {
130  for (size_t j = 0; j < 3; ++j) {
131  normalization_factor +=
132  inverse_spatial_metric.get(i, j) * collocation_points.get(i) *
133  collocation_points.get(j) /
134  square(extraction_radius *
135  (1.0 + amplitude * sin(frequency * time)));
136  }
137  }
138  normalization_factor = sqrt(normalization_factor);
139  }
140 
141  Scalar<DataVector> dr_lapse{number_of_angular_points};
142  get(dr_lapse) = (get<0>(collocation_points) * get<0>(d_lapse) +
143  get<1>(collocation_points) * get<1>(d_lapse) +
144  get<2>(collocation_points) * get<2>(d_lapse)) /
145  (extraction_radius * normalization_factor);
146  tnsr::I<DataVector, 3> dr_shift{number_of_angular_points};
147  for (size_t i = 0; i < 3; ++i) {
148  dr_shift.get(i) = (get<0>(collocation_points) * d_shift.get(0, i) +
149  get<1>(collocation_points) * d_shift.get(1, i) +
150  get<2>(collocation_points) * d_shift.get(2, i)) /
151  (extraction_radius * normalization_factor);
152  }
153  tnsr::ii<DataVector, 3> dr_spatial_metric{number_of_angular_points};
154  for (size_t i = 0; i < 3; ++i) {
155  for (size_t j = i; j < 3; ++j) {
156  dr_spatial_metric.get(i, j) =
157  (get<0>(collocation_points) * d_spatial_metric.get(0, i, j) +
158  get<1>(collocation_points) * d_spatial_metric.get(1, i, j) +
159  get<2>(collocation_points) * d_spatial_metric.get(2, i, j)) /
160  (extraction_radius * normalization_factor);
161  }
162  }
163 
164  if (convert_to_goldberg) {
165  *lapse_coefficients =
166  TestHelpers::tensor_to_goldberg_coefficients(lapse, l_max);
167  *dt_lapse_coefficients =
168  TestHelpers::tensor_to_goldberg_coefficients(dt_lapse, l_max);
169  *dr_lapse_coefficients =
170  TestHelpers::tensor_to_goldberg_coefficients(dr_lapse, l_max);
171 
172  *shift_coefficients =
173  TestHelpers::tensor_to_goldberg_coefficients(shift, l_max);
174  *dt_shift_coefficients =
175  TestHelpers::tensor_to_goldberg_coefficients(dt_shift, l_max);
176  *dr_shift_coefficients =
177  TestHelpers::tensor_to_goldberg_coefficients(dr_shift, l_max);
178 
179  *spatial_metric_coefficients =
180  TestHelpers::tensor_to_goldberg_coefficients(spatial_metric, l_max);
181  *dt_spatial_metric_coefficients =
182  TestHelpers::tensor_to_goldberg_coefficients(dt_spatial_metric, l_max);
183  *dr_spatial_metric_coefficients =
184  TestHelpers::tensor_to_goldberg_coefficients(dr_spatial_metric, l_max);
185  } else {
186  *lapse_coefficients =
187  TestHelpers::tensor_to_libsharp_coefficients(lapse, l_max);
188  *dt_lapse_coefficients =
189  TestHelpers::tensor_to_libsharp_coefficients(dt_lapse, l_max);
190  *dr_lapse_coefficients =
191  TestHelpers::tensor_to_libsharp_coefficients(dr_lapse, l_max);
192 
193  *shift_coefficients =
194  TestHelpers::tensor_to_libsharp_coefficients(shift, l_max);
195  *dt_shift_coefficients =
196  TestHelpers::tensor_to_libsharp_coefficients(dt_shift, l_max);
197  *dr_shift_coefficients =
198  TestHelpers::tensor_to_libsharp_coefficients(dr_shift, l_max);
199 
200  *spatial_metric_coefficients =
201  TestHelpers::tensor_to_libsharp_coefficients(spatial_metric, l_max);
202  *dt_spatial_metric_coefficients =
203  TestHelpers::tensor_to_libsharp_coefficients(dt_spatial_metric, l_max);
204  *dr_spatial_metric_coefficients =
205  TestHelpers::tensor_to_libsharp_coefficients(dr_spatial_metric, l_max);
206  }
207 }
208 
209 template <typename AnalyticSolution>
210 void write_test_file(const AnalyticSolution& solution,
211  const std::string& filename, const double target_time,
212  const double extraction_radius, const double frequency,
213  const double amplitude, const size_t l_max) noexcept {
214  const size_t goldberg_size = square(l_max + 1);
215  tnsr::ii<ComplexModalVector, 3> spatial_metric_coefficients{goldberg_size};
216  tnsr::ii<ComplexModalVector, 3> dt_spatial_metric_coefficients{goldberg_size};
217  tnsr::ii<ComplexModalVector, 3> dr_spatial_metric_coefficients{goldberg_size};
218  tnsr::I<ComplexModalVector, 3> shift_coefficients{goldberg_size};
219  tnsr::I<ComplexModalVector, 3> dt_shift_coefficients{goldberg_size};
220  tnsr::I<ComplexModalVector, 3> dr_shift_coefficients{goldberg_size};
221  Scalar<ComplexModalVector> lapse_coefficients{goldberg_size};
222  Scalar<ComplexModalVector> dt_lapse_coefficients{goldberg_size};
223  Scalar<ComplexModalVector> dr_lapse_coefficients{goldberg_size};
224 
225  // write times to file for several steps before and after the target time
226  if (file_system::check_if_file_exists(filename)) {
227  file_system::rm(filename, true);
228  }
229  // scoped to close the file
230  {
231  TestHelpers::WorldtubeModeRecorder recorder{filename, l_max};
232  for (size_t t = 0; t < 30; ++t) {
233  const double time = 0.1 * t + target_time - 1.5;
234  TestHelpers::create_fake_time_varying_modal_data(
235  make_not_null(&spatial_metric_coefficients),
236  make_not_null(&dt_spatial_metric_coefficients),
237  make_not_null(&dr_spatial_metric_coefficients),
238  make_not_null(&shift_coefficients),
239  make_not_null(&dt_shift_coefficients),
240  make_not_null(&dr_shift_coefficients),
241  make_not_null(&lapse_coefficients),
242  make_not_null(&dt_lapse_coefficients),
243  make_not_null(&dr_lapse_coefficients), solution, extraction_radius,
244  amplitude, frequency, time, l_max);
245  for (size_t i = 0; i < 3; ++i) {
246  for (size_t j = i; j < 3; ++j) {
247  recorder.append_worldtube_mode_data(
248  detail::dataset_name_for_component("/g", i, j), time,
249  spatial_metric_coefficients.get(i, j), l_max);
250  recorder.append_worldtube_mode_data(
251  detail::dataset_name_for_component("/Drg", i, j), time,
252  dr_spatial_metric_coefficients.get(i, j), l_max);
253  recorder.append_worldtube_mode_data(
254  detail::dataset_name_for_component("/Dtg", i, j), time,
255  dt_spatial_metric_coefficients.get(i, j), l_max);
256  }
257  recorder.append_worldtube_mode_data(
258  detail::dataset_name_for_component("/Shift", i), time,
259  shift_coefficients.get(i), l_max);
260  recorder.append_worldtube_mode_data(
261  detail::dataset_name_for_component("/DrShift", i), time,
262  dr_shift_coefficients.get(i), l_max);
263  recorder.append_worldtube_mode_data(
264  detail::dataset_name_for_component("/DtShift", i), time,
265  dt_shift_coefficients.get(i), l_max);
266  }
267  recorder.append_worldtube_mode_data(
268  detail::dataset_name_for_component("/Lapse"), time,
269  get(lapse_coefficients), l_max);
270  recorder.append_worldtube_mode_data(
271  detail::dataset_name_for_component("/DrLapse"), time,
272  get(dr_lapse_coefficients), l_max);
273  recorder.append_worldtube_mode_data(
274  detail::dataset_name_for_component("/DtLapse"), time,
275  get(dt_lapse_coefficients), l_max);
276  }
277  }
278 }
279 } // namespace TestHelpers
280 } // namespace Cce
constexpr size_t size_of_libsharp_coefficient_vector(const size_t l_max) noexcept
Convenience function for determining the number of spin-weighted spherical harmonics coefficients tha...
Definition: SwshCoefficients.hpp:34
bool check_if_file_exists(const std::string &file)
Returns true if the regular file or link to the regular file exists.
Definition: FileSystem.cpp:159
Definition: TestCreation.hpp:14
The set of utilities for performing Cauchy characteristic evolution and Cauchy characteristic matchin...
Definition: BoundaryComputeAndSendToEvolution.hpp:24
tnsr::ii< DataType, SpatialDim, Frame > spatial_metric(const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric) noexcept
Compute spatial metric from spacetime metric.
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
void libsharp_to_goldberg_modes(gsl::not_null< SpinWeighted< ComplexModalVector, Spin > *> goldberg_modes, const SpinWeighted< ComplexModalVector, Spin > &libsharp_modes, size_t l_max) noexcept
Compute the set of Goldberg Spin-weighted spherical harmonic modes (in the convention of ) from a lib...
const CollocationMetadata< Representation > & cached_collocation_metadata(const size_t l_max) noexcept
precomputation function for those collocation grids that are requested
Definition: SwshCollocation.cpp:66
Defines Functions for calculating spacetime tensors from 3+1 quantities.
Prefix indicating a time derivative.
Definition: Prefixes.hpp:32
T sin(T... args)
void rm(const std::string &path, bool recursive)
Deletes a file or directory.
Definition: FileSystem.cpp:234
SpinWeighted< ComplexModalVector, Spin > swsh_transform(const size_t l_max, const size_t number_of_radial_points, const SpinWeighted< ComplexDataVector, Spin > &collocation) noexcept
Perform a forward libsharp spin-weighted spherical harmonic transform on a single supplied SpinWeight...
Definition: SwshTransform.cpp:107
Scalar< DataType > lapse(const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric) noexcept
Compute lapse from shift and spacetime metric.
Make a spin-weighted type T with spin-weight Spin. Mathematical operators are restricted to addition...
Definition: SpinWeighted.hpp:24
auto determinant_and_inverse(const Tensor< T, Symm, tmpl::list< Index0, Index1 >> &tensor) noexcept -> std::pair< Scalar< T >, Tensor< T, Symm, tmpl::list< change_index_up_lo< Index1 >, change_index_up_lo< Index0 >>>>
Computes the determinant and inverse of a rank-2 Tensor.
Definition: DeterminantAndInverse.hpp:357
T cos(T... args)
A class for storing complex spectral coefficients on a spectral grid.
Definition: ComplexModalVector.hpp:39
Defines classes for Tensor.
Defines a list of useful type aliases for tensors.
const DataVector & collocation_points(const Mesh< 1 > &mesh) noexcept
Collocation points for a one-dimensional mesh.
tnsr::I< DataType, SpatialDim, Frame > shift(const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric, const tnsr::II< DataType, SpatialDim, Frame > &inverse_spatial_metric) noexcept
Compute shift from spacetime metric and inverse spatial metric.
Stores a collection of function values.
Definition: DataVector.hpp:42
Declares functions to do file system manipulations.
Defines functions and classes from the GSL.
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, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:879
T sqrt(T... args)
Defines Functions for calculating spacetime tensors from 3+1 quantities.
decltype(auto) constexpr square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:54
Tensor< T, Symmetry<>, index_list<> > Scalar
Scalar type.
Definition: TypeAliases.hpp:21
Require a pointer to not be a nullptr
Definition: Gsl.hpp:182