SpecBoundaryData.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"
10 #include "DataStructures/DataVector.hpp"
11 #include "DataStructures/SpinWeighted.hpp"
14 #include "Evolution/Systems/Cce/BoundaryData.hpp"
15 #include "Evolution/Systems/Cce/BoundaryDataTags.hpp"
16 #include "Evolution/Systems/Cce/Tags.hpp"
17 #include "NumericalAlgorithms/Spectral/SwshCollocation.hpp"
18 #include "NumericalAlgorithms/Spectral/SwshDerivatives.hpp"
21 #include "Utilities/Gsl.hpp"
22 
23 namespace Cce {
24 /*
25  * \brief Compute \f$\gamma_{i j}\f$, \f$\gamma^{i j}\f$,
26  * \f$\partial_i \gamma_{j k}\f$, and
27  * \f$\partial_t g_{i j}\f$ from input libsharp-compatible modal spatial
28  * metric quantities.
29  *
30  * \details This function will apply a correction factor associated with a SpEC
31  * bug.
32  */
33 void cartesian_spatial_metric_and_derivatives_from_unnormalized_spec_modes(
34  gsl::not_null<tnsr::ii<DataVector, 3>*> cartesian_spatial_metric,
35  gsl::not_null<tnsr::II<DataVector, 3>*> inverse_cartesian_spatial_metric,
36  gsl::not_null<tnsr::ijj<DataVector, 3>*> d_cartesian_spatial_metric,
37  gsl::not_null<tnsr::ii<DataVector, 3>*> dt_cartesian_spatial_metric,
39  interpolation_modal_buffer,
41  interpolation_buffer,
43  gsl::not_null<Scalar<DataVector>*> radial_correction_factor,
44  const tnsr::ii<ComplexModalVector, 3>& spatial_metric_coefficients,
45  const tnsr::ii<ComplexModalVector, 3>& dr_spatial_metric_coefficients,
46  const tnsr::ii<ComplexModalVector, 3>& dt_spatial_metric_coefficients,
47  const CartesianiSphericalJ& inverse_cartesian_to_spherical_jacobian,
48  const tnsr::I<DataVector, 3>& unit_cartesian_coords, size_t l_max) noexcept;
49 
50 /*!
51  * \brief Compute \f$\beta^{i}\f$, \f$\partial_i \beta^{j}\f$, and
52  * \f$\partial_t \beta^i\f$ from input libsharp-compatible modal spatial
53  * metric quantities.
54  *
55  * \details This function will apply a correction factor associated with a SpEC
56  * bug.
57  */
58 void cartesian_shift_and_derivatives_from_unnormalized_spec_modes(
59  gsl::not_null<tnsr::I<DataVector, 3>*> cartesian_shift,
60  gsl::not_null<tnsr::iJ<DataVector, 3>*> d_cartesian_shift,
61  gsl::not_null<tnsr::I<DataVector, 3>*> dt_cartesian_shift,
63  interpolation_modal_buffer,
65  interpolation_buffer,
67  const tnsr::I<ComplexModalVector, 3>& shift_coefficients,
68  const tnsr::I<ComplexModalVector, 3>& dr_shift_coefficients,
69  const tnsr::I<ComplexModalVector, 3>& dt_shift_coefficients,
70  const CartesianiSphericalJ& inverse_cartesian_to_spherical_jacobian,
71  const Scalar<DataVector>& radial_derivative_correction_factor,
72  size_t l_max) noexcept;
73 
74 /*!
75  * \brief Compute \f$\alpha\f$, \f$\partial_i \alpha\f$, and
76  * \f$\partial_t \beta^i\f$ from input libsharp-compatible modal spatial
77  * metric quantities.
78  *
79  * \details This function will apply a correction factor associated with a SpEC
80  * bug.
81  */
82 void cartesian_lapse_and_derivatives_from_unnormalized_spec_modes(
83  gsl::not_null<Scalar<DataVector>*> cartesian_lapse,
84  gsl::not_null<tnsr::i<DataVector, 3>*> d_cartesian_lapse,
85  gsl::not_null<Scalar<DataVector>*> dt_cartesian_lapse,
87  interpolation_modal_buffer,
89  interpolation_buffer,
91  const Scalar<ComplexModalVector>& lapse_coefficients,
92  const Scalar<ComplexModalVector>& dr_lapse_coefficients,
93  const Scalar<ComplexModalVector>& dt_lapse_coefficients,
94  const CartesianiSphericalJ& inverse_cartesian_to_spherical_jacobian,
95  const Scalar<DataVector>& radial_derivative_correction_factor,
96  size_t l_max) noexcept;
97 
98 /*!
99  * \brief Process the worldtube data from modal metric components and
100  * derivatives with incorrectly normalized radial derivatives from an old
101  * version of SpEC to desired Bondi quantities, placing the result in the passed
102  * \ref DataBoxGroup.
103  *
104  * \details
105  * The mathematics are a bit complicated for all of the coordinate
106  * transformations that are necessary to obtain the Bondi gauge quantities.
107  * For full mathematical details, see the documentation for functions in
108  * `BoundaryData.hpp` and \cite Barkett2019uae \cite Bishop1998uk.
109  *
110  * This function takes as input the full set of ADM metric data and its radial
111  * and time derivatives on a two-dimensional surface of constant \f$r\f$ and
112  * \f$t\f$ in numerical coordinates. This data must be provided as spherical
113  * harmonic coefficients in the libsharp format. This data is provided in nine
114  * `Tensor`s.
115  *
116  * Sufficient tags to provide full worldtube boundary data at a particular
117  * time are set in `bondi_boundary_data`. In particular, the set of tags in
118  * `Tags::characteristic_worldtube_boundary_tags` in the provided \ref
119  * DataBoxGroup are assigned to the worldtube boundary values associated with
120  * the input metric components.
121  *
122  * The majority of the mathematical transformations are implemented as a set of
123  * individual cascaded functions below. The details of the manipulations that
124  * are performed to the input data may be found in the individual functions
125  * themselves, which are called in the following order:
126  * - `trigonometric_functions_on_swsh_collocation()`
127  * - `cartesian_to_spherical_coordinates_and_jacobians()`
128  * - `cartesian_spatial_metric_and_derivatives_from_unnormalized_spec_modes()`
129  * - `cartesian_shift_and_derivatives_from_unnormalized_spec_modes()`
130  * - `cartesian_lapse_and_derivatives_from_unnormalized_spec_modes()`
131  * - `GeneralizedHarmonic::phi()`
132  * - `gr::time_derivative_of_spacetime_metric`
133  * - `gr::spacetime_metric`
134  * - `generalized_harmonic_quantities()`
135  * - `worldtube_normal_and_derivatives()`
136  * - `null_vector_l_and_derivatives()`
137  * - `null_metric_and_derivative()`
138  * - `dlambda_null_metric_and_inverse()`
139  * - `bondi_r()`
140  * - `d_bondi_r()`
141  * - `dyads()`
142  * - `beta_worldtube_data()`
143  * - `bondi_u_worldtube_data()`
144  * - `bondi_w_worldtube_data()`
145  * - `bondi_j_worldtube_data()`
146  * - `dr_bondi_j()`
147  * - `d2lambda_bondi_r()`
148  * - `bondi_q_worldtube_data()`
149  * - `bondi_h_worldtube_data()`
150  * - `du_j_worldtube_data()`
151  */
152 template <typename DataBoxTagList>
154  const gsl::not_null<db::DataBox<DataBoxTagList>*> bondi_boundary_data,
155  const tnsr::ii<ComplexModalVector, 3>& spatial_metric_coefficients,
156  const tnsr::ii<ComplexModalVector, 3>& dt_spatial_metric_coefficients,
157  const tnsr::ii<ComplexModalVector, 3>& dr_spatial_metric_coefficients,
158  const tnsr::I<ComplexModalVector, 3>& shift_coefficients,
159  const tnsr::I<ComplexModalVector, 3>& dt_shift_coefficients,
160  const tnsr::I<ComplexModalVector, 3>& dr_shift_coefficients,
161  const Scalar<ComplexModalVector>& lapse_coefficients,
162  const Scalar<ComplexModalVector>& dt_lapse_coefficients,
163  const Scalar<ComplexModalVector>& dr_lapse_coefficients,
164  const double extraction_radius, const size_t l_max) noexcept {
165  const size_t size = Spectral::Swsh::number_of_swsh_collocation_points(l_max);
166 
167  // Most allocations required for the full boundary computation are merged into
168  // a single, large Variables allocation. There remain a handful of cases in
169  // the computational functions called where an intermediate quantity that is
170  // not re-used is allocated rather than taking a buffer. These cases are
171  // marked with code comments 'Allocation'; In future, allocations are
172  // identified as a point to optimize, those buffers may be allocated here and
173  // passed as function arguments
174  Variables<tmpl::list<
175  Tags::detail::CosPhi, Tags::detail::CosTheta, Tags::detail::SinPhi,
176  Tags::detail::SinTheta, Tags::detail::CartesianCoordinates,
177  Tags::detail::CartesianToSphericalJacobian,
178  Tags::detail::InverseCartesianToSphericalJacobian,
181  ::Tags::deriv<gr::Tags::SpatialMetric<3, ::Frame::Inertial, DataVector>,
182  tmpl::size_t<3>, ::Frame::Inertial>,
185  ::Tags::deriv<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>,
186  tmpl::size_t<3>, ::Frame::Inertial>,
189  ::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
190  ::Frame::Inertial>,
195  Tags::detail::WorldtubeNormal, ::Tags::dt<Tags::detail::WorldtubeNormal>,
196  Tags::detail::NullL, ::Tags::dt<Tags::detail::NullL>,
197  // for the detail function called at the end
201  Tags::detail::AngularDNullL,
202  Tags::detail::DLambda<
203  gr::Tags::SpacetimeMetric<3, Frame::RadialNull, DataVector>>,
204  Tags::detail::DLambda<
205  gr::Tags::InverseSpacetimeMetric<3, Frame::RadialNull, DataVector>>,
206  ::Tags::spacetime_deriv<Tags::detail::RealBondiR, tmpl::size_t<3>,
208  Tags::detail::DLambda<Tags::detail::DLambda<Tags::detail::RealBondiR>>,
209  ::Tags::deriv<Tags::detail::DLambda<Tags::detail::RealBondiR>,
210  tmpl::size_t<2>, Frame::RadialNull>,
212  computation_variables{size};
213 
214  Variables<
215  tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
219  derivative_buffers{size};
220  auto& cos_phi = get<Tags::detail::CosPhi>(computation_variables);
221  auto& cos_theta = get<Tags::detail::CosTheta>(computation_variables);
222  auto& sin_phi = get<Tags::detail::SinPhi>(computation_variables);
223  auto& sin_theta = get<Tags::detail::SinTheta>(computation_variables);
224  trigonometric_functions_on_swsh_collocation(
225  make_not_null(&cos_phi), make_not_null(&cos_theta),
226  make_not_null(&sin_phi), make_not_null(&sin_theta), l_max);
227 
228  // NOTE: to handle the singular values of polar coordinates, the phi
229  // components of all tensors are scaled according to their sin(theta)
230  // prefactors.
231  // so, any down-index component get<2>(A) represents 1/sin(theta) A_\phi,
232  // and any up-index component get<2>(A) represents sin(theta) A^\phi.
233  // This holds for Jacobians, and so direct application of the Jacobians
234  // brings the factors through.
235  auto& cartesian_coords =
236  get<Tags::detail::CartesianCoordinates>(computation_variables);
237  auto& cartesian_to_spherical_jacobian =
238  get<Tags::detail::CartesianToSphericalJacobian>(computation_variables);
239  auto& inverse_cartesian_to_spherical_jacobian =
240  get<Tags::detail::InverseCartesianToSphericalJacobian>(
241  computation_variables);
242  cartesian_to_spherical_coordinates_and_jacobians(
243  make_not_null(&cartesian_coords),
244  make_not_null(&cartesian_to_spherical_jacobian),
245  make_not_null(&inverse_cartesian_to_spherical_jacobian), cos_phi,
246  cos_theta, sin_phi, sin_theta, extraction_radius);
247 
248  auto& cartesian_spatial_metric =
249  get<gr::Tags::SpatialMetric<3, ::Frame::Inertial, DataVector>>(
250  computation_variables);
251  auto& inverse_spatial_metric =
252  get<gr::Tags::InverseSpatialMetric<3, ::Frame::Inertial, DataVector>>(
253  computation_variables);
254  auto& d_cartesian_spatial_metric = get<
255  ::Tags::deriv<gr::Tags::SpatialMetric<3, ::Frame::Inertial, DataVector>,
256  tmpl::size_t<3>, ::Frame::Inertial>>(computation_variables);
257  auto& dt_cartesian_spatial_metric = get<
258  ::Tags::dt<gr::Tags::SpatialMetric<3, ::Frame::Inertial, DataVector>>>(
259  computation_variables);
260  auto& interpolation_buffer =
261  get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
262  std::integral_constant<int, 0>>>(
263  derivative_buffers);
264  Scalar<SpinWeighted<ComplexModalVector, 0>> interpolation_modal_buffer{size};
265  auto& eth_buffer =
266  get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
268  derivative_buffers);
269  auto& radial_correction_factor =
270  get<::Tags::TempScalar<0, DataVector>>(computation_variables);
271  cartesian_spatial_metric_and_derivatives_from_unnormalized_spec_modes(
272  make_not_null(&cartesian_spatial_metric),
273  make_not_null(&inverse_spatial_metric),
274  make_not_null(&d_cartesian_spatial_metric),
275  make_not_null(&dt_cartesian_spatial_metric),
276  make_not_null(&interpolation_modal_buffer),
277  make_not_null(&interpolation_buffer), make_not_null(&eth_buffer),
278  make_not_null(&radial_correction_factor), spatial_metric_coefficients,
279  dr_spatial_metric_coefficients, dt_spatial_metric_coefficients,
280  inverse_cartesian_to_spherical_jacobian, cartesian_coords, l_max);
281 
282  auto& cartesian_shift =
283  get<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>(
284  computation_variables);
285  auto& d_cartesian_shift =
286  get<::Tags::deriv<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>,
287  tmpl::size_t<3>, ::Frame::Inertial>>(
288  computation_variables);
289  auto& dt_cartesian_shift =
290  get<::Tags::dt<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>>(
291  computation_variables);
292 
293  cartesian_shift_and_derivatives_from_unnormalized_spec_modes(
294  make_not_null(&cartesian_shift), make_not_null(&d_cartesian_shift),
295  make_not_null(&dt_cartesian_shift),
296  make_not_null(&interpolation_modal_buffer),
297  make_not_null(&interpolation_buffer), make_not_null(&eth_buffer),
298  shift_coefficients, dr_shift_coefficients, dt_shift_coefficients,
299  inverse_cartesian_to_spherical_jacobian, radial_correction_factor, l_max);
300 
301  auto& cartesian_lapse =
302  get<gr::Tags::Lapse<DataVector>>(computation_variables);
303  auto& d_cartesian_lapse =
304  get<::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
305  ::Frame::Inertial>>(computation_variables);
306  auto& dt_cartesian_lapse =
307  get<::Tags::dt<gr::Tags::Lapse<DataVector>>>(computation_variables);
308  cartesian_lapse_and_derivatives_from_unnormalized_spec_modes(
309  make_not_null(&cartesian_lapse), make_not_null(&d_cartesian_lapse),
310  make_not_null(&dt_cartesian_lapse),
311  make_not_null(&interpolation_modal_buffer),
312  make_not_null(&interpolation_buffer), make_not_null(&eth_buffer),
313  lapse_coefficients, dr_lapse_coefficients, dt_lapse_coefficients,
314  inverse_cartesian_to_spherical_jacobian, radial_correction_factor, l_max);
315 
316  auto& phi = get<GeneralizedHarmonic::Tags::Phi<3, ::Frame::Inertial>>(
317  computation_variables);
318  auto& dt_spacetime_metric = get<
319  ::Tags::dt<gr::Tags::SpacetimeMetric<3, ::Frame::Inertial, DataVector>>>(
320  computation_variables);
321  auto& spacetime_metric =
322  get<gr::Tags::SpacetimeMetric<3, ::Frame::Inertial, DataVector>>(
323  computation_variables);
324  GeneralizedHarmonic::phi(
325  make_not_null(&phi), cartesian_lapse, d_cartesian_lapse, cartesian_shift,
326  d_cartesian_shift, cartesian_spatial_metric, d_cartesian_spatial_metric);
327  gr::time_derivative_of_spacetime_metric(
328  make_not_null(&dt_spacetime_metric), cartesian_lapse, dt_cartesian_lapse,
329  cartesian_shift, dt_cartesian_shift, cartesian_spatial_metric,
330  dt_cartesian_spatial_metric);
331  gr::spacetime_metric(make_not_null(&spacetime_metric), cartesian_lapse,
332  cartesian_shift, cartesian_spatial_metric);
333 
334  auto& dt_worldtube_normal =
335  get<::Tags::dt<Tags::detail::WorldtubeNormal>>(computation_variables);
336  auto& worldtube_normal =
337  get<Tags::detail::WorldtubeNormal>(computation_variables);
338  worldtube_normal_and_derivatives(
339  make_not_null(&worldtube_normal), make_not_null(&dt_worldtube_normal),
340  cos_phi, cos_theta, spacetime_metric, dt_spacetime_metric, sin_phi,
341  sin_theta, inverse_spatial_metric);
342 
343  auto& du_null_l = get<::Tags::dt<Tags::detail::NullL>>(computation_variables);
344  auto& null_l = get<Tags::detail::NullL>(computation_variables);
345  null_vector_l_and_derivatives(
346  make_not_null(&du_null_l), make_not_null(&null_l), dt_worldtube_normal,
347  dt_cartesian_lapse, dt_spacetime_metric, dt_cartesian_shift,
348  cartesian_lapse, spacetime_metric, cartesian_shift, worldtube_normal);
349 
350  // pass to the next step that is common between the 'modal' input and 'GH'
351  // input strategies
352  detail::create_bondi_boundary_data(
353  bondi_boundary_data, make_not_null(&computation_variables),
354  make_not_null(&derivative_buffers), dt_spacetime_metric, phi,
355  spacetime_metric, null_l, du_null_l, cartesian_to_spherical_jacobian,
356  l_max, extraction_radius);
357 }
358 } // namespace Cce
Auxiliary variable which is analytically the spatial derivative of the spacetime metric.
Definition: Tags.hpp:41
Defines function computing the determinant and inverse of a tensor.
Definition: Tags.hpp:17
Contains functionality for Cauchy Characteristic Extraction.
Definition: BoundaryData.cpp:21
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
Definition: Tags.hpp:59
Inverse of the spatial metric.
Definition: Tags.hpp:36
Defines Functions for calculating spacetime tensors from 3+1 quantities.
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.
Prefix indicating a time derivative.
Definition: Prefixes.hpp:30
Definition: Variables.hpp:909
Given a Tag with a type of Tensor<VectorType, ...>, acts as a new version of the tag with type of Ten...
Definition: Tags.hpp:58
Make a spin-weighted type T with spin-weight Spin. Mathematical operators are restricted to addition...
Definition: SpinWeighted.hpp:25
Defines classes and functions used for manipulating DataBox&#39;s.
Definition: Tags.hpp:28
The frame for the spherical metric in which the radial coordinate is an affine parameter along outwar...
Definition: BoundaryDataTags.hpp:12
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&#39;s equations...
Definition: ComputeGhQuantities.cpp:59
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
Definition: Tags.hpp:54
Defines a list of useful type aliases for tensors.
Definition: Tags.hpp:22
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
Definition: IndexType.hpp:44
void create_bondi_boundary_data_from_unnormalized_spec_modes(const gsl::not_null< db::DataBox< DataBoxTagList > *> bondi_boundary_data, const tnsr::ii< ComplexModalVector, 3 > &spatial_metric_coefficients, const tnsr::ii< ComplexModalVector, 3 > &dt_spatial_metric_coefficients, const tnsr::ii< ComplexModalVector, 3 > &dr_spatial_metric_coefficients, const tnsr::I< ComplexModalVector, 3 > &shift_coefficients, const tnsr::I< ComplexModalVector, 3 > &dt_shift_coefficients, const tnsr::I< ComplexModalVector, 3 > &dr_shift_coefficients, const Scalar< ComplexModalVector > &lapse_coefficients, const Scalar< ComplexModalVector > &dt_lapse_coefficients, const Scalar< ComplexModalVector > &dr_lapse_coefficients, const double extraction_radius, const size_t l_max) noexcept
Process the worldtube data from modal metric components and derivatives with incorrectly normalized r...
Definition: SpecBoundaryData.hpp:153
Defines Functions for calculating spacetime tensors from 3+1 quantities.
Tensor< T, Symmetry<>, index_list<> > Scalar
Scalar type.
Definition: TypeAliases.hpp:21
Require a pointer to not be a nullptr
Definition: Gsl.hpp:182