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