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