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/TaggedTuple.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              : 
      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