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 "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
28 : #include "Options/Auto.hpp"
29 : #include "Options/String.hpp"
30 : #include "ParallelAlgorithms/EventsAndTriggers/Trigger.hpp"
31 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp"
32 : #include "PointwiseFunctions/GeneralRelativity/Surfaces/Tags.hpp"
33 : #include "Utilities/EqualWithinRoundoff.hpp"
34 : #include "Utilities/Gsl.hpp"
35 : #include "Utilities/Serialization/Serialize.hpp"
36 : #include "Utilities/TaggedTuple.hpp"
37 :
38 : /// \cond
39 : namespace Tags {
40 : struct Time;
41 : } // namespace Tags
42 : namespace OptionTags {
43 : struct InitialTime;
44 : } // namespace OptionTags
45 : /// \endcond
46 :
47 : namespace CurvedScalarWave::Worldtube {
48 : /*!
49 : * \brief Option tags for the worldtube
50 : */
51 1 : namespace OptionTags {
52 : /*!
53 : * \brief Options for the worldtube
54 : */
55 1 : struct Worldtube {
56 0 : static constexpr Options::String help = {"Options for the Worldtube"};
57 : };
58 :
59 : /*!
60 : * \brief The value of the scalar charge in units of the black hole mass M.
61 : */
62 1 : struct Charge {
63 0 : using type = double;
64 0 : static constexpr Options::String help{
65 : "The value of the scalar charge in units of the black hole mass M."};
66 0 : using group = Worldtube;
67 : };
68 :
69 : /*!
70 : * \brief Options for the scalar self-force. Select `None` for a purely geodesic
71 : * evolution
72 : */
73 1 : struct SelfForceOptions {
74 0 : static constexpr Options::String help = {
75 : "Options for the scalar self-force. Select `None` for a purely geodesic "
76 : "evolution"};
77 0 : using group = Worldtube;
78 0 : using type = Options::Auto<SelfForceOptions, Options::AutoLabel::None>;
79 :
80 0 : struct Mass {
81 0 : using type = double;
82 0 : static constexpr Options::String help{
83 : "The mass of the scalar particle in units of the black hole mass M."};
84 0 : static double lower_bound() { return 0.; }
85 : };
86 :
87 0 : struct Iterations {
88 0 : using type = size_t;
89 0 : static constexpr Options::String help{
90 : "The number of iterations used to compute the particle acceleration. "
91 : "Must be at least 1 as 0 iterations corresponds to the geodesic "
92 : "acceleration."};
93 0 : static size_t lower_bound() { return 1; }
94 : };
95 0 : void pup(PUP::er& p) {
96 : p | mass;
97 : p | iterations;
98 : }
99 :
100 0 : using options = tmpl::list<Mass, Iterations>;
101 :
102 0 : double mass{};
103 0 : size_t iterations{};
104 : };
105 :
106 : /*!
107 : * \brief Name of the excision sphere designated to act as a worldtube
108 : */
109 1 : struct ExcisionSphere {
110 0 : using type = std::string;
111 0 : static constexpr Options::String help{
112 : "The name of the excision sphere as returned by the domain."};
113 0 : using group = Worldtube;
114 : };
115 :
116 : /*!
117 : * \brief Triggers at which to write the coefficients of the worldtube's
118 : * internal Taylor series to file.
119 : */
120 1 : struct ObserveCoefficientsTrigger {
121 0 : using type = std::unique_ptr<Trigger>;
122 0 : static constexpr Options::String help{
123 : "Specifies a non-dense trigger in which the coefficients of the internal "
124 : "regular field expansion are written to file."};
125 0 : using group = Worldtube;
126 : };
127 :
128 : /*!
129 : * \brief The internal expansion order of the worldtube solution.
130 : */
131 1 : struct ExpansionOrder {
132 0 : using type = size_t;
133 0 : static constexpr Options::String help{
134 : "The internal expansion order of the worldtube solution. Currently "
135 : "orders 0 and 1 are implemented"};
136 0 : static size_t upper_bound() { return 1; }
137 0 : using group = Worldtube;
138 : };
139 : } // namespace OptionTags
140 :
141 : /*!
142 : * \brief Tags related to the worldtube
143 : */
144 : namespace Tags {
145 : /*!
146 : * \brief Dummy tag that throws an error if the input file does not describe a
147 : * circular orbit.
148 : */
149 : template <size_t Dim, typename BackgroundType>
150 1 : struct CheckInputFile : db::SimpleTag {
151 0 : using type = bool;
152 0 : using option_tags = tmpl::list<
153 : domain::OptionTags::DomainCreator<Dim>, OptionTags::ExcisionSphere,
154 : CurvedScalarWave::OptionTags::BackgroundSpacetime<BackgroundType>>;
155 :
156 : // puncture field is specialised on Kerr-Schild bakckground.
157 : static_assert(std::is_same_v<BackgroundType, gr::Solutions::KerrSchild>);
158 0 : static constexpr bool pass_metavariables = false;
159 0 : static bool create_from_options(
160 : const std::unique_ptr<::DomainCreator<Dim>>& domain_creator,
161 : const std::string& excision_sphere_name,
162 : const BackgroundType& kerr_schild_background) {
163 : if (not kerr_schild_background.zero_spin()) {
164 : ERROR(
165 : "Black hole spin is not implemented yet but you requested non-zero "
166 : "spin.");
167 : }
168 : if (not equal_within_roundoff(kerr_schild_background.center(),
169 : make_array(0., 0., 0.))) {
170 : ERROR("The central black hole must be centered at [0., 0., 0.].");
171 : }
172 : if (not equal_within_roundoff(kerr_schild_background.mass(), 1.)) {
173 : ERROR("The central black hole must have mass 1.");
174 : }
175 : const auto domain = domain_creator->create_domain();
176 : const auto& excision_spheres = domain.excision_spheres();
177 : const auto& excision_sphere = excision_spheres.at(excision_sphere_name);
178 : const double orbital_radius = get<0>(excision_sphere.center());
179 : const auto& functions_of_time = domain_creator->functions_of_time();
180 : if (not functions_of_time.count("Rotation")) {
181 : ERROR("Expected functions of time to contain 'Rotation'.");
182 : }
183 : // dynamic cast to access `angle_func_and_deriv` method
184 : const auto* rotation_function_of_time =
185 : dynamic_cast<domain::FunctionsOfTime::QuaternionFunctionOfTime<3>*>(
186 : &*functions_of_time.at("Rotation"));
187 : if (rotation_function_of_time == nullptr) {
188 : ERROR("Failed dynamic cast to QuaternionFunctionOfTime.");
189 : }
190 : const auto angular_velocity =
191 : rotation_function_of_time->angle_func_and_deriv(0.).at(1);
192 : if (equal_within_roundoff(orbital_radius, 0.)) {
193 : ERROR("The orbital radius was set to 0.");
194 : }
195 : if (not equal_within_roundoff(
196 : angular_velocity,
197 : DataVector{0.0, 0.0, pow(orbital_radius, -1.5)})) {
198 : ERROR(
199 : "Only circular orbits are implemented at the moment so the "
200 : "angular velocity should be [0., 0., orbital_radius^(-3/2)] = "
201 : << "[0., 0., " << pow(orbital_radius, -1.5) << "]");
202 : }
203 : return true;
204 : }
205 : };
206 :
207 : /*!
208 : * \brief The excision sphere corresponding to the worldtube
209 : */
210 : template <size_t Dim>
211 1 : struct ExcisionSphere : db::SimpleTag {
212 0 : using type = ::ExcisionSphere<Dim>;
213 0 : using option_tags = tmpl::list<domain::OptionTags::DomainCreator<Dim>,
214 : OptionTags::ExcisionSphere>;
215 0 : static constexpr bool pass_metavariables = false;
216 0 : static ::ExcisionSphere<Dim> create_from_options(
217 : const std::unique_ptr<::DomainCreator<Dim>>& domain_creator,
218 : const std::string& excision_sphere) {
219 : const auto domain = domain_creator->create_domain();
220 : const auto& excision_spheres = domain.excision_spheres();
221 : if (excision_spheres.count(excision_sphere) == 0) {
222 : ERROR("Specified excision sphere '"
223 : << excision_sphere
224 : << "' not available. Available excision spheres are: "
225 : << keys_of(excision_spheres));
226 : }
227 : return excision_spheres.at(excision_sphere);
228 : }
229 : };
230 :
231 : /*!
232 : * \brief Triggers at which to write the coefficients of the worldtube's
233 : * internal Taylor series to file.
234 : */
235 1 : struct ObserveCoefficientsTrigger : db::SimpleTag {
236 0 : using type = std::unique_ptr<Trigger>;
237 0 : using option_tags = tmpl::list<OptionTags::ObserveCoefficientsTrigger>;
238 0 : static constexpr bool pass_metavariables = false;
239 0 : static std::unique_ptr<Trigger> create_from_options(
240 : const std::unique_ptr<Trigger>& trigger) {
241 : return deserialize<type>(serialize<type>(trigger).data());
242 : }
243 : };
244 :
245 : /*!
246 : * \brief The value of the scalar charge
247 : */
248 1 : struct Charge : db::SimpleTag {
249 0 : using type = double;
250 0 : using option_tags = tmpl::list<OptionTags::Charge>;
251 0 : static constexpr bool pass_metavariables = false;
252 0 : static double create_from_options(const double charge) { return charge; };
253 : };
254 :
255 : /*!
256 : * \brief The mass of the scalar charge. Only has a value if the scalar self
257 : * force is applied.
258 : */
259 1 : struct Mass : db::SimpleTag {
260 0 : using type = std::optional<double>;
261 0 : using option_tags = tmpl::list<OptionTags::SelfForceOptions>;
262 0 : static constexpr bool pass_metavariables = false;
263 0 : static std::optional<double> create_from_options(
264 : const std::optional<OptionTags::SelfForceOptions>& self_force_options) {
265 : return self_force_options.has_value()
266 : ? std::make_optional(self_force_options->mass)
267 : : std::nullopt;
268 : }
269 : };
270 :
271 : /*!
272 : * \brief The maximum number of iterations that will be applied to the
273 : * acceleration of the particle.
274 : */
275 1 : struct MaxIterations : db::SimpleTag {
276 0 : using type = size_t;
277 0 : using option_tags = tmpl::list<OptionTags::SelfForceOptions>;
278 0 : static constexpr bool pass_metavariables = false;
279 0 : static size_t create_from_options(
280 : const std::optional<OptionTags::SelfForceOptions>& self_force_options) {
281 : return self_force_options.has_value() ? self_force_options->iterations : 0;
282 : }
283 : };
284 :
285 : /*!
286 : * \brief The current number of iterations that has been applied to the
287 : * acceleration of the particle.
288 : */
289 1 : struct CurrentIteration : db::SimpleTag {
290 0 : using type = size_t;
291 : };
292 :
293 : /*!
294 : * \brief The initial position and velocity of the scalar charge in inertial
295 : * coordinates.
296 : */
297 1 : struct InitialPositionAndVelocity : db::SimpleTag {
298 0 : using type = std::array<tnsr::I<double, 3, Frame::Inertial>, 2>;
299 0 : using option_tags =
300 : tmpl::list<domain::OptionTags::DomainCreator<3>,
301 : OptionTags::ExcisionSphere, ::OptionTags::InitialTime>;
302 0 : static constexpr bool pass_metavariables = false;
303 0 : static type create_from_options(
304 : const std::unique_ptr<::DomainCreator<3>>& domain_creator,
305 : const std::string& excision_sphere_name, const double initial_time) {
306 : // only evaluated at initial time, so expiration times don't matter
307 : const auto initial_fot = domain_creator->functions_of_time();
308 : const auto domain = domain_creator->create_domain();
309 : const auto& excision_sphere =
310 : domain.excision_spheres().at(excision_sphere_name);
311 : ASSERT(excision_sphere.is_time_dependent(),
312 : "excision_sphere not time dependent");
313 : const auto& maps = excision_sphere.moving_mesh_grid_to_inertial_map();
314 : const auto mapped_tuple = maps.coords_frame_velocity_jacobians(
315 : excision_sphere.center(), initial_time, initial_fot);
316 : return {std::get<0>(mapped_tuple), std::get<3>(mapped_tuple)};
317 : }
318 : };
319 :
320 : /// @{
321 : /*!
322 : * \brief The position and velocity of the scalar charge particle orbiting a
323 : * central black hole given in inertial coordinates. This compute tag is meant
324 : * to be used by the elements.
325 : */
326 : template <size_t Dim>
327 1 : struct ParticlePositionVelocity : db::SimpleTag {
328 0 : using type = std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>;
329 : };
330 :
331 : template <size_t Dim>
332 0 : struct ParticlePositionVelocityCompute : ParticlePositionVelocity<Dim>,
333 : db::ComputeTag {
334 0 : using base = ParticlePositionVelocity<Dim>;
335 0 : using return_type = std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>;
336 0 : using argument_tags = tmpl::list<ExcisionSphere<Dim>, ::Tags::Time,
337 : domain::Tags::FunctionsOfTime>;
338 0 : static void function(
339 : gsl::not_null<std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>*>
340 : position_velocity,
341 : const ::ExcisionSphere<Dim>& excision_sphere, double time,
342 : const std::unordered_map<
343 : std::string,
344 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
345 : functions_of_time);
346 : };
347 : /// @}
348 :
349 : /*!
350 : * \brief The position of the scalar charge evolved by the worldtube singleton.
351 : * This tag is meant to be used by the worldtube singleton to evolve the orbit.
352 : */
353 : template <size_t Dim>
354 1 : struct EvolvedPosition : db::SimpleTag {
355 0 : using type = tnsr::I<DataVector, Dim>;
356 : };
357 :
358 : /*!
359 : * \brief The velocity of the scalar charge evolved by the worldtube singleton.
360 : * This tag is meant to be used by the worldtube singleton to evolve the orbit.
361 : */
362 : template <size_t Dim>
363 1 : struct EvolvedVelocity : db::SimpleTag {
364 0 : using type = tnsr::I<DataVector, Dim>;
365 : };
366 :
367 : /*!
368 : * \brief The position and velocity of the scalar charge particle orbiting a
369 : * central black hole given in inertial coordinates. This compute tag is meant
370 : * to be used by the worldtube singleton which evolves the position and velocity
371 : * according to an ODE along with the DG evolution.
372 : */
373 : template <size_t Dim>
374 1 : struct EvolvedParticlePositionVelocityCompute : ParticlePositionVelocity<Dim>,
375 : db::ComputeTag {
376 0 : using base = ParticlePositionVelocity<Dim>;
377 0 : using return_type = std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>;
378 0 : using argument_tags = tmpl::list<EvolvedPosition<Dim>, EvolvedVelocity<Dim>>;
379 0 : static void function(
380 : gsl::not_null<std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>*>
381 : position_velocity,
382 : const tnsr::I<DataVector, Dim>& evolved_position,
383 : const tnsr::I<DataVector, Dim>& evolved_velocity);
384 : };
385 :
386 : /// @{
387 : /*!
388 : * \brief Computes the coordinate geodesic acceleration of the particle in the
389 : * inertial frame in Kerr-Schild coordinates.
390 : */
391 : template <size_t Dim>
392 1 : struct GeodesicAcceleration : db::SimpleTag {
393 0 : using type = tnsr::I<double, Dim, Frame::Inertial>;
394 : };
395 :
396 : template <size_t Dim>
397 0 : struct GeodesicAccelerationCompute : GeodesicAcceleration<Dim>, db::ComputeTag {
398 0 : using base = GeodesicAcceleration<Dim>;
399 0 : using return_type = tnsr::I<double, Dim, Frame::Inertial>;
400 0 : using argument_tags = tmpl::list<
401 : ParticlePositionVelocity<Dim>,
402 : CurvedScalarWave::Tags::BackgroundSpacetime<gr::Solutions::KerrSchild>>;
403 0 : static void function(
404 : gsl::not_null<tnsr::I<double, Dim, Frame::Inertial>*> acceleration,
405 : const std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>&
406 : position_velocity,
407 : const gr::Solutions::KerrSchild& background_spacetime);
408 : };
409 : /// @}
410 :
411 : /*!
412 : * \brief The coordinate time dilation factor of the scalar charge, i.e. the 0th
413 : * component of its 4-velocity.
414 : */
415 1 : struct TimeDilationFactor : db::SimpleTag {
416 0 : using type = Scalar<double>;
417 : };
418 :
419 : /// @{
420 : /*!
421 : * \brief A tuple of Tensors evaluated at the charge depending only the
422 : * background and the particle's position and velocity. These values are
423 : * effectively cached between different iterations of the worldtube scheme.
424 : */
425 : template <size_t Dim>
426 1 : struct BackgroundQuantities : db::SimpleTag {
427 0 : using type =
428 : tuples::TaggedTuple<gr::Tags::SpacetimeMetric<double, Dim>,
429 : gr::Tags::InverseSpacetimeMetric<double, Dim>,
430 : Tags::TimeDilationFactor>;
431 : };
432 :
433 : template <size_t Dim>
434 0 : struct BackgroundQuantitiesCompute : BackgroundQuantities<Dim>, db::ComputeTag {
435 0 : using base = BackgroundQuantities<Dim>;
436 0 : using return_type =
437 : tuples::TaggedTuple<gr::Tags::SpacetimeMetric<double, Dim>,
438 : gr::Tags::InverseSpacetimeMetric<double, Dim>,
439 : Tags::TimeDilationFactor>;
440 :
441 0 : using argument_tags = tmpl::list<
442 : ParticlePositionVelocity<Dim>,
443 : CurvedScalarWave::Tags::BackgroundSpacetime<gr::Solutions::KerrSchild>>;
444 0 : static void function(gsl::not_null<return_type*> result,
445 : const std::array<tnsr::I<double, Dim, Frame::Inertial>,
446 : 2>& position_velocity,
447 : const gr::Solutions::KerrSchild& background_spacetime);
448 : };
449 : /// @}
450 :
451 : /// @{
452 : /*!
453 : * \brief An optional that holds the coordinates of an element face abutting the
454 : * worldtube excision sphere. If the element does not abut the worldtube, this
455 : * holds std::nullopt. This tag should be in the databox of element chares. The
456 : * available frames are Grid and Inertial. The Centered template tag can be
457 : * turned on to center the coordinates around the position of the scalar
458 : * charge.
459 : */
460 : template <size_t Dim, typename Frame, bool Centered>
461 1 : struct FaceCoordinates : db::SimpleTag {
462 0 : using type = std::optional<tnsr::I<DataVector, Dim, Frame>>;
463 : };
464 :
465 : template <size_t Dim, typename Frame, bool Centered>
466 0 : struct FaceCoordinatesCompute : FaceCoordinates<Dim, Frame, Centered>,
467 : db::ComputeTag {
468 0 : using base = FaceCoordinates<Dim, Frame, Centered>;
469 0 : static constexpr bool needs_inertial_wt_coords =
470 : (Centered and std::is_same_v<Frame, ::Frame::Inertial>);
471 0 : using argument_tags = tmpl::flatten<
472 : tmpl::list<ExcisionSphere<Dim>, domain::Tags::Element<Dim>,
473 : domain::Tags::Coordinates<Dim, Frame>, domain::Tags::Mesh<Dim>,
474 : tmpl::conditional_t<needs_inertial_wt_coords,
475 : tmpl::list<ParticlePositionVelocity<Dim>>,
476 : tmpl::list<>>>>;
477 :
478 0 : using return_type = std::optional<tnsr::I<DataVector, Dim, Frame>>;
479 0 : static void function(
480 : const gsl::not_null<std::optional<tnsr::I<DataVector, Dim, Frame>>*>
481 : result,
482 : const ::ExcisionSphere<Dim>& excision_sphere, const Element<Dim>& element,
483 : const tnsr::I<DataVector, Dim, Frame>& coords, const Mesh<Dim>& mesh);
484 :
485 0 : static void function(
486 : const gsl::not_null<
487 : std::optional<tnsr::I<DataVector, Dim, ::Frame::Inertial>>*>
488 : result,
489 : const ::ExcisionSphere<Dim>& excision_sphere, const Element<Dim>& element,
490 : const tnsr::I<DataVector, Dim, ::Frame::Inertial>& coords,
491 : const Mesh<Dim>& mesh,
492 : const std::array<tnsr::I<double, Dim, ::Frame::Inertial>, 2>&
493 : particle_position_velocity);
494 : };
495 : /// @}
496 :
497 : /// @{
498 : /*!
499 : * \brief The value of the scalar field and its time derivative on element faces
500 : * forming the worldtube boundary, as well as the Euclidean area element of the
501 : * face.
502 : *
503 : * \details If the element does not abut the worldtube, this will be
504 : * `std::nullopt`.
505 : */
506 1 : struct FaceQuantities : db::SimpleTag {
507 0 : using type = std::optional<Variables<tmpl::list<
508 : CurvedScalarWave::Tags::Psi, ::Tags::dt<CurvedScalarWave::Tags::Psi>,
509 : gr::surfaces::Tags::AreaElement<DataVector>>>>;
510 : };
511 :
512 0 : struct FaceQuantitiesCompute : FaceQuantities, db::ComputeTag {
513 0 : static constexpr size_t Dim = 3;
514 0 : using base = FaceQuantities;
515 0 : using return_type = std::optional<Variables<tmpl::list<
516 : CurvedScalarWave::Tags::Psi, ::Tags::dt<CurvedScalarWave::Tags::Psi>,
517 : gr::surfaces::Tags::AreaElement<DataVector>>>>;
518 0 : using tags_to_slice_to_face =
519 : tmpl::list<CurvedScalarWave::Tags::Psi, CurvedScalarWave::Tags::Pi,
520 : CurvedScalarWave::Tags::Phi<Dim>,
521 : gr::Tags::Shift<DataVector, Dim>, gr::Tags::Lapse<DataVector>,
522 : domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
523 : Frame::Inertial>>;
524 0 : using argument_tags = tmpl::flatten<
525 : tmpl::list<tags_to_slice_to_face, ExcisionSphere<Dim>,
526 : domain::Tags::Element<Dim>, domain::Tags::Mesh<Dim>>>;
527 :
528 0 : static void function(
529 : gsl::not_null<return_type*> result, const Scalar<DataVector>& psi,
530 : const Scalar<DataVector>& pi, const tnsr::i<DataVector, Dim>& phi,
531 : const tnsr::I<DataVector, Dim>& shift, const Scalar<DataVector>& lapse,
532 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
533 : Frame::Inertial>& inv_jacobian,
534 : const ::ExcisionSphere<Dim>& excision_sphere, const Element<Dim>& element,
535 : const Mesh<Dim>& mesh);
536 : };
537 : /// @}
538 :
539 : /*!
540 : * \brief The internal expansion order of the worldtube solution.
541 : */
542 1 : struct ExpansionOrder : db::SimpleTag {
543 0 : using type = size_t;
544 0 : static constexpr bool pass_metavariables = false;
545 0 : using option_tags = tmpl::list<OptionTags::ExpansionOrder>;
546 0 : static size_t create_from_options(const size_t order) { return order; }
547 : };
548 :
549 : /// @{
550 : /*!
551 : * Computes the puncture field on an element face abutting the worldtube
552 : * assuming geodesic acceleration. If the current element does not abut the
553 : * worldtube this holds a std::nullopt.
554 : */
555 : template <size_t Dim>
556 1 : struct PunctureField : db::SimpleTag {
557 0 : using type = std::optional<Variables<tmpl::list<
558 : CurvedScalarWave::Tags::Psi, ::Tags::dt<CurvedScalarWave::Tags::Psi>,
559 : ::Tags::deriv<CurvedScalarWave::Tags::Psi, tmpl::size_t<3>,
560 : Frame::Inertial>>>>;
561 : };
562 :
563 : template <size_t Dim>
564 0 : struct PunctureFieldCompute : PunctureField<Dim>, db::ComputeTag {
565 0 : using base = PunctureField<Dim>;
566 0 : using argument_tags =
567 : tmpl::list<FaceCoordinates<Dim, Frame::Inertial, true>,
568 : ParticlePositionVelocity<Dim>, GeodesicAcceleration<Dim>,
569 : Charge, ExpansionOrder>;
570 0 : using return_type = std::optional<Variables<tmpl::list<
571 : CurvedScalarWave::Tags::Psi, ::Tags::dt<CurvedScalarWave::Tags::Psi>,
572 : ::Tags::deriv<CurvedScalarWave::Tags::Psi, tmpl::size_t<3>,
573 : Frame::Inertial>>>>;
574 0 : static void function(
575 : const gsl::not_null<return_type*> result,
576 : const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
577 : inertial_face_coords_centered,
578 : const std::array<tnsr::I<double, Dim, ::Frame::Inertial>, 2>&
579 : particle_position_velocity,
580 : const tnsr::I<double, Dim>& particle_acceleration, double charge,
581 : const size_t expansion_order);
582 : };
583 : /// @}
584 :
585 : /*!
586 : * \brief A map that holds the grid coordinates centered on the worldtube of
587 : * all element faces abutting the worldtube with the corresponding ElementIds.
588 : */
589 : template <size_t Dim>
590 1 : struct ElementFacesGridCoordinates : db::SimpleTag {
591 0 : using type =
592 : std::unordered_map<ElementId<Dim>, tnsr::I<DataVector, Dim, Frame::Grid>>;
593 : };
594 :
595 : /*!
596 : * \brief The solution inside the worldtube, evaluated at the face coordinates
597 : * of an abutting element. This tag is used to provide boundary conditions to
598 : * the element in \ref CurvedScalarWave::BoundaryConditions::Worldtube .
599 : */
600 : template <size_t Dim>
601 1 : struct WorldtubeSolution : db::SimpleTag {
602 0 : using type = Variables<
603 : tmpl::list<::CurvedScalarWave::Tags::Psi, ::CurvedScalarWave::Tags::Pi,
604 : ::CurvedScalarWave::Tags::Phi<Dim>>>;
605 : };
606 :
607 : /*!
608 : * \brief The scalar field inside the worldtube.
609 : *
610 : * \details This tag is used as a base tag for Stf::Tags::StfTensor
611 : */
612 1 : struct PsiWorldtube : db::SimpleTag {
613 0 : using type = Scalar<double>;
614 : };
615 :
616 : /*!
617 : * \brief Holds the constant coefficient of the regular field inside the
618 : * worldtube.
619 : *
620 : * \details At orders n = 0 or 1 this is just equal to the monopole, but at n =
621 : * 2, the monopole gets an additional contribution from the trace of the second
622 : * order coefficient. At this point, this tag is used to solve an ODE based on
623 : * the expanded Klein-Gordon equation. It is implemented as a `Scalar` of size 1
624 : * because the evolution system does not work with doubles.
625 : */
626 1 : struct Psi0 : db::SimpleTag {
627 0 : using type = Scalar<DataVector>;
628 : };
629 :
630 : /*!
631 : * \brief Holds the time derivative of Psi0 which is used as a reduction
632 : * variable.
633 : */
634 1 : struct dtPsi0 : db::SimpleTag {
635 0 : using type = Scalar<DataVector>;
636 : };
637 :
638 : } // namespace Tags
639 : } // namespace CurvedScalarWave::Worldtube
|
