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/DataVectorHelpers.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 "DataStructures/Tensor/Tensor.hpp"
      16             : #include "Domain/Creators/Tags/ObjectCenter.hpp"
      17             : #include "Domain/Structure/ObjectLabel.hpp"
      18             : #include "Options/String.hpp"
      19             : #include "Parallel/GlobalCache.hpp"
      20             : #include "Utilities/ErrorHandling/Assert.hpp"
      21             : #include "Utilities/ProtocolHelpers.hpp"
      22             : #include "Utilities/TMPL.hpp"
      23             : #include "Utilities/TaggedTuple.hpp"
      24             : 
      25             : /// \cond
      26             : namespace domain::Tags {
      27             : template <size_t Dim>
      28             : struct Domain;
      29             : struct FunctionsOfTime;
      30             : }  // namespace domain::Tags
      31             : /// \endcond
      32             : 
      33             : namespace control_system {
      34             : namespace ControlErrors {
      35             : /*!
      36             :  * \brief Control error in the 3D
      37             :  * \link domain::CoordinateMaps::TimeDependent::Rotation Rotation \endlink
      38             :  * coordinate map
      39             :  *
      40             :  * \details Computes the error in the angle rotated by the system using a
      41             :  * slightly modified version of Eq. (41) from \cite Ossokine2013zga. The
      42             :  * equation is
      43             :  *
      44             :  * \f[ \delta\vec{q} = \frac{\vec{C}\times\vec{X}}{\vec{C}\cdot\vec{X}} \f]
      45             :  *
      46             :  * where \f$\vec{X} = \vec{x}_A - \vec{x}_B\f$ and \f$\vec{C} = \vec{c}_A -
      47             :  * \vec{c}_B\f$. Here, object A is located on the positive x-axis and object B
      48             :  * is located on the negative x-axis, \f$\vec{X}\f$ is the difference in
      49             :  * positions of the centers of the mapped objects, and
      50             :  * \f$\vec{C}\f$ is the difference of the centers of the excision spheres, all
      51             :  * in the grid frame.
      52             :  *
      53             :  * Requirements:
      54             :  * - This control error requires that there be exactly two objects in the
      55             :  *   simulation
      56             :  * - Currently both these objects must be black holes
      57             :  * - Currently this control system can only be used with the \link
      58             :  *   control_system::Systems::Rotation Rotation \endlink control system
      59             :  */
      60           1 : struct Rotation : tt::ConformsTo<protocols::ControlError> {
      61           0 :   static constexpr size_t expected_number_of_excisions = 2;
      62             : 
      63           0 :   using object_centers =
      64             :       domain::object_list<domain::ObjectLabel::A, domain::ObjectLabel::B>;
      65             : 
      66           0 :   using options = tmpl::list<>;
      67           0 :   static constexpr Options::String help{
      68             :       "Computes the control error for rotation control. This should not "
      69             :       "take any options."};
      70             : 
      71           0 :   void pup(PUP::er& /*p*/) {}
      72             : 
      73             :   template <typename Metavariables, typename... TupleTags>
      74           0 :   DataVector operator()(const ::TimescaleTuner<true>& /*unused*/,
      75             :                         const Parallel::GlobalCache<Metavariables>& cache,
      76             :                         const double /*time*/,
      77             :                         const std::string& /*function_of_time_name*/,
      78             :                         const tuples::TaggedTuple<TupleTags...>& measurements) {
      79             :     using center_A =
      80             :         control_system::QueueTags::Center<::domain::ObjectLabel::A>;
      81             :     using center_B =
      82             :         control_system::QueueTags::Center<::domain::ObjectLabel::B>;
      83             : 
      84             :     const tnsr::I<double, 3, Frame::Grid>& grid_position_of_A =
      85             :         Parallel::get<domain::Tags::ObjectCenter<domain::ObjectLabel::A>>(
      86             :             cache);
      87             :     const tnsr::I<double, 3, Frame::Grid>& grid_position_of_B =
      88             :         Parallel::get<domain::Tags::ObjectCenter<domain::ObjectLabel::B>>(
      89             :             cache);
      90             :     const DataVector& current_position_of_A = get<center_A>(measurements);
      91             :     const DataVector& current_position_of_B = get<center_B>(measurements);
      92             : 
      93             :     // A is to the right of B in grid frame. To get positive differences,
      94             :     // take A - B
      95             :     const auto grid_diff_tnsr = tenex::evaluate<ti::I>(
      96             :         grid_position_of_A(ti::I) - grid_position_of_B(ti::I));
      97             :     const DataVector grid_diff{
      98             :         {grid_diff_tnsr[0], grid_diff_tnsr[1], grid_diff_tnsr[2]}};
      99             :     const DataVector current_diff =
     100             :         current_position_of_A - current_position_of_B;
     101             : 
     102             :     const double grid_dot_current = dot(grid_diff, current_diff);
     103             :     const DataVector grid_cross_current = cross(grid_diff, current_diff);
     104             : 
     105             :     return grid_cross_current / grid_dot_current;
     106             :   }
     107             : };
     108             : }  // namespace ControlErrors
     109             : }  // namespace control_system