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
|
