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 :
8 : #include "ControlSystem/Averager.hpp"
9 : #include "ControlSystem/CalculateMeasurementTimescales.hpp"
10 : #include "ControlSystem/Component.hpp"
11 : #include "ControlSystem/ControlErrors/Size/Update.hpp"
12 : #include "ControlSystem/Controller.hpp"
13 : #include "ControlSystem/ExpirationTimes.hpp"
14 : #include "ControlSystem/IsSize.hpp"
15 : #include "ControlSystem/Metafunctions.hpp"
16 : #include "ControlSystem/Tags/IsActiveMap.hpp"
17 : #include "ControlSystem/Tags/MeasurementTimescales.hpp"
18 : #include "ControlSystem/Tags/SystemTags.hpp"
19 : #include "ControlSystem/TimescaleTuner.hpp"
20 : #include "ControlSystem/UpdateFunctionOfTime.hpp"
21 : #include "ControlSystem/UpdateTimescaleTuner.hpp"
22 : #include "ControlSystem/WriteData.hpp"
23 : #include "DataStructures/DataBox/DataBox.hpp"
24 : #include "DataStructures/DataVector.hpp"
25 : #include "DataStructures/LinkedMessageId.hpp"
26 : #include "DataStructures/TaggedTuple.hpp"
27 : #include "Domain/FunctionsOfTime/Tags.hpp"
28 : #include "Domain/Tags.hpp"
29 : #include "IO/Logging/Verbosity.hpp"
30 : #include "Parallel/GlobalCache.hpp"
31 : #include "Parallel/Printf/Printf.hpp"
32 : #include "Utilities/Gsl.hpp"
33 : #include "Utilities/MakeString.hpp"
34 : #include "Utilities/StdHelpers.hpp"
35 :
36 : /// \cond
37 : namespace domain::Tags {
38 : struct FunctionsOfTime;
39 : } // namespace domain::Tags
40 : /// \endcond
41 :
42 : namespace control_system {
43 : /*!
44 : * \brief Functor for updating control systems when they are ready.
45 : *
46 : * \details The apply operator of this struct is meant to be called by the
47 : * UpdateMessageQueue action once an entire measurement has been made.
48 : *
49 : * Requires a few tags to be in the DataBox of the ControlComponent that this
50 : * is running on:
51 : * - \link Tags::Averager Averager \endlink
52 : * - \link Tags::Controller Controller \endlink
53 : * - \link Tags::TimescaleTuner TimescaleTuner \endlink
54 : * - \link Tags::ControlError ControlError \endlink
55 : * - \link Tags::WriteDataToDisk WriteDataToDisk \endlink
56 : * - \link Tags::CurrentNumberOfMeasurements CurrentNumberOfMeasurements
57 : * \endlink
58 : *
59 : * And the \link control_system::Tags::MeasurementsPerUpdate \endlink must be in
60 : * the GlobalCache. If these tags are not present, a build error will occur.
61 : *
62 : * The algorithm to determine whether or not to update the functions of time is
63 : * as follows:
64 : * 1. Ensure this control system is active. If it isn't, end here and don't
65 : * process the measurements.
66 : * 2. Increment the current measurement stored by
67 : * `control_system::Tags::CurrentNumberOfMeasurements`. This keeps track of
68 : * which measurement we are on out of
69 : * `control_system::Tags::MeasurementsPerUpdate`.
70 : * 3. Calculate the control error and update the averager (store the current
71 : * measurement). If the averager doesn't have enough information to determine
72 : * the derivatives of the control error, exit now (this usually only happens
73 : * for the first couple measurements of a simulation). If this is \link
74 : * control_system::Systems::Size size \endlink control, an extra step happens
75 : * after we calculate the control error, but before we update the averager.
76 : * See `control_system::size::update_averager` for this step.
77 : * 4. If the `control_system::Tags::WriteDataToDisk` tag is set to `true`, write
78 : * the time, function of time values, and the control error and its
79 : * derivative to disk.
80 : * 5. Determine if we need to update. We only want to update when the current
81 : * measurement is equal to the number of measurements per update. If it's not
82 : * time to update, exit now. Once we determine that it is time to update, set
83 : * the current measurement back to 0.
84 : * 6. Compute the control signal using the control error and its derivatives and
85 : * update the damping timescales in the TimescaleTuner. If this control
86 : * system can suggest a damping timescale, there is an extra step after we
87 : * update the damping timescale. See `control_system::update_timescale_tuner`
88 : * for this step.
89 : * 7. Calculate the new measurement timescale based off the updated damping
90 : * timescales and the number of measurements per update.
91 : * 8. Determine the new expiration times for the
92 : * `::domain::Tags::FunctionsOfTime` and
93 : * `control_system::Tags::MeasurementTimescales`. Call the
94 : * `control_system::AggregateUpdate` simple action on the first control
95 : * system in the `component_list` of the metavariables. This simple action
96 : * will mutate the global cache tags when it has enough data.
97 : */
98 : template <typename ControlSystem>
99 1 : struct UpdateControlSystem {
100 0 : static constexpr size_t deriv_order = ControlSystem::deriv_order;
101 :
102 : template <typename DbTags, typename Metavariables, typename ArrayIndex,
103 : typename... TupleTags>
104 0 : static void apply(const gsl::not_null<db::DataBox<DbTags>*> box,
105 : Parallel::GlobalCache<Metavariables>& cache,
106 : const ArrayIndex& /*array_index*/, const double time,
107 : tuples::TaggedTuple<TupleTags...> data) {
108 : const std::string& function_of_time_name = ControlSystem::name();
109 :
110 : // Begin step 1
111 : // If this control system isn't active, don't do anything
112 : ASSERT(Parallel::get<control_system::Tags::IsActiveMap>(cache).contains(
113 : function_of_time_name),
114 : "Couldn't find function of time '"
115 : << function_of_time_name
116 : << "' in the control_system::Tags::IsActiveMap tag");
117 : if (not Parallel::get<control_system::Tags::IsActiveMap>(cache).at(
118 : function_of_time_name)) {
119 : return;
120 : }
121 :
122 : // Begin step 2
123 : const int measurements_per_update =
124 : get<control_system::Tags::MeasurementsPerUpdate>(cache);
125 : const bool delay_update = get<control_system::Tags::DelayUpdate>(cache);
126 : int& current_measurement = db::get_mutable_reference<
127 : control_system::Tags::CurrentNumberOfMeasurements>(box);
128 :
129 : ++current_measurement;
130 :
131 : const auto& functions_of_time =
132 : Parallel::get<::domain::Tags::FunctionsOfTime>(cache);
133 : const auto& function_of_time = functions_of_time.at(function_of_time_name);
134 :
135 : // Begin step 3
136 : // Get the averager, controller, tuner, and control error from the box
137 : auto& averager = db::get_mutable_reference<
138 : control_system::Tags::Averager<ControlSystem>>(box);
139 : auto& controller = db::get_mutable_reference<
140 : control_system::Tags::Controller<ControlSystem>>(box);
141 : auto& tuner = db::get_mutable_reference<
142 : control_system::Tags::TimescaleTuner<ControlSystem>>(box);
143 : auto& control_error = db::get_mutable_reference<
144 : control_system::Tags::ControlError<ControlSystem>>(box);
145 : const DataVector& current_timescale = tuner.current_timescale();
146 :
147 : // Compute control error
148 : const DataVector Q =
149 : control_error(tuner, cache, time, function_of_time_name, data);
150 :
151 : if (Parallel::get<Tags::Verbosity>(cache) >= ::Verbosity::Verbose) {
152 : Parallel::printf(
153 : "%s, time = %.16f: current measurement = %d, control error = %s\n",
154 : function_of_time_name, time, current_measurement, Q);
155 : }
156 :
157 : if constexpr (size::is_size_v<ControlSystem>) {
158 : // This function must be called after we update the control error because
159 : // it uses the new state in its logic (to repopulate the averager).
160 : size::update_averager(make_not_null(&averager),
161 : make_not_null(&control_error), cache, time,
162 : tuner.current_timescale(), function_of_time_name,
163 : current_measurement);
164 : }
165 :
166 : // Update the averager. We do this for every measurement because we still
167 : // want the averager to be up-to-date even if we aren't updating at this
168 : // time
169 : averager.update(time, Q, current_timescale);
170 :
171 : // Get the averaged values of the control error and its derivatives
172 : const auto& opt_avg_values = averager(time);
173 :
174 : if (not opt_avg_values.has_value()) {
175 : if (Parallel::get<Tags::Verbosity>(cache) >= ::Verbosity::Verbose) {
176 : Parallel::printf(
177 : "%s, time = %.16f: current measurement = %d, Averager does not "
178 : "have enough data.\n",
179 : function_of_time_name, time, current_measurement);
180 : }
181 : return;
182 : }
183 :
184 : // Begin step 4
185 : // Write data for every measurement
186 : std::array<DataVector, 2> q_and_dtq{{Q, {Q.size(), 0.0}}};
187 : q_and_dtq[0] = (*opt_avg_values)[0];
188 : if constexpr (deriv_order > 1) {
189 : q_and_dtq[1] = (*opt_avg_values)[1];
190 : }
191 :
192 : if (Parallel::get<control_system::Tags::WriteDataToDisk>(cache)) {
193 : // LCOV_EXCL_START
194 : write_components_to_disk<ControlSystem>(time, cache, function_of_time,
195 : q_and_dtq, current_timescale);
196 : // LCOV_EXCL_STOP
197 : }
198 :
199 : // Begin step 5
200 : // Check if it is time to update
201 : if (current_measurement != measurements_per_update) {
202 : if (Parallel::get<Tags::Verbosity>(cache) >= ::Verbosity::Verbose) {
203 : Parallel::printf(
204 : "%s, time = %.16f: current measurement = %d is not at "
205 : "measurements_per_update = %d. Waiting for more data.\n",
206 : function_of_time_name, time, current_measurement,
207 : measurements_per_update);
208 : }
209 : return;
210 : }
211 :
212 : // Set current measurement back to 0 because we're updating now
213 : current_measurement = 0;
214 :
215 : // Begin step 6
216 : const double time_offset =
217 : averager.last_time_updated() - averager.average_time(time);
218 : const double time_offset_0th =
219 : averager.using_average_0th_deriv_of_q() ? time_offset : 0.0;
220 :
221 : // Calculate the control signal which will be used to update the highest
222 : // derivative of the FunctionOfTime
223 : const DataVector control_signal =
224 : controller(time, current_timescale, opt_avg_values.value(),
225 : time_offset_0th, time_offset);
226 :
227 : tuner.update_timescale(q_and_dtq);
228 :
229 : update_timescale_tuner(make_not_null(&tuner), make_not_null(&control_error),
230 : Parallel::get<Tags::Verbosity>(cache), time,
231 : function_of_time_name);
232 :
233 : // Begin step 7
234 : // Calculate new measurement timescales with updated damping timescales
235 : const DataVector new_measurement_timescale =
236 : calculate_measurement_timescales(controller, tuner,
237 : measurements_per_update);
238 :
239 : const auto& measurement_timescales =
240 : Parallel::get<Tags::MeasurementTimescales>(cache);
241 : const auto& system_to_combined_names =
242 : Parallel::get<Tags::SystemToCombinedNames>(cache);
243 : const auto& measurement_timescale = measurement_timescales.at(
244 : system_to_combined_names.at(function_of_time_name));
245 : const double current_fot_expiration_time =
246 : function_of_time->time_bounds()[1];
247 : const double current_measurement_expiration_time =
248 : measurement_timescale->time_bounds()[1];
249 : const DataVector old_measurement_timescale =
250 : measurement_timescale->func(current_measurement_expiration_time)[0];
251 :
252 : // Begin step 8
253 : // Calculate the next expiration times for both the functions of time and
254 : // the measurement timescales based on the current time. Then, actually
255 : // update the functions of time and measurement timescales
256 : const double new_fot_expiration_time = function_of_time_expiration_time(
257 : time, old_measurement_timescale, new_measurement_timescale,
258 : measurements_per_update, delay_update);
259 :
260 : const double new_measurement_expiration_time = measurement_expiration_time(
261 : time, old_measurement_timescale, new_measurement_timescale,
262 : measurements_per_update, delay_update);
263 :
264 : if (Parallel::get<Tags::Verbosity>(cache) >= ::Verbosity::Verbose) {
265 : Parallel::printf("%s, time = %.16f: Control signal = %s\n",
266 : function_of_time_name, time, control_signal);
267 : Parallel::printf(
268 : "%s, time = %.16f: Updating the functions of time.\n"
269 : " min(old_measure_timescale) = %.16f\n"
270 : " min(new_measure_timescale) = %.16f\n"
271 : " old_measure_expr_time = %.16f\n"
272 : " new_measure_expr_time = %.16f\n"
273 : " old_function_expr_time = %.16f\n"
274 : " new_function_expr_time = %.16f\n",
275 : function_of_time_name, time, min(old_measurement_timescale),
276 : min(new_measurement_timescale), current_measurement_expiration_time,
277 : new_measurement_expiration_time, current_fot_expiration_time,
278 : new_fot_expiration_time);
279 : }
280 :
281 : using first_control_component =
282 : tmpl::front<metafunctions::all_control_components<Metavariables>>;
283 :
284 : auto& first_control_system_proxy =
285 : Parallel::get_parallel_component<first_control_component>(cache);
286 :
287 : Parallel::simple_action<AggregateUpdate<ControlSystem>>(
288 : first_control_system_proxy, new_measurement_timescale,
289 : current_measurement_expiration_time, new_measurement_expiration_time,
290 : control_signal, current_fot_expiration_time, new_fot_expiration_time);
291 : }
292 : };
293 : } // namespace control_system
|
