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 <unordered_map>
9 : #include <utility>
10 :
11 : #include "ControlSystem/UpdateFunctionOfTime.hpp"
12 : #include "DataStructures/DataBox/DataBox.hpp"
13 : #include "DataStructures/DataBox/Prefixes.hpp"
14 : #include "Domain/Creators/Tags/Domain.hpp"
15 : #include "Evolution/Systems/CurvedScalarWave/Worldtube/Inboxes.hpp"
16 : #include "Evolution/Systems/CurvedScalarWave/Worldtube/RadiusFunctions.hpp"
17 : #include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp"
18 : #include "Evolution/Systems/CurvedScalarWave/Worldtube/Worldtube.hpp"
19 : #include "Parallel/AlgorithmExecution.hpp"
20 : #include "Parallel/GlobalCache.hpp"
21 : #include "ParallelAlgorithms/EventsAndTriggers/Tags.hpp"
22 : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
23 : #include "Time/Tags/TimeStepId.hpp"
24 : #include "Time/TimeStepId.hpp"
25 : #include "Time/TimeSteppers/TimeStepper.hpp"
26 : #include "Utilities/Gsl.hpp"
27 : #include "Utilities/TMPL.hpp"
28 :
29 : namespace control_system {
30 : // forward declare
31 : struct UpdateMultipleFunctionsOfTime;
32 : } // namespace control_system
33 :
34 : namespace CurvedScalarWave::Worldtube::Actions {
35 :
36 : /*!
37 : * \brief Updates the functions of time according to the motion of the
38 : * worldtube.
39 : *
40 : * \details We demand that the scalar charge is always in the center of the
41 : * worldtube and therefore deform the grid to the track the particle's motion.
42 : * In addition, the worldtube and black hole excision sphere radii are adjusted
43 : * according to smooth_broken_power_law.
44 : */
45 1 : struct UpdateFunctionsOfTime {
46 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
47 : typename ArrayIndex, typename ActionList,
48 : typename ParallelComponent>
49 0 : static Parallel::iterable_action_return_t apply(
50 : db::DataBox<DbTagsList>& box,
51 : tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
52 : Parallel::GlobalCache<Metavariables>& cache,
53 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
54 : const ParallelComponent* const /*meta*/) {
55 : const double current_expiration_time = db::get<Tags::ExpirationTime>(box);
56 : const double time = db::get<::Tags::Time>(box);
57 : const auto& particle_pos_vel =
58 : db::get<Tags::ParticlePositionVelocity<3>>(box);
59 : const double& x = get<0>(particle_pos_vel[0]);
60 : const double& y = get<1>(particle_pos_vel[0]);
61 : const double& xdot = get<0>(particle_pos_vel[1]);
62 : const double& ydot = get<1>(particle_pos_vel[1]);
63 : const double r = hypot(x, y);
64 : const double radial_vel = (xdot * x + ydot * y) / r;
65 : const auto& excision_sphere = db::get<Tags::ExcisionSphere<3>>(box);
66 : const double grid_radius_particle =
67 : get(magnitude(excision_sphere.center()));
68 :
69 : const DataVector angular_update{atan2(y, x),
70 : (x * ydot - y * xdot) / square(r)};
71 : const DataVector expansion_update{r / grid_radius_particle,
72 : radial_vel / grid_radius_particle};
73 :
74 : const auto& wt_radius_params =
75 : db::get<Tags::WorldtubeRadiusParameters>(box);
76 : const auto& bh_radius_params =
77 : db::get<Tags::BlackHoleRadiusParameters>(box);
78 : const double wt_radius_grid = excision_sphere.radius();
79 : const double wt_radius_inertial =
80 : smooth_broken_power_law(r, wt_radius_params[0], wt_radius_params[1],
81 : wt_radius_params[2], wt_radius_params[3]);
82 : const double wt_radius_derivative = smooth_broken_power_law_derivative(
83 : r, wt_radius_params[0], wt_radius_params[1], wt_radius_params[2],
84 : wt_radius_params[3]);
85 : const double wt_radius_time_derivative = wt_radius_derivative * radial_vel;
86 :
87 : const auto& bh_excision_sphere =
88 : db::get<domain::Tags::Domain<3>>(box).excision_spheres().at(
89 : "ExcisionSphereB");
90 : const double bh_excision_radius_grid = bh_excision_sphere.radius();
91 : const double bh_excision_radius_inertial =
92 : smooth_broken_power_law(r, bh_radius_params[0], bh_radius_params[1],
93 : bh_radius_params[2], bh_radius_params[3]);
94 : const double bh_excision_radius_derivative =
95 : smooth_broken_power_law_derivative(
96 : r, bh_radius_params[0], bh_radius_params[1], bh_radius_params[2],
97 : bh_radius_params[3]);
98 : const double bh_excision_radius_time_derivative =
99 : bh_excision_radius_derivative * radial_vel;
100 :
101 : const double sqrt_4_pi = sqrt(4. * M_PI);
102 : // The expansion map has already stretched the excision spheres and we need
103 : // to account for that.
104 : const DataVector size_a_update{
105 : sqrt_4_pi * (wt_radius_grid - wt_radius_inertial / expansion_update[0]),
106 : sqrt_4_pi *
107 : (wt_radius_inertial * expansion_update[1] -
108 : wt_radius_time_derivative * expansion_update[0]) /
109 : square(expansion_update[0])};
110 :
111 : const DataVector size_b_update{
112 : sqrt_4_pi * (bh_excision_radius_grid -
113 : bh_excision_radius_inertial / expansion_update[0]),
114 : sqrt_4_pi *
115 : (bh_excision_radius_inertial * expansion_update[1] -
116 : bh_excision_radius_time_derivative * expansion_update[0]) /
117 : square(expansion_update[0])};
118 :
119 : // we set the new expiration time to 1/100th of a time step next to the
120 : // current time. This is small enough that it can handle rapid time step
121 : // decreases but large enough to avoid floating point precision issues.
122 : const double new_expiration_time =
123 : time +
124 : 0.01 * (db::get<::Tags::Next<::Tags::TimeStepId>>(box).substep_time() -
125 : time);
126 :
127 : db::mutate<Tags::ExpirationTime>(
128 : [&new_expiration_time](const auto expiration_time) {
129 : *expiration_time = new_expiration_time;
130 : },
131 : make_not_null(&box));
132 : db::mutate<Tags::WorldtubeRadius>(
133 : [&wt_radius_inertial](const auto wt_radius) {
134 : *wt_radius = wt_radius_inertial;
135 : },
136 : make_not_null(&box));
137 : std::unordered_map<std::string, std::pair<DataVector, double>>
138 : all_updates{};
139 : all_updates["Rotation"] =
140 : std::make_pair(angular_update, new_expiration_time);
141 : all_updates["Expansion"] =
142 : std::make_pair(expansion_update, new_expiration_time);
143 : all_updates["SizeA"] = std::make_pair(size_a_update, new_expiration_time);
144 : all_updates["SizeB"] = std::make_pair(size_b_update, new_expiration_time);
145 :
146 : Parallel::mutate<::domain::Tags::FunctionsOfTime,
147 : control_system::UpdateMultipleFunctionsOfTime>(
148 : cache, current_expiration_time, all_updates);
149 :
150 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
151 : }
152 : };
153 : } // namespace CurvedScalarWave::Worldtube::Actions
|
