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 <pup.h>
       8             : 
       9             : #include "ControlSystem/Measurements/BothHorizons.hpp"
      10             : #include "ControlSystem/Protocols/ControlError.hpp"
      11             : #include "ControlSystem/Tags/QueueTags.hpp"
      12             : #include "ControlSystem/Tags/SystemTags.hpp"
      13             : #include "DataStructures/DataBox/DataBox.hpp"
      14             : #include "DataStructures/DataVector.hpp"
      15             : #include "Domain/Creators/Tags/ObjectCenter.hpp"
      16             : #include "Domain/Structure/ObjectLabel.hpp"
      17             : #include "Options/String.hpp"
      18             : #include "Parallel/GlobalCache.hpp"
      19             : #include "Utilities/ErrorHandling/Assert.hpp"
      20             : #include "Utilities/ProtocolHelpers.hpp"
      21             : #include "Utilities/TMPL.hpp"
      22             : #include "Utilities/TaggedTuple.hpp"
      23             : 
      24             : /// \cond
      25             : namespace domain::Tags {
      26             : template <size_t Dim>
      27             : struct Domain;
      28             : struct FunctionsOfTime;
      29             : }  // namespace domain::Tags
      30             : /// \endcond
      31             : 
      32             : namespace control_system {
      33           1 : namespace ControlErrors {
      34             : /*!
      35             :  * \brief Control error in the 3D
      36             :  * \link domain::CoordinateMaps::TimeDependent::CubicScale CubicScale \endlink
      37             :  * coordinate map
      38             :  *
      39             :  * \details Computes the error in the map parameter \f$a(t)\f$ using Eq. (40)
      40             :  * from \cite Ossokine2013zga (see \link
      41             :  * domain::CoordinateMaps::TimeDependent::CubicScale CubicScale \endlink for a
      42             :  * definition of \f$a(t)\f$). The equation is
      43             :  *
      44             :  * \f{align}{
      45             :  *   \delta a &= a\left( \frac{\vec{X}\cdot\vec{C}}{||\vec{C}||^2} - 1 \right)
      46             :  * \\ \delta a &= a\left( \frac{X_0}{C_0} - 1 \right) \f}
      47             :  *
      48             :  * where \f$\vec{X} = \vec{x}_A - \vec{x}_B\f$ and \f$\vec{C} = \vec{c}_A -
      49             :  * \vec{c}_B\f$. Here, object A is located on the positive x-axis and object B
      50             :  * is located on the negative x-axis, \f$\vec{X}\f$ is the difference in
      51             :  * positions of the centers of the mapped objects, and
      52             :  * \f$\vec{C}\f$ is the difference of the centers of the excision spheres, all
      53             :  * in the grid frame. It is assumed that the positions of the excision spheres
      54             :  * are exactly along the x-axis, which is why we were able to make the
      55             :  * simplification in the second line above.
      56             :  *
      57             :  * Requirements:
      58             :  * - This control error requires that there be exactly two objects in the
      59             :  *   simulation
      60             :  * - Currently both these objects must be black holes
      61             :  * - Currently this control system can only be used with the \link
      62             :  *   control_system::Systems::Expansion Expansion \endlink control system
      63             :  */
      64           1 : struct Expansion : tt::ConformsTo<protocols::ControlError> {
      65           0 :   static constexpr size_t expected_number_of_excisions = 2;
      66             : 
      67           0 :   using object_centers =
      68             :       domain::object_list<domain::ObjectLabel::A, domain::ObjectLabel::B>;
      69             : 
      70           0 :   using options = tmpl::list<>;
      71           0 :   static constexpr Options::String help{
      72             :       "Computes the control error for expansion control. This should not "
      73             :       "take any options."};
      74             : 
      75           0 :   void pup(PUP::er& /*p*/) {}
      76             : 
      77             :   template <typename Metavariables, typename... TupleTags>
      78           0 :   DataVector operator()(const ::TimescaleTuner<true>& /*unused*/,
      79             :                         const Parallel::GlobalCache<Metavariables>& cache,
      80             :                         const double time,
      81             :                         const std::string& function_of_time_name,
      82             :                         const tuples::TaggedTuple<TupleTags...>& measurements) {
      83             :     const auto& functions_of_time = get<domain::Tags::FunctionsOfTime>(cache);
      84             : 
      85             :     const double current_expansion_factor =
      86             :         functions_of_time.at(function_of_time_name)->func(time)[0][0];
      87             : 
      88             :     using center_A =
      89             :         control_system::QueueTags::Center<::domain::ObjectLabel::A>;
      90             :     using center_B =
      91             :         control_system::QueueTags::Center<::domain::ObjectLabel::B>;
      92             : 
      93             :     const double grid_position_of_A =
      94             :         Parallel::get<domain::Tags::ObjectCenter<domain::ObjectLabel::A>>(
      95             :             cache)[0];
      96             :     const double grid_position_of_B =
      97             :         Parallel::get<domain::Tags::ObjectCenter<domain::ObjectLabel::B>>(
      98             :             cache)[0];
      99             :     const double current_position_of_A = get<center_A>(measurements)[0];
     100             :     const double current_position_of_B = get<center_B>(measurements)[0];
     101             : 
     102             :     // A is to the right of B in grid frame. To get positive differences,
     103             :     // take A - B
     104             :     const double expected_expansion_factor =
     105             :         current_expansion_factor *
     106             :         (current_position_of_A - current_position_of_B) /
     107             :         (grid_position_of_A - grid_position_of_B);
     108             : 
     109             :     return {expected_expansion_factor - current_expansion_factor};
     110             :   }
     111             : };
     112             : }  // namespace ControlErrors
     113             : }  // namespace control_system