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