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