Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 :
8 : #include "DataStructures/DataBox/DataBox.hpp"
9 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
10 : #include "DataStructures/DataBox/Prefixes.hpp"
11 : #include "DataStructures/DataVector.hpp"
12 : #include "DataStructures/SpinWeighted.hpp"
13 : #include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 : #include "Evolution/Systems/Cce/BoundaryDataTags.hpp"
16 : #include "Evolution/Systems/Cce/Tags.hpp"
17 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshCollocation.hpp"
18 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshDerivatives.hpp"
19 : #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/Phi.hpp"
20 : #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/TimeDerivOfLapse.hpp"
21 : #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/TimeDerivOfShift.hpp"
22 : #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/TimeDerivativeOfSpacetimeMetric.hpp"
23 : #include "PointwiseFunctions/GeneralRelativity/Lapse.hpp"
24 : #include "PointwiseFunctions/GeneralRelativity/Shift.hpp"
25 : #include "PointwiseFunctions/GeneralRelativity/SpacetimeMetric.hpp"
26 : #include "PointwiseFunctions/GeneralRelativity/SpacetimeNormalVector.hpp"
27 : #include "PointwiseFunctions/GeneralRelativity/SpatialMetric.hpp"
28 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
29 : #include "PointwiseFunctions/GeneralRelativity/TimeDerivativeOfSpacetimeMetric.hpp"
30 : #include "Utilities/Gsl.hpp"
31 :
32 : /// \cond
33 : class DataVector;
34 : class ComplexDataVector;
35 : /// \endcond
36 :
37 : namespace Cce {
38 :
39 : /*!
40 : * \brief Constructs the collocation values for \f$\cos(\phi)\f$,
41 : * \f$\cos(\theta)\f$, \f$\sin(\phi)\f$, and \f$\sin(\theta)\f$, returned by
42 : * `not_null` pointer in that order.
43 : *
44 : * \details These are needed for coordinate transformations from the input
45 : * Cartesian-like coordinates.
46 : */
47 1 : void trigonometric_functions_on_swsh_collocation(
48 : gsl::not_null<Scalar<DataVector>*> cos_phi,
49 : gsl::not_null<Scalar<DataVector>*> cos_theta,
50 : gsl::not_null<Scalar<DataVector>*> sin_phi,
51 : gsl::not_null<Scalar<DataVector>*> sin_theta, size_t l_max);
52 :
53 : /*!
54 : * \brief Creates both the Jacobian and inverse Jacobian between Cartesian and
55 : * spherical coordinates, and the coordinates themselves
56 : *
57 : * \details The `cartesian_to_spherical_jacobian` is
58 : * \f$dx^i/d\tilde{x}^{\tilde j}\f$,
59 : * where the Cartesian components are in order \f$x^i = \{x, y, z\}\f$
60 : * and the spherical coordinates are
61 : * \f$\tilde{x}^{\tilde j} = \{r, \theta, \phi\}\f$.
62 : * The Cartesian coordinates given are the standard unit sphere coordinates:
63 : *
64 : * \f{align*}{
65 : * x &= \cos(\phi) \sin(\theta)\\
66 : * y &= \sin(\phi) \sin(\theta)\\
67 : * z &= \cos(\theta)
68 : * \f}
69 : *
70 : * \note These Jacobians are adjusted to improve regularity near the pole, in
71 : * particular the \f$\partial \phi / \partial x^i\f$ components have been scaled
72 : * by \f$\sin \theta\f$ (omitting a \f$1/\sin(\theta)\f$) and the
73 : * \f$\partial x^i/\partial \phi\f$ components have been scaled by
74 : * \f$1/\sin(\theta)\f$ (omitting a \f$\sin(\theta)\f$). The reason is that in
75 : * most careful calculations, these problematic sin factors can actually be
76 : * omitted because they cancel. In cases where they are actually required, they
77 : * must be put in by hand.
78 : */
79 1 : void cartesian_to_spherical_coordinates_and_jacobians(
80 : gsl::not_null<tnsr::I<DataVector, 3>*> unit_cartesian_coords,
81 : gsl::not_null<SphericaliCartesianJ*> cartesian_to_spherical_jacobian,
82 : gsl::not_null<CartesianiSphericalJ*>
83 : inverse_cartesian_to_spherical_jacobian,
84 : const Scalar<DataVector>& cos_phi, const Scalar<DataVector>& cos_theta,
85 : const Scalar<DataVector>& sin_phi, const Scalar<DataVector>& sin_theta,
86 : double extraction_radius);
87 :
88 : /*
89 : * \brief Compute \f$g_{i j}\f$, \f$g^{i j}\f$, \f$\partial_i g_{j k}\f$, and
90 : * \f$\partial_t g_{i j}\f$ from input libsharp-compatible modal spatial
91 : * metric quantities.
92 : *
93 : * \details This function interpolates the modes of
94 : * input \f$g_{ij}\f$, \f$\partial_r g_{i j}\f$, and \f$\partial_r g_{i j}\f$ to
95 : * the libsharp-compatible grid. This function then applies the necessary
96 : * jacobian factors and angular derivatives to determine the full \f$\partial_i
97 : * g_{j k}\f$.
98 : */
99 0 : void cartesian_spatial_metric_and_derivatives_from_modes(
100 : gsl::not_null<tnsr::ii<DataVector, 3>*> cartesian_spatial_metric,
101 : gsl::not_null<tnsr::II<DataVector, 3>*> inverse_cartesian_spatial_metric,
102 : gsl::not_null<tnsr::ijj<DataVector, 3>*> d_cartesian_spatial_metric,
103 : gsl::not_null<tnsr::ii<DataVector, 3>*> dt_cartesian_spatial_metric,
104 : gsl::not_null<Scalar<SpinWeighted<ComplexModalVector, 0>>*>
105 : interpolation_modal_buffer,
106 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
107 : interpolation_buffer,
108 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> eth_buffer,
109 : const tnsr::ii<ComplexModalVector, 3>& spatial_metric_coefficients,
110 : const tnsr::ii<ComplexModalVector, 3>& dr_spatial_metric_coefficients,
111 : const tnsr::ii<ComplexModalVector, 3>& dt_spatial_metric_coefficients,
112 : const CartesianiSphericalJ& inverse_cartesian_to_spherical_jacobian,
113 : size_t l_max);
114 :
115 : /*!
116 : * \brief Compute \f$\beta^{i}\f$, \f$\partial_i \beta^{j}\f$, and
117 : * \f$\partial_t \beta^i\f$ from input libsharp-compatible modal spatial
118 : * metric quantities.
119 : *
120 : * \details This function interpolates the modes of
121 : * input \f$\beta^i\f$, \f$\partial_r \beta^i\f$, and \f$\partial_r \beta^i\f$
122 : * to the libsharp-compatible grid. This function then applies the necessary
123 : * jacobian factors and angular derivatives to determine the full \f$\partial_i
124 : * \beta^i\f$.
125 : */
126 1 : void cartesian_shift_and_derivatives_from_modes(
127 : gsl::not_null<tnsr::I<DataVector, 3>*> cartesian_shift,
128 : gsl::not_null<tnsr::iJ<DataVector, 3>*> d_cartesian_shift,
129 : gsl::not_null<tnsr::I<DataVector, 3>*> dt_cartesian_shift,
130 : gsl::not_null<Scalar<SpinWeighted<ComplexModalVector, 0>>*>
131 : interpolation_modal_buffer,
132 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
133 : interpolation_buffer,
134 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> eth_buffer,
135 : const tnsr::I<ComplexModalVector, 3>& shift_coefficients,
136 : const tnsr::I<ComplexModalVector, 3>& dr_shift_coefficients,
137 : const tnsr::I<ComplexModalVector, 3>& dt_shift_coefficients,
138 : const CartesianiSphericalJ& inverse_cartesian_to_spherical_jacobian,
139 : size_t l_max);
140 :
141 : /*!
142 : * \brief Compute \f$\alpha\f$, \f$\partial_i \alpha\f$, and
143 : * \f$\partial_t \beta^i\f$ from input libsharp-compatible modal spatial
144 : * metric quantities.
145 : *
146 : * \details This function interpolates the modes of input \f$\alpha\f$,
147 : * \f$\partial_r \alpha\f$, and \f$\partial_r \alpha\f$ to the
148 : * libsharp-compatible grid. This function then applies the necessary jacobian
149 : * factors and angular derivatives to determine the full \f$\partial_i
150 : * \alpha\f$.
151 : */
152 1 : void cartesian_lapse_and_derivatives_from_modes(
153 : gsl::not_null<Scalar<DataVector>*> cartesian_lapse,
154 : gsl::not_null<tnsr::i<DataVector, 3>*> d_cartesian_lapse,
155 : gsl::not_null<Scalar<DataVector>*> dt_cartesian_lapse,
156 : gsl::not_null<Scalar<SpinWeighted<ComplexModalVector, 0>>*>
157 : interpolation_modal_buffer,
158 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
159 : interpolation_buffer,
160 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> eth_buffer,
161 : const Scalar<ComplexModalVector>& lapse_coefficients,
162 : const Scalar<ComplexModalVector>& dr_lapse_coefficients,
163 : const Scalar<ComplexModalVector>& dt_lapse_coefficients,
164 : const CartesianiSphericalJ& inverse_cartesian_to_spherical_jacobian,
165 : size_t l_max);
166 :
167 : /*!
168 : * \brief Computes spatial derivatives of cartesian metric, shift, and lapse
169 : * from nodal metric quantities on a spherical worldtube.
170 : *
171 : * \details See the details for
172 : * `cartesian_spatial_metric_and_derivatives_from_modes`,
173 : * `cartesian_shift_and_derivatives_from_modes`, and
174 : * `cartesian_lapse_and_derivatives_from_modes` ignoring the transformation from
175 : * modal to nodal (since the metric quantities are already in nodal form).
176 : */
177 1 : void deriv_cartesian_metric_lapse_shift_from_nodes(
178 : gsl::not_null<tnsr::ijj<DataVector, 3>*> d_cartesian_spatial_metric,
179 : gsl::not_null<tnsr::iJ<DataVector, 3>*> d_cartesian_shift,
180 : gsl::not_null<tnsr::i<DataVector, 3>*> d_cartesian_lapse,
181 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*> buffer,
182 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> eth_buffer,
183 : const tnsr::ii<DataVector, 3>& cartesian_spatial_metric,
184 : const tnsr::ii<DataVector, 3>& dr_cartesian_spatial_metric,
185 : const tnsr::I<DataVector, 3>& cartesian_shift,
186 : const tnsr::I<DataVector, 3>& dr_cartesian_shift,
187 : const Scalar<DataVector>& cartesian_lapse,
188 : const Scalar<DataVector>& dr_cartesian_lapse,
189 : const CartesianiSphericalJ& inverse_cartesian_to_spherical_jacobian,
190 : size_t l_max);
191 :
192 : /*!
193 : * \brief Computes the spacetime metric and its first derivative in the
194 : * intermediate radial null coordinates
195 : *
196 : * \details These components are obtained by the steps in
197 : * Section II-A of \cite Barkett2019uae, which is based on the computation from
198 : * Section 4.3 of \cite Bishop1998uk. The most direct comparison is to be made
199 : * with equation (31) of \cite Barkett2019uae, which gives the null metric
200 : * components explicitly. The time derivative is then (using notation from
201 : * equation (31) of \cite Barkett2019uae):
202 : *
203 : * \f{align}{
204 : * \partial_{\bar u} g_{\bar u \bar \lambda} =
205 : * \partial_{\bar u} g_{\bar \lambda \bar \lambda} =
206 : * \partial_{\bar u} g_{\bar \lambda \bar A} &= 0 \\
207 : * \partial_{\bar u} g_{\bar u \bar u} &=
208 : * \partial_{\breve t} g_{\breve t \breve t} \\
209 : * \partial_{\bar u} g_{\bar u \bar A} &=
210 : * \frac{\partial \breve x^{\breve i}}{\partial \bar x^{\bar A}}\\
211 : * g_{\breve i \breve t}
212 : * \partial_{\bar u} g_{\bar A \bar B}
213 : * &= \frac{\partial \breve x^{\breve i}}{\partial \bar x^{\bar A}}
214 : * \frac{\partial \breve x^{\breve j}}{\partial \bar x^{\bar B}}
215 : * g_{\breve i \breve j}
216 : * \f}
217 : */
218 1 : void null_metric_and_derivative(
219 : gsl::not_null<tnsr::aa<DataVector, 3, Frame::RadialNull>*> du_null_metric,
220 : gsl::not_null<tnsr::aa<DataVector, 3, Frame::RadialNull>*> null_metric,
221 : const SphericaliCartesianJ& cartesian_to_spherical_jacobian,
222 : const tnsr::aa<DataVector, 3>& dt_spacetime_metric,
223 : const tnsr::aa<DataVector, 3>& spacetime_metric);
224 :
225 : /*!
226 : * \brief Computes the spatial unit normal vector \f$s^i\f$ to the spherical
227 : * worldtube surface and its first time derivative.
228 : *
229 : * \details Refer to equation (20) of \cite Barkett2019uae for the expression of
230 : * the spatial unit normal vector, and equation (23) of \cite Barkett2019uae for
231 : * the first time derivative. Refer to \cite Bishop1998uk for more exposition
232 : * about the overall construction of the coordinate transformations used for the
233 : * intermediate null coordinates.
234 : */
235 1 : void worldtube_normal_and_derivatives(
236 : gsl::not_null<tnsr::I<DataVector, 3>*> worldtube_normal,
237 : gsl::not_null<tnsr::I<DataVector, 3>*> dt_worldtube_normal,
238 : const Scalar<DataVector>& cos_phi, const Scalar<DataVector>& cos_theta,
239 : const tnsr::aa<DataVector, 3>& spacetime_metric,
240 : const tnsr::aa<DataVector, 3>& dt_spacetime_metric,
241 : const Scalar<DataVector>& sin_phi, const Scalar<DataVector>& sin_theta,
242 : const tnsr::II<DataVector, 3>& inverse_spatial_metric);
243 :
244 : /*!
245 : * \brief Computes the null 4-vector \f$l^\mu\f$ on the worldtube surface that
246 : * is to be used as the CCE hypersurface generator, and the first time
247 : * derivative \f$\partial_u l^\mu\f$.
248 : *
249 : * \details For mathematical description of our choice of the null generator,
250 : * refer to equation (22) of \cite Barkett2019uae, and for the first time
251 : * derivative see (25) of \cite Barkett2019uae. Refer to \cite Bishop1998uk for
252 : * more exposition about the overall construction of the coordinate
253 : * transformations used for the intermediate null coordinates.
254 : */
255 1 : void null_vector_l_and_derivatives(
256 : gsl::not_null<tnsr::A<DataVector, 3>*> du_null_l,
257 : gsl::not_null<tnsr::A<DataVector, 3>*> null_l,
258 : const tnsr::I<DataVector, 3>& dt_worldtube_normal,
259 : const Scalar<DataVector>& dt_lapse,
260 : const tnsr::aa<DataVector, 3>& dt_spacetime_metric,
261 : const tnsr::I<DataVector, 3>& dt_shift, const Scalar<DataVector>& lapse,
262 : const tnsr::aa<DataVector, 3>& spacetime_metric,
263 : const tnsr::I<DataVector, 3>& shift,
264 : const tnsr::I<DataVector, 3>& worldtube_normal);
265 :
266 : /*!
267 : * \brief Computes the partial derivative of the spacetime metric and inverse
268 : * spacetime metric in the intermediate null radial coordinates with respect to
269 : * the null generator \f$l^\mu\f$
270 : *
271 : * \details For full expressions of the \f$l^\mu \partial_\mu g_{a b}\f$ and
272 : * \f$l^\mu \partial_\mu g^{a b}\f$ computed in this function, see equation (31)
273 : * and (32) of \cite Barkett2019uae. Refer to \cite Bishop1998uk for more
274 : * exposition about the overall construction of the coordinate transformations
275 : * used for the intermediate null coordinates.
276 : */
277 1 : void dlambda_null_metric_and_inverse(
278 : gsl::not_null<tnsr::aa<DataVector, 3, Frame::RadialNull>*>
279 : dlambda_null_metric,
280 : gsl::not_null<tnsr::AA<DataVector, 3, Frame::RadialNull>*>
281 : dlambda_inverse_null_metric,
282 : const AngulariCartesianA& angular_d_null_l,
283 : const SphericaliCartesianJ& cartesian_to_spherical_jacobian,
284 : const tnsr::iaa<DataVector, 3>& phi,
285 : const tnsr::aa<DataVector, 3>& dt_spacetime_metric,
286 : const tnsr::A<DataVector, 3>& du_null_l,
287 : const tnsr::AA<DataVector, 3, Frame::RadialNull>& inverse_null_metric,
288 : const tnsr::A<DataVector, 3>& null_l,
289 : const tnsr::aa<DataVector, 3>& spacetime_metric);
290 :
291 : /*!
292 : * \brief Computes the Bondi radius of the worldtube.
293 : *
294 : * \details Note that unlike the Cauchy coordinate radius, the Bondi radius is
295 : * not constant over the worldtube. Instead, it is obtained by the determinant
296 : * of the angular part of the metric in the intermediate null coordinates (see
297 : * \cite Barkett2019uae).
298 : *
299 : * \f[
300 : * r = \left(\frac{\det g_{A B}}{ q_{A B}}\right)^{1/4},
301 : * \f]
302 : *
303 : * where \f$q_{A B}\f$ is the unit sphere metric.
304 : */
305 1 : void bondi_r(gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*> bondi_r,
306 : const tnsr::aa<DataVector, 3, Frame::RadialNull>& null_metric);
307 :
308 : /*!
309 : * \brief Computes the full 4-dimensional partial of the Bondi radius with
310 : * respect to the intermediate null coordinates.
311 : *
312 : * \details The expression evaluated is obtained from differentiating the
313 : * determinant equation for `bondi_r`, from (35) of \cite Barkett2019uae :
314 : *
315 : * \f[
316 : * \partial_\alpha r = \frac{r}{4} \left(g^{A B} \partial_\alpha g_{A B}
317 : * - \frac{\partial_\alpha \det q_{A B}}{\det q_{A B}}\right)
318 : * \f]
319 : *
320 : * Note that for the angular derivatives, we just numerically differentiate
321 : * using the utilities in `Spectral::Swsh::angular_derivative()`. For the time
322 : * and radial derivatives, the second term in the above equation vanishes.
323 : */
324 1 : void d_bondi_r(
325 : gsl::not_null<tnsr::a<DataVector, 3, Frame::RadialNull>*> d_bondi_r,
326 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
327 : const tnsr::aa<DataVector, 3, Frame::RadialNull>& dlambda_null_metric,
328 : const tnsr::aa<DataVector, 3, Frame::RadialNull>& du_null_metric,
329 : const tnsr::AA<DataVector, 3, Frame::RadialNull>& inverse_null_metric,
330 : size_t l_max);
331 :
332 : /*!
333 : * \brief Compute the complex angular dyads used to define the spin-weighted
334 : * scalars in the CCE system.
335 : *
336 : * \details We use the typically chosen angular dyads in CCE
337 : * \cite Barkett2019uae \cite Bishop1997ik :
338 : *
339 : * \f{align*}{
340 : * q_A &= \{-1, -i \sin(\theta)\}\\
341 : * q^A &= \left\{-1, -i \frac{1}{\sin \theta}\right\}
342 : * \f}
343 : *
344 : * However, to maintain regularity and for compatibility with the more regular
345 : * Jacobians from `Cce::cartesian_to_spherical_coordinates_and_jacobians()`, in
346 : * the code we omit the factors of \f$\sin \theta\f$ from the above equations.
347 : */
348 1 : void dyads(
349 : gsl::not_null<tnsr::i<ComplexDataVector, 2, Frame::RadialNull>*> down_dyad,
350 : gsl::not_null<tnsr::I<ComplexDataVector, 2, Frame::RadialNull>*> up_dyad);
351 :
352 : /*!
353 : * \brief Compute the \f$\beta\f$ (lapse) function for the CCE Bondi-like
354 : * metric.
355 : *
356 : * \details The Bondi-like metric has \f$g^{u r} = - e^{2 \beta}\f$, and the
357 : * value of \f$\beta\f$ is obtained from the intermediate null metric by (see
358 : * equation (51) of \cite Barkett2019uae) using:
359 : *
360 : * \f[
361 : * \beta = -\frac{1}{2} \ln \partial_{\lambda} r
362 : * \f]
363 : */
364 1 : void beta_worldtube_data(
365 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*> beta,
366 : const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r);
367 :
368 : /*!
369 : * \brief Compute the \f$U\f$ (shift) function for the CCE Bondi-like metric.
370 : *
371 : * \details The Bondi-like metric has \f$g^{r A} = -e^{-2 \beta} U^A\f$, and the
372 : * spin-weighted vector \f$U = U^A q_A\f$. The value of \f$U^A\f$ can be
373 : * computed from the intermediate null metric quantities (see equation (54) of
374 : * \cite Barkett2019uae) using:
375 : *
376 : * \f[
377 : * U = -(\partial_\lambda r g^{\lambda A} + \partial_B r g^{A B}) q_A
378 : * / \partial_\lambda r \f]
379 : *
380 : */
381 1 : void bondi_u_worldtube_data(
382 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> bondi_u,
383 : const tnsr::i<ComplexDataVector, 2, Frame::RadialNull>& dyad,
384 : const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
385 : const tnsr::AA<DataVector, 3, Frame::RadialNull>& inverse_null_metric);
386 :
387 : /*!
388 : * \brief Compute the \f$W\f$ (mass aspect) function for the CCE Bondi-like
389 : * metric.
390 : *
391 : * \details The Bondi-like metric has \f$g^{rr} = e^{-2 \beta}(1 + r W)\f$. The
392 : * value of \f$W\f$ can be computed from the null metric quantities (see
393 : * equation (55) of \cite Barkett2019uae) using:
394 : *
395 : * \f[
396 : * W = \frac{1}{r} \left(-1
397 : * + \frac{g^{\lambda \lambda} (\partial_\lambda r)^2
398 : * + 2 \partial_\lambda r \left(\partial_A r g^{\lambda A}
399 : * - \partial_u r\right) + \partial_A r \partial_B r g^{A B}}
400 : * {\partial_\lambda r}\right) \f]
401 : */
402 1 : void bondi_w_worldtube_data(
403 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*> bondi_w,
404 : const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
405 : const tnsr::AA<DataVector, 3, Frame::RadialNull>& inverse_null_metric,
406 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r);
407 :
408 : /*!
409 : * \brief Compute the \f$J\f$ (intuitively similar to the transverse-traceless
410 : * part of the angular metric) function for the CCE Bondi-like metric.
411 : *
412 : * \details The Bondi-like metric has \f$J = \frac{1}{2 r^2} q^A q^B g_{A B}\f$.
413 : * This expression holds both for the right-hand side in the Bondi coordinates
414 : * and for the right-hand side in the intermediate null coordinates (see
415 : * equation (45) of \cite Barkett2019uae).
416 : */
417 1 : void bondi_j_worldtube_data(
418 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> bondi_j,
419 : const tnsr::aa<DataVector, 3, Frame::RadialNull>& null_metric,
420 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
421 : const tnsr::I<ComplexDataVector, 2, Frame::RadialNull>& dyad);
422 :
423 : /*!
424 : * \brief Compute the radial derivative of the angular metric spin-weighted
425 : * scalar \f$\partial_r J\f$ in the CCE Bondi-like metric.
426 : *
427 : * \details The radial derivative of the angular spin-weighted scalar \f$J\f$
428 : * can be computed from the null metric components by (c.f. equation (47) of
429 : * \cite Barkett2019uae):
430 : *
431 : * \f[
432 : * \partial_r J = \frac{\partial_\lambda J}{\partial_\lambda r} =
433 : * \frac{q^A q^B \partial_\lambda g_{A B} / (2 r^2)
434 : * - 2 \partial_\lambda r J / r}{\partial_\lambda r}
435 : * \f]
436 : */
437 1 : void dr_bondi_j(
438 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> dr_bondi_j,
439 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
440 : denominator_buffer,
441 : const tnsr::aa<DataVector, 3, Frame::RadialNull>& dlambda_null_metric,
442 : const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
443 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& bondi_j,
444 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
445 : const tnsr::I<ComplexDataVector, 2, Frame::RadialNull>& dyad);
446 :
447 : /*!
448 : * \brief Compute the second derivative of the Bondi radius with respect to the
449 : * intermediate null coordinate radius \f$\partial_\lambda^2 r\f$.
450 : *
451 : * \details To determine this second derivative quantity without resorting to
452 : * depending on second-derivative metric inputs, we need to take advantage of
453 : * one of the Einstein field equations. Combining equations (53) and (52) of
454 : * \cite Barkett2019uae, we have:
455 : *
456 : * \f[
457 : * \partial_\lambda^2 r = \frac{-r}{4} \left(
458 : * \partial_\lambda J \partial_\lambda \bar J - (\partial_\lambda K)^2\right)
459 : * \f],
460 : *
461 : * where the first derivative of \f$K\f$ can be obtained from \f$K = \sqrt{1 + J
462 : * \bar J}\f$ and the first derivative of \f$J\f$ can be obtained from (47) of
463 : * \cite Barkett2019uae
464 : */
465 1 : void d2lambda_bondi_r(
466 : gsl::not_null<Scalar<DataVector>*> d2lambda_bondi_r,
467 : const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
468 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& dr_bondi_j,
469 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& bondi_j,
470 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r);
471 :
472 : /*!
473 : * \brief Compute the Bondi metric contribution \f$Q\f$ (radial derivative of
474 : * shift).
475 : *
476 : * \details The definition of \f$Q\f$ in terms of the Bondi metric components is
477 : *
478 : * \f[
479 : * Q = q^A e^{-2 \beta} g_{A B} \partial_r U^B.
480 : * \f]
481 : *
482 : * $Q$ can be derived from the intermediate null metric quantities via (see
483 : * equations (56) and (57) of \cite Barkett2019uae)
484 : *
485 : * \f[
486 : * \partial_\lambda U = - \left(\partial_\lambda g^{\lambda A}
487 : * + \frac{\partial_A \partial_\lambda r}{\partial_\lambda r} g^{A B}
488 : * + \frac{\partial_B r}{\partial_\lambda r} \partial_\lambda g^{A B}\right) q_A
489 : * + 2 \partial_\lambda \beta (U + g^{\lambda A} q_A)
490 : * \f]
491 : *
492 : * and
493 : *
494 : * \f[
495 : * Q = r^2 (J \partial_\lambda \bar U + K \partial_\lambda U)
496 : * \f]
497 : *
498 : * also provided is \f$\partial_r U\f$, which is separately useful to cache for
499 : * other intermediate steps in the CCE computation.
500 : */
501 1 : void bondi_q_worldtube_data(
502 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> bondi_q,
503 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> dr_bondi_u,
504 : const Scalar<DataVector>& d2lambda_r,
505 : const tnsr::AA<DataVector, 3, Frame::RadialNull>&
506 : dlambda_inverse_null_metric,
507 : const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
508 : const tnsr::i<ComplexDataVector, 2, Frame::RadialNull>& dyad,
509 : const tnsr::i<DataVector, 2, Frame::RadialNull>& angular_d_dlambda_r,
510 : const tnsr::AA<DataVector, 3, Frame::RadialNull>& inverse_null_metric,
511 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& bondi_j,
512 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
513 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& bondi_u);
514 :
515 : /*!
516 : * \brief Compute the Bondi metric contribution \f$(\partial_u J)_{y} \equiv
517 : * H\f$ (the retarded time derivative evaluated at fixed $y$ coordinate) on the
518 : * worldtube boundary.
519 : *
520 : * \details The numerical time derivative (along the worldtube, rather than
521 : * along the surface of constant Bondi \f$r\f$) is computed by (see equation
522 : * (48) of \cite Barkett2019uae)
523 : *
524 : * \f[
525 : * (\partial_u J)_y = \frac{1}{2 r^2} q^A q^B \partial_u g_{A B}
526 : * - \frac{2 \partial_u r}{r} J
527 : * \f]
528 : *
529 : * \note There is the regrettable notation difference with the primary reference
530 : * for these formulas \cite Barkett2019uae in that we denote with \f$H\f$ the
531 : * time derivative at constant numerical radius, where \cite Barkett2019uae uses
532 : * \f$H\f$ to denote the time derivative at constant Bondi radius.
533 : */
534 1 : void bondi_h_worldtube_data(
535 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> bondi_h,
536 : const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
537 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& bondi_j,
538 : const tnsr::aa<DataVector, 3, Frame::RadialNull>& du_null_metric,
539 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
540 : const tnsr::I<ComplexDataVector, 2, Frame::RadialNull>& dyad);
541 :
542 : /*!
543 : * \brief Compute the Bondi metric contribution \f$(\partial_u J)_r\f$ (the
544 : * retarded time derivative at fixed coordinate $r$) on the worldtube boundary.
545 : *
546 : * \details The numerical time derivative (along the surface of constant r, not
547 : * along the worldtube) is computed by (see equation (50) of
548 : * \cite Barkett2019uae)
549 : *
550 : * \f[
551 : * \partial_u J = \frac{1}{2 r^2} q^A q^B \left(\partial_u g_{A B}
552 : * - \frac{ \partial_u r}{ \partial_\lambda r} \partial_\lambda g_{A B}\right)
553 : * \f]
554 : *
555 : * \note There is the regrettable notation difference with the primary reference
556 : * for these formulas \cite Barkett2019uae in that we denote with \f$H\f$ the
557 : * time derivative at constant numerical radius, where \cite Barkett2019uae uses
558 : * \f$H\f$ to denote the time derivative at constant Bondi radius.
559 : */
560 1 : void du_j_worldtube_data(
561 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> du_bondi_j,
562 : const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
563 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& bondi_j,
564 : const tnsr::aa<DataVector, 3, Frame::RadialNull>& du_null_metric,
565 : const tnsr::aa<DataVector, 3, Frame::RadialNull>& dlambda_null_metric,
566 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
567 : const tnsr::I<ComplexDataVector, 2, Frame::RadialNull>& dyad);
568 :
569 : namespace Tags {
570 : /*!
571 : * \brief The collection of tags mutated by `create_bondi_boundary_data`
572 : *
573 : * \details This list is used in the evolution of `CharacteristicExtract`
574 : */
575 : template <template <typename> class BoundaryPrefix>
576 1 : using characteristic_worldtube_boundary_tags = db::wrap_tags_in<
577 : BoundaryPrefix,
578 : tmpl::list<Tags::BondiBeta, Tags::BondiU, Tags::Dr<Tags::BondiU>,
579 : Tags::BondiQ, Tags::BondiW, Tags::BondiJ, Tags::Dr<Tags::BondiJ>,
580 : Tags::BondiH, Tags::Du<Tags::BondiJ>, Tags::BondiR,
581 : Tags::Du<Tags::BondiR>, Tags::DuRDividedByR>>;
582 :
583 : /*!
584 : * \brief The collection of tags for worldtube quantities that need to be
585 : * written to disk during the Cauchy evolution that the `CharacateristicExtract`
586 : * need.
587 : *
588 : * \details This list is used when writing Bondi quantities to disk.
589 : */
590 : template <template <typename> class BoundaryPrefix = Cce::Tags::BoundaryValue>
591 1 : using worldtube_boundary_tags_for_writing = db::wrap_tags_in<
592 : BoundaryPrefix,
593 : tmpl::list<Cce::Tags::BondiBeta, Cce::Tags::Dr<Cce::Tags::BondiJ>,
594 : Cce::Tags::Du<Cce::Tags::BondiR>, Cce::Tags::BondiJ,
595 : Cce::Tags::Du<Cce::Tags::BondiJ>, Cce::Tags::BondiQ,
596 : Cce::Tags::BondiR, Cce::Tags::BondiU, Cce::Tags::BondiW>>;
597 :
598 0 : using klein_gordon_worldtube_boundary_tags =
599 : tmpl::list<Tags::BoundaryValue<Tags::KleinGordonPsi>,
600 : Tags::BoundaryValue<Tags::KleinGordonPi>>;
601 : } // namespace Tags
602 :
603 : namespace detail {
604 : // the common step between the modal input and the Generalized harmonic input
605 : // that performs the final gauge processing to Bondi scalars and places them in
606 : // the Variables.
607 : template <typename BoundaryTagList, typename BufferTagList,
608 : typename ComplexBufferTagList>
609 : void create_bondi_boundary_data(
610 : const gsl::not_null<Variables<BoundaryTagList>*> bondi_boundary_data,
611 : const gsl::not_null<Variables<BufferTagList>*> computation_variables,
612 : const gsl::not_null<Variables<ComplexBufferTagList>*> derivative_buffers,
613 : const tnsr::aa<DataVector, 3>& dt_spacetime_metric,
614 : const tnsr::iaa<DataVector, 3>& phi,
615 : const tnsr::aa<DataVector, 3>& spacetime_metric,
616 : const tnsr::A<DataVector, 3>& null_l,
617 : const tnsr::A<DataVector, 3>& du_null_l,
618 : const SphericaliCartesianJ& cartesian_to_spherical_jacobian,
619 : const size_t l_max, const double extraction_radius) {
620 : const size_t size = Spectral::Swsh::number_of_swsh_collocation_points(l_max);
621 :
622 : // unfortunately, because the dyads are not themselves spin-weighted, they
623 : // need a separate Variables
624 : Variables<tmpl::list<Tags::detail::DownDyad, Tags::detail::UpDyad>>
625 : dyad_variables{size};
626 :
627 : auto& null_metric =
628 : get<gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>(
629 : *computation_variables);
630 : auto& du_null_metric = get<
631 : ::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>>(
632 : *computation_variables);
633 : null_metric_and_derivative(
634 : make_not_null(&du_null_metric), make_not_null(&null_metric),
635 : cartesian_to_spherical_jacobian, dt_spacetime_metric, spacetime_metric);
636 :
637 : auto& inverse_null_metric =
638 : get<gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>>(
639 : *computation_variables);
640 :
641 : // the below scaling process is used to reduce accumulation of numerical
642 : // error in the determinant evaluation
643 :
644 : // buffer reuse because the scaled null metric is only needed until the
645 : // `determinant_and_inverse` call
646 : auto& scaled_null_metric =
647 : get<gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>>(
648 : *computation_variables);
649 : for (size_t i = 0; i < 4; ++i) {
650 : for (size_t j = i; j < 4; ++j) {
651 : if (i > 1 and j > 1) {
652 : scaled_null_metric.get(i, j) =
653 : null_metric.get(i, j) / square(extraction_radius);
654 : } else if (i > 1 or j > 1) {
655 : scaled_null_metric.get(i, j) =
656 : null_metric.get(i, j) / extraction_radius;
657 : } else {
658 : scaled_null_metric.get(i, j) = null_metric.get(i, j);
659 : }
660 : }
661 : }
662 : // Allocation
663 : const auto scaled_inverse_null_metric =
664 : determinant_and_inverse(scaled_null_metric).second;
665 : for (size_t i = 0; i < 4; ++i) {
666 : for (size_t j = i; j < 4; ++j) {
667 : if (i > 1 and j > 1) {
668 : inverse_null_metric.get(i, j) =
669 : scaled_inverse_null_metric.get(i, j) / square(extraction_radius);
670 : } else if (i > 1 or j > 1) {
671 : inverse_null_metric.get(i, j) =
672 : scaled_inverse_null_metric.get(i, j) / extraction_radius;
673 : } else {
674 : inverse_null_metric.get(i, j) = scaled_inverse_null_metric.get(i, j);
675 : }
676 : }
677 : }
678 :
679 : auto& angular_d_null_l =
680 : get<Tags::detail::AngularDNullL>(*computation_variables);
681 : auto& buffer_for_derivatives =
682 : get(get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
683 : std::integral_constant<int, 0>>>(
684 : *derivative_buffers));
685 : auto& eth_buffer =
686 : get(get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
687 : std::integral_constant<int, 1>>>(
688 : *derivative_buffers));
689 : for (size_t a = 0; a < 4; ++a) {
690 : buffer_for_derivatives.data() =
691 : std::complex<double>(1.0, 0.0) * null_l.get(a);
692 : Spectral::Swsh::angular_derivatives<tmpl::list<Spectral::Swsh::Tags::Eth>>(
693 : l_max, 1, make_not_null(ð_buffer), buffer_for_derivatives);
694 : angular_d_null_l.get(0, a) = -real(eth_buffer.data());
695 : angular_d_null_l.get(1, a) = -imag(eth_buffer.data());
696 : }
697 :
698 : auto& dlambda_null_metric = get<Tags::detail::DLambda<
699 : gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>>(
700 : *computation_variables);
701 : auto& dlambda_inverse_null_metric = get<Tags::detail::DLambda<
702 : gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>>>(
703 : *computation_variables);
704 : dlambda_null_metric_and_inverse(
705 : make_not_null(&dlambda_null_metric),
706 : make_not_null(&dlambda_inverse_null_metric), angular_d_null_l,
707 : cartesian_to_spherical_jacobian, phi, dt_spacetime_metric, du_null_l,
708 : inverse_null_metric, null_l, spacetime_metric);
709 :
710 : auto& r = get<Tags::BoundaryValue<Tags::BondiR>>(*bondi_boundary_data);
711 : bondi_r(make_not_null(&r), null_metric);
712 :
713 : auto& d_r =
714 : get<::Tags::spacetime_deriv<Tags::detail::RealBondiR, tmpl::size_t<3>,
715 : Frame::RadialNull>>(*computation_variables);
716 : d_bondi_r(make_not_null(&d_r), r, dlambda_null_metric, du_null_metric,
717 : inverse_null_metric, l_max);
718 : get(get<Tags::BoundaryValue<Tags::DuRDividedByR>>(*bondi_boundary_data))
719 : .data() = std::complex<double>(1.0, 0.0) * get<0>(d_r) / get(r).data();
720 : get(get<Tags::BoundaryValue<Tags::Du<Tags::BondiR>>>(*bondi_boundary_data))
721 : .data() = std::complex<double>(1.0, 0.0) * get<0>(d_r);
722 :
723 : auto& down_dyad = get<Tags::detail::DownDyad>(dyad_variables);
724 : auto& up_dyad = get<Tags::detail::UpDyad>(dyad_variables);
725 : dyads(make_not_null(&down_dyad), make_not_null(&up_dyad));
726 :
727 : beta_worldtube_data(make_not_null(&get<Tags::BoundaryValue<Tags::BondiBeta>>(
728 : *bondi_boundary_data)),
729 : d_r);
730 :
731 : auto& bondi_u = get<Tags::BoundaryValue<Tags::BondiU>>(*bondi_boundary_data);
732 : bondi_u_worldtube_data(make_not_null(&bondi_u), down_dyad, d_r,
733 : inverse_null_metric);
734 :
735 : bondi_w_worldtube_data(make_not_null(&get<Tags::BoundaryValue<Tags::BondiW>>(
736 : *bondi_boundary_data)),
737 : d_r, inverse_null_metric, r);
738 :
739 : auto& bondi_j = get<Tags::BoundaryValue<Tags::BondiJ>>(*bondi_boundary_data);
740 : bondi_j_worldtube_data(make_not_null(&bondi_j), null_metric, r, up_dyad);
741 :
742 : auto& dr_j =
743 : get<Tags::BoundaryValue<Tags::Dr<Tags::BondiJ>>>(*bondi_boundary_data);
744 : auto& denominator_buffer =
745 : get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
746 : std::integral_constant<int, 0>>>(
747 : *derivative_buffers);
748 : dr_bondi_j(make_not_null(&get<Tags::BoundaryValue<Tags::Dr<Tags::BondiJ>>>(
749 : *bondi_boundary_data)),
750 : make_not_null(&denominator_buffer), dlambda_null_metric, d_r,
751 : bondi_j, r, up_dyad);
752 :
753 : auto& d2lambda_r = get<
754 : Tags::detail::DLambda<Tags::detail::DLambda<Tags::detail::RealBondiR>>>(
755 : *computation_variables);
756 : d2lambda_bondi_r(make_not_null(&d2lambda_r), d_r, dr_j, bondi_j, r);
757 :
758 : auto& angular_d_dlambda_r =
759 : get<::Tags::deriv<Tags::detail::DLambda<Tags::detail::RealBondiR>,
760 : tmpl::size_t<2>, Frame::RadialNull>>(
761 : *computation_variables);
762 : buffer_for_derivatives.data() = std::complex<double>(1.0, 0.0) * get<1>(d_r);
763 : Spectral::Swsh::angular_derivatives<tmpl::list<Spectral::Swsh::Tags::Eth>>(
764 : l_max, 1, make_not_null(ð_buffer), buffer_for_derivatives);
765 : angular_d_dlambda_r.get(0) = -real(eth_buffer.data());
766 : angular_d_dlambda_r.get(1) = -imag(eth_buffer.data());
767 :
768 : bondi_q_worldtube_data(
769 : make_not_null(
770 : &get<Tags::BoundaryValue<Tags::BondiQ>>(*bondi_boundary_data)),
771 : make_not_null(&get<Tags::BoundaryValue<Tags::Dr<Tags::BondiU>>>(
772 : *bondi_boundary_data)),
773 : d2lambda_r, dlambda_inverse_null_metric, d_r, down_dyad,
774 : angular_d_dlambda_r, inverse_null_metric, bondi_j, r, bondi_u);
775 :
776 : bondi_h_worldtube_data(make_not_null(&get<Tags::BoundaryValue<Tags::BondiH>>(
777 : *bondi_boundary_data)),
778 : d_r, bondi_j, du_null_metric, r, up_dyad);
779 :
780 : du_j_worldtube_data(
781 : make_not_null(&get<Tags::BoundaryValue<Tags::Du<Tags::BondiJ>>>(
782 : *bondi_boundary_data)),
783 : d_r, bondi_j, du_null_metric, dlambda_null_metric, r, up_dyad);
784 : }
785 : } // namespace detail
786 :
787 : /*!
788 : * \brief Process the worldtube data from generalized harmonic quantities
789 : * to desired Bondi quantities, placing the result in the passed
790 : * `Variables`.
791 : *
792 : * \details
793 : * The mathematics are a bit complicated for all of the coordinate
794 : * transformations that are necessary to obtain the Bondi gauge quantities.
795 : * For full mathematical details, see the documentation for functions in
796 : * `BoundaryData.hpp` and \cite Barkett2019uae \cite Bishop1998uk.
797 : *
798 : * This function takes as input the full set of Generalized harmonic metric data
799 : * on a two-dimensional surface of constant \f$r\f$ and \f$t\f$ in numerical
800 : * coordinates.
801 : *
802 : * Sufficient tags to provide full worldtube boundary data at a particular
803 : * time are set in `bondi_boundary_data`. In particular, the set of tags in
804 : * `Tags::characteristic_worldtube_boundary_tags` in the provided `Variables`
805 : * are assigned to the worldtube boundary values associated with the input
806 : * metric components.
807 : *
808 : * The majority of the mathematical transformations are implemented as a set of
809 : * individual cascaded functions below. The details of the manipulations that
810 : * are performed to the input data may be found in the individual functions
811 : * themselves, which are called in the following order:
812 : * - `trigonometric_functions_on_swsh_collocation()`
813 : * - `gr::shift()`
814 : * - `gr::lapse()`
815 : * - `worldtube_normal_and_derivatives()`
816 : * - `gr::spacetime_normal_vector()`
817 : * - `gh::time_deriv_of_lapse()`
818 : * - `gh::time_deriv_of_shift()`
819 : * - `null_vector_l_and_derivatives()`
820 : * - `cartesian_to_spherical_coordinates_and_jacobians()`
821 : * - `null_metric_and_derivative()`
822 : * - `dlambda_null_metric_and_inverse()`
823 : * - `bondi_r()`
824 : * - `d_bondi_r()`
825 : * - `dyads()`
826 : * - `beta_worldtube_data()`
827 : * - `bondi_u_worldtube_data()`
828 : * - `bondi_w_worldtube_data()`
829 : * - `bondi_j_worldtube_data()`
830 : * - `dr_bondi_j()`
831 : * - `d2lambda_bondi_r()`
832 : * - `bondi_q_worldtube_data()`
833 : * - `bondi_h_worldtube_data()`
834 : * - `du_j_worldtube_data()`
835 : */
836 : template <typename BoundaryTagList>
837 1 : void create_bondi_boundary_data(
838 : const gsl::not_null<Variables<BoundaryTagList>*> bondi_boundary_data,
839 : const tnsr::iaa<DataVector, 3>& phi, const tnsr::aa<DataVector, 3>& pi,
840 : const tnsr::aa<DataVector, 3>& spacetime_metric,
841 : const double extraction_radius, const size_t l_max) {
842 : const size_t size = Spectral::Swsh::number_of_swsh_collocation_points(l_max);
843 : // Most allocations required for the full boundary computation are merged into
844 : // a single, large Variables allocation. There remain a handful of cases in
845 : // the computational functions called where an intermediate quantity that is
846 : // not re-used is allocated rather than taking a buffer. These cases are
847 : // marked with code comments 'Allocation'; In the future, if allocations are
848 : // identified as a point to optimize, those buffers may be allocated here and
849 : // passed as function arguments
850 : Variables<tmpl::list<
851 : Tags::detail::CosPhi, Tags::detail::CosTheta, Tags::detail::SinPhi,
852 : Tags::detail::SinTheta, Tags::detail::CartesianCoordinates,
853 : Tags::detail::CartesianToSphericalJacobian,
854 : Tags::detail::InverseCartesianToSphericalJacobian,
855 : gr::Tags::SpatialMetric<DataVector, 3>,
856 : gr::Tags::InverseSpatialMetric<DataVector, 3>,
857 : gr::Tags::Shift<DataVector, 3>,
858 : ::Tags::dt<gr::Tags::Shift<DataVector, 3>>, gr::Tags::Lapse<DataVector>,
859 : ::Tags::dt<gr::Tags::Lapse<DataVector>>,
860 : ::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3>>,
861 : Tags::detail::WorldtubeNormal, ::Tags::dt<Tags::detail::WorldtubeNormal>,
862 : gr::Tags::SpacetimeNormalVector<DataVector, 3>, Tags::detail::NullL,
863 : ::Tags::dt<Tags::detail::NullL>,
864 : // for the detail function called at the end
865 : gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>,
866 : ::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>,
867 : gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>,
868 : Tags::detail::AngularDNullL,
869 : Tags::detail::DLambda<
870 : gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>,
871 : Tags::detail::DLambda<
872 : gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>>,
873 : ::Tags::spacetime_deriv<Tags::detail::RealBondiR, tmpl::size_t<3>,
874 : Frame::RadialNull>,
875 : Tags::detail::DLambda<Tags::detail::DLambda<Tags::detail::RealBondiR>>,
876 : ::Tags::deriv<Tags::detail::DLambda<Tags::detail::RealBondiR>,
877 : tmpl::size_t<2>, Frame::RadialNull>>>
878 : computation_variables{size};
879 :
880 : Variables<
881 : tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
882 : std::integral_constant<int, 0>>,
883 : ::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
884 : std::integral_constant<int, 1>>>>
885 : derivative_buffers{size};
886 :
887 : auto& cos_phi = get<Tags::detail::CosPhi>(computation_variables);
888 : auto& cos_theta = get<Tags::detail::CosTheta>(computation_variables);
889 : auto& sin_phi = get<Tags::detail::SinPhi>(computation_variables);
890 : auto& sin_theta = get<Tags::detail::SinTheta>(computation_variables);
891 : trigonometric_functions_on_swsh_collocation(
892 : make_not_null(&cos_phi), make_not_null(&cos_theta),
893 : make_not_null(&sin_phi), make_not_null(&sin_theta), l_max);
894 :
895 : // NOTE: to handle the singular values of polar coordinates, the phi
896 : // components of all tensors are scaled according to their sin(theta)
897 : // prefactors.
898 : // so, any down-index component get<2>(A) represents 1/sin(theta) A_\phi,
899 : // and any up-index component get<2>(A) represents sin(theta) A^\phi.
900 : // This holds for Jacobians, and so direct application of the Jacobians
901 : // brings the factors through.
902 : auto& cartesian_coords =
903 : get<Tags::detail::CartesianCoordinates>(computation_variables);
904 : auto& cartesian_to_spherical_jacobian =
905 : get<Tags::detail::CartesianToSphericalJacobian>(computation_variables);
906 : auto& inverse_cartesian_to_spherical_jacobian =
907 : get<Tags::detail::InverseCartesianToSphericalJacobian>(
908 : computation_variables);
909 : cartesian_to_spherical_coordinates_and_jacobians(
910 : make_not_null(&cartesian_coords),
911 : make_not_null(&cartesian_to_spherical_jacobian),
912 : make_not_null(&inverse_cartesian_to_spherical_jacobian), cos_phi,
913 : cos_theta, sin_phi, sin_theta, extraction_radius);
914 :
915 : auto& spatial_metric =
916 : get<gr::Tags::SpatialMetric<DataVector, 3>>(computation_variables);
917 : gr::spatial_metric(make_not_null(&spatial_metric), spacetime_metric);
918 :
919 : auto& inverse_spatial_metric =
920 : get<gr::Tags::InverseSpatialMetric<DataVector, 3>>(computation_variables);
921 : // Allocation
922 : inverse_spatial_metric = determinant_and_inverse(spatial_metric).second;
923 :
924 : auto& shift = get<gr::Tags::Shift<DataVector, 3>>(computation_variables);
925 : gr::shift(make_not_null(&shift), spacetime_metric, inverse_spatial_metric);
926 :
927 : auto& lapse = get<gr::Tags::Lapse<DataVector>>(computation_variables);
928 : gr::lapse(make_not_null(&lapse), shift, spacetime_metric);
929 :
930 : auto& dt_spacetime_metric =
931 : get<::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3>>>(
932 : computation_variables);
933 :
934 : gh::time_derivative_of_spacetime_metric(make_not_null(&dt_spacetime_metric),
935 : lapse, shift, pi, phi);
936 :
937 : auto& dt_worldtube_normal =
938 : get<::Tags::dt<Tags::detail::WorldtubeNormal>>(computation_variables);
939 : auto& worldtube_normal =
940 : get<Tags::detail::WorldtubeNormal>(computation_variables);
941 : worldtube_normal_and_derivatives(
942 : make_not_null(&worldtube_normal), make_not_null(&dt_worldtube_normal),
943 : cos_phi, cos_theta, spacetime_metric, dt_spacetime_metric, sin_phi,
944 : sin_theta, inverse_spatial_metric);
945 : auto& spacetime_unit_normal =
946 : get<gr::Tags::SpacetimeNormalVector<DataVector, 3>>(
947 : computation_variables);
948 : gr::spacetime_normal_vector(make_not_null(&spacetime_unit_normal), lapse,
949 : shift);
950 : auto& dt_lapse =
951 : get<::Tags::dt<gr::Tags::Lapse<DataVector>>>(computation_variables);
952 : gh::time_deriv_of_lapse(make_not_null(&dt_lapse), lapse, shift,
953 : spacetime_unit_normal, phi, pi);
954 : auto& dt_shift =
955 : get<::Tags::dt<gr::Tags::Shift<DataVector, 3>>>(computation_variables);
956 : gh::time_deriv_of_shift(make_not_null(&dt_shift), lapse, shift,
957 : inverse_spatial_metric, spacetime_unit_normal, phi,
958 : pi);
959 :
960 : auto& du_null_l = get<::Tags::dt<Tags::detail::NullL>>(computation_variables);
961 : auto& null_l = get<Tags::detail::NullL>(computation_variables);
962 : null_vector_l_and_derivatives(make_not_null(&du_null_l),
963 : make_not_null(&null_l), dt_worldtube_normal,
964 : dt_lapse, dt_spacetime_metric, dt_shift, lapse,
965 : spacetime_metric, shift, worldtube_normal);
966 :
967 : // pass to the next step that is common between the 'modal' input and 'GH'
968 : // input strategies
969 : detail::create_bondi_boundary_data(
970 : bondi_boundary_data, make_not_null(&computation_variables),
971 : make_not_null(&derivative_buffers), dt_spacetime_metric, phi,
972 : spacetime_metric, null_l, du_null_l, cartesian_to_spherical_jacobian,
973 : l_max, extraction_radius);
974 : }
975 :
976 : /*!
977 : * \brief Process the worldtube data from modal metric components and
978 : * derivatives to desired Bondi quantities, placing the result in the passed
979 : * `Variables`.
980 : *
981 : * \details
982 : * The mathematics are a bit complicated for all of the coordinate
983 : * transformations that are necessary to obtain the Bondi gauge quantities.
984 : * For full mathematical details, see the documentation for functions in
985 : * `BoundaryData.hpp` and \cite Barkett2019uae \cite Bishop1998uk.
986 : *
987 : * This function takes as input the full set of ADM metric data and its radial
988 : * and time derivatives on a two-dimensional surface of constant \f$r\f$ and
989 : * \f$t\f$ in numerical coordinates. This data must be provided as spherical
990 : * harmonic coefficients in the libsharp format. This data is provided in nine
991 : * `Tensor`s.
992 : *
993 : * Sufficient tags to provide full worldtube boundary data at a particular
994 : * time are set in `bondi_boundary_data`. In particular, the set of tags in
995 : * `Tags::characteristic_worldtube_boundary_tags` in the provided `Variables`
996 : * are assigned to the worldtube boundary values associated with the input
997 : * metric components.
998 : *
999 : * The majority of the mathematical transformations are implemented as a set of
1000 : * individual cascaded functions below. The details of the manipulations that
1001 : * are performed to the input data may be found in the individual functions
1002 : * themselves, which are called in the following order:
1003 : * - `trigonometric_functions_on_swsh_collocation()`
1004 : * - `cartesian_to_spherical_coordinates_and_jacobians()`
1005 : * - `cartesian_spatial_metric_and_derivatives_from_modes()`
1006 : * - `cartesian_shift_and_derivatives_from_modes()`
1007 : * - `cartesian_lapse_and_derivatives_from_modes()`
1008 : * - `gh::phi()`
1009 : * - `gr::time_derivative_of_spacetime_metric`
1010 : * - `gr::spacetime_metric`
1011 : * - `generalized_harmonic_quantities()`
1012 : * - `worldtube_normal_and_derivatives()`
1013 : * - `null_vector_l_and_derivatives()`
1014 : * - `null_metric_and_derivative()`
1015 : * - `dlambda_null_metric_and_inverse()`
1016 : * - `bondi_r()`
1017 : * - `d_bondi_r()`
1018 : * - `dyads()`
1019 : * - `beta_worldtube_data()`
1020 : * - `bondi_u_worldtube_data()`
1021 : * - `bondi_w_worldtube_data()`
1022 : * - `bondi_j_worldtube_data()`
1023 : * - `dr_bondi_j()`
1024 : * - `d2lambda_bondi_r()`
1025 : * - `bondi_q_worldtube_data()`
1026 : * - `bondi_h_worldtube_data()`
1027 : * - `du_j_worldtube_data()`
1028 : */
1029 : template <typename BoundaryTagList>
1030 1 : void create_bondi_boundary_data(
1031 : const gsl::not_null<Variables<BoundaryTagList>*> bondi_boundary_data,
1032 : const tnsr::ii<ComplexModalVector, 3>& spatial_metric_coefficients,
1033 : const tnsr::ii<ComplexModalVector, 3>& dt_spatial_metric_coefficients,
1034 : const tnsr::ii<ComplexModalVector, 3>& dr_spatial_metric_coefficients,
1035 : const tnsr::I<ComplexModalVector, 3>& shift_coefficients,
1036 : const tnsr::I<ComplexModalVector, 3>& dt_shift_coefficients,
1037 : const tnsr::I<ComplexModalVector, 3>& dr_shift_coefficients,
1038 : const Scalar<ComplexModalVector>& lapse_coefficients,
1039 : const Scalar<ComplexModalVector>& dt_lapse_coefficients,
1040 : const Scalar<ComplexModalVector>& dr_lapse_coefficients,
1041 : const double extraction_radius, const size_t l_max) {
1042 : const size_t size = Spectral::Swsh::number_of_swsh_collocation_points(l_max);
1043 :
1044 : // Most allocations required for the full boundary computation are merged into
1045 : // a single, large Variables allocation. There remain a handful of cases in
1046 : // the computational functions called where an intermediate quantity that is
1047 : // not re-used is allocated rather than taking a buffer. These cases are
1048 : // marked with code comments 'Allocation'; In the future, if allocations are
1049 : // identified as a point to optimize, those buffers may be allocated here and
1050 : // passed as function arguments
1051 : Variables<tmpl::list<
1052 : Tags::detail::CosPhi, Tags::detail::CosTheta, Tags::detail::SinPhi,
1053 : Tags::detail::SinTheta, Tags::detail::CartesianCoordinates,
1054 : Tags::detail::CartesianToSphericalJacobian,
1055 : Tags::detail::InverseCartesianToSphericalJacobian,
1056 : gr::Tags::SpatialMetric<DataVector, 3>,
1057 : gr::Tags::InverseSpatialMetric<DataVector, 3>,
1058 : ::Tags::deriv<gr::Tags::SpatialMetric<DataVector, 3>, tmpl::size_t<3>,
1059 : ::Frame::Inertial>,
1060 : ::Tags::dt<gr::Tags::SpatialMetric<DataVector, 3>>,
1061 : gr::Tags::Shift<DataVector, 3>,
1062 : ::Tags::deriv<gr::Tags::Shift<DataVector, 3>, tmpl::size_t<3>,
1063 : ::Frame::Inertial>,
1064 : ::Tags::dt<gr::Tags::Shift<DataVector, 3>>, gr::Tags::Lapse<DataVector>,
1065 : ::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
1066 : ::Frame::Inertial>,
1067 : ::Tags::dt<gr::Tags::Lapse<DataVector>>,
1068 : gr::Tags::SpacetimeMetric<DataVector, 3>,
1069 : ::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3>>,
1070 : gh::Tags::Phi<DataVector, 3>, Tags::detail::WorldtubeNormal,
1071 : ::Tags::dt<Tags::detail::WorldtubeNormal>, Tags::detail::NullL,
1072 : ::Tags::dt<Tags::detail::NullL>,
1073 : // for the detail function called at the end
1074 : gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>,
1075 : ::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>,
1076 : gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>,
1077 : Tags::detail::AngularDNullL,
1078 : Tags::detail::DLambda<
1079 : gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>,
1080 : Tags::detail::DLambda<
1081 : gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>>,
1082 : ::Tags::spacetime_deriv<Tags::detail::RealBondiR, tmpl::size_t<3>,
1083 : Frame::RadialNull>,
1084 : Tags::detail::DLambda<Tags::detail::DLambda<Tags::detail::RealBondiR>>,
1085 : ::Tags::deriv<Tags::detail::DLambda<Tags::detail::RealBondiR>,
1086 : tmpl::size_t<2>, Frame::RadialNull>>>
1087 : computation_variables{size};
1088 :
1089 : Variables<
1090 : tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
1091 : std::integral_constant<int, 0>>,
1092 : ::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
1093 : std::integral_constant<int, 1>>>>
1094 : derivative_buffers{size};
1095 : auto& cos_phi = get<Tags::detail::CosPhi>(computation_variables);
1096 : auto& cos_theta = get<Tags::detail::CosTheta>(computation_variables);
1097 : auto& sin_phi = get<Tags::detail::SinPhi>(computation_variables);
1098 : auto& sin_theta = get<Tags::detail::SinTheta>(computation_variables);
1099 : trigonometric_functions_on_swsh_collocation(
1100 : make_not_null(&cos_phi), make_not_null(&cos_theta),
1101 : make_not_null(&sin_phi), make_not_null(&sin_theta), l_max);
1102 :
1103 : // NOTE: to handle the singular values of polar coordinates, the phi
1104 : // components of all tensors are scaled according to their sin(theta)
1105 : // prefactors.
1106 : // so, any down-index component get<2>(A) represents 1/sin(theta) A_\phi,
1107 : // and any up-index component get<2>(A) represents sin(theta) A^\phi.
1108 : // This holds for Jacobians, and so direct application of the Jacobians
1109 : // brings the factors through.
1110 : auto& cartesian_coords =
1111 : get<Tags::detail::CartesianCoordinates>(computation_variables);
1112 : auto& cartesian_to_spherical_jacobian =
1113 : get<Tags::detail::CartesianToSphericalJacobian>(computation_variables);
1114 : auto& inverse_cartesian_to_spherical_jacobian =
1115 : get<Tags::detail::InverseCartesianToSphericalJacobian>(
1116 : computation_variables);
1117 : cartesian_to_spherical_coordinates_and_jacobians(
1118 : make_not_null(&cartesian_coords),
1119 : make_not_null(&cartesian_to_spherical_jacobian),
1120 : make_not_null(&inverse_cartesian_to_spherical_jacobian), cos_phi,
1121 : cos_theta, sin_phi, sin_theta, extraction_radius);
1122 :
1123 : auto& cartesian_spatial_metric =
1124 : get<gr::Tags::SpatialMetric<DataVector, 3>>(computation_variables);
1125 : auto& inverse_spatial_metric =
1126 : get<gr::Tags::InverseSpatialMetric<DataVector, 3>>(computation_variables);
1127 : auto& d_cartesian_spatial_metric =
1128 : get<::Tags::deriv<gr::Tags::SpatialMetric<DataVector, 3>, tmpl::size_t<3>,
1129 : ::Frame::Inertial>>(computation_variables);
1130 : auto& dt_cartesian_spatial_metric =
1131 : get<::Tags::dt<gr::Tags::SpatialMetric<DataVector, 3>>>(
1132 : computation_variables);
1133 : auto& interpolation_buffer =
1134 : get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
1135 : std::integral_constant<int, 0>>>(
1136 : derivative_buffers);
1137 : Scalar<SpinWeighted<ComplexModalVector, 0>> interpolation_modal_buffer{size};
1138 : auto& eth_buffer =
1139 : get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
1140 : std::integral_constant<int, 1>>>(
1141 : derivative_buffers);
1142 : cartesian_spatial_metric_and_derivatives_from_modes(
1143 : make_not_null(&cartesian_spatial_metric),
1144 : make_not_null(&inverse_spatial_metric),
1145 : make_not_null(&d_cartesian_spatial_metric),
1146 : make_not_null(&dt_cartesian_spatial_metric),
1147 : make_not_null(&interpolation_modal_buffer),
1148 : make_not_null(&interpolation_buffer), make_not_null(ð_buffer),
1149 : spatial_metric_coefficients, dr_spatial_metric_coefficients,
1150 : dt_spatial_metric_coefficients, inverse_cartesian_to_spherical_jacobian,
1151 : l_max);
1152 :
1153 : auto& cartesian_shift =
1154 : get<gr::Tags::Shift<DataVector, 3>>(computation_variables);
1155 : auto& d_cartesian_shift =
1156 : get<::Tags::deriv<gr::Tags::Shift<DataVector, 3>, tmpl::size_t<3>,
1157 : ::Frame::Inertial>>(computation_variables);
1158 : auto& dt_cartesian_shift =
1159 : get<::Tags::dt<gr::Tags::Shift<DataVector, 3>>>(computation_variables);
1160 :
1161 : cartesian_shift_and_derivatives_from_modes(
1162 : make_not_null(&cartesian_shift), make_not_null(&d_cartesian_shift),
1163 : make_not_null(&dt_cartesian_shift),
1164 : make_not_null(&interpolation_modal_buffer),
1165 : make_not_null(&interpolation_buffer), make_not_null(ð_buffer),
1166 : shift_coefficients, dr_shift_coefficients, dt_shift_coefficients,
1167 : inverse_cartesian_to_spherical_jacobian, l_max);
1168 :
1169 : auto& cartesian_lapse =
1170 : get<gr::Tags::Lapse<DataVector>>(computation_variables);
1171 : auto& d_cartesian_lapse =
1172 : get<::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
1173 : ::Frame::Inertial>>(computation_variables);
1174 : auto& dt_cartesian_lapse =
1175 : get<::Tags::dt<gr::Tags::Lapse<DataVector>>>(computation_variables);
1176 : cartesian_lapse_and_derivatives_from_modes(
1177 : make_not_null(&cartesian_lapse), make_not_null(&d_cartesian_lapse),
1178 : make_not_null(&dt_cartesian_lapse),
1179 : make_not_null(&interpolation_modal_buffer),
1180 : make_not_null(&interpolation_buffer), make_not_null(ð_buffer),
1181 : lapse_coefficients, dr_lapse_coefficients, dt_lapse_coefficients,
1182 : inverse_cartesian_to_spherical_jacobian, l_max);
1183 :
1184 : auto& phi = get<gh::Tags::Phi<DataVector, 3>>(computation_variables);
1185 : auto& dt_spacetime_metric =
1186 : get<::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3>>>(
1187 : computation_variables);
1188 : auto& spacetime_metric =
1189 : get<gr::Tags::SpacetimeMetric<DataVector, 3>>(computation_variables);
1190 : gh::phi(make_not_null(&phi), cartesian_lapse, d_cartesian_lapse,
1191 : cartesian_shift, d_cartesian_shift, cartesian_spatial_metric,
1192 : d_cartesian_spatial_metric);
1193 : gr::time_derivative_of_spacetime_metric(
1194 : make_not_null(&dt_spacetime_metric), cartesian_lapse, dt_cartesian_lapse,
1195 : cartesian_shift, dt_cartesian_shift, cartesian_spatial_metric,
1196 : dt_cartesian_spatial_metric);
1197 : gr::spacetime_metric(make_not_null(&spacetime_metric), cartesian_lapse,
1198 : cartesian_shift, cartesian_spatial_metric);
1199 :
1200 : auto& dt_worldtube_normal =
1201 : get<::Tags::dt<Tags::detail::WorldtubeNormal>>(computation_variables);
1202 : auto& worldtube_normal =
1203 : get<Tags::detail::WorldtubeNormal>(computation_variables);
1204 : worldtube_normal_and_derivatives(
1205 : make_not_null(&worldtube_normal), make_not_null(&dt_worldtube_normal),
1206 : cos_phi, cos_theta, spacetime_metric, dt_spacetime_metric, sin_phi,
1207 : sin_theta, inverse_spatial_metric);
1208 :
1209 : auto& du_null_l = get<::Tags::dt<Tags::detail::NullL>>(computation_variables);
1210 : auto& null_l = get<Tags::detail::NullL>(computation_variables);
1211 : null_vector_l_and_derivatives(
1212 : make_not_null(&du_null_l), make_not_null(&null_l), dt_worldtube_normal,
1213 : dt_cartesian_lapse, dt_spacetime_metric, dt_cartesian_shift,
1214 : cartesian_lapse, spacetime_metric, cartesian_shift, worldtube_normal);
1215 :
1216 : // pass to the next step that is common between the 'modal' input and 'GH'
1217 : // input strategies
1218 : detail::create_bondi_boundary_data(
1219 : bondi_boundary_data, make_not_null(&computation_variables),
1220 : make_not_null(&derivative_buffers), dt_spacetime_metric, phi,
1221 : spacetime_metric, null_l, du_null_l, cartesian_to_spherical_jacobian,
1222 : l_max, extraction_radius);
1223 : }
1224 :
1225 : /*!
1226 : * \brief Process the worldtube data from nodal metric components and
1227 : * derivatives to desired Bondi quantities, placing the result in the passed
1228 : * `Variables`.
1229 : *
1230 : * \details
1231 : * The mathematics are a bit complicated for all of the coordinate
1232 : * transformations that are necessary to obtain the Bondi gauge quantities.
1233 : * For full mathematical details, see the documentation for functions in
1234 : * `BoundaryData.hpp` and \cite Barkett2019uae \cite Bishop1998uk.
1235 : *
1236 : * This function takes as input the full set of ADM metric data and its radial
1237 : * and time derivatives on a two-dimensional surface of constant \f$r\f$ and
1238 : * \f$t\f$ in numerical coordinates. This data must be provided as values at
1239 : * the collocation points compatible with the libsharp format. This data is
1240 : * provided in nine `Tensor`s.
1241 : *
1242 : * Sufficient tags to provide full worldtube boundary data at a particular
1243 : * time are set in `bondi_boundary_data`. In particular, the set of tags in
1244 : * `Tags::characteristic_worldtube_boundary_tags` in the provided `Variables`
1245 : * are assigned to the worldtube boundary values associated with the input
1246 : * metric components.
1247 : *
1248 : * The majority of the mathematical transformations are implemented as a set of
1249 : * individual cascaded functions below. The details of the manipulations that
1250 : * are performed to the input data may be found in the individual functions
1251 : * themselves, which are called in the following order:
1252 : * - `trigonometric_functions_on_swsh_collocation()`
1253 : * - `gh::phi()`
1254 : * - `gr::time_derivative_of_spacetime_metric`
1255 : * - `gr::spacetime_metric`
1256 : * - `worldtube_normal_and_derivatives()`
1257 : * - `null_vector_l_and_derivatives()`
1258 : * - `cartesian_to_spherical_coordinates_and_jacobians()`
1259 : * - `null_metric_and_derivative()`
1260 : * - `dlambda_null_metric_and_inverse()`
1261 : * - `bondi_r()`
1262 : * - `d_bondi_r()`
1263 : * - `dyads()`
1264 : * - `beta_worldtube_data()`
1265 : * - `bondi_u_worldtube_data()`
1266 : * - `bondi_w_worldtube_data()`
1267 : * - `bondi_j_worldtube_data()`
1268 : * - `dr_bondi_j()`
1269 : * - `d2lambda_bondi_r()`
1270 : * - `bondi_q_worldtube_data()`
1271 : * - `bondi_h_worldtube_data()`
1272 : * - `du_j_worldtube_data()`
1273 : */
1274 : template <typename BoundaryTagList>
1275 1 : void create_bondi_boundary_data(
1276 : const gsl::not_null<Variables<BoundaryTagList>*> bondi_boundary_data,
1277 : const tnsr::ii<DataVector, 3>& cartesian_spatial_metric,
1278 : const tnsr::ii<DataVector, 3>& cartesian_dt_spatial_metric,
1279 : const tnsr::ii<DataVector, 3>& cartesian_dr_spatial_metric,
1280 : const tnsr::I<DataVector, 3>& cartesian_shift,
1281 : const tnsr::I<DataVector, 3>& cartesian_dt_shift,
1282 : const tnsr::I<DataVector, 3>& cartesian_dr_shift,
1283 : const Scalar<DataVector>& cartesian_lapse,
1284 : const Scalar<DataVector>& cartesian_dt_lapse,
1285 : const Scalar<DataVector>& cartesian_dr_lapse,
1286 : const double extraction_radius, const size_t l_max) {
1287 : const size_t size = Spectral::Swsh::number_of_swsh_collocation_points(l_max);
1288 :
1289 : Variables<tmpl::list<
1290 : Tags::detail::CosPhi, Tags::detail::CosTheta, Tags::detail::SinPhi,
1291 : Tags::detail::SinTheta, Tags::detail::CartesianCoordinates,
1292 : Tags::detail::CartesianToSphericalJacobian,
1293 : Tags::detail::InverseCartesianToSphericalJacobian,
1294 : ::Tags::deriv<gr::Tags::SpatialMetric<DataVector, 3>, tmpl::size_t<3>,
1295 : ::Frame::Inertial>,
1296 : ::Tags::deriv<gr::Tags::Shift<DataVector, 3>, tmpl::size_t<3>,
1297 : ::Frame::Inertial>,
1298 : ::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
1299 : ::Frame::Inertial>,
1300 : gr::Tags::InverseSpatialMetric<DataVector, 3>,
1301 : gr::Tags::SpacetimeMetric<DataVector, 3>,
1302 : ::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3>>,
1303 : gh::Tags::Phi<DataVector, 3>, Tags::detail::WorldtubeNormal,
1304 : ::Tags::dt<Tags::detail::WorldtubeNormal>, Tags::detail::NullL,
1305 : ::Tags::dt<Tags::detail::NullL>,
1306 : // for the detail function called at the end
1307 : gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>,
1308 : ::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>,
1309 : gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>,
1310 : Tags::detail::AngularDNullL,
1311 : Tags::detail::DLambda<
1312 : gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>,
1313 : Tags::detail::DLambda<
1314 : gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>>,
1315 : ::Tags::spacetime_deriv<Tags::detail::RealBondiR, tmpl::size_t<3>,
1316 : Frame::RadialNull>,
1317 : Tags::detail::DLambda<Tags::detail::DLambda<Tags::detail::RealBondiR>>,
1318 : ::Tags::deriv<Tags::detail::DLambda<Tags::detail::RealBondiR>,
1319 : tmpl::size_t<2>, Frame::RadialNull>>>
1320 : computation_variables{size};
1321 :
1322 : Variables<
1323 : tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
1324 : std::integral_constant<int, 0>>,
1325 : ::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
1326 : std::integral_constant<int, 1>>>>
1327 : derivative_buffers{size};
1328 :
1329 : auto& cos_phi = get<Tags::detail::CosPhi>(computation_variables);
1330 : auto& cos_theta = get<Tags::detail::CosTheta>(computation_variables);
1331 : auto& sin_phi = get<Tags::detail::SinPhi>(computation_variables);
1332 : auto& sin_theta = get<Tags::detail::SinTheta>(computation_variables);
1333 : trigonometric_functions_on_swsh_collocation(
1334 : make_not_null(&cos_phi), make_not_null(&cos_theta),
1335 : make_not_null(&sin_phi), make_not_null(&sin_theta), l_max);
1336 :
1337 : // NOTE: to handle the singular values of polar coordinates, the phi
1338 : // components of all tensors are scaled according to their sin(theta)
1339 : // prefactors.
1340 : // so, any down-index component get<2>(A) represents 1/sin(theta) A_\phi,
1341 : // and any up-index component get<2>(A) represents sin(theta) A^\phi.
1342 : // This holds for Jacobians, and so direct application of the Jacobians
1343 : // brings the factors through.
1344 : auto& unused_cartesian_coords =
1345 : get<Tags::detail::CartesianCoordinates>(computation_variables);
1346 : auto& inverse_cartesian_to_spherical_jacobian =
1347 : get<Tags::detail::InverseCartesianToSphericalJacobian>(
1348 : computation_variables);
1349 : auto& cartesian_to_spherical_jacobian =
1350 : get<Tags::detail::CartesianToSphericalJacobian>(computation_variables);
1351 : cartesian_to_spherical_coordinates_and_jacobians(
1352 : make_not_null(&unused_cartesian_coords),
1353 : make_not_null(&cartesian_to_spherical_jacobian),
1354 : make_not_null(&inverse_cartesian_to_spherical_jacobian), cos_phi,
1355 : cos_theta, sin_phi, sin_theta, extraction_radius);
1356 :
1357 : auto& phi = get<gh::Tags::Phi<DataVector, 3>>(computation_variables);
1358 : auto& dt_spacetime_metric =
1359 : get<::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3>>>(
1360 : computation_variables);
1361 : auto& spacetime_metric =
1362 : get<gr::Tags::SpacetimeMetric<DataVector, 3>>(computation_variables);
1363 :
1364 : auto& interpolation_buffer =
1365 : get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
1366 : std::integral_constant<int, 0>>>(
1367 : derivative_buffers);
1368 : auto& eth_buffer =
1369 : get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
1370 : std::integral_constant<int, 1>>>(
1371 : derivative_buffers);
1372 :
1373 : auto& d_cartesian_spatial_metric =
1374 : get<::Tags::deriv<gr::Tags::SpatialMetric<DataVector, 3>, tmpl::size_t<3>,
1375 : ::Frame::Inertial>>(computation_variables);
1376 : auto& d_cartesian_shift =
1377 : get<::Tags::deriv<gr::Tags::Shift<DataVector, 3>, tmpl::size_t<3>,
1378 : ::Frame::Inertial>>(computation_variables);
1379 : auto& d_cartesian_lapse =
1380 : get<::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
1381 : ::Frame::Inertial>>(computation_variables);
1382 :
1383 : deriv_cartesian_metric_lapse_shift_from_nodes(
1384 : make_not_null(&d_cartesian_spatial_metric),
1385 : make_not_null(&d_cartesian_shift), make_not_null(&d_cartesian_lapse),
1386 : make_not_null(&interpolation_buffer), make_not_null(ð_buffer),
1387 : cartesian_spatial_metric, cartesian_dr_spatial_metric, cartesian_shift,
1388 : cartesian_dr_shift, cartesian_lapse, cartesian_dr_lapse,
1389 : inverse_cartesian_to_spherical_jacobian, l_max);
1390 :
1391 : gh::phi(make_not_null(&phi), cartesian_lapse, d_cartesian_lapse,
1392 : cartesian_shift, d_cartesian_shift, cartesian_spatial_metric,
1393 : d_cartesian_spatial_metric);
1394 : gr::time_derivative_of_spacetime_metric(
1395 : make_not_null(&dt_spacetime_metric), cartesian_lapse, cartesian_dt_lapse,
1396 : cartesian_shift, cartesian_dt_shift, cartesian_spatial_metric,
1397 : cartesian_dt_spatial_metric);
1398 : gr::spacetime_metric(make_not_null(&spacetime_metric), cartesian_lapse,
1399 : cartesian_shift, cartesian_spatial_metric);
1400 : auto& inverse_spatial_metric =
1401 : get<gr::Tags::InverseSpatialMetric<DataVector, 3>>(computation_variables);
1402 : inverse_spatial_metric =
1403 : determinant_and_inverse(cartesian_spatial_metric).second;
1404 :
1405 : auto& dt_worldtube_normal =
1406 : get<::Tags::dt<Tags::detail::WorldtubeNormal>>(computation_variables);
1407 : auto& worldtube_normal =
1408 : get<Tags::detail::WorldtubeNormal>(computation_variables);
1409 : worldtube_normal_and_derivatives(
1410 : make_not_null(&worldtube_normal), make_not_null(&dt_worldtube_normal),
1411 : cos_phi, cos_theta, spacetime_metric, dt_spacetime_metric, sin_phi,
1412 : sin_theta, inverse_spatial_metric);
1413 :
1414 : auto& du_null_l = get<::Tags::dt<Tags::detail::NullL>>(computation_variables);
1415 : auto& null_l = get<Tags::detail::NullL>(computation_variables);
1416 : null_vector_l_and_derivatives(
1417 : make_not_null(&du_null_l), make_not_null(&null_l), dt_worldtube_normal,
1418 : cartesian_dt_lapse, dt_spacetime_metric, cartesian_dt_shift,
1419 : cartesian_lapse, spacetime_metric, cartesian_shift, worldtube_normal);
1420 :
1421 : // pass to the next step that is common between the 'modal' input and 'GH'
1422 : // input strategies
1423 : detail::create_bondi_boundary_data(
1424 : bondi_boundary_data, make_not_null(&computation_variables),
1425 : make_not_null(&derivative_buffers), dt_spacetime_metric, phi,
1426 : spacetime_metric, null_l, du_null_l, cartesian_to_spherical_jacobian,
1427 : l_max, extraction_radius);
1428 : }
1429 : } // namespace Cce
|
