Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cmath>
       7             : #include <limits>
       8             : #include <memory>
       9             : 
      10             : #include "ControlSystem/Tags/IsActiveMap.hpp"
      11             : #include "ControlSystem/Tags/OptionTags.hpp"
      12             : #include "DataStructures/DataBox/DataBox.hpp"
      13             : #include "DataStructures/DataBox/Tag.hpp"
      14             : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
      15             : #include "Domain/FunctionsOfTime/SettleToConstantQuaternion.hpp"
      16             : #include "Domain/FunctionsOfTime/Tags.hpp"
      17             : #include "Domain/Structure/ElementId.hpp"
      18             : #include "Domain/Tags.hpp"
      19             : #include "Parallel/AlgorithmExecution.hpp"
      20             : #include "Parallel/GlobalCache.hpp"
      21             : #include "ParallelAlgorithms/EventsAndTriggers/EventsAndTriggers.hpp"
      22             : #include "ParallelAlgorithms/EventsAndTriggers/Tags.hpp"
      23             : #include "ParallelAlgorithms/EventsAndTriggers/WhenToCheck.hpp"
      24             : #include "Time/Tags/Time.hpp"
      25             : #include "Time/Tags/TimeStepId.hpp"
      26             : #include "Utilities/ConstantExpressions.hpp"
      27             : #include "Utilities/ErrorHandling/Error.hpp"
      28             : 
      29             : /// \cond
      30             : namespace Parallel {
      31             : template <typename Metavariables>
      32             : class GlobalCache;
      33             : }  // namespace Parallel
      34             : namespace Tags {
      35             : struct TimeStep;
      36             : struct TimeStepId;
      37             : template <typename StepperInterface>
      38             : struct TimeStepper;
      39             : }  // namespace Tags
      40             : namespace tuples {
      41             : template <class... Tags>
      42             : class TaggedTuple;
      43             : }  // namespace tuples
      44             : /// \endcond
      45             : 
      46           1 : namespace control_system {
      47             : /*!
      48             :  * \brief Holds options used to control when to start disabling the rotation
      49             :  * control system in a BNS simulation.
      50             :  */
      51           1 : struct DisableRotationWhen {
      52             :   /// The separation at which we start turning off rotation.
      53           1 :   struct DisableAtSeparation {
      54           0 :     using type = double;
      55           0 :     static type lower_bound() { return 2.0; }
      56           0 :     static constexpr Options::String help{
      57             :         "The separation at which we start turning off rotation."};
      58             :   };
      59             : 
      60             :   /// The timescale in code units over which grid rotation is disabled. A
      61             :   /// reasonable value is 30M-60M.
      62           1 :   struct RotationDecayTimescale {
      63           0 :     using type = double;
      64           0 :     static type lower_bound() { return 30.0; }
      65           0 :     static type upper_bound() { return 200.0; }
      66           0 :     static constexpr Options::String help{
      67             :         "The timescale in code units over which grid rotation is disabled. A "
      68             :         "reasonable value is 40M-60M."};
      69             :   };
      70             : 
      71           0 :   using options = tmpl::list<DisableAtSeparation, RotationDecayTimescale>;
      72             : 
      73           0 :   static constexpr Options::String help{
      74             :       "Constrols the separation at which the rotation control system is "
      75             :       "disabled and the rotation function of time starts settling to a "
      76             :       "constant."};
      77             : 
      78           0 :   void pup(PUP::er& p);
      79             : 
      80           0 :   double disable_at_separation{std::numeric_limits<double>::signaling_NaN()};
      81           0 :   double rotation_decay_timescale{std::numeric_limits<double>::signaling_NaN()};
      82             : };
      83             : 
      84           1 : namespace OptionTags {
      85             : /// Option tag for controlling when and how the rotation map is disabled.
      86           1 : struct DisableRotationWhen {
      87           0 :   static constexpr Options::String help =
      88             :       "Options for controlling how the rotation control system stops as the "
      89             :       "two neutron stars merge.";
      90           0 :   using type = control_system::DisableRotationWhen;
      91           0 :   using group = control_system::OptionTags::ControlSystemGroup;
      92             : };
      93             : }  // namespace OptionTags
      94             : 
      95           1 : namespace Tags {
      96             : /// Tag for controlling when and how the rotation map is disabled.
      97           1 : struct DisableRotationWhen : db::SimpleTag {
      98           0 :   using type = control_system::DisableRotationWhen;
      99           0 :   using option_tags = tmpl::list<OptionTags::DisableRotationWhen>;
     100             : 
     101           0 :   static constexpr bool pass_metavariables = false;
     102           0 :   static control_system::DisableRotationWhen create_from_options(
     103             :       const control_system::DisableRotationWhen& disable_rotation_when);
     104             : };
     105             : }  // namespace Tags
     106             : 
     107           1 : namespace Actions {
     108             : /*!
     109             :  * \brief Checks if the Rotation function of time has been updated because the
     110             :  * separation between the neutron star grid centers is small enough.
     111             :  *
     112             :  * \note This is an iterable action that is to be run on the elements. It is
     113             :  * only run on the element that satisfies `is_zeroth_element()`.
     114             :  *
     115             :  * \warning This action should only ever be run in the
     116             :  * `Parallel::Phase::DisableRotationControl` phase.
     117             :  *
     118             :  * The desired separation is controlled via the
     119             :  * control_systems::Tags::DisableRotationWhen tag. The main use for this
     120             :  * functionality is to disable the rotation control system in binary neutron
     121             :  * star mergers when the stars are sufficiently close because as the stars
     122             :  * merge, the dual control system starts to lose anything to lock on to and can
     123             :  * even start counter rotating.
     124             :  *
     125             :  * Tags used and modified:
     126             :  * - DataBox:
     127             :  *   - `domain::Tags::Element<3>`
     128             :  *   - `::Tags::TimeStepId`
     129             :  *   - `::Tags::Time`
     130             :  * - MutableGlobalCache:
     131             :  *  - Uses:
     132             :  *   - `domain::Tags::FunctionsOfTime` ("GridCenters" and "Rotation")
     133             :  *  - Modifies:
     134             :  *    - `domain::Tags::FunctionsOfTime` ("Rotation")
     135             :  *    - `control_system::Tags::IsActiveMap` ("Rotation")
     136             :  */
     137           1 : struct SwitchGridRotationToSettle {
     138             :   /// Invokable that changes the rotation function of time to a
     139             :   /// QuaternionSettleToConstant matching the current function of time values
     140             :   /// and derivatives for the rotation and decaying over a specified timescale.
     141             :   ///
     142             :   /// Note that the final orientation of the domain is arbitrary. If we want
     143             :   /// to rotate to a specific angle (e.g. aligning logical and inertial
     144             :   /// coordinate axes), we need a (new) QuaternionSettleToSpecifiedValue
     145             :   /// function of time.
     146           1 :   struct UpdateRotationToSettle {
     147           0 :     static void apply(
     148             :         gsl::not_null<std::unordered_map<
     149             :             std::string,
     150             :             std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>*>
     151             :             f_of_t_list,
     152             :         const std::string& function_of_time_name,
     153             :         const std::array<DataVector, 3>& initial_func_and_derivs,
     154             :         double match_time, double decay_time);
     155             :   };
     156             : 
     157             :   /// Invokable used to disable the rotation control system via a call to
     158             :   /// Parallel::mutate.
     159           1 :   struct DisableControlSystem {
     160           0 :     static void apply(
     161             :         gsl::not_null<std::unordered_map<std::string, bool>*> is_active_map,
     162             :         const std::string& control_system_name);
     163             :   };
     164             : 
     165           0 :   using const_global_cache_tags =
     166             :       tmpl::list<control_system::Tags::DisableRotationWhen>;
     167             : 
     168             :   template <typename DbTags, typename... InboxTags, typename Metavariables,
     169             :             typename ActionList, typename ParallelComponent>
     170           0 :   static Parallel::iterable_action_return_t apply(
     171             :       db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
     172             :       Parallel::GlobalCache<Metavariables>& cache,
     173             :       const ElementId<3>& element_id, const ActionList /*meta*/,
     174             :       const ParallelComponent* const /*meta*/) {
     175             :     if (not is_zeroth_element(element_id)) {
     176             :       // Only one element needs to modify the control system and function of
     177             :       // time.
     178             :       return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     179             :     }
     180             :     const auto& time_step_id = db::get<::Tags::TimeStepId>(box);
     181             :     if (not time_step_id.is_at_slab_boundary()) {
     182             :       ERROR(
     183             :           "Expected to be at a Slab boundary when changing the Rotation "
     184             :           "function of time to a SettleToConstant. Current TimeStepId is "
     185             :           << time_step_id);
     186             :     }
     187             : 
     188             :     const auto& bns_rotation_control =
     189             :         db::get<control_system::Tags::DisableRotationWhen>(box);
     190             : 
     191             :     const double time = db::get<::Tags::Time>(box);
     192             :     if (not get<domain::Tags::FunctionsOfTime>(cache).contains("GridCenters")) {
     193             :       ERROR(
     194             :           "There is no function of time named 'GridCenters', which is required "
     195             :           "when disabling the rotation control system since in a binary "
     196             :           "neutron star simulation we need to track the grid centers of the "
     197             :           "stars to decide when to disable rotation control.");
     198             :     }
     199             :     const domain::FunctionsOfTime::FunctionOfTime& grid_centers_fot =
     200             :         *get<domain::Tags::FunctionsOfTime>(cache).at("GridCenters");
     201             :     const DataVector grid_centers = grid_centers_fot.func(time)[0];
     202             :     const double separation = sqrt(square(grid_centers[0] - grid_centers[3]) +
     203             :                                    square(grid_centers[1] - grid_centers[4]) +
     204             :                                    square(grid_centers[2] - grid_centers[5]));
     205             :     if (separation > bns_rotation_control.disable_at_separation) {
     206             :       ERROR(
     207             :           "Disabling the rotation control system should happen when the "
     208             :           "separation is less than or equal to "
     209             :           << bns_rotation_control.disable_at_separation
     210             :           << " but the separation at time " << time << " is calculated to be "
     211             :           << separation);
     212             :     }
     213             : 
     214             :     if (not get<domain::Tags::FunctionsOfTime>(cache).contains("Rotation")) {
     215             :       ERROR(
     216             :           "There is no function of time named 'Rotation', which means that it "
     217             :           "cannot be disabled.");
     218             :     }
     219             :     const domain::FunctionsOfTime::FunctionOfTime& rotation_fot =
     220             :         *get<domain::Tags::FunctionsOfTime>(cache).at("Rotation");
     221             :     const std::array<DataVector, 3> current_func_and_derivs =
     222             :         rotation_fot.func_and_2_derivs(time);
     223             :     Parallel::mutate<domain::Tags::FunctionsOfTime, UpdateRotationToSettle>(
     224             :         cache, std::string{"Rotation"}, current_func_and_derivs, time,
     225             :         bns_rotation_control.rotation_decay_timescale);
     226             :     Parallel::mutate<control_system::Tags::IsActiveMap, DisableControlSystem>(
     227             :         cache, std::string{"Rotation"});
     228             :     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     229             :   }
     230             : };
     231             : }  // namespace Actions
     232             : }  // namespace control_system