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/SpinWeighted.hpp"
9 : #include "DataStructures/Tensor/TypeAliases.hpp"
10 : #include "Utilities/Gsl.hpp"
11 :
12 : /// \cond
13 : class ComplexDataVector;
14 : /// \endcond
15 :
16 : namespace Cce::InitializeJ::CauchySecondOrder_detail {
17 :
18 : /*!
19 : * \brief Solve the H-hypersurface equation at the worldtube for
20 : * \f$\partial_y^2 J\f$.
21 : *
22 : * \details Evaluating the H-hypersurface equation at \f$y=-1\f$ and isolating
23 : * the contributions coming from the swsh Jacobians at the worldtube produces
24 : * the complex identity
25 : *
26 : * \f[
27 : * c_1 \, \partial_y^2 J + c_2 \, \partial_y^2 \bar J = -c_3,
28 : * \f]
29 : *
30 : * where \f$c_1\f$, \f$c_2\f$, and \f$c_3\f$ depend only on quantities known at
31 : * the worldtube. The conjugate of this identity eliminates
32 : * \f$\partial_y^2 \bar J\f$ and yields \f$\partial_y^2 J\f$.
33 : *
34 : * The whole computation reuses the volume CCE machinery rather than any
35 : * hand-written angular-Jacobian expressions. The worldtube right-minus-left
36 : * side of the H-hypersurface equation is assembled as a function
37 : * \f$F(\partial_y^2 J)\f$: the angular derivatives come from
38 : * `Spectral::Swsh::angular_derivatives` and are converted from the numerical
39 : * (constant \f$y\f$) to the physical (constant \f$r\f$) coordinate with
40 : * `Cce::ApplySwshJacobianInplace`, and the hypersurface right-hand sides use
41 : * the same `Cce::ComputeBondiIntegrand` specializations as the evolution, with
42 : * \f$\partial_y^2 J\f$ threaded through every dependence. Because \f$F\f$ is
43 : * affine in \f$(\partial_y^2 J, \partial_y^2 \bar J)\f$, evaluating it at
44 : * \f$\partial_y^2 J = 0, 1, i\f$ fixes \f$c_3 = F(0)\f$ and the linear
45 : * coefficients \f$c_1\f$, \f$c_2\f$ exactly.
46 : *
47 : * The caller supplies the physical worldtube data; `compute_dy_dy_j` converts
48 : * it to the numerical (constant \f$y\f$) coordinate that
49 : * `evaluate_worldtube_h_residual` works in, using the worldtube Jacobian
50 : * \f$\partial_y J = (R / 2) \partial_r J\f$,
51 : *
52 : * \f{align*}{
53 : * \partial_y J &= \tfrac{1}{2} R \, \partial_r J, \\
54 : * \breve{H} &= H + \partial_u R \, \partial_r J, \\
55 : * \partial_y \breve{H} &= \tfrac{1}{2}\left(\partial_u R \, \partial_r J
56 : * + R \, \partial_{\breve u} \partial_r J\right),
57 : * \f}
58 : *
59 : * where \f$\breve{H} = \partial_{\breve u} J = (\partial_u J)_y\f$ is the
60 : * numerical-coordinate \f$H\f$. The time derivative
61 : * enters only through the primitive radial quantity
62 : * `du_dr_j` \f$= \partial_{\breve u} \partial_r J\f$.
63 : *
64 : * \note The coordinate held fixed by a time derivative is written two
65 : * equivalent ways: a breve accent marks "at constant numerical coordinate
66 : * \f$y\f$" (following the worldtube), while an unaccented symbol means "at
67 : * constant Bondi \f$r\f$"; `BoundaryData.hpp` writes the same distinction with
68 : * an explicit \f$(\,\cdot\,)_y\f$ / \f$(\,\cdot\,)_r\f$ subscript. So, for
69 : * \f$H\f$,
70 : * - \f$\breve{H} = (\partial_u J)_y\f$ (constant \f$y\f$, "numerical")
71 : * \f$\;\leftrightarrow\;\f$ `Cce::Tags::BondiH` (`= ::Tags::dt<BondiJ>`),
72 : * - \f$H = (\partial_u J)_r\f$ (constant Bondi \f$r\f$)
73 : * \f$\;\leftrightarrow\;\f$ `Cce::Tags::Du<BondiJ>`.
74 : *
75 : * The tag \f$\breve{H}\f$ drops the accent because the evolution only ever uses
76 : * the constant-\f$y\f$ \f$H\f$, so there is nothing to distinguish it from.
77 : *
78 : * \see evaluate_worldtube_h_residual for the function \f$F\f$ that is probed.
79 : */
80 : void compute_dy_dy_j(
81 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> dy_dy_j,
82 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& j,
83 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& u,
84 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& w,
85 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& beta,
86 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& q,
87 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& du_j,
88 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& dr_j,
89 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& du_dr_j,
90 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& du_r,
91 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& r, size_t l_max);
92 :
93 : /*!
94 : * \brief The worldtube residual \f$F(\partial_y^2 J)\f$ of the H-hypersurface
95 : * equation whose root `compute_dy_dy_j` returns.
96 : *
97 : * \details Returns the worldtube (\f$y = -1\f$) right-minus-left side of the
98 : * H-hypersurface equation,
99 : *
100 : * \f[
101 : * F(\partial_y^2 J) = \mathrm{rhs}(H) - \mathrm{lhs}(H),
102 : * \f]
103 : *
104 : * evaluated with the supplied trial value `dy_dy_j_value` for
105 : * \f$\partial_y^2 J\f$ (with \f$\partial_y^2 \bar J\f$ taken to be its complex
106 : * conjugate). All worldtube inputs are supplied in the numerical (constant
107 : * \f$y\f$) coordinate: `dy_j` \f$= \partial_y J\f$, `h`
108 : * \f$= \breve{H} = (\partial_u J)_y\f$, and `dy_h`
109 : * \f$= \partial_y \breve{H}\f$; no physical radial-derivative quantity is
110 : * passed. Every angular derivative is converted from the numerical to the
111 : * physical coordinate with `Cce::ApplySwshJacobianInplace`, and every
112 : * hypersurface right-hand side comes from the `Cce::ComputeBondiIntegrand`
113 : * specializations, so the residual is assembled entirely from the same volume
114 : * machinery the evolution uses. `compute_dy_dy_j` returns the value of
115 : * \f$\partial_y^2 J\f$ for which this residual vanishes, so feeding that value
116 : * back in reproduces the H-hypersurface equation to machine precision.
117 : */
118 : Scalar<SpinWeighted<ComplexDataVector, 2>> evaluate_worldtube_h_residual(
119 : const ComplexDataVector& dy_dy_j_value,
120 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& j,
121 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& u,
122 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& w,
123 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& beta,
124 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& q,
125 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& dy_j,
126 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& h,
127 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& dy_h,
128 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& du_r,
129 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& r, size_t l_max);
130 :
131 : } // namespace Cce::InitializeJ::CauchySecondOrder_detail
|
