AnalyticDataHelpers.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
7 
8 #include <complex>
9 #include <cstddef>
10 
11 #include "DataStructures/ComplexDataVector.hpp"
12 #include "DataStructures/DataVector.hpp"
13 #include "DataStructures/SpinWeighted.hpp"
15 #include "Evolution/Systems/Cce/AnalyticSolutions/SphericalMetricData.hpp"
16 #include "Evolution/Systems/Cce/Tags.hpp"
17 #include "Framework/Pypp.hpp"
18 #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/SpatialDerivOfLapse.hpp"
19 #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/SpatialDerivOfShift.hpp"
20 #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/TimeDerivOfLapse.hpp"
21 #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/TimeDerivOfShift.hpp"
22 #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/TimeDerivOfSpatialMetric.hpp"
23 #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/TimeDerivativeOfSpacetimeMetric.hpp"
24 #include "PointwiseFunctions/GeneralRelativity/InverseSpacetimeMetric.hpp"
25 #include "PointwiseFunctions/GeneralRelativity/Lapse.hpp"
26 #include "PointwiseFunctions/GeneralRelativity/Shift.hpp"
27 #include "PointwiseFunctions/GeneralRelativity/SpacetimeNormalVector.hpp"
28 #include "PointwiseFunctions/GeneralRelativity/SpatialMetric.hpp"
29 #include "Utilities/TaggedTuple.hpp"
30 
31 namespace Cce {
32 namespace Solutions {
33 namespace TestHelpers {
34 
35 template <typename SphericalSolution>
36 struct SphericalSolutionWrapper : public SphericalSolution {
37  using taglist =
38  tmpl::list<gr::Tags::SpacetimeMetric<
44  Tags::News>;
45  using SphericalSolution::SphericalSolution;
46 
47  template <typename... Args>
48  void test_spherical_metric(const std::string python_file, const size_t l_max,
49  const double time,
50  const Args... args) const noexcept {
51  const size_t size =
53  Scalar<DataVector> sin_theta{size};
54  Scalar<DataVector> cos_theta{size};
55  const auto& collocation = Spectral::Swsh::cached_collocation_metadata<
56  Spectral::Swsh::ComplexRepresentation::Interleaved>(l_max);
57  for (const auto& collocation_point : collocation) {
58  get(sin_theta)[collocation_point.offset] = sin(collocation_point.theta);
59  get(cos_theta)[collocation_point.offset] = cos(collocation_point.theta);
60  }
61  tnsr::aa<DataVector, 3, ::Frame::Spherical<::Frame::Inertial>>
62  local_spherical_metric{size};
63  tnsr::aa<DataVector, 3, ::Frame::Spherical<::Frame::Inertial>>
64  local_dr_spherical_metric{size};
65  tnsr::aa<DataVector, 3, ::Frame::Spherical<::Frame::Inertial>>
66  local_dt_spherical_metric{size};
67  Scalar<SpinWeighted<ComplexDataVector, -2>> local_news{size};
68 
69  this->spherical_metric(make_not_null(&local_spherical_metric), l_max, time);
70  this->dr_spherical_metric(make_not_null(&local_dr_spherical_metric), l_max,
71  time);
72  this->dt_spherical_metric(make_not_null(&local_dt_spherical_metric), l_max,
73  time);
74  this->variables_impl(make_not_null(&local_news), l_max, time,
75  tmpl::type_<Tags::News>{});
76 
77  // Pypp call expects all of the objects to be the same category -- here we
78  // need to use tensors, so we must pack up the `double` arguments into
79  // tensors.
80  Scalar<DataVector> time_vector;
81  get(time_vector) = DataVector{size, time};
82 
83  const auto py_spherical_metric = pypp::call<
84  tnsr::aa<DataVector, 3, ::Frame::Spherical<::Frame::Inertial>>>(
85  python_file, "spherical_metric", sin_theta, cos_theta, time_vector,
86  Scalar<DataVector>{DataVector{size, args}}...);
87  const auto py_dt_spherical_metric = pypp::call<
88  tnsr::aa<DataVector, 3, ::Frame::Spherical<::Frame::Inertial>>>(
89  python_file, "dt_spherical_metric", sin_theta, cos_theta, time_vector,
90  Scalar<DataVector>{DataVector{size, args}}...);
91  const auto py_dr_spherical_metric = pypp::call<
92  tnsr::aa<DataVector, 3, ::Frame::Spherical<::Frame::Inertial>>>(
93  python_file, "dr_spherical_metric", sin_theta, cos_theta, time_vector,
94  Scalar<DataVector>{DataVector{size, args}}...);
95  const auto py_news = pypp::call<Scalar<SpinWeighted<ComplexDataVector, 2>>>(
96  python_file, "news", sin_theta, time_vector,
97  Scalar<DataVector>{DataVector{size, args}}...);
98 
99  for (size_t a = 0; a < 4; ++a) {
100  for (size_t b = a; b < 4; ++b) {
101  CAPTURE(a);
102  CAPTURE(b);
103  const auto& lhs = local_spherical_metric.get(a, b);
104  const auto& rhs = py_spherical_metric.get(a, b);
105  CHECK_ITERABLE_APPROX(lhs, rhs);
106  const auto& dr_lhs = local_dr_spherical_metric.get(a, b);
107  const auto& dr_rhs = py_dr_spherical_metric.get(a, b);
108  CHECK_ITERABLE_APPROX(dr_lhs, dr_rhs);
109  const auto& dt_lhs = local_dt_spherical_metric.get(a, b);
110  const auto& dt_rhs = py_dt_spherical_metric.get(a, b);
111  CHECK_ITERABLE_APPROX(dt_lhs, dt_rhs);
112  }
113  }
114  }
115 
116  protected:
117  using SphericalSolution::extraction_radius_;
118 };
119 
120 // This function determines the Bondi-Sachs scalars from a Cartesian spacetime
121 // metric, assuming that the metric is already in null form, so the spatial
122 // coordinates are related to standard Bondi-Sachs coordinates by just the
123 // standard Cartesian to spherical Jacobian.
125 extract_bondi_scalars_from_cartesian_metric(
126  const tnsr::aa<DataVector, 3>& spacetime_metric,
127  const CartesianiSphericalJ& inverse_jacobian,
128  double extraction_radius) noexcept;
129 
130 // This function determines the time derivative of the Bondi-Sachs scalars
131 // from the time derivative of a Cartesian spacetime metric, the Cartesian
132 // metric, and Jacobian factors. This procedure assumes that the metric is
133 // already in null form, so the spatial coordinates are related to standard
134 // Bondi-Sachs coordinates by just the standard cartesian to spherical Jacobian.
137 extract_dt_bondi_scalars_from_cartesian_metric(
138  const tnsr::aa<DataVector, 3>& dt_spacetime_metric,
139  const tnsr::aa<DataVector, 3>& spacetime_metric,
140  const CartesianiSphericalJ& inverse_jacobian,
141  double extraction_radius) noexcept;
142 
143 // This function determines the radial derivative of the Bondi-Sachs scalars
144 // from the radial derivative of a Cartesian spacetime metric, the Cartesian
145 // metric, and Jacobian factors. This procedure assumes that the metric is
146 // already in null form, so the spatial coordinates are related to standard
147 // Bondi-Sachs coordinates by just the standard cartesian to spherical Jacobian.
150 extract_dr_bondi_scalars_from_cartesian_metric(
151  const tnsr::aa<DataVector, 3>& dr_spacetime_metric,
152  const tnsr::aa<DataVector, 3>& spacetime_metric,
153  const CartesianiSphericalJ& inverse_jacobian,
154  const CartesianiSphericalJ& dr_inverse_jacobian,
155  double extraction_radius) noexcept;
156 
157 // This function checks the consistency of the 3+1 ADM quantities in the
158 // `boundary_tuple` with quantities computed from `expected_spacetime_metric`,
159 // `expected_dt_spacetime_metric`, and `expected_d_spacetime_metric`. If the
160 // expected quantities are also extracted from the tuple, this simply checks
161 // that the boundary computation has produced a consistent set of quantities
162 // and has not generated NaNs or other pathological values (e.g. a degenerate
163 // spacetime metric) in the process.
164 template <typename... TupleTags>
165 void check_adm_metric_quantities(
166  const tuples::TaggedTuple<TupleTags...>& boundary_tuple,
167  const tnsr::aa<DataVector, 3>& expected_spacetime_metric,
168  const tnsr::aa<DataVector, 3>& expected_dt_spacetime_metric,
169  const tnsr::iaa<DataVector, 3>& expected_d_spacetime_metric) noexcept {
170  // check the 3+1 quantities are computed correctly in the abstract base class
171  // `WorldtubeData`
172  const auto& dr_cartesian_coordinates =
173  get<Tags::Dr<Tags::CauchyCartesianCoords>>(boundary_tuple);
174 
175  const auto& lapse = get<gr::Tags::Lapse<DataVector>>(boundary_tuple);
176  const auto& shift =
177  get<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>(boundary_tuple);
178  const auto& spatial_metric =
179  get<gr::Tags::SpatialMetric<3, ::Frame::Inertial, DataVector>>(
180  boundary_tuple);
181 
182  const auto expected_spatial_metric =
183  gr::spatial_metric(expected_spacetime_metric);
184  const auto expected_inverse_spatial_metric =
185  determinant_and_inverse(expected_spatial_metric).second;
186  const auto expected_shift =
187  gr::shift(expected_spacetime_metric, expected_inverse_spatial_metric);
188  const auto expected_lapse =
189  gr::lapse(expected_shift, expected_spacetime_metric);
190  CHECK_ITERABLE_APPROX(spatial_metric, expected_spatial_metric);
191  CHECK_ITERABLE_APPROX(shift, expected_shift);
192  CHECK_ITERABLE_APPROX(lapse, expected_lapse);
193 
194  const auto& pi =
195  get<GeneralizedHarmonic::Tags::Pi<3, ::Frame::Inertial>>(boundary_tuple);
196  const auto dt_spacetime_metric_from_pi =
198  expected_lapse, expected_shift, pi, expected_d_spacetime_metric);
199  CHECK_ITERABLE_APPROX(expected_dt_spacetime_metric,
200  dt_spacetime_metric_from_pi);
201  // Check that the time derivative values are consistent with the Generalized
202  // Harmonic `pi` -- these are redundant if the boundary calculation uses `pi`
203  // to derive these, but it is not required to do so.
204  const auto& dt_lapse =
205  get<::Tags::dt<gr::Tags::Lapse<DataVector>>>(boundary_tuple);
206  const auto& dt_shift =
207  get<::Tags::dt<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>>(
208  boundary_tuple);
209  const auto& dt_spatial_metric = get<
211  boundary_tuple);
212  const auto expected_dt_spatial_metric =
214  expected_lapse, expected_shift, expected_d_spacetime_metric, pi);
215  const auto expected_spacetime_unit_normal =
216  gr::spacetime_normal_vector(expected_lapse, expected_shift);
217  const auto expected_dt_lapse = GeneralizedHarmonic::time_deriv_of_lapse(
218  expected_lapse, expected_shift, expected_spacetime_unit_normal,
219  expected_d_spacetime_metric, pi);
220  const auto expected_dt_shift = GeneralizedHarmonic::time_deriv_of_shift(
221  expected_lapse, expected_shift, expected_inverse_spatial_metric,
222  expected_spacetime_unit_normal, expected_d_spacetime_metric, pi);
223  CHECK_ITERABLE_APPROX(dt_lapse, expected_dt_lapse);
224  CHECK_ITERABLE_APPROX(dt_shift, expected_dt_shift);
225  CHECK_ITERABLE_APPROX(dt_spatial_metric, expected_dt_spatial_metric);
226 
227  const auto& dr_lapse =
228  get<Tags::Dr<gr::Tags::Lapse<DataVector>>>(boundary_tuple);
229  const auto& dr_shift =
230  get<Tags::Dr<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>>(
231  boundary_tuple);
232  const auto& dr_spatial_metric =
233  get<Tags::Dr<gr::Tags::SpatialMetric<3, ::Frame::Inertial, DataVector>>>(
234  boundary_tuple);
235  const auto expected_spatial_derivative_of_lapse =
237  expected_lapse, expected_spacetime_unit_normal,
238  expected_d_spacetime_metric);
239  const auto expected_inverse_spacetime_metric = gr::inverse_spacetime_metric(
240  expected_lapse, expected_shift, expected_inverse_spatial_metric);
241  const auto expected_spatial_derivative_of_shift =
243  expected_lapse, expected_inverse_spacetime_metric,
244  expected_spacetime_unit_normal, expected_d_spacetime_metric);
245  DataVector expected_buffer =
246  get<0>(dr_cartesian_coordinates) *
247  get<0>(expected_spatial_derivative_of_lapse) +
248  get<1>(dr_cartesian_coordinates) *
249  get<1>(expected_spatial_derivative_of_lapse) +
250  get<2>(dr_cartesian_coordinates) *
251  get<2>(expected_spatial_derivative_of_lapse);
252  CHECK_ITERABLE_APPROX(expected_buffer, get(dr_lapse));
253  for (size_t i = 0; i < 3; ++i) {
254  expected_buffer = get<0>(dr_cartesian_coordinates) *
255  expected_spatial_derivative_of_shift.get(0, i) +
256  get<1>(dr_cartesian_coordinates) *
257  expected_spatial_derivative_of_shift.get(1, i) +
258  get<2>(dr_cartesian_coordinates) *
259  expected_spatial_derivative_of_shift.get(2, i);
260  CHECK_ITERABLE_APPROX(expected_buffer, dr_shift.get(i));
261  for (size_t j = i; j < 3; ++j) {
262  expected_buffer = get<0>(dr_cartesian_coordinates) *
263  expected_d_spacetime_metric.get(0, i + 1, j + 1) +
264  get<1>(dr_cartesian_coordinates) *
265  expected_d_spacetime_metric.get(1, i + 1, j + 1) +
266  get<2>(dr_cartesian_coordinates) *
267  expected_d_spacetime_metric.get(2, i + 1, j + 1);
268  CHECK_ITERABLE_APPROX(expected_buffer, dr_spatial_metric.get(i, j));
269  }
270  }
271 }
272 
273 } // namespace TestHelpers
274 } // namespace Solutions
275 } // 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
GeneralizedHarmonic::spatial_deriv_of_shift
void spatial_deriv_of_shift(gsl::not_null< tnsr::iJ< DataType, SpatialDim, Frame > * > deriv_shift, const Scalar< DataType > &lapse, const tnsr::AA< DataType, SpatialDim, Frame > &inverse_spacetime_metric, const tnsr::A< DataType, SpatialDim, Frame > &spacetime_unit_normal, const tnsr::iaa< DataType, SpatialDim, Frame > &phi) noexcept
Computes spatial derivatives of the shift vector from the generalized harmonic and geometric variable...
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:639
TestingFramework.hpp
GeneralizedHarmonic::time_deriv_of_spatial_metric
void time_deriv_of_spatial_metric(gsl::not_null< tnsr::ii< DataType, SpatialDim, Frame > * > dt_spatial_metric, const Scalar< DataType > &lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::iaa< DataType, SpatialDim, Frame > &phi, const tnsr::aa< DataType, SpatialDim, Frame > &pi) noexcept
Computes time derivative of the spatial metric.
pypp::call
R call(const std::string &module_name, const std::string &function_name, const Args &... t)
Calls a Python function from a module/file with given parameters.
Definition: Pypp.hpp:336
CHECK_ITERABLE_APPROX
#define CHECK_ITERABLE_APPROX(a, b)
A wrapper around Catch's CHECK macro that checks approximate equality of entries in iterable containe...
Definition: TestingFramework.hpp:139
GeneralizedHarmonic::spatial_deriv_of_lapse
void spatial_deriv_of_lapse(gsl::not_null< tnsr::i< DataType, SpatialDim, Frame > * > deriv_lapse, const Scalar< DataType > &lapse, const tnsr::A< DataType, SpatialDim, Frame > &spacetime_unit_normal, const tnsr::iaa< DataType, SpatialDim, Frame > &phi) noexcept
Computes spatial derivatives of lapse (N) from the generalized harmonic variables and spacetime unit ...
SpinWeighted
Make a spin-weighted type T with spin-weight Spin. Mathematical operators are restricted to addition,...
Definition: SpinWeighted.hpp:24
GeneralizedHarmonic::time_deriv_of_shift
void time_deriv_of_shift(gsl::not_null< tnsr::I< DataType, SpatialDim, Frame > * > dt_shift, const Scalar< DataType > &lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::II< DataType, SpatialDim, Frame > &inverse_spatial_metric, const tnsr::A< DataType, SpatialDim, Frame > &spacetime_unit_normal, const tnsr::iaa< DataType, SpatialDim, Frame > &phi, const tnsr::aa< DataType, SpatialDim, Frame > &pi) noexcept
Computes time derivative of the shift vector from the generalized harmonic and geometric variables.
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
Pypp.hpp
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
gr::Tags::SpacetimeMetric
Definition: Tags.hpp:17
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:42
GeneralizedHarmonic::time_derivative_of_spacetime_metric
void time_derivative_of_spacetime_metric(gsl::not_null< tnsr::aa< DataType, SpatialDim, Frame > * > dt_spacetime_metric, const Scalar< DataType > &lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::aa< DataType, SpatialDim, Frame > &pi, const tnsr::iaa< DataType, SpatialDim, Frame > &phi) noexcept
Computes the time derivative of the spacetime metric from the generalized harmonic quantities ,...
Cce::Solutions::TestHelpers::SphericalSolutionWrapper
Definition: AnalyticDataHelpers.hpp:36
Cce::Tags::News
Definition: Tags.hpp:151
Frame::Spherical
Represents a spherical-coordinate frame that is associated with a Cartesian frame,...
Definition: IndexType.hpp:54
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:29
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
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
GeneralizedHarmonic::time_deriv_of_lapse
void time_deriv_of_lapse(gsl::not_null< Scalar< DataType > * > dt_lapse, const Scalar< DataType > &lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::A< DataType, SpatialDim, Frame > &spacetime_unit_normal, const tnsr::iaa< DataType, SpatialDim, Frame > &phi, const tnsr::aa< DataType, SpatialDim, Frame > &pi) noexcept
Computes time derivative of lapse (N) from the generalized harmonic variables, lapse,...
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.
complex
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
gr::inverse_spacetime_metric
void inverse_spacetime_metric(gsl::not_null< tnsr::AA< DataType, SpatialDim, Frame > * > inverse_spacetime_metric, const Scalar< DataType > &lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::II< DataType, SpatialDim, Frame > &inverse_spatial_metric) noexcept
Compute inverse spacetime metric from inverse spatial metric, lapse and shift.
gr::spacetime_normal_vector
tnsr::A< DataType, SpatialDim, Frame > spacetime_normal_vector(const Scalar< DataType > &lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift) noexcept
Computes spacetime normal vector from lapse and shift.
Cce::Tags::Dr
The derivative with respect to Bondi .
Definition: Tags.hpp:124
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