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"
18 #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/Phi.hpp"
19 #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/Pi.hpp"
20 #include "PointwiseFunctions/GeneralRelativity/SpacetimeMetric.hpp"
21 #include "Utilities/FileSystem.hpp"
22 #include "Utilities/Gsl.hpp"
23 
24 namespace Cce {
25 namespace TestHelpers {
26 template <typename... Structure>
27 Tensor<ComplexModalVector, Structure...> tensor_to_goldberg_coefficients(
28  const Tensor<DataVector, Structure...>& nodal_data, size_t l_max) noexcept {
29  Tensor<ComplexModalVector, Structure...> goldberg_modal_data{
30  square(l_max + 1)};
31  SpinWeighted<ComplexDataVector, 0> transform_buffer{
33  for (size_t i = 0; i < nodal_data.size(); ++i) {
34  transform_buffer.data() = std::complex<double>(1.0, 0.0) * nodal_data[i];
35  goldberg_modal_data[i] =
37  Spectral::Swsh::swsh_transform(l_max, 1, transform_buffer), l_max)
38  .data();
39  }
40  return goldberg_modal_data;
41 }
42 
43 template <typename... Structure>
44 Tensor<ComplexModalVector, Structure...> tensor_to_libsharp_coefficients(
45  const Tensor<DataVector, Structure...>& nodal_data,
46  const size_t l_max) // NOLINT(readability-avoid-const-params-in-decls)
47  noexcept {
48  Tensor<ComplexModalVector, Structure...> libsharp_modal_data{
50  SpinWeighted<ComplexDataVector, 0> transform_buffer{
52  for (size_t i = 0; i < nodal_data.size(); ++i) {
53  transform_buffer.data() = std::complex<double>(1.0, 0.0) * nodal_data[i];
54  libsharp_modal_data[i] =
55  Spectral::Swsh::swsh_transform(l_max, 1, transform_buffer).data();
56  }
57  return libsharp_modal_data;
58 }
59 
60 template <typename AnalyticSolution>
61 void create_fake_time_varying_gh_nodal_data(
62  const gsl::not_null<tnsr::aa<DataVector, 3>*> spacetime_metric,
63  const gsl::not_null<tnsr::iaa<DataVector, 3>*> phi,
64  const gsl::not_null<tnsr::aa<DataVector, 3>*> pi,
65  const AnalyticSolution& solution, const double extraction_radius,
66  const double amplitude, const double frequency, const double time,
67  const size_t l_max) noexcept {
68  const size_t number_of_angular_points =
70  // create the vector of collocation points that we want to interpolate to
71 
72  tnsr::I<DataVector, 3> collocation_points{number_of_angular_points};
73  const auto& collocation = Spectral::Swsh::cached_collocation_metadata<
74  Spectral::Swsh::ComplexRepresentation::Interleaved>(l_max);
75  for (const auto& collocation_point : collocation) {
76  get<0>(collocation_points)[collocation_point.offset] =
77  extraction_radius * (1.0 + amplitude * sin(frequency * time)) *
78  sin(collocation_point.theta) * cos(collocation_point.phi);
79  get<1>(collocation_points)[collocation_point.offset] =
80  extraction_radius * (1.0 + amplitude * sin(frequency * time)) *
81  sin(collocation_point.theta) * sin(collocation_point.phi);
82  get<2>(collocation_points)[collocation_point.offset] =
83  extraction_radius * (1.0 + amplitude * sin(frequency * time)) *
84  cos(collocation_point.theta);
85  }
86 
87  const auto kerr_schild_variables = solution.variables(
88  collocation_points, 0.0, gr::Solutions::KerrSchild::tags<DataVector>{});
89 
90  const Scalar<DataVector>& lapse =
91  get<gr::Tags::Lapse<DataVector>>(kerr_schild_variables);
92  const Scalar<DataVector>& dt_lapse =
93  get<::Tags::dt<gr::Tags::Lapse<DataVector>>>(kerr_schild_variables);
94  const auto& d_lapse = get<gr::Solutions::KerrSchild::DerivLapse<DataVector>>(
95  kerr_schild_variables);
96 
97  const auto& shift = get<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>(
98  kerr_schild_variables);
99  const auto& dt_shift =
100  get<::Tags::dt<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>>(
101  kerr_schild_variables);
102  const auto& d_shift = get<gr::Solutions::KerrSchild::DerivShift<DataVector>>(
103  kerr_schild_variables);
104 
105  const auto& spatial_metric =
106  get<gr::Tags::SpatialMetric<3, ::Frame::Inertial, DataVector>>(
107  kerr_schild_variables);
108  const auto& dt_spatial_metric = get<
110  kerr_schild_variables);
111  const auto& d_spatial_metric =
112  get<gr::Solutions::KerrSchild::DerivSpatialMetric<DataVector>>(
113  kerr_schild_variables);
114 
117  d_spatial_metric);
118  GeneralizedHarmonic::pi(pi, lapse, dt_lapse, shift, dt_shift, spatial_metric,
119  dt_spatial_metric, *phi);
120 }
121 
122 template <typename AnalyticSolution>
123 void create_fake_time_varying_modal_data(
124  const gsl::not_null<tnsr::ii<ComplexModalVector, 3>*>
125  spatial_metric_coefficients,
126  const gsl::not_null<tnsr::ii<ComplexModalVector, 3>*>
127  dt_spatial_metric_coefficients,
128  const gsl::not_null<tnsr::ii<ComplexModalVector, 3>*>
129  dr_spatial_metric_coefficients,
130  const gsl::not_null<tnsr::I<ComplexModalVector, 3>*> shift_coefficients,
131  const gsl::not_null<tnsr::I<ComplexModalVector, 3>*> dt_shift_coefficients,
132  const gsl::not_null<tnsr::I<ComplexModalVector, 3>*> dr_shift_coefficients,
133  const gsl::not_null<Scalar<ComplexModalVector>*> lapse_coefficients,
134  const gsl::not_null<Scalar<ComplexModalVector>*> dt_lapse_coefficients,
135  const gsl::not_null<Scalar<ComplexModalVector>*> dr_lapse_coefficients,
136  const AnalyticSolution& solution, const double extraction_radius,
137  const double amplitude, const double frequency, const double time,
138  const size_t l_max, const bool convert_to_goldberg = true,
139  const bool apply_normalization_bug = false) noexcept {
140  const size_t number_of_angular_points =
142  // create the vector of collocation points that we want to interpolate to
143 
144  tnsr::I<DataVector, 3> collocation_points{number_of_angular_points};
145  const auto& collocation = Spectral::Swsh::cached_collocation_metadata<
146  Spectral::Swsh::ComplexRepresentation::Interleaved>(l_max);
147  for (const auto& collocation_point : collocation) {
148  get<0>(collocation_points)[collocation_point.offset] =
149  extraction_radius * (1.0 + amplitude * sin(frequency * time)) *
150  sin(collocation_point.theta) * cos(collocation_point.phi);
151  get<1>(collocation_points)[collocation_point.offset] =
152  extraction_radius * (1.0 + amplitude * sin(frequency * time)) *
153  sin(collocation_point.theta) * sin(collocation_point.phi);
154  get<2>(collocation_points)[collocation_point.offset] =
155  extraction_radius * (1.0 + amplitude * sin(frequency * time)) *
156  cos(collocation_point.theta);
157  }
158 
159  const auto kerr_schild_variables = solution.variables(
160  collocation_points, 0.0, gr::Solutions::KerrSchild::tags<DataVector>{});
161 
162  const Scalar<DataVector>& lapse =
163  get<gr::Tags::Lapse<DataVector>>(kerr_schild_variables);
164  const Scalar<DataVector>& dt_lapse =
165  get<::Tags::dt<gr::Tags::Lapse<DataVector>>>(kerr_schild_variables);
166  const auto& d_lapse = get<gr::Solutions::KerrSchild::DerivLapse<DataVector>>(
167  kerr_schild_variables);
168 
169  const auto& shift = get<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>(
170  kerr_schild_variables);
171  const auto& dt_shift =
172  get<::Tags::dt<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>>(
173  kerr_schild_variables);
174  const auto& d_shift = get<gr::Solutions::KerrSchild::DerivShift<DataVector>>(
175  kerr_schild_variables);
176 
177  const auto& spatial_metric =
178  get<gr::Tags::SpatialMetric<3, ::Frame::Inertial, DataVector>>(
179  kerr_schild_variables);
180  const auto& dt_spatial_metric = get<
182  kerr_schild_variables);
183  const auto& d_spatial_metric =
184  get<gr::Solutions::KerrSchild::DerivSpatialMetric<DataVector>>(
185  kerr_schild_variables);
186 
187  DataVector normalization_factor{number_of_angular_points, 1.0};
188  if (apply_normalization_bug) {
189  normalization_factor = 0.0;
190  const auto inverse_spatial_metric =
192  for (size_t i = 0; i < 3; ++i) {
193  for (size_t j = 0; j < 3; ++j) {
194  normalization_factor +=
195  inverse_spatial_metric.get(i, j) * collocation_points.get(i) *
196  collocation_points.get(j) /
197  square(extraction_radius *
198  (1.0 + amplitude * sin(frequency * time)));
199  }
200  }
201  normalization_factor = sqrt(normalization_factor);
202  }
203 
204  Scalar<DataVector> dr_lapse{number_of_angular_points};
205  get(dr_lapse) = (get<0>(collocation_points) * get<0>(d_lapse) +
206  get<1>(collocation_points) * get<1>(d_lapse) +
207  get<2>(collocation_points) * get<2>(d_lapse)) /
208  (extraction_radius * normalization_factor);
209  tnsr::I<DataVector, 3> dr_shift{number_of_angular_points};
210  for (size_t i = 0; i < 3; ++i) {
211  dr_shift.get(i) = (get<0>(collocation_points) * d_shift.get(0, i) +
212  get<1>(collocation_points) * d_shift.get(1, i) +
213  get<2>(collocation_points) * d_shift.get(2, i)) /
214  (extraction_radius * normalization_factor);
215  }
216  tnsr::ii<DataVector, 3> dr_spatial_metric{number_of_angular_points};
217  for (size_t i = 0; i < 3; ++i) {
218  for (size_t j = i; j < 3; ++j) {
219  dr_spatial_metric.get(i, j) =
220  (get<0>(collocation_points) * d_spatial_metric.get(0, i, j) +
221  get<1>(collocation_points) * d_spatial_metric.get(1, i, j) +
222  get<2>(collocation_points) * d_spatial_metric.get(2, i, j)) /
223  (extraction_radius * normalization_factor);
224  }
225  }
226 
227  if (convert_to_goldberg) {
228  *lapse_coefficients =
229  TestHelpers::tensor_to_goldberg_coefficients(lapse, l_max);
230  *dt_lapse_coefficients =
231  TestHelpers::tensor_to_goldberg_coefficients(dt_lapse, l_max);
232  *dr_lapse_coefficients =
233  TestHelpers::tensor_to_goldberg_coefficients(dr_lapse, l_max);
234 
235  *shift_coefficients =
236  TestHelpers::tensor_to_goldberg_coefficients(shift, l_max);
237  *dt_shift_coefficients =
238  TestHelpers::tensor_to_goldberg_coefficients(dt_shift, l_max);
239  *dr_shift_coefficients =
240  TestHelpers::tensor_to_goldberg_coefficients(dr_shift, l_max);
241 
242  *spatial_metric_coefficients =
243  TestHelpers::tensor_to_goldberg_coefficients(spatial_metric, l_max);
244  *dt_spatial_metric_coefficients =
245  TestHelpers::tensor_to_goldberg_coefficients(dt_spatial_metric, l_max);
246  *dr_spatial_metric_coefficients =
247  TestHelpers::tensor_to_goldberg_coefficients(dr_spatial_metric, l_max);
248  } else {
249  *lapse_coefficients =
250  TestHelpers::tensor_to_libsharp_coefficients(lapse, l_max);
251  *dt_lapse_coefficients =
252  TestHelpers::tensor_to_libsharp_coefficients(dt_lapse, l_max);
253  *dr_lapse_coefficients =
254  TestHelpers::tensor_to_libsharp_coefficients(dr_lapse, l_max);
255 
256  *shift_coefficients =
257  TestHelpers::tensor_to_libsharp_coefficients(shift, l_max);
258  *dt_shift_coefficients =
259  TestHelpers::tensor_to_libsharp_coefficients(dt_shift, l_max);
260  *dr_shift_coefficients =
261  TestHelpers::tensor_to_libsharp_coefficients(dr_shift, l_max);
262 
263  *spatial_metric_coefficients =
264  TestHelpers::tensor_to_libsharp_coefficients(spatial_metric, l_max);
265  *dt_spatial_metric_coefficients =
266  TestHelpers::tensor_to_libsharp_coefficients(dt_spatial_metric, l_max);
267  *dr_spatial_metric_coefficients =
268  TestHelpers::tensor_to_libsharp_coefficients(dr_spatial_metric, l_max);
269  }
270 }
271 
272 template <typename AnalyticSolution>
273 void write_test_file(const AnalyticSolution& solution,
274  const std::string& filename, const double target_time,
275  const double extraction_radius, const double frequency,
276  const double amplitude, const size_t l_max) noexcept {
277  const size_t goldberg_size = square(l_max + 1);
278  tnsr::ii<ComplexModalVector, 3> spatial_metric_coefficients{goldberg_size};
279  tnsr::ii<ComplexModalVector, 3> dt_spatial_metric_coefficients{goldberg_size};
280  tnsr::ii<ComplexModalVector, 3> dr_spatial_metric_coefficients{goldberg_size};
281  tnsr::I<ComplexModalVector, 3> shift_coefficients{goldberg_size};
282  tnsr::I<ComplexModalVector, 3> dt_shift_coefficients{goldberg_size};
283  tnsr::I<ComplexModalVector, 3> dr_shift_coefficients{goldberg_size};
284  Scalar<ComplexModalVector> lapse_coefficients{goldberg_size};
285  Scalar<ComplexModalVector> dt_lapse_coefficients{goldberg_size};
286  Scalar<ComplexModalVector> dr_lapse_coefficients{goldberg_size};
287 
288  // write times to file for several steps before and after the target time
289  if (file_system::check_if_file_exists(filename)) {
290  file_system::rm(filename, true);
291  }
292  // scoped to close the file
293  {
294  TestHelpers::WorldtubeModeRecorder recorder{filename, l_max};
295  for (size_t t = 0; t < 30; ++t) {
296  const double time = 0.1 * t + target_time - 1.5;
297  TestHelpers::create_fake_time_varying_modal_data(
298  make_not_null(&spatial_metric_coefficients),
299  make_not_null(&dt_spatial_metric_coefficients),
300  make_not_null(&dr_spatial_metric_coefficients),
301  make_not_null(&shift_coefficients),
302  make_not_null(&dt_shift_coefficients),
303  make_not_null(&dr_shift_coefficients),
304  make_not_null(&lapse_coefficients),
305  make_not_null(&dt_lapse_coefficients),
306  make_not_null(&dr_lapse_coefficients), solution, extraction_radius,
307  amplitude, frequency, time, l_max);
308  for (size_t i = 0; i < 3; ++i) {
309  for (size_t j = i; j < 3; ++j) {
310  recorder.append_worldtube_mode_data(
311  detail::dataset_name_for_component("/g", i, j), time,
312  spatial_metric_coefficients.get(i, j), l_max);
313  recorder.append_worldtube_mode_data(
314  detail::dataset_name_for_component("/Drg", i, j), time,
315  dr_spatial_metric_coefficients.get(i, j), l_max);
316  recorder.append_worldtube_mode_data(
317  detail::dataset_name_for_component("/Dtg", i, j), time,
318  dt_spatial_metric_coefficients.get(i, j), l_max);
319  }
320  recorder.append_worldtube_mode_data(
321  detail::dataset_name_for_component("/Shift", i), time,
322  shift_coefficients.get(i), l_max);
323  recorder.append_worldtube_mode_data(
324  detail::dataset_name_for_component("/DrShift", i), time,
325  dr_shift_coefficients.get(i), l_max);
326  recorder.append_worldtube_mode_data(
327  detail::dataset_name_for_component("/DtShift", i), time,
328  dt_shift_coefficients.get(i), l_max);
329  }
330  recorder.append_worldtube_mode_data(
331  detail::dataset_name_for_component("/Lapse"), time,
332  get(lapse_coefficients), l_max);
333  recorder.append_worldtube_mode_data(
334  detail::dataset_name_for_component("/DrLapse"), time,
335  get(dr_lapse_coefficients), l_max);
336  recorder.append_worldtube_mode_data(
337  detail::dataset_name_for_component("/DtLapse"), time,
338  get(dt_lapse_coefficients), l_max);
339  }
340  }
341 }
342 } // namespace TestHelpers
343 } // namespace Cce
GeneralizedHarmonic::pi
void pi(gsl::not_null< tnsr::aa< DataType, SpatialDim, Frame > * > pi, const Scalar< DataType > &lapse, const Scalar< DataType > &dt_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::I< DataType, SpatialDim, Frame > &dt_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ii< DataType, SpatialDim, Frame > &dt_spatial_metric, const tnsr::iaa< DataType, SpatialDim, Frame > &phi) noexcept
Computes the conjugate momentum of the spacetime metric .
std::string
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:689
std::cos
T cos(T... args)
square
constexpr decltype(auto) square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:55
Spectral::Swsh::swsh_transform
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
Spectral::Swsh::size_of_libsharp_coefficient_vector
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
SpinWeighted
Make a spin-weighted type T with spin-weight Spin. Mathematical operators are restricted to addition,...
Definition: SpinWeighted.hpp:24
FileSystem.hpp
std::sqrt
T sqrt(T... args)
file_system::check_if_file_exists
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
gr::lapse
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.
determinant_and_inverse
void determinant_and_inverse(const gsl::not_null< Scalar< T > * > det, const gsl::not_null< Tensor< T, Symm, tmpl::list< change_index_up_lo< Index1 >, change_index_up_lo< Index0 >>> * > inv, const Tensor< T, Symm, tmpl::list< Index0, Index1 >> &tensor) noexcept
Computes the determinant and inverse of a rank-2 Tensor.
Definition: DeterminantAndInverse.hpp:371
cstddef
Spectral::Swsh::cached_collocation_metadata
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
Spectral::Swsh::libsharp_to_goldberg_modes
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...
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:42
GeneralizedHarmonic::phi
void phi(gsl::not_null< tnsr::iaa< DataType, SpatialDim, Frame > * > phi, const Scalar< DataType > &lapse, const tnsr::i< DataType, SpatialDim, Frame > &deriv_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::iJ< DataType, SpatialDim, Frame > &deriv_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ijj< DataType, SpatialDim, Frame > &deriv_spatial_metric) noexcept
Computes the auxiliary variable used by the generalized harmonic formulation of Einstein's equations...
std::sin
T sin(T... args)
Spectral::collocation_points
const DataVector & collocation_points(const Mesh< 1 > &mesh) noexcept
Collocation points for a one-dimensional mesh.
ComplexModalVector
A class for storing complex spectral coefficients on a spectral grid.
Definition: ComplexModalVector.hpp:39
file_system::rm
void rm(const std::string &path, bool recursive)
Deletes a file or directory.
Definition: FileSystem.cpp:234
gr::shift
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.
Tags::dt
Prefix indicating a time derivative.
Definition: Prefixes.hpp:32
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
TypeAliases.hpp
gr::spacetime_metric
void spacetime_metric(gsl::not_null< tnsr::aa< DataType, Dim, Frame > * > spacetime_metric, const Scalar< DataType > &lapse, const tnsr::I< DataType, Dim, Frame > &shift, const tnsr::ii< DataType, Dim, Frame > &spatial_metric) noexcept
Computes the spacetime metric from the spatial metric, lapse, and shift.
Tensor.hpp
gr::spatial_metric
tnsr::ii< DataType, SpatialDim, Frame > spatial_metric(const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric) noexcept
Compute spatial metric from spacetime metric.
std::complex< double >
std::time
T time(T... args)
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
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183
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