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/Prefixes.hpp"
9 : #include "DataStructures/DataVector.hpp"
10 : #include "DataStructures/Tensor/Tensor.hpp"
11 : #include "DataStructures/Variables.hpp"
12 : #include "Evolution/Systems/CurvedScalarWave/Tags.hpp"
13 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
14 : #include "Options/Context.hpp"
15 : #include "Options/Options.hpp"
16 : #include "Options/String.hpp"
17 : #include "Utilities/Gsl.hpp"
18 :
19 : /// \cond
20 : namespace PUP {
21 : class er;
22 : } // namespace PUP
23 : /// \endcond
24 :
25 : namespace CurvedScalarWave::Worldtube {
26 : /*!
27 : * \brief Dispatcher to compute the puncture/singular field for a scalar charge
28 : * on a generic orbit in Schwarzschild or Kerr spacetime. The Kerr puncture
29 : * reduces to Schwarzsschild for zero spin but is faster to evaluate.
30 : */
31 1 : class PunctureField {
32 : public:
33 0 : enum class Type { Schwarzschild, Kerr };
34 :
35 : /*!
36 : * \brief Use the Schwarzschild puncture field model.
37 : */
38 1 : struct Schwarzschild {
39 0 : using type = Schwarzschild;
40 0 : static constexpr Options::String help = {
41 : "Use the Schwarzschild puncture field model."};
42 :
43 : /*!
44 : * \brief Puncture field expansion order. Currently orders 0 and 1 are
45 : * implemented.
46 : */
47 1 : struct ExpansionOrder {
48 0 : using type = size_t;
49 0 : static constexpr Options::String help{
50 : "Puncture field expansion order. Currently orders 0 and 1 are "
51 : "implemented."};
52 0 : static size_t upper_bound() { return 1; }
53 : };
54 :
55 : /*!
56 : * \brief The mass of the central black hole.
57 : */
58 1 : struct BlackHoleMass {
59 0 : using type = double;
60 0 : static constexpr Options::String help{
61 : "The mass of the central black hole."};
62 0 : static double lower_bound() { return 0.; }
63 : };
64 :
65 0 : using options = tmpl::list<ExpansionOrder, BlackHoleMass>;
66 :
67 0 : Schwarzschild() = default;
68 0 : Schwarzschild(size_t expansion_order_in, double black_hole_mass_in,
69 : const Options::Context& context = {});
70 :
71 0 : size_t expansion_order{};
72 0 : double black_hole_mass{};
73 : };
74 :
75 : /*!
76 : * \brief Use the Kerr puncture field model. This option is currently parsed
77 : * but not yet implemented at runtime.
78 : */
79 1 : struct Kerr {
80 0 : using type = Kerr;
81 0 : static constexpr Options::String help = {
82 : "Use the Kerr puncture field model. This option is currently parsed "
83 : "but not yet implemented at runtime."};
84 :
85 : /*!
86 : * \brief Puncture field expansion order. Currently orders 0 and 1 are
87 : * implemented.
88 : */
89 1 : struct ExpansionOrder {
90 0 : using type = size_t;
91 0 : static constexpr Options::String help{
92 : "Puncture field expansion order. Currently orders 0 and 1 are "
93 : "implemented."};
94 0 : static size_t upper_bound() { return 1; }
95 : };
96 :
97 : /*!
98 : * \brief The mass of the central black hole.
99 : */
100 1 : struct BlackHoleMass {
101 0 : using type = double;
102 0 : static constexpr Options::String help{
103 : "The mass of the central black hole."};
104 0 : static double lower_bound() { return 0.; }
105 : };
106 :
107 : /*!
108 : * \brief The dimensionless z-component of the black-hole spin.
109 : */
110 1 : struct SpinAlongZAxis {
111 0 : using type = double;
112 0 : static constexpr Options::String help{
113 : "The dimensionless z-component of the black-hole spin."};
114 : };
115 :
116 0 : using options = tmpl::list<ExpansionOrder, BlackHoleMass, SpinAlongZAxis>;
117 :
118 0 : Kerr() = default;
119 0 : Kerr(size_t expansion_order_in, double black_hole_mass_in,
120 : double spin_along_z_axis_in, const Options::Context& context = {});
121 :
122 0 : size_t expansion_order{};
123 0 : double black_hole_mass{};
124 0 : double spin_along_z_axis{};
125 : };
126 :
127 0 : using options = tmpl::list<
128 : Options::Alternatives<tmpl::list<Schwarzschild>, tmpl::list<Kerr>>>;
129 :
130 0 : static constexpr Options::String help = {
131 : "Configuration and dispatcher for puncture-field expressions. Choose "
132 : "either Schwarzschild or Kerr."};
133 :
134 0 : PunctureField() = default;
135 0 : explicit PunctureField(const Schwarzschild& schwarzschild,
136 : const Options::Context& context = {});
137 0 : explicit PunctureField(const Kerr& kerr,
138 : const Options::Context& context = {});
139 :
140 0 : void pup(PUP::er& p);
141 :
142 0 : Type type() const;
143 0 : size_t expansion_order() const;
144 0 : double black_hole_mass() const;
145 0 : double spin_along_z_axis() const;
146 :
147 : /*!
148 : * \brief Compute and write the puncture field and its derivatives for the
149 : * configured puncture model and expansion order.
150 : */
151 1 : void apply_puncture(
152 : gsl::not_null<Variables<tmpl::list<
153 : CurvedScalarWave::Tags::Psi, ::Tags::dt<CurvedScalarWave::Tags::Psi>,
154 : ::Tags::deriv<CurvedScalarWave::Tags::Psi, tmpl::size_t<3>,
155 : Frame::Inertial>>>*>
156 : result,
157 : const tnsr::I<DataVector, 3, Frame::Inertial>& centered_coords,
158 : const tnsr::I<double, 3>& particle_position,
159 : const tnsr::I<double, 3>& particle_velocity,
160 : const tnsr::I<double, 3>& particle_acceleration) const;
161 :
162 : /*!
163 : * \brief Compute and write the corrections to the
164 : * puncture field for the configured puncture model and expansion order. These
165 : * terms arise at non-geodesic accelerations such as the self-force.
166 : */
167 1 : void apply_acceleration_terms(
168 : gsl::not_null<Variables<tmpl::list<
169 : CurvedScalarWave::Tags::Psi, ::Tags::dt<CurvedScalarWave::Tags::Psi>,
170 : ::Tags::deriv<CurvedScalarWave::Tags::Psi, tmpl::size_t<3>,
171 : Frame::Inertial>>>*>
172 : result,
173 : const tnsr::I<DataVector, 3, Frame::Inertial>& centered_coords,
174 : const tnsr::I<double, 3>& particle_position,
175 : const tnsr::I<double, 3>& particle_velocity,
176 : const tnsr::I<double, 3>& particle_acceleration, double ft, double fx,
177 : double fy, double dt_ft, double dt_fx, double dt_fy, double Du_ft,
178 : double Du_fx, double Du_fy, double dt_Du_ft, double dt_Du_fx,
179 : double dt_Du_fy) const;
180 :
181 : private:
182 0 : Type type_{Type::Schwarzschild};
183 0 : size_t expansion_order_{0};
184 0 : double black_hole_mass_{1.};
185 0 : double spin_along_z_axis_{0.};
186 : };
187 : /*!
188 : * \brief Computes the puncture/singular field \f$\Psi^\mathcal{P}\f$ of a
189 : * scalar charge on a generic orbit in Schwarzschild spacetime.
190 : * described in \cite Detweiler2003.
191 : *
192 : * \details The field is computed using a Detweiler-Whiting singular
193 : * Green's function and perturbatively expanded in the geodesic distance from
194 : * the particle. It solves the inhomogeneous wave equation
195 : *
196 : * \f{align*}{
197 : * \Box \Psi^\mathcal{P} = -4 \pi q \int \sqrt{-g} \delta^4(x^i, z(\tau)) d \tau
198 : * \f}
199 : *
200 : * where \f$q\f$ is the scalar charge and \f$z(\tau)\f$ is the worldline of the
201 : * particle. The expression is expanded up to a certain order in geodesic
202 : * distance and transformed to Kerr-Schild coordinates.
203 : *
204 : * The function given here assumes that the particle has scalar charge \f$q=1\f$
205 : * and is on a fixed geodesic orbit. It returns the
206 : * singular field at the requested coordinates as well as its time and spatial
207 : * derivative. For non-geodesic orbits, corresponding acceleration terms have to
208 : * be added to the puncture field.
209 : *
210 : * \note The expressions were computed with Mathematica and optimized by
211 : * applying common subexpression elimination with sympy. The memory allocations
212 : * of temporaries were optimized manually.
213 : */
214 1 : void puncture_field_0(
215 : gsl::not_null<Variables<tmpl::list<
216 : CurvedScalarWave::Tags::Psi, ::Tags::dt<CurvedScalarWave::Tags::Psi>,
217 : ::Tags::deriv<CurvedScalarWave::Tags::Psi, tmpl::size_t<3>,
218 : Frame::Inertial>>>*>
219 : result,
220 : const tnsr::I<DataVector, 3, Frame::Inertial>& centered_coords,
221 : const tnsr::I<double, 3>& particle_position,
222 : const tnsr::I<double, 3>& particle_velocity,
223 : const tnsr::I<double, 3>& particle_acceleration, double bh_mass);
224 :
225 : /*!
226 : * \brief Computes the puncture/singular field \f$\Psi^\mathcal{P}\f$ of a
227 : * scalar charge on a generic orbit in Schwarzschild spacetime.
228 : * described in \cite Detweiler2003.
229 : *
230 : * \details For non-geodesic orbits, there are additional contributions, see
231 : * `acceleration_terms_0`.
232 : */
233 1 : void puncture_field_1(
234 : gsl::not_null<Variables<tmpl::list<
235 : CurvedScalarWave::Tags::Psi, ::Tags::dt<CurvedScalarWave::Tags::Psi>,
236 : ::Tags::deriv<CurvedScalarWave::Tags::Psi, tmpl::size_t<3>,
237 : Frame::Inertial>>>*>
238 : result,
239 : const tnsr::I<DataVector, 3, Frame::Inertial>& centered_coords,
240 : const tnsr::I<double, 3>& particle_position,
241 : const tnsr::I<double, 3>& particle_velocity,
242 : const tnsr::I<double, 3>& particle_acceleration, double bh_mass);
243 :
244 : /*!
245 : * \brief Computes the acceleration terms of a puncture/singular field
246 : * \f$\Psi^\mathcal{P}\f$ of a scalar charge on a generic orbit in Schwarzschild
247 : * spacetime up to zeroth order in coordinate distance.
248 : * \details The appropriate expression can be found in Eq. (37) of
249 : * \cite Wittek:2024gxn. The values ft, fx, fy are the time, x and y component
250 : * of the self force per unit mass evaluated at the position of the particle;
251 : * dt_ft, dt_fx, dt_fy are the respective total time derivatives. The code in
252 : * this function was auto-generated by generating the full expressions with
253 : * Mathematica and employing common subexpression elimination with sympy. The
254 : * mathematica file and generating script can be found at
255 : * https://github.com/nikwit/puncture-field.
256 : */
257 1 : void acceleration_terms_0(
258 : gsl::not_null<Variables<tmpl::list<
259 : CurvedScalarWave::Tags::Psi, ::Tags::dt<CurvedScalarWave::Tags::Psi>,
260 : ::Tags::deriv<CurvedScalarWave::Tags::Psi, tmpl::size_t<3>,
261 : Frame::Inertial>>>*>
262 : result,
263 : const tnsr::I<DataVector, 3, Frame::Inertial>& centered_coords,
264 : const tnsr::I<double, 3>& particle_position,
265 : const tnsr::I<double, 3>& particle_velocity,
266 : const tnsr::I<double, 3>& particle_acceleration, double ft, double fx,
267 : double fy, double dt_ft, double dt_fx, double dt_fy, double bh_mass);
268 :
269 : /*!
270 : * \brief Computes the acceleration terms of a puncture/singular field
271 : * \f$\Psi^\mathcal{P}\f$ of a scalar charge on a generic orbit in Schwarzschild
272 : * spacetime up to first order in coordinate distance (i.e. zeroth and first
273 : * order).
274 : * \details The appropriate expression can be found in Eq. (37) of
275 : * \cite Wittek:2024gxn. The values ft, fx, fy are the time, x and y component
276 : * of the self force per unit mass evaluated at the position of the particle;
277 : * dt_ft, dt_fx, dt_fy are the respective total time derivatives. The code in
278 : * this function was auto-generated by generating the full expressions with
279 : * Mathematica and employing common subexpression elimination with sympy. The
280 : * mathematica file and generating script can be found at
281 : * https://github.com/nikwit/puncture-field.
282 : *
283 : */
284 1 : void acceleration_terms_1(
285 : gsl::not_null<Variables<tmpl::list<
286 : CurvedScalarWave::Tags::Psi, ::Tags::dt<CurvedScalarWave::Tags::Psi>,
287 : ::Tags::deriv<CurvedScalarWave::Tags::Psi, tmpl::size_t<3>,
288 : Frame::Inertial>>>*>
289 : result,
290 : const tnsr::I<DataVector, 3, Frame::Inertial>& centered_coords,
291 : const tnsr::I<double, 3>& particle_position,
292 : const tnsr::I<double, 3>& particle_velocity,
293 : const tnsr::I<double, 3>& particle_acceleration, double ft, double fx,
294 : double fy, double dt_ft, double dt_fx, double dt_fy, double Du_ft,
295 : double Du_fx, double Du_fy, double dt_Du_ft, double dt_Du_fx,
296 : double dt_Du_fy, double bh_mass);
297 : } // namespace CurvedScalarWave::Worldtube
|
