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