Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <algorithm>
7 : #include <memory>
8 : #include <optional>
9 :
10 : #include "ControlSystem/CombinedName.hpp"
11 : #include "ControlSystem/Metafunctions.hpp"
12 : #include "ControlSystem/Tags/IsActiveMap.hpp"
13 : #include "DataStructures/DataBox/DataBox.hpp"
14 : #include "Parallel/AlgorithmExecution.hpp"
15 : #include "Parallel/ArrayCollection/IsDgElementCollection.hpp"
16 : #include "Parallel/ArrayCollection/PerformAlgorithmOnElement.hpp"
17 : #include "Parallel/ArrayCollection/Tags/ElementLocations.hpp"
18 : #include "Parallel/ArrayComponentId.hpp"
19 : #include "Parallel/Callback.hpp"
20 : #include "Parallel/GlobalCache.hpp"
21 : #include "ParallelAlgorithms/Actions/GetItemFromDistributedObject.hpp"
22 : #include "Time/ChangeSlabSize/ChangeSlabSize.hpp"
23 : #include "Time/Tags/HistoryEvolvedVariables.hpp"
24 : #include "Time/Tags/MinimumTimeStep.hpp"
25 : #include "Utilities/ErrorHandling/Assert.hpp"
26 : #include "Utilities/ErrorHandling/Error.hpp"
27 : #include "Utilities/Gsl.hpp"
28 : #include "Utilities/TMPL.hpp"
29 :
30 : /// \cond
31 : namespace Tags {
32 : struct TimeStep;
33 : struct TimeStepId;
34 : template <typename StepperInterface>
35 : struct TimeStepper;
36 : } // namespace Tags
37 : class TimeStepper;
38 : namespace control_system::Tags {
39 : template <typename ControlSystems>
40 : struct FutureMeasurements;
41 : struct MeasurementTimescales;
42 : } // namespace control_system::Tags
43 : namespace domain::Tags {
44 : struct FunctionsOfTime;
45 : } // namespace domain::Tags
46 : namespace tuples {
47 : template <typename... Tags>
48 : class TaggedTuple;
49 : } // namespace tuples
50 : /// \endcond
51 :
52 : namespace control_system::Actions {
53 : /// \ingroup ControlSystemGroup
54 : /// \brief Limit the step size in a GTS evolution to prevent deadlocks from
55 : /// control system measurements.
56 : ///
57 : /// \details Most time steppers require evaluations of the coordinates
58 : /// at several times during the step before they can produce dense
59 : /// output. If any of those evaluations require a function-of-time
60 : /// update depending on a measurement within the step, the evolution
61 : /// will deadlock. This action reduces the step size if necessary to
62 : /// prevent that from happening.
63 : ///
64 : /// Specifically:
65 : /// 1. The chosen step will never be longer than the unmodified step,
66 : /// and will be short enough to avoid relevant function-of-time
67 : /// expirations.
68 : /// 2. Given the previous, the step will cover as many control-system
69 : /// updates as possible.
70 : /// 3. If the next step is likely to be limited by this action, adjust
71 : /// the length of the current step so that this step and the next
72 : /// step will be as close as possible to the same size.
73 : template <typename ControlSystems>
74 1 : struct LimitTimeStep {
75 : private:
76 0 : using control_system_groups =
77 : tmpl::transform<metafunctions::measurements_t<ControlSystems>,
78 : metafunctions::control_systems_with_measurement<
79 : tmpl::pin<ControlSystems>, tmpl::_1>>;
80 :
81 : template <typename Group>
82 0 : struct GroupExpiration {
83 0 : using type = double;
84 : };
85 :
86 : public:
87 0 : using const_global_cache_tags = tmpl::list<::Tags::MinimumTimeStep>;
88 :
89 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
90 : size_t Dim, typename ActionList, typename ParallelComponent>
91 0 : static Parallel::iterable_action_return_t apply(
92 : db::DataBox<DbTagsList>& box,
93 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
94 : Parallel::GlobalCache<Metavariables>& cache,
95 : const ElementId<Dim>& array_index, ActionList /*meta*/,
96 : const ParallelComponent* const /*meta*/) {
97 : const auto& time_step_id = db::get<::Tags::TimeStepId>(box);
98 : if (time_step_id.substep() != 0) {
99 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
100 : }
101 :
102 : const auto& time_stepper = db::get<::Tags::TimeStepper<TimeStepper>>(box);
103 : if (time_stepper.monotonic()) {
104 : // Monotonic steppers order operations in the same manner at the
105 : // control system, so they cannot introduce deadlocks.
106 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
107 : }
108 :
109 : auto& proxy = ::Parallel::get_parallel_component<ParallelComponent>(cache);
110 :
111 : // Minimum expiration time for any FoT in the measurement group.
112 : tmpl::wrap<tmpl::transform<control_system_groups,
113 : tmpl::bind<GroupExpiration, tmpl::_1>>,
114 : tuples::TaggedTuple>
115 : group_expiration_times{};
116 :
117 : bool ready = true;
118 : // Calculate group_expiration_times
119 : tmpl::for_each<control_system_groups>([&](auto group_v) {
120 : if (not ready) {
121 : return;
122 : }
123 : using group = tmpl::type_from<decltype(group_v)>;
124 :
125 : auto& future_measurements =
126 : db::get_mutable_reference<Tags::FutureMeasurements<group>>(
127 : make_not_null(&box));
128 :
129 : std::optional<double> group_update = future_measurements.next_update();
130 : if (not group_update.has_value()) {
131 : Parallel::mutable_cache_item_is_ready<
132 : control_system::Tags::MeasurementTimescales>(
133 : cache,
134 : Parallel::make_array_component_id<ParallelComponent>(array_index),
135 : [&](const auto& measurement_timescales) {
136 : const auto& group_timescale =
137 : *measurement_timescales.at(combined_name<group>());
138 : future_measurements.update(group_timescale);
139 : group_update = future_measurements.next_update();
140 : ready = group_update.has_value();
141 : if constexpr (Parallel::is_dg_element_collection_v<
142 : ParallelComponent>) {
143 : // Note: The ArrayComponentId is still created with the
144 : // array_index (ElementId) because we only support 1 callback
145 : // per ArrayComponentId. This would mean for nodegroups we would
146 : // discard a lot of callbacks that we need. Alternatively, the
147 : // callback could do a broadcast to all elements on the
148 : // nodegroup.
149 : const auto element_location = static_cast<int>(
150 : Parallel::local_synchronous_action<
151 : ::Parallel::Actions::GetItemFromDistributedOject<
152 : Parallel::Tags::ElementLocations<Dim>>>(proxy)
153 : ->at(array_index));
154 : return ready ? std::unique_ptr<Parallel::Callback>{}
155 : : std::unique_ptr<Parallel::Callback>(
156 : new Parallel::ThreadedActionCallback<
157 : Parallel::Actions::
158 : PerformAlgorithmOnElement<false>,
159 : decltype(proxy[element_location]),
160 : std::decay_t<decltype(array_index)>>{
161 : proxy[element_location], array_index});
162 : } else {
163 : return ready ? std::unique_ptr<Parallel::Callback>{}
164 : : std::unique_ptr<Parallel::Callback>(
165 : new Parallel::PerformAlgorithmCallback(
166 : proxy[array_index]));
167 : }
168 : });
169 : if (not ready) {
170 : return;
171 : }
172 : }
173 :
174 : auto& group_expiration =
175 : get<GroupExpiration<group>>(group_expiration_times);
176 : group_expiration = std::numeric_limits<double>::infinity();
177 :
178 : if (*group_update == std::numeric_limits<double>::infinity()) {
179 : // Control measurement is not active.
180 : return;
181 : }
182 :
183 : const auto& is_active_map =
184 : Parallel::get<control_system::Tags::IsActiveMap>(cache);
185 :
186 : // Calculate group_expiration
187 : Parallel::mutable_cache_item_is_ready<domain::Tags::FunctionsOfTime>(
188 : cache,
189 : Parallel::make_array_component_id<ParallelComponent>(array_index),
190 : [&](const auto& functions_of_time) {
191 : // Reset default readiness state to true because
192 : // mutable_cache_item_is_ready() calls this lambda twice. The first
193 : // call may have set ready = false, but FoT could have been updated
194 : // before the second call, so we default to true and then do a fresh
195 : // readiness check from the fresh FoT below.
196 : ready = true;
197 : tmpl::for_each<group>([&](auto system) {
198 : using System = tmpl::type_from<decltype(system)>;
199 : if (not ready or not is_active_map.at(System::name())) {
200 : return;
201 : }
202 : const auto& fot = *functions_of_time.at(System::name());
203 : ready = fot.time_bounds()[1] > *group_update;
204 : if (ready) {
205 : group_expiration = std::min(
206 : group_expiration, fot.expiration_after(*group_update));
207 : }
208 : });
209 : if constexpr (Parallel::is_dg_element_collection_v<
210 : ParallelComponent>) {
211 : // Note: The ArrayComponentId is still created with the
212 : // array_index (ElementId) because we only support 1 callback
213 : // per ArrayComponentId. This would mean for nodegroups we would
214 : // discard a lot of callbacks that we need. Alternatively, the
215 : // callback could do a broadcast to all elements on the
216 : // nodegroup.
217 : const auto element_location = static_cast<int>(
218 : Parallel::local_synchronous_action<
219 : ::Parallel::Actions::GetItemFromDistributedOject<
220 : Parallel::Tags::ElementLocations<Dim>>>(proxy)
221 : ->at(array_index));
222 : return ready
223 : ? std::unique_ptr<Parallel::Callback>{}
224 : : std::unique_ptr<Parallel::Callback>(
225 : new Parallel::ThreadedActionCallback<
226 : Parallel::Actions::PerformAlgorithmOnElement<
227 : false>,
228 : decltype(proxy[element_location]),
229 : std::decay_t<decltype(array_index)>>{
230 : proxy[element_location], array_index});
231 : } else {
232 : return ready ? std::unique_ptr<Parallel::Callback>{}
233 : : std::unique_ptr<Parallel::Callback>(
234 : new Parallel::PerformAlgorithmCallback(
235 : proxy[array_index]));
236 : }
237 : });
238 : });
239 : if (not ready) {
240 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
241 : }
242 :
243 : const double orig_step_start = time_step_id.step_time().value();
244 : const double orig_step_end =
245 : (time_step_id.step_time() + db::get<::Tags::TimeStep>(box)).value();
246 :
247 : // Smallest of the current step end and the FoT expirations. We
248 : // can't step any farther than this.
249 : const double latest_valid_step =
250 : tmpl::as_pack<control_system_groups>([&](auto... groups) {
251 : return std::min(
252 : {orig_step_end,
253 : get<GroupExpiration<tmpl::type_from<decltype(groups)>>>(
254 : group_expiration_times)...});
255 : });
256 :
257 : if (not tmpl::as_pack<::Tags::get_all_history_tags<DbTagsList>>(
258 : [&]<typename... HistoryTags>(tmpl::type_<HistoryTags>... /*meta*/) {
259 : return (... and time_stepper.can_change_step_size(
260 : time_step_id, db::get<HistoryTags>(box)));
261 : })) {
262 : if (orig_step_end > latest_valid_step) {
263 : ERROR(
264 : "Step must be decreased to avoid control-system deadlock, but "
265 : "time-stepper requires a fixed step size.");
266 : }
267 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
268 : }
269 : ASSERT(db::get<::Tags::TimeStep>(box).fraction() == 1,
270 : "Trying to change GTS step, but it isn't a full slab. Non-slab "
271 : "steps should only happen during self-start, but the preceding "
272 : "check should have ended the action if this is self-start.");
273 :
274 : // The last update that we can perform on the next step. Don't
275 : // shrink the step past this time since that will force another
276 : // step to take the measurement.
277 : double last_update_time = orig_step_start;
278 : // Step time that produces a balanced step with the following
279 : // step, ignoring the restrictions on the current step.
280 : double preferred_step_time = orig_step_end;
281 :
282 : tmpl::for_each<control_system_groups>([&](auto group_v) {
283 : using group = tmpl::type_from<decltype(group_v)>;
284 :
285 : // This was used above, so it is not nullopt.
286 : const double group_update =
287 : db::get<Tags::FutureMeasurements<group>>(box).next_update().value();
288 : if (group_update <= latest_valid_step) {
289 : // We've satisfied this measurement.
290 : last_update_time = std::max(last_update_time, group_update);
291 : } else {
292 : // We can't make it far enough to do the final measurement.
293 : // Try to avoid a small step by choosing two equal-sized steps
294 : // to the expiration time.
295 : const double equal_step_time =
296 : 0.5 * (orig_step_start +
297 : get<GroupExpiration<group>>(group_expiration_times));
298 : preferred_step_time = std::min(preferred_step_time, equal_step_time);
299 : }
300 : });
301 :
302 : const double new_step_end =
303 : std::clamp(preferred_step_time, last_update_time, latest_valid_step);
304 :
305 : change_slab_size(make_not_null(&box), new_step_end);
306 :
307 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
308 : }
309 : };
310 :
311 : /// \ingroup ControlSystemGroup
312 : /// \brief No-control-system specialization that does nothing.
313 : ///
314 : /// Exists just to avoid having to put conditionals in metavariables
315 : /// that may or may not have control systems.
316 : template <>
317 1 : struct LimitTimeStep<tmpl::list<>> {
318 : public:
319 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
320 : size_t Dim, typename ActionList, typename ParallelComponent>
321 0 : static Parallel::iterable_action_return_t apply(
322 : db::DataBox<DbTagsList>& /*box*/,
323 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
324 : Parallel::GlobalCache<Metavariables>& /*cache*/,
325 : const ElementId<Dim>& /*array_index*/, ActionList /*meta*/,
326 : const ParallelComponent* const /*meta*/) {
327 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
328 : }
329 : };
330 : } // namespace control_system::Actions
|
