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 : #include <memory>
8 : #include <optional>
9 : #include <string>
10 : #include <type_traits>
11 : #include <unordered_map>
12 :
13 : #include "DataStructures/DataBox/Tag.hpp"
14 : #include "DataStructures/TaggedTuple.hpp"
15 : #include "DataStructures/Tensor/Tensor.hpp"
16 : #include "DataStructures/Variables.hpp"
17 : #include "Domain/Creators/DomainCreator.hpp"
18 : #include "Domain/Creators/OptionTags.hpp"
19 : #include "Domain/Domain.hpp"
20 : #include "Domain/ExcisionSphere.hpp"
21 : #include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp"
22 : #include "Domain/FunctionsOfTime/Tags.hpp"
23 : #include "Domain/Structure/Element.hpp"
24 : #include "Domain/Structure/ElementId.hpp"
25 : #include "Domain/Tags.hpp"
26 : #include "Evolution/Systems/CurvedScalarWave/BackgroundSpacetime.hpp"
27 : #include "Evolution/Systems/CurvedScalarWave/Tags.hpp"
28 : #include "Evolution/Systems/CurvedScalarWave/Worldtube/PunctureField.hpp"
29 : #include "IO/Logging/Verbosity.hpp"
30 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
31 : #include "Options/Auto.hpp"
32 : #include "Options/String.hpp"
33 : #include "ParallelAlgorithms/EventsAndTriggers/Trigger.hpp"
34 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp"
35 : #include "PointwiseFunctions/GeneralRelativity/Surfaces/Tags.hpp"
36 : #include "Utilities/EqualWithinRoundoff.hpp"
37 : #include "Utilities/Gsl.hpp"
38 : #include "Utilities/Serialization/Serialize.hpp"
39 :
40 : /// \cond
41 : namespace Tags {
42 : struct Time;
43 : } // namespace Tags
44 : namespace OptionTags {
45 : struct InitialTime;
46 : } // namespace OptionTags
47 : /// \endcond
48 :
49 : namespace CurvedScalarWave::Worldtube {
50 : /*!
51 : * \brief Option tags for the worldtube
52 : */
53 1 : namespace OptionTags {
54 : /*!
55 : * \brief Options for the worldtube
56 : */
57 1 : struct Worldtube {
58 0 : static constexpr Options::String help = {"Options for the Worldtube"};
59 : };
60 :
61 : /*!
62 : * \brief The value of the scalar charge in units of the black hole mass M.
63 : */
64 1 : struct Charge {
65 0 : using type = double;
66 0 : static constexpr Options::String help{
67 : "The value of the scalar charge in units of the black hole mass M."};
68 0 : using group = Worldtube;
69 : };
70 :
71 : /*!
72 : * \brief The dimensionless spin vector of the central black hole.
73 : */
74 1 : struct Spin {
75 0 : using type = std::array<double, 3>;
76 0 : static constexpr Options::String help{
77 : "The value of the scalar charge in units of the black hole mass M."};
78 0 : using group = Worldtube;
79 : };
80 :
81 : /*!
82 : * \brief Options for the scalar self-force. Select `None` for a purely geodesic
83 : * evolution
84 : *
85 : *\details The self force is turned on using the smooth transition function
86 : *
87 : * \begin{equation}
88 : * w(t) = 1 - \exp{ \left(- \left(\frac{t - t_1}{\sigma} \right)^4 \right)}.
89 : * \end{equation}
90 : *
91 : * The turn on time is given by \f$t_1\f$ and the turn on interval is given by
92 : *\f$sigma\f$.
93 : */
94 1 : struct SelfForceOptions {
95 0 : static constexpr Options::String help = {
96 : "Options for the scalar self-force. Select `None` for a purely geodesic "
97 : "evolution"};
98 0 : using group = Worldtube;
99 0 : using type = Options::Auto<SelfForceOptions, Options::AutoLabel::None>;
100 :
101 0 : struct Mass {
102 0 : using type = double;
103 0 : static constexpr Options::String help{
104 : "The mass of the scalar particle in units of the black hole mass M."};
105 0 : static double lower_bound() { return 0.; }
106 : };
107 :
108 0 : struct Iterations {
109 0 : using type = size_t;
110 0 : static constexpr Options::String help{
111 : "The number of iterations used to compute the particle acceleration. "
112 : "Must be at least 1 as 0 iterations corresponds to the geodesic "
113 : "acceleration."};
114 0 : static size_t lower_bound() { return 1; }
115 : };
116 :
117 0 : struct TurnOnTime {
118 0 : using type = double;
119 0 : static constexpr Options::String help{
120 : "The time at which the scalar self force is turned on."};
121 0 : static double lower_bound() { return 0.; }
122 : };
123 :
124 0 : struct TurnOnInterval {
125 0 : using type = double;
126 0 : static constexpr Options::String help{
127 : "The interval over which the scalar self force is smoothly turned on. "
128 : "We require a minimum of 1 M for the interval."};
129 0 : static double lower_bound() { return 1.; }
130 : };
131 :
132 0 : SelfForceOptions();
133 0 : SelfForceOptions(double mass_in, size_t iterations_in, double turn_on_time_in,
134 : double turn_on_interval_in);
135 0 : void pup(PUP::er& p);
136 :
137 0 : using options = tmpl::list<Mass, Iterations, TurnOnTime, TurnOnInterval>;
138 :
139 0 : double mass{};
140 0 : size_t iterations{};
141 0 : double turn_on_time{};
142 0 : double turn_on_interval{};
143 : };
144 :
145 : /*!
146 : * \brief Options for the excision sphere radii which are adjusted according to
147 : * `smooth_broken_power_law`. If `IsWorldtube` is true, these options control
148 : * the worldtube growth around the scalar charge. Else, they control the growth
149 : * of the excision sphere within the central black hole.
150 : */
151 : template <bool IsWorldtube>
152 1 : struct RadiusOptions {
153 0 : static constexpr Options::String help = {
154 : "Options for the radii of the excision spheres"};
155 0 : using group = Worldtube;
156 0 : using type = RadiusOptions;
157 :
158 0 : static std::string name() {
159 : if constexpr (IsWorldtube) {
160 : return "WorldtubeRadiusOptions";
161 : } else {
162 : return "BlackHoleRadiusOptions";
163 : }
164 : }
165 :
166 0 : struct Exponent {
167 0 : using type = double;
168 0 : static constexpr Options::String help{
169 : "The exponent alpha according to which the excision sphere grows with "
170 : "orbital radius until the transition radius."};
171 0 : static double lower_bound() { return 0.; }
172 0 : static double upper_bound() { return 4.; }
173 : };
174 :
175 0 : struct Amplitude {
176 0 : using type = double;
177 0 : static constexpr Options::String help{
178 : "The amplitude A of the smoothly broken power law."};
179 0 : static double lower_bound() { return 0.; }
180 : };
181 :
182 0 : struct TransitionRadius {
183 0 : using type = double;
184 0 : static constexpr Options::String help{
185 : "The transition radius rb of the smoothly broken power law. At this "
186 : "point the radius transitions to a constant value."};
187 0 : static double lower_bound() { return 0.; }
188 : };
189 :
190 0 : struct TransitionWidth {
191 0 : using type = double;
192 0 : static constexpr Options::String help{
193 : "The width delta of the transition region."};
194 0 : static double lower_bound() { return 1e-3; }
195 : };
196 :
197 0 : RadiusOptions();
198 0 : RadiusOptions(double exponent_in, double amplitude_in,
199 : double transition_radius_in, double transition_width_in);
200 0 : void pup(PUP::er& p);
201 :
202 0 : using options =
203 : tmpl::list<Exponent, Amplitude, TransitionRadius, TransitionWidth>;
204 :
205 0 : double exponent{};
206 0 : double amplitude{};
207 0 : double transition_radius{};
208 0 : double transition_width{};
209 : };
210 :
211 : /*!
212 : * \brief Name of the excision sphere designated to act as a worldtube
213 : */
214 1 : struct ExcisionSphere {
215 0 : using type = std::string;
216 0 : static constexpr Options::String help{
217 : "The name of the excision sphere as returned by the domain."};
218 0 : using group = Worldtube;
219 : };
220 :
221 : /*!
222 : * \brief Triggers at which to write the coefficients of the worldtube's
223 : * internal Taylor series to file.
224 : */
225 1 : struct ObserveCoefficientsTrigger {
226 0 : using type = std::unique_ptr<Trigger>;
227 0 : static constexpr Options::String help{
228 : "Specifies a non-dense trigger in which the coefficients of the internal "
229 : "regular field expansion are written to file."};
230 0 : using group = Worldtube;
231 : };
232 :
233 0 : struct PunctureField {
234 0 : using type = CurvedScalarWave::Worldtube::PunctureField;
235 0 : static constexpr Options::String help{
236 : "Options controlling puncture-field evaluation."};
237 0 : using group = Worldtube;
238 : };
239 :
240 : /*!
241 : * \brief The verbosity of the worldtube executable.
242 : */
243 1 : struct Verbosity {
244 0 : using type = ::Verbosity;
245 0 : static constexpr Options::String help{
246 : "Gives the verbosity of the worldtube executable."};
247 0 : using group = Worldtube;
248 : };
249 : } // namespace OptionTags
250 :
251 : /*!
252 : * \brief Tags related to the worldtube
253 : */
254 : namespace Tags {
255 : /*!
256 : * \brief The excision sphere corresponding to the worldtube
257 : */
258 : template <size_t Dim>
259 1 : struct ExcisionSphere : db::SimpleTag {
260 0 : using type = ::ExcisionSphere<Dim>;
261 0 : using option_tags = tmpl::list<domain::OptionTags::DomainCreator<Dim>,
262 : OptionTags::ExcisionSphere>;
263 0 : static constexpr bool pass_metavariables = false;
264 0 : static ::ExcisionSphere<Dim> create_from_options(
265 : const std::unique_ptr<::DomainCreator<Dim>>& domain_creator,
266 : const std::string& excision_sphere) {
267 : const auto domain = domain_creator->create_domain();
268 : const auto& excision_spheres = domain.excision_spheres();
269 : if (excision_spheres.count(excision_sphere) == 0) {
270 : ERROR("Specified excision sphere '"
271 : << excision_sphere
272 : << "' not available. Available excision spheres are: "
273 : << keys_of(excision_spheres));
274 : }
275 : return excision_spheres.at(excision_sphere);
276 : }
277 : };
278 :
279 : /*!
280 : * \brief Triggers at which to write the coefficients of the worldtube's
281 : * internal Taylor series to file.
282 : */
283 1 : struct ObserveCoefficientsTrigger : db::SimpleTag {
284 0 : using type = std::unique_ptr<Trigger>;
285 0 : using option_tags = tmpl::list<OptionTags::ObserveCoefficientsTrigger>;
286 0 : static constexpr bool pass_metavariables = false;
287 0 : static std::unique_ptr<Trigger> create_from_options(
288 : const std::unique_ptr<Trigger>& trigger) {
289 : return serialize_and_deserialize<type>(trigger);
290 : }
291 : };
292 :
293 : /*!
294 : * \brief The value of the scalar charge
295 : */
296 1 : struct Charge : db::SimpleTag {
297 0 : using type = double;
298 0 : using option_tags = tmpl::list<OptionTags::Charge>;
299 0 : static constexpr bool pass_metavariables = false;
300 0 : static double create_from_options(const double charge) { return charge; };
301 : };
302 :
303 : /*!
304 : * \brief The dimensionless spin vector of the central black hole.
305 : */
306 1 : struct Spin : db::SimpleTag {
307 0 : using type = std::array<double, 3>;
308 0 : using option_tags = tmpl::list<OptionTags::Spin>;
309 0 : static constexpr bool pass_metavariables = false;
310 0 : static std::array<double, 3> create_from_options(
311 : const std::array<double, 3> spin) {
312 : return spin;
313 : };
314 : };
315 :
316 : /*!
317 : * \brief The time at which the self-force is smoothly turned on.
318 : *
319 : * \details The self force is turned on using the smooth transition function
320 : *
321 : * \begin{equation}
322 : * w(t) = 1 - \exp{ \left(- \left(\frac{t - t_1}{\sigma} \right)^4 \right)}.
323 : * \end{equation}
324 : *
325 : * The turn on time is given by \f$t_1\f$.
326 : */
327 1 : struct SelfForceTurnOnTime : db::SimpleTag {
328 0 : using type = std::optional<double>;
329 0 : using option_tags = tmpl::list<OptionTags::SelfForceOptions>;
330 0 : static constexpr bool pass_metavariables = false;
331 0 : static std::optional<double> create_from_options(
332 : const std::optional<OptionTags::SelfForceOptions>& self_force_options) {
333 : return self_force_options.has_value()
334 : ? std::make_optional(self_force_options->turn_on_time)
335 : : std::nullopt;
336 : };
337 : };
338 :
339 : /*!
340 : * \brief The interval over which the self-force is smoothly turned on.
341 : *
342 : * \details The self force is turned on using the smooth transition function
343 : *
344 : * \begin{equation}
345 : * w(t) = 1 - \exp{ \left(- \left(\frac{t - t_1}{\sigma} \right)^4 \right)}.
346 : * \end{equation}
347 : *
348 : * The turn on interval is given by \f$\sigma\f$.
349 : */
350 1 : struct SelfForceTurnOnInterval : db::SimpleTag {
351 0 : using type = std::optional<double>;
352 0 : using option_tags = tmpl::list<OptionTags::SelfForceOptions>;
353 0 : static constexpr bool pass_metavariables = false;
354 0 : static std::optional<double> create_from_options(
355 : const std::optional<OptionTags::SelfForceOptions>& self_force_options) {
356 : return self_force_options.has_value()
357 : ? std::make_optional(self_force_options->turn_on_interval)
358 : : std::nullopt;
359 : };
360 : };
361 :
362 : /*!
363 : * \brief The mass of the scalar charge. Only has a value if the scalar self
364 : * force is applied.
365 : */
366 1 : struct Mass : db::SimpleTag {
367 0 : using type = std::optional<double>;
368 0 : using option_tags = tmpl::list<OptionTags::SelfForceOptions>;
369 0 : static constexpr bool pass_metavariables = false;
370 0 : static std::optional<double> create_from_options(
371 : const std::optional<OptionTags::SelfForceOptions>& self_force_options) {
372 : return self_force_options.has_value()
373 : ? std::make_optional(self_force_options->mass)
374 : : std::nullopt;
375 : }
376 : };
377 :
378 : /*!
379 : * \brief The maximum number of iterations that will be applied to the
380 : * acceleration of the particle.
381 : */
382 1 : struct MaxIterations : db::SimpleTag {
383 0 : using type = size_t;
384 0 : using option_tags = tmpl::list<OptionTags::SelfForceOptions>;
385 0 : static constexpr bool pass_metavariables = false;
386 0 : static size_t create_from_options(
387 : const std::optional<OptionTags::SelfForceOptions>& self_force_options) {
388 : return self_force_options.has_value() ? self_force_options->iterations : 0;
389 : }
390 : };
391 :
392 : /*!
393 : * \brief The verbosity of the worldtube executable.
394 : */
395 1 : struct Verbosity : db::SimpleTag {
396 0 : using type = ::Verbosity;
397 0 : using option_tags = tmpl::list<OptionTags::Verbosity>;
398 0 : static constexpr bool pass_metavariables = false;
399 0 : static ::Verbosity create_from_options(const ::Verbosity& verbosity) {
400 : return verbosity;
401 : }
402 : };
403 :
404 : /*!
405 : * \brief The current number of iterations that has been applied to the
406 : * acceleration of the particle.
407 : */
408 1 : struct CurrentIteration : db::SimpleTag {
409 0 : using type = size_t;
410 : };
411 :
412 : /*!
413 : * \brief The current expiration time of the functions of time which are
414 : * controlled by the worldtube singleton.
415 : */
416 1 : struct ExpirationTime : db::SimpleTag {
417 0 : using type = double;
418 : };
419 :
420 : /*!
421 : * \brief The current worldtube radius held by the singleton.
422 : */
423 1 : struct WorldtubeRadius : db::SimpleTag {
424 0 : using type = double;
425 : };
426 :
427 : /*!
428 : * \brief The initial position and velocity of the scalar charge in inertial
429 : * coordinates.
430 : */
431 1 : struct InitialPositionAndVelocity : db::SimpleTag {
432 0 : using type = std::array<tnsr::I<double, 3, Frame::Inertial>, 2>;
433 0 : using option_tags =
434 : tmpl::list<domain::OptionTags::DomainCreator<3>,
435 : OptionTags::ExcisionSphere, ::OptionTags::InitialTime>;
436 0 : static constexpr bool pass_metavariables = false;
437 0 : static type create_from_options(
438 : const std::unique_ptr<::DomainCreator<3>>& domain_creator,
439 : const std::string& excision_sphere_name, const double initial_time) {
440 : // only evaluated at initial time, so expiration times don't matter
441 : const auto initial_fot = domain_creator->functions_of_time();
442 : const auto domain = domain_creator->create_domain();
443 : const auto& excision_sphere =
444 : domain.excision_spheres().at(excision_sphere_name);
445 : ASSERT(excision_sphere.is_time_dependent(),
446 : "excision_sphere not time dependent");
447 : const auto& maps = excision_sphere.moving_mesh_grid_to_inertial_map();
448 : const auto mapped_tuple = maps.coords_frame_velocity_jacobians(
449 : excision_sphere.center(), initial_time, initial_fot);
450 : return {std::get<0>(mapped_tuple), std::get<3>(mapped_tuple)};
451 : }
452 : };
453 :
454 : /// @{
455 : /*!
456 : * \brief The position and velocity of the scalar charge particle orbiting a
457 : * central black hole given in inertial coordinates. This compute tag is meant
458 : * to be used by the elements.
459 : */
460 : template <size_t Dim>
461 1 : struct ParticlePositionVelocity : db::SimpleTag {
462 0 : using type = std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>;
463 : };
464 :
465 : template <size_t Dim>
466 0 : struct ParticlePositionVelocityCompute : ParticlePositionVelocity<Dim>,
467 : db::ComputeTag {
468 0 : using base = ParticlePositionVelocity<Dim>;
469 0 : using return_type = std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>;
470 0 : using argument_tags = tmpl::list<ExcisionSphere<Dim>, ::Tags::Time,
471 : domain::Tags::FunctionsOfTime>;
472 0 : static void function(
473 : gsl::not_null<std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>*>
474 : position_velocity,
475 : const ::ExcisionSphere<Dim>& excision_sphere, double time,
476 : const std::unordered_map<
477 : std::string,
478 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
479 : functions_of_time);
480 : };
481 : /// @}
482 :
483 : /*!
484 : * \brief The position of the scalar charge evolved by the worldtube singleton.
485 : * This tag is meant to be used by the worldtube singleton to evolve the orbit.
486 : */
487 : template <size_t Dim>
488 1 : struct EvolvedPosition : db::SimpleTag {
489 0 : using type = tnsr::I<DataVector, Dim>;
490 : };
491 :
492 : /*!
493 : * \brief The velocity of the scalar charge evolved by the worldtube singleton.
494 : * This tag is meant to be used by the worldtube singleton to evolve the orbit.
495 : */
496 : template <size_t Dim>
497 1 : struct EvolvedVelocity : db::SimpleTag {
498 0 : using type = tnsr::I<DataVector, Dim>;
499 : };
500 :
501 : /*!
502 : * \brief The position and velocity of the scalar charge particle orbiting a
503 : * central black hole given in inertial coordinates. This compute tag is meant
504 : * to be used by the worldtube singleton which evolves the position and velocity
505 : * according to an ODE along with the DG evolution.
506 : */
507 : template <size_t Dim>
508 1 : struct EvolvedParticlePositionVelocityCompute : ParticlePositionVelocity<Dim>,
509 : db::ComputeTag {
510 0 : using base = ParticlePositionVelocity<Dim>;
511 0 : using return_type = std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>;
512 0 : using argument_tags = tmpl::list<EvolvedPosition<Dim>, EvolvedVelocity<Dim>>;
513 0 : static void function(
514 : gsl::not_null<std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>*>
515 : position_velocity,
516 : const tnsr::I<DataVector, Dim>& evolved_position,
517 : const tnsr::I<DataVector, Dim>& evolved_velocity);
518 : };
519 :
520 : /// @{
521 : /*!
522 : * \brief Computes the coordinate geodesic acceleration of the particle in the
523 : * inertial frame in Kerr-Schild coordinates.
524 : */
525 : template <size_t Dim>
526 1 : struct GeodesicAcceleration : db::SimpleTag {
527 0 : using type = tnsr::I<double, Dim, Frame::Inertial>;
528 : };
529 :
530 : template <size_t Dim>
531 0 : struct GeodesicAccelerationCompute : GeodesicAcceleration<Dim>, db::ComputeTag {
532 0 : using base = GeodesicAcceleration<Dim>;
533 0 : using return_type = tnsr::I<double, Dim, Frame::Inertial>;
534 0 : using argument_tags = tmpl::list<
535 : ParticlePositionVelocity<Dim>,
536 : CurvedScalarWave::Tags::BackgroundSpacetime<gr::Solutions::KerrSchild>>;
537 0 : static void function(
538 : gsl::not_null<tnsr::I<double, Dim, Frame::Inertial>*> acceleration,
539 : const std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>&
540 : position_velocity,
541 : const gr::Solutions::KerrSchild& background_spacetime);
542 : };
543 : /// @}
544 :
545 : /*!
546 : * \brief The coordinate time dilation factor of the scalar charge, i.e. the 0th
547 : * component of its 4-velocity.
548 : */
549 1 : struct TimeDilationFactor : db::SimpleTag {
550 0 : using type = Scalar<double>;
551 : };
552 :
553 : /*!
554 : * \brief The parameters controlling the growth of the worldtube excision
555 : * sphere, see smooth_broken_power_law. The parameters here are, in order, the
556 : * amplitude, the transition radius, the transition width and the exponent.
557 : */
558 1 : struct WorldtubeRadiusParameters : db::SimpleTag {
559 0 : using type = std::array<double, 4>;
560 0 : using option_tags = tmpl::list<OptionTags::RadiusOptions<true>>;
561 0 : static constexpr bool pass_metavariables = false;
562 0 : static std::array<double, 4> create_from_options(
563 : const OptionTags::RadiusOptions<true>& params) {
564 : return {{params.exponent, params.amplitude, params.transition_radius,
565 : params.transition_width}};
566 : }
567 : };
568 :
569 : /*!
570 : * \brief The parameters controlling the growth of the black holes excision
571 : * sphere, see smooth_broken_power_law. The parameters here are, in order, the
572 : * amplitude, the transition radius, the transition width and the exponent.
573 : */
574 1 : struct BlackHoleRadiusParameters : db::SimpleTag {
575 0 : using type = std::array<double, 4>;
576 0 : using option_tags = tmpl::list<OptionTags::RadiusOptions<false>>;
577 0 : static constexpr bool pass_metavariables = false;
578 0 : static std::array<double, 4> create_from_options(
579 : const OptionTags::RadiusOptions<false>& params) {
580 : return {{params.exponent, params.amplitude, params.transition_radius,
581 : params.transition_width}};
582 : }
583 : };
584 :
585 : /// @{
586 : /*!
587 : * \brief A tuple of Tensors evaluated at the charge depending only the
588 : * background and the particle's position and velocity. These values are
589 : * effectively cached between different iterations of the worldtube scheme.
590 : */
591 : template <size_t Dim>
592 1 : struct BackgroundQuantities : db::SimpleTag {
593 0 : using type = tuples::TaggedTuple<
594 : gr::Tags::SpacetimeMetric<double, Dim>,
595 : gr::Tags::InverseSpacetimeMetric<double, Dim>,
596 : gr::Tags::SpacetimeChristoffelSecondKind<double, Dim>,
597 : gr::Tags::TraceSpacetimeChristoffelSecondKind<double, Dim>,
598 : Tags::TimeDilationFactor>;
599 : };
600 :
601 : template <size_t Dim>
602 0 : struct BackgroundQuantitiesCompute : BackgroundQuantities<Dim>, db::ComputeTag {
603 0 : using base = BackgroundQuantities<Dim>;
604 0 : using return_type = tuples::TaggedTuple<
605 : gr::Tags::SpacetimeMetric<double, Dim>,
606 : gr::Tags::InverseSpacetimeMetric<double, Dim>,
607 : gr::Tags::SpacetimeChristoffelSecondKind<double, Dim>,
608 : gr::Tags::TraceSpacetimeChristoffelSecondKind<double, Dim>,
609 : Tags::TimeDilationFactor>;
610 :
611 0 : using argument_tags = tmpl::list<
612 : ParticlePositionVelocity<Dim>,
613 : CurvedScalarWave::Tags::BackgroundSpacetime<gr::Solutions::KerrSchild>>;
614 0 : static void function(gsl::not_null<return_type*> result,
615 : const std::array<tnsr::I<double, Dim, Frame::Inertial>,
616 : 2>& position_velocity,
617 : const gr::Solutions::KerrSchild& background_spacetime);
618 : };
619 : /// @}
620 :
621 : /// @{
622 : /*!
623 : * \brief An optional that holds the coordinates of an element face abutting the
624 : * worldtube excision sphere. If the element does not abut the worldtube, this
625 : * holds std::nullopt. This tag should be in the databox of element chares. The
626 : * available frames are Grid and Inertial. The Centered template tag can be
627 : * turned on to center the coordinates around the position of the scalar
628 : * charge.
629 : */
630 : template <size_t Dim, typename Frame, bool Centered>
631 1 : struct FaceCoordinates : db::SimpleTag {
632 0 : using type = std::optional<tnsr::I<DataVector, Dim, Frame>>;
633 : };
634 :
635 : template <size_t Dim, typename Frame, bool Centered>
636 0 : struct FaceCoordinatesCompute : FaceCoordinates<Dim, Frame, Centered>,
637 : db::ComputeTag {
638 0 : using base = FaceCoordinates<Dim, Frame, Centered>;
639 0 : static constexpr bool needs_inertial_wt_coords =
640 : (Centered and std::is_same_v<Frame, ::Frame::Inertial>);
641 0 : using argument_tags = tmpl::flatten<
642 : tmpl::list<ExcisionSphere<Dim>, domain::Tags::Element<Dim>,
643 : domain::Tags::Coordinates<Dim, Frame>, domain::Tags::Mesh<Dim>,
644 : tmpl::conditional_t<needs_inertial_wt_coords,
645 : tmpl::list<ParticlePositionVelocity<Dim>>,
646 : tmpl::list<>>>>;
647 :
648 0 : using return_type = std::optional<tnsr::I<DataVector, Dim, Frame>>;
649 0 : static void function(
650 : const gsl::not_null<std::optional<tnsr::I<DataVector, Dim, Frame>>*>
651 : result,
652 : const ::ExcisionSphere<Dim>& excision_sphere, const Element<Dim>& element,
653 : const tnsr::I<DataVector, Dim, Frame>& coords, const Mesh<Dim>& mesh);
654 :
655 0 : static void function(
656 : const gsl::not_null<
657 : std::optional<tnsr::I<DataVector, Dim, ::Frame::Inertial>>*>
658 : result,
659 : const ::ExcisionSphere<Dim>& excision_sphere, const Element<Dim>& element,
660 : const tnsr::I<DataVector, Dim, ::Frame::Inertial>& coords,
661 : const Mesh<Dim>& mesh,
662 : const std::array<tnsr::I<double, Dim, ::Frame::Inertial>, 2>&
663 : particle_position_velocity);
664 : };
665 : /// @}
666 :
667 : /// @{
668 : /*!
669 : * \brief The value of the scalar field and its time derivative on element faces
670 : * forming the worldtube boundary, as well as the Euclidean area element of the
671 : * face.
672 : *
673 : * \details If the element does not abut the worldtube, this will be
674 : * `std::nullopt`.
675 : */
676 1 : struct FaceQuantities : db::SimpleTag {
677 0 : using type = std::optional<Variables<tmpl::list<
678 : CurvedScalarWave::Tags::Psi, ::Tags::dt<CurvedScalarWave::Tags::Psi>,
679 : gr::surfaces::Tags::AreaElement<DataVector>>>>;
680 : };
681 :
682 0 : struct FaceQuantitiesCompute : FaceQuantities, db::ComputeTag {
683 0 : static constexpr size_t Dim = 3;
684 0 : using base = FaceQuantities;
685 0 : using return_type = std::optional<Variables<tmpl::list<
686 : CurvedScalarWave::Tags::Psi, ::Tags::dt<CurvedScalarWave::Tags::Psi>,
687 : gr::surfaces::Tags::AreaElement<DataVector>>>>;
688 0 : using tags_to_slice_to_face =
689 : tmpl::list<CurvedScalarWave::Tags::Psi, CurvedScalarWave::Tags::Pi,
690 : CurvedScalarWave::Tags::Phi<Dim>,
691 : gr::Tags::Shift<DataVector, Dim>, gr::Tags::Lapse<DataVector>,
692 : domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
693 : Frame::Inertial>>;
694 0 : using argument_tags = tmpl::flatten<
695 : tmpl::list<tags_to_slice_to_face, ExcisionSphere<Dim>,
696 : domain::Tags::Element<Dim>, domain::Tags::Mesh<Dim>>>;
697 :
698 0 : static void function(
699 : gsl::not_null<return_type*> result, const Scalar<DataVector>& psi,
700 : const Scalar<DataVector>& pi, const tnsr::i<DataVector, Dim>& phi,
701 : const tnsr::I<DataVector, Dim>& shift, const Scalar<DataVector>& lapse,
702 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
703 : Frame::Inertial>& inv_jacobian,
704 : const ::ExcisionSphere<Dim>& excision_sphere, const Element<Dim>& element,
705 : const Mesh<Dim>& mesh);
706 : };
707 : /// @}
708 :
709 : /*!
710 : * \brief Configuration and dispatch object for puncture-field expressions.
711 : */
712 1 : struct PunctureFieldConfig : db::SimpleTag {
713 0 : using type = CurvedScalarWave::Worldtube::PunctureField;
714 0 : static constexpr bool pass_metavariables = false;
715 0 : using option_tags = tmpl::list<OptionTags::PunctureField>;
716 0 : static type create_from_options(
717 : const CurvedScalarWave::Worldtube::PunctureField& puncture_field) {
718 : return puncture_field;
719 : }
720 : };
721 :
722 : /*!
723 : * \brief The internal expansion order of the worldtube solution.
724 : */
725 1 : struct ExpansionOrder : db::SimpleTag {
726 0 : using type = size_t;
727 0 : static constexpr bool pass_metavariables = false;
728 0 : using option_tags = tmpl::list<OptionTags::PunctureField>;
729 0 : static size_t create_from_options(
730 : const CurvedScalarWave::Worldtube::PunctureField& puncture_field) {
731 : return puncture_field.expansion_order();
732 : }
733 : };
734 :
735 : /// @{
736 : /*!
737 : * Computes the puncture field on an element face abutting the worldtube
738 : * assuming geodesic acceleration. If the current element does not abut the
739 : * worldtube this holds a std::nullopt.
740 : */
741 : template <size_t Dim>
742 1 : struct GeodesicPunctureField : db::SimpleTag {
743 0 : using type = std::optional<Variables<tmpl::list<
744 : CurvedScalarWave::Tags::Psi, ::Tags::dt<CurvedScalarWave::Tags::Psi>,
745 : ::Tags::deriv<CurvedScalarWave::Tags::Psi, tmpl::size_t<3>,
746 : Frame::Inertial>>>>;
747 : };
748 :
749 : template <size_t Dim>
750 0 : struct GeodesicPunctureFieldCompute : GeodesicPunctureField<Dim>,
751 : db::ComputeTag {
752 0 : using base = GeodesicPunctureField<Dim>;
753 0 : using argument_tags =
754 : tmpl::list<FaceCoordinates<Dim, Frame::Inertial, true>,
755 : ParticlePositionVelocity<Dim>, GeodesicAcceleration<Dim>,
756 : Charge, PunctureFieldConfig>;
757 0 : using return_type = std::optional<Variables<tmpl::list<
758 : CurvedScalarWave::Tags::Psi, ::Tags::dt<CurvedScalarWave::Tags::Psi>,
759 : ::Tags::deriv<CurvedScalarWave::Tags::Psi, tmpl::size_t<3>,
760 : Frame::Inertial>>>>;
761 0 : static void function(
762 : const gsl::not_null<return_type*> result,
763 : const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
764 : inertial_face_coords_centered,
765 : const std::array<tnsr::I<double, Dim, ::Frame::Inertial>, 2>&
766 : particle_position_velocity,
767 : const tnsr::I<double, Dim>& particle_acceleration, double charge,
768 : const CurvedScalarWave::Worldtube::PunctureField& puncture_field);
769 : };
770 : /// @}
771 :
772 : /*!
773 : * \brief Holds the current iteration of the puncture field computed with the
774 : * current iteration of the acceleration which includes the scalar self-force.
775 : * It is computed in `Actions::IteratePunctureField`.
776 : */
777 : template <size_t Dim>
778 1 : struct IteratedPunctureField : db::SimpleTag {
779 0 : using type = std::optional<Variables<tmpl::list<
780 : CurvedScalarWave::Tags::Psi, ::Tags::dt<CurvedScalarWave::Tags::Psi>,
781 : ::Tags::deriv<CurvedScalarWave::Tags::Psi, tmpl::size_t<3>,
782 : Frame::Inertial>>>>;
783 : };
784 :
785 : /*!
786 : * The acceleration terms computed by the worldtube singleton and sent to the
787 : * neighboring elements used to compute the puncture field.
788 : */
789 1 : struct AccelerationTerms : db::SimpleTag {
790 0 : using type = Scalar<DataVector>;
791 : };
792 :
793 : /*!
794 : * \brief A map that holds the grid coordinates centered on the worldtube of
795 : * all element faces abutting the worldtube with the corresponding ElementIds.
796 : */
797 : template <size_t Dim>
798 1 : struct ElementFacesGridCoordinates : db::SimpleTag {
799 0 : using type =
800 : std::unordered_map<ElementId<Dim>, tnsr::I<DataVector, Dim, Frame::Grid>>;
801 : };
802 :
803 : /*!
804 : * \brief The solution inside the worldtube, evaluated at the face coordinates
805 : * of an abutting element. This tag is used to provide boundary conditions to
806 : * the element in \ref CurvedScalarWave::BoundaryConditions::Worldtube .
807 : */
808 : template <size_t Dim>
809 1 : struct WorldtubeSolution : db::SimpleTag {
810 0 : using type = Variables<
811 : tmpl::list<::CurvedScalarWave::Tags::Psi, ::CurvedScalarWave::Tags::Pi,
812 : ::CurvedScalarWave::Tags::Phi<Dim>>>;
813 : };
814 :
815 : /*!
816 : * \brief The scalar field inside the worldtube.
817 : *
818 : * \details This tag is used as a base tag for Stf::Tags::StfTensor
819 : */
820 1 : struct PsiWorldtube : db::SimpleTag {
821 0 : using type = Scalar<double>;
822 : };
823 :
824 : /*!
825 : * \brief Holds the constant coefficient of the regular field inside the
826 : * worldtube.
827 : *
828 : * \details At orders n = 0 or 1 this is just equal to the monopole, but at n =
829 : * 2, the monopole gets an additional contribution from the trace of the second
830 : * order coefficient. At this point, this tag is used to solve an ODE based on
831 : * the expanded Klein-Gordon equation. It is implemented as a `Scalar` of size 1
832 : * because the evolution system does not work with doubles.
833 : */
834 1 : struct Psi0 : db::SimpleTag {
835 0 : using type = Scalar<DataVector>;
836 : };
837 :
838 : /*!
839 : * \brief Holds the time derivative of Psi0 which is used as a reduction
840 : * variable.
841 : */
842 1 : struct dtPsi0 : db::SimpleTag {
843 0 : using type = Scalar<DataVector>;
844 : };
845 :
846 : /*!
847 : * \brief Sets Gamma1 to zero throughout the domain. The equations are given in
848 : * Initialization::InitializeConstraintDampingGammas.
849 : */
850 1 : struct ConstraintGamma1Compute : CurvedScalarWave::Tags::ConstraintGamma1,
851 : db::ComputeTag {
852 0 : static constexpr size_t Dim = 3;
853 0 : using base = CurvedScalarWave::Tags::ConstraintGamma1;
854 0 : using return_type = Scalar<DataVector>;
855 0 : using argument_tags =
856 : tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>>;
857 0 : static void function(gsl::not_null<Scalar<DataVector>*> gamma1,
858 : const tnsr::I<DataVector, Dim, Frame::Inertial>& coords);
859 : };
860 :
861 : /*!
862 : * \brief Sets Gamma2 to a Gaussian that falls off to a constant value centered
863 : * on the position of the particle. This was found to be necessary for a stable
864 : * evolution. The equations are given in
865 : * Initialization::InitializeConstraintDampingGammas.
866 : */
867 1 : struct ConstraintGamma2Compute : CurvedScalarWave::Tags::ConstraintGamma2,
868 : db::ComputeTag {
869 0 : static constexpr size_t Dim = 3;
870 0 : using base = CurvedScalarWave::Tags::ConstraintGamma2;
871 0 : using return_type = Scalar<DataVector>;
872 0 : using argument_tags =
873 : tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>,
874 : ParticlePositionVelocity<Dim>>;
875 0 : static void function(
876 : gsl::not_null<Scalar<DataVector>*> gamma2,
877 : const tnsr::I<DataVector, Dim, Frame::Inertial>& coords,
878 : const std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>& pos_vel);
879 : };
880 :
881 : } // namespace Tags
882 : } // namespace CurvedScalarWave::Worldtube
|
