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