SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions - UpdateFunctionsOfTime.hpp Hit Total Coverage
Commit: bb4d647fd2cc14353206e6e15cae7930cab6f182 Lines: 1 3 33.3 %
Date: 2025-02-05 08:19:47
Legend: Lines: hit not hit

          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

Generated by: LCOV version 1.14