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 "Domain/FunctionsOfTime/Tags.hpp"
27 : #include "Domain/Tags.hpp"
28 : #include "IO/Logging/Verbosity.hpp"
29 : #include "Parallel/GlobalCache.hpp"
30 : #include "Parallel/Printf/Printf.hpp"
31 : #include "Utilities/Gsl.hpp"
32 : #include "Utilities/MakeString.hpp"
33 : #include "Utilities/StdHelpers.hpp"
34 : #include "Utilities/TaggedTuple.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 : int& current_measurement = db::get_mutable_reference<
126 : control_system::Tags::CurrentNumberOfMeasurements>(box);
127 :
128 : ++current_measurement;
129 :
130 : const auto& functions_of_time =
131 : Parallel::get<::domain::Tags::FunctionsOfTime>(cache);
132 : const auto& function_of_time = functions_of_time.at(function_of_time_name);
133 :
134 : // Begin step 3
135 : // Get the averager, controller, tuner, and control error from the box
136 : auto& averager = db::get_mutable_reference<
137 : control_system::Tags::Averager<ControlSystem>>(box);
138 : auto& controller = db::get_mutable_reference<
139 : control_system::Tags::Controller<ControlSystem>>(box);
140 : auto& tuner = db::get_mutable_reference<
141 : control_system::Tags::TimescaleTuner<ControlSystem>>(box);
142 : auto& control_error = db::get_mutable_reference<
143 : control_system::Tags::ControlError<ControlSystem>>(box);
144 : const DataVector& current_timescale = tuner.current_timescale();
145 :
146 : // Compute control error
147 : const DataVector Q =
148 : control_error(tuner, cache, time, function_of_time_name, data);
149 :
150 : if (Parallel::get<Tags::Verbosity>(cache) >= ::Verbosity::Verbose) {
151 : Parallel::printf(
152 : "%s, time = %.16f: current measurement = %d, control error = %s\n",
153 : function_of_time_name, time, current_measurement, Q);
154 : }
155 :
156 : if constexpr (size::is_size_v<ControlSystem>) {
157 : // This function must be called after we update the control error because
158 : // it uses the new state in its logic (to repopulate the averager).
159 : size::update_averager(make_not_null(&averager),
160 : make_not_null(&control_error), cache, time,
161 : tuner.current_timescale(), function_of_time_name,
162 : current_measurement);
163 : }
164 :
165 : // Update the averager. We do this for every measurement because we still
166 : // want the averager to be up-to-date even if we aren't updating at this
167 : // time
168 : averager.update(time, Q, current_timescale);
169 :
170 : // Get the averaged values of the control error and its derivatives
171 : const auto& opt_avg_values = averager(time);
172 :
173 : if (not opt_avg_values.has_value()) {
174 : if (Parallel::get<Tags::Verbosity>(cache) >= ::Verbosity::Verbose) {
175 : Parallel::printf(
176 : "%s, time = %.16f: current measurement = %d, Averager does not "
177 : "have enough data.\n",
178 : function_of_time_name, time, current_measurement);
179 : }
180 : return;
181 : }
182 :
183 : // Begin step 4
184 : // Write data for every measurement
185 : std::array<DataVector, 2> q_and_dtq{{Q, {Q.size(), 0.0}}};
186 : q_and_dtq[0] = (*opt_avg_values)[0];
187 : if constexpr (deriv_order > 1) {
188 : q_and_dtq[1] = (*opt_avg_values)[1];
189 : }
190 :
191 : if (Parallel::get<control_system::Tags::WriteDataToDisk>(cache)) {
192 : // LCOV_EXCL_START
193 : write_components_to_disk<ControlSystem>(time, cache, function_of_time,
194 : q_and_dtq, current_timescale);
195 : // LCOV_EXCL_STOP
196 : }
197 :
198 : // Begin step 5
199 : // Check if it is time to update
200 : if (current_measurement != measurements_per_update) {
201 : if (Parallel::get<Tags::Verbosity>(cache) >= ::Verbosity::Verbose) {
202 : Parallel::printf(
203 : "%s, time = %.16f: current measurement = %d is not at "
204 : "measurements_per_update = %d. Waiting for more data.\n",
205 : function_of_time_name, time, current_measurement,
206 : measurements_per_update);
207 : }
208 : return;
209 : }
210 :
211 : // Set current measurement back to 0 because we're updating now
212 : current_measurement = 0;
213 :
214 : // Begin step 6
215 : const double time_offset =
216 : averager.last_time_updated() - averager.average_time(time);
217 : const double time_offset_0th =
218 : averager.using_average_0th_deriv_of_q() ? time_offset : 0.0;
219 :
220 : // Calculate the control signal which will be used to update the highest
221 : // derivative of the FunctionOfTime
222 : const DataVector control_signal =
223 : controller(time, current_timescale, opt_avg_values.value(),
224 : time_offset_0th, time_offset);
225 :
226 : tuner.update_timescale(q_and_dtq);
227 :
228 : update_timescale_tuner(make_not_null(&tuner), make_not_null(&control_error),
229 : Parallel::get<Tags::Verbosity>(cache), time,
230 : function_of_time_name);
231 :
232 : // Begin step 7
233 : // Calculate new measurement timescales with updated damping timescales
234 : const DataVector new_measurement_timescale =
235 : calculate_measurement_timescales(controller, tuner,
236 : measurements_per_update);
237 :
238 : const auto& measurement_timescales =
239 : Parallel::get<Tags::MeasurementTimescales>(cache);
240 : const auto& system_to_combined_names =
241 : Parallel::get<Tags::SystemToCombinedNames>(cache);
242 : const auto& measurement_timescale = measurement_timescales.at(
243 : system_to_combined_names.at(function_of_time_name));
244 : const double current_fot_expiration_time =
245 : function_of_time->time_bounds()[1];
246 : const double current_measurement_expiration_time =
247 : measurement_timescale->time_bounds()[1];
248 : // This call is ok because the measurement timescales are still valid
249 : // because the measurement timescales expire half a measurement after this
250 : // time.
251 : const DataVector old_measurement_timescale =
252 : measurement_timescale->func(time)[0];
253 :
254 : // Begin step 8
255 : // Calculate the next expiration times for both the functions of time and
256 : // the measurement timescales based on the current time. Then, actually
257 : // update the functions of time and measurement timescales
258 : const double new_fot_expiration_time = function_of_time_expiration_time(
259 : time, old_measurement_timescale, new_measurement_timescale,
260 : measurements_per_update);
261 :
262 : const double new_measurement_expiration_time = measurement_expiration_time(
263 : time, old_measurement_timescale, new_measurement_timescale,
264 : measurements_per_update);
265 :
266 : if (Parallel::get<Tags::Verbosity>(cache) >= ::Verbosity::Verbose) {
267 : Parallel::printf("%s, time = %.16f: Control signal = %s\n",
268 : function_of_time_name, time, control_signal);
269 : Parallel::printf(
270 : "%s, time = %.16f: Updating the functions of time.\n"
271 : " min(old_measure_timescale) = %.16f\n"
272 : " min(new_measure_timescale) = %.16f\n"
273 : " old_measure_expr_time = %.16f\n"
274 : " new_measure_expr_time = %.16f\n"
275 : " old_function_expr_time = %.16f\n"
276 : " new_function_expr_time = %.16f\n",
277 : function_of_time_name, time, min(old_measurement_timescale),
278 : min(new_measurement_timescale), current_measurement_expiration_time,
279 : new_measurement_expiration_time, current_fot_expiration_time,
280 : new_fot_expiration_time);
281 : }
282 :
283 : using first_control_component =
284 : tmpl::front<metafunctions::all_control_components<Metavariables>>;
285 :
286 : auto& first_control_system_proxy =
287 : Parallel::get_parallel_component<first_control_component>(cache);
288 :
289 : Parallel::simple_action<AggregateUpdate<ControlSystem>>(
290 : first_control_system_proxy, new_measurement_timescale,
291 : current_measurement_expiration_time, new_measurement_expiration_time,
292 : control_signal, current_fot_expiration_time, new_fot_expiration_time);
293 : }
294 : };
295 : } // namespace control_system
|
