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 "DataStructures/TaggedTuple.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            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 :   using object_centers =
      66              :       domain::object_list<domain::ObjectLabel::A, domain::ObjectLabel::B>;
      67              : 
      68            0 :   using options = tmpl::list<>;
      69            0 :   static constexpr Options::String help{
      70              :       "Computes the control error for expansion control. This should not "
      71              :       "take any options."};
      72              : 
      73              :   // NOLINTNEXTLINE(readability-convert-member-functions-to-static)
      74            0 :   std::optional<double> get_suggested_timescale() const { return std::nullopt; }
      75              : 
      76            0 :   void reset() {}
      77              : 
      78            0 :   void pup(PUP::er& /*p*/) {}
      79              : 
      80              :   template <typename Metavariables, typename... TupleTags>
      81            0 :   DataVector operator()(const ::TimescaleTuner<true>& /*unused*/,
      82              :                         const Parallel::GlobalCache<Metavariables>& cache,
      83              :                         const double time,
      84              :                         const std::string& function_of_time_name,
      85              :                         const tuples::TaggedTuple<TupleTags...>& measurements) {
      86              :     const auto& functions_of_time = get<domain::Tags::FunctionsOfTime>(cache);
      87              : 
      88              :     const double current_expansion_factor =
      89              :         functions_of_time.at(function_of_time_name)->func(time)[0][0];
      90              : 
      91              :     using center_A = control_system::QueueTags::Center<::domain::ObjectLabel::A,
      92              :                                                        Frame::Grid>;
      93              :     using center_B = control_system::QueueTags::Center<::domain::ObjectLabel::B,
      94              :                                                        Frame::Grid>;
      95              : 
      96              :     const double grid_position_of_A =
      97              :         Parallel::get<domain::Tags::ObjectCenter<domain::ObjectLabel::A>>(
      98              :             cache)[0];
      99              :     const double grid_position_of_B =
     100              :         Parallel::get<domain::Tags::ObjectCenter<domain::ObjectLabel::B>>(
     101              :             cache)[0];
     102              :     const double current_position_of_A = get<center_A>(measurements)[0];
     103              :     const double current_position_of_B = get<center_B>(measurements)[0];
     104              : 
     105              :     // A is to the right of B in grid frame. To get positive differences,
     106              :     // take A - B
     107              :     const double expected_expansion_factor =
     108              :         current_expansion_factor *
     109              :         (current_position_of_A - current_position_of_B) /
     110              :         (grid_position_of_A - grid_position_of_B);
     111              : 
     112              :     return {expected_expansion_factor - current_expansion_factor};
     113              :   }
     114              : };
     115              : }  // namespace ControlErrors
     116              : }  // namespace control_system