Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <map>
8 : #include <optional>
9 : #include <set>
10 : #include <utility>
11 : #include <vector>
12 :
13 : #include "DataStructures/DataBox/DataBox.hpp"
14 : #include "DataStructures/DataBox/Protocols/Mutator.hpp"
15 : #include "DataStructures/LinkedMessageId.hpp"
16 : #include "Evolution/Systems/GeneralizedHarmonic/Bbh/CompletionCriteria.hpp"
17 : #include "Parallel/Algorithms/AlgorithmSingletonDeclarations.hpp"
18 : #include "Parallel/ArrayCollection/IsDgElementCollection.hpp"
19 : #include "Parallel/ArrayCollection/SimpleActionOnElement.hpp"
20 : #include "Parallel/GlobalCache.hpp"
21 : #include "Parallel/Info.hpp"
22 : #include "Parallel/Invoke.hpp"
23 : #include "Parallel/Local.hpp"
24 : #include "Parallel/Phase.hpp"
25 : #include "Parallel/PhaseDependentActionList.hpp"
26 : #include "Parallel/Printf/Printf.hpp"
27 : #include "ParallelAlgorithms/Actions/AddSimpleTags.hpp"
28 : #include "ParallelAlgorithms/Actions/TerminatePhase.hpp"
29 : #include "Utilities/ErrorHandling/Error.hpp"
30 : #include "Utilities/Gsl.hpp"
31 : #include "Utilities/ProtocolHelpers.hpp"
32 : #include "Utilities/TMPL.hpp"
33 :
34 : namespace gh::bbh {
35 : namespace Tags {
36 : /// Successful AhC finds keyed by temporal id, storing the corresponding `LMax`.
37 1 : struct CommonHorizonSuccessRecords : db::SimpleTag {
38 0 : using type = std::map<LinkedMessageId<double>, size_t>;
39 0 : using option_tags = tmpl::list<>;
40 :
41 0 : static constexpr bool pass_metavariables = false;
42 0 : static type create_from_options() { return {}; }
43 : };
44 :
45 : /// Reduced constraint maxima keyed by time.
46 1 : struct ConstraintCheckRecords : db::SimpleTag {
47 0 : using key_type = double;
48 0 : using mapped_type = std::pair<double, double>;
49 0 : using type = std::map<key_type, mapped_type>;
50 0 : using option_tags = tmpl::list<>;
51 :
52 0 : static constexpr bool pass_metavariables = false;
53 0 : static type create_from_options() { return {}; }
54 : };
55 :
56 : /// Constraint checks already reported at verbose level.
57 1 : struct ReportedConstraintCheckRecords : db::SimpleTag {
58 0 : using type = std::set<gh::bbh::Tags::ConstraintCheckRecords::key_type>;
59 0 : using option_tags = tmpl::list<>;
60 :
61 0 : static constexpr bool pass_metavariables = false;
62 0 : static type create_from_options() { return {}; }
63 : };
64 : } // namespace Tags
65 :
66 0 : namespace Actions {
67 : /// Element simple action that latches completion-request state in the DataBox.
68 1 : struct SetElementCompletionRequested {
69 : template <typename ParallelComponent, typename DbTags, typename Metavariables,
70 : typename ArrayIndex>
71 0 : static void apply(db::DataBox<DbTags>& box,
72 : Parallel::GlobalCache<Metavariables>& /*cache*/,
73 : const ArrayIndex& /*array_index*/) {
74 : db::mutate<gh::bbh::Tags::ElementCompletionRequested>(
75 : [](const gsl::not_null<bool*> element_completion_requested) {
76 : *element_completion_requested = true;
77 : },
78 : make_not_null(&box));
79 : }
80 : };
81 : } // namespace Actions
82 :
83 : namespace detail {
84 : template <typename Metavariables>
85 : struct BroadcastCompletionRequestToElements {
86 : static void apply(Parallel::GlobalCache<Metavariables>& cache) {
87 : using dg_array = typename Metavariables::gh_dg_element_array;
88 : if constexpr (Parallel::is_dg_element_collection_v<dg_array>) {
89 : Parallel::threaded_action<Parallel::Actions::SimpleActionOnElement<
90 : gh::bbh::Actions::SetElementCompletionRequested, true>>(
91 : Parallel::get_parallel_component<dg_array>(cache));
92 : } else {
93 : Parallel::simple_action<gh::bbh::Actions::SetElementCompletionRequested>(
94 : Parallel::get_parallel_component<dg_array>(cache));
95 : }
96 : }
97 : };
98 :
99 : template <typename DbTags, typename Metavariables>
100 : void recompute_completion_state(const gsl::not_null<db::DataBox<DbTags>*> box,
101 : Parallel::GlobalCache<Metavariables>& cache) {
102 : const size_t min_successes =
103 : Parallel::get<gh::bbh::Tags::MinCommonHorizonSuccessesBeforeChecks>(
104 : cache);
105 : const size_t max_successes =
106 : Parallel::get<gh::bbh::Tags::MaxCommonHorizonSuccesses>(cache);
107 : const size_t l_max_threshold =
108 : Parallel::get<gh::bbh::Tags::CommonHorizonLMaxThreshold>(cache);
109 : const double gauge_constraint_threshold =
110 : Parallel::get<gh::bbh::Tags::GaugeConstraintLinfThreshold>(cache);
111 : const double three_index_constraint_threshold =
112 : Parallel::get<gh::bbh::Tags::ThreeIndexConstraintLinfThreshold>(cache);
113 : const bool verbose =
114 : Parallel::get<gh::bbh::Tags::ConstraintCheckVerbose>(cache);
115 : if (max_successes < min_successes) {
116 : ERROR_NO_TRACE("MaxCommonHorizonSuccesses ("
117 : << max_successes << ") must be >= "
118 : << "MinCommonHorizonSuccessesBeforeChecks (" << min_successes
119 : << ").");
120 : }
121 :
122 : const auto& horizon_records =
123 : db::get<gh::bbh::Tags::CommonHorizonSuccessRecords>(*box);
124 : const size_t success_count = horizon_records.size();
125 : bool lmax_latched = false;
126 : std::optional<std::pair<double, size_t>> first_lmax_match{};
127 : for (const auto& [time_id, l_max] : horizon_records) {
128 : if (l_max <= l_max_threshold) {
129 : lmax_latched = true;
130 : first_lmax_match = std::pair{time_id.id, l_max};
131 : break;
132 : }
133 : }
134 :
135 : std::optional<double> first_gauge_match_time{};
136 : std::optional<double> first_three_index_match_time{};
137 : size_t successes_up_to_constraint_time = 0;
138 : auto horizon_it = horizon_records.begin();
139 : const auto& constraint_records =
140 : db::get<gh::bbh::Tags::ConstraintCheckRecords>(*box);
141 : const auto& reported_constraint_records =
142 : db::get<gh::bbh::Tags::ReportedConstraintCheckRecords>(*box);
143 : std::vector<gh::bbh::Tags::ConstraintCheckRecords::key_type>
144 : newly_reported_constraint_records{};
145 : for (const auto& [constraint_time, maxima] : constraint_records) {
146 : while (horizon_it != horizon_records.end() and
147 : horizon_it->first.id <= constraint_time) {
148 : ++successes_up_to_constraint_time;
149 : ++horizon_it;
150 : }
151 : if (successes_up_to_constraint_time < min_successes) {
152 : continue;
153 : }
154 : if (verbose and not reported_constraint_records.contains(constraint_time)) {
155 : Parallel::printf(
156 : "BBH completion constraint check at t=%.16f: "
157 : "Linf(GaugeConstraint)=%.16e (threshold %.16e), "
158 : "Linf(ThreeIndexConstraint)=%.16e (threshold %.16e).\n",
159 : constraint_time, maxima.first, gauge_constraint_threshold,
160 : maxima.second, three_index_constraint_threshold);
161 : newly_reported_constraint_records.push_back(constraint_time);
162 : }
163 : if (not first_gauge_match_time.has_value() and
164 : maxima.first >= gauge_constraint_threshold) {
165 : first_gauge_match_time = constraint_time;
166 : }
167 : if (not first_three_index_match_time.has_value() and
168 : maxima.second >= three_index_constraint_threshold) {
169 : first_three_index_match_time = constraint_time;
170 : }
171 : }
172 : const bool horizon_completion_requested =
173 : success_count >= min_successes and
174 : (success_count >= max_successes or lmax_latched);
175 : const bool completion_requested = horizon_completion_requested or
176 : first_gauge_match_time.has_value() or
177 : first_three_index_match_time.has_value();
178 :
179 : const bool old_gauge_exceeded =
180 : db::get<gh::bbh::Tags::GaugeConstraintExceeded>(*box);
181 : const bool old_three_index_exceeded =
182 : db::get<gh::bbh::Tags::ThreeIndexConstraintExceeded>(*box);
183 : const bool old_lmax_latched =
184 : db::get<gh::bbh::Tags::CommonHorizonLMaxBelowOrEqualThreshold>(*box);
185 : const size_t old_success_count =
186 : db::get<gh::bbh::Tags::CommonHorizonSuccessCount>(*box);
187 : const bool old_completion_requested =
188 : db::get<gh::bbh::Tags::CompletionRequested>(*box);
189 : db::mutate<gh::bbh::Tags::GaugeConstraintExceeded,
190 : gh::bbh::Tags::ThreeIndexConstraintExceeded,
191 : gh::bbh::Tags::CommonHorizonLMaxBelowOrEqualThreshold,
192 : gh::bbh::Tags::CommonHorizonSuccessCount,
193 : gh::bbh::Tags::CompletionRequested>(
194 : [&first_gauge_match_time, &first_three_index_match_time, &lmax_latched,
195 : &success_count, &completion_requested](
196 : const gsl::not_null<bool*> gauge_constraint_exceeded_flag,
197 : const gsl::not_null<bool*> three_index_constraint_exceeded_flag,
198 : const gsl::not_null<bool*> lmax_latched_flag,
199 : const gsl::not_null<size_t*> common_horizon_success_count,
200 : const gsl::not_null<bool*> completion_requested_flag) {
201 : *gauge_constraint_exceeded_flag = first_gauge_match_time.has_value();
202 : *three_index_constraint_exceeded_flag =
203 : first_three_index_match_time.has_value();
204 : *lmax_latched_flag = lmax_latched;
205 : *common_horizon_success_count = success_count;
206 : *completion_requested_flag = completion_requested;
207 : },
208 : box);
209 : if (not newly_reported_constraint_records.empty()) {
210 : db::mutate<gh::bbh::Tags::ReportedConstraintCheckRecords>(
211 : [&newly_reported_constraint_records](
212 : const gsl::not_null<
213 : gh::bbh::Tags::ReportedConstraintCheckRecords::type*>
214 : records) {
215 : for (const auto& key : newly_reported_constraint_records) {
216 : records->insert(key);
217 : }
218 : },
219 : box);
220 : }
221 : if (old_success_count < min_successes and success_count >= min_successes) {
222 : Parallel::printf(
223 : "BBH completion criterion armed: AhC successes reached %zu (minimum "
224 : "required: %zu).\n",
225 : success_count, min_successes);
226 : }
227 : if (not old_lmax_latched and lmax_latched and first_lmax_match.has_value()) {
228 : Parallel::printf(
229 : "BBH completion criterion met at t=%.16f: AhC Lmax=%zu <= %zu.\n",
230 : first_lmax_match->first, first_lmax_match->second, l_max_threshold);
231 : }
232 : if (not old_gauge_exceeded and first_gauge_match_time.has_value()) {
233 : Parallel::printf(
234 : "BBH completion criterion met at t=%.16f: "
235 : "Linf(GaugeConstraint) >= %.16e.\n",
236 : *first_gauge_match_time, gauge_constraint_threshold);
237 : }
238 : if (not old_three_index_exceeded and
239 : first_three_index_match_time.has_value()) {
240 : Parallel::printf(
241 : "BBH completion criterion met at t=%.16f: "
242 : "Linf(ThreeIndexConstraint) >= %.16e.\n",
243 : *first_three_index_match_time, three_index_constraint_threshold);
244 : }
245 : if (not old_completion_requested and completion_requested) {
246 : Parallel::printf("BBH completion criteria request latched.\n");
247 : BroadcastCompletionRequestToElements<Metavariables>::apply(cache);
248 : }
249 : }
250 :
251 : struct InitializeCompletionState : tt::ConformsTo<db::protocols::Mutator> {
252 : using return_tags =
253 : tmpl::list<gh::bbh::Tags::GaugeConstraintExceeded,
254 : gh::bbh::Tags::ThreeIndexConstraintExceeded,
255 : gh::bbh::Tags::CommonHorizonLMaxBelowOrEqualThreshold,
256 : gh::bbh::Tags::CommonHorizonSuccessCount,
257 : gh::bbh::Tags::CompletionRequested,
258 : gh::bbh::Tags::CommonHorizonSuccessRecords,
259 : gh::bbh::Tags::ConstraintCheckRecords,
260 : gh::bbh::Tags::ReportedConstraintCheckRecords>;
261 : using argument_tags = tmpl::list<>;
262 :
263 : static void apply(
264 : const gsl::not_null<bool*> gauge_constraint_exceeded,
265 : const gsl::not_null<bool*> three_index_constraint_exceeded,
266 : const gsl::not_null<bool*> common_horizon_lmax_below_or_equal_threshold,
267 : const gsl::not_null<size_t*> common_horizon_success_count,
268 : const gsl::not_null<bool*> completion_requested,
269 : const gsl::not_null<gh::bbh::Tags::CommonHorizonSuccessRecords::type*>
270 : common_horizon_success_records,
271 : const gsl::not_null<gh::bbh::Tags::ConstraintCheckRecords::type*>
272 : constraint_check_records,
273 : const gsl::not_null<gh::bbh::Tags::ReportedConstraintCheckRecords::type*>
274 : reported_constraint_check_records) {
275 : *gauge_constraint_exceeded = false;
276 : *three_index_constraint_exceeded = false;
277 : *common_horizon_lmax_below_or_equal_threshold = false;
278 : *common_horizon_success_count = 0;
279 : *completion_requested = false;
280 : common_horizon_success_records->clear();
281 : constraint_check_records->clear();
282 : reported_constraint_check_records->clear();
283 : }
284 : };
285 : } // namespace detail
286 :
287 : namespace Actions {
288 : /// Records a successful AhC find and updates BBH completion state.
289 1 : struct RecordCommonHorizonSuccess {
290 : template <typename ParallelComponent, typename DbTags, typename Metavariables,
291 : typename ArrayIndex>
292 0 : static void apply(db::DataBox<DbTags>& box,
293 : Parallel::GlobalCache<Metavariables>& cache,
294 : const ArrayIndex& /*array_index*/,
295 : const LinkedMessageId<double>& temporal_id,
296 : const size_t l_max) {
297 : auto& common_horizon_success_records =
298 : db::get_mutable_reference<gh::bbh::Tags::CommonHorizonSuccessRecords>(
299 : make_not_null(&box));
300 : const bool inserted =
301 : common_horizon_success_records.emplace(temporal_id, l_max).second;
302 : if (not inserted) {
303 : ERROR("Duplicate common-horizon completion record for temporal id "
304 : << temporal_id << ".");
305 : }
306 : detail::recompute_completion_state(make_not_null(&box), cache);
307 : }
308 : };
309 :
310 : /// Reduction target callback that records reduced constraint maxima and updates
311 : /// BBH completion state.
312 1 : struct ProcessConstraintMaxima {
313 : template <typename ParallelComponent, typename DbTags, typename Metavariables,
314 : typename ArrayIndex>
315 0 : static void apply(db::DataBox<DbTags>& box,
316 : Parallel::GlobalCache<Metavariables>& cache,
317 : const ArrayIndex& /*array_index*/, const double time,
318 : const double max_gauge_linf,
319 : const double max_three_index_linf) {
320 : auto& constraint_check_records =
321 : db::get_mutable_reference<gh::bbh::Tags::ConstraintCheckRecords>(
322 : make_not_null(&box));
323 : const bool inserted =
324 : constraint_check_records
325 : .emplace(time,
326 : gh::bbh::Tags::ConstraintCheckRecords::mapped_type{
327 : max_gauge_linf, max_three_index_linf})
328 : .second;
329 : if (not inserted) {
330 : ERROR("Duplicate BBH completion constraint-max record at t=" << time
331 : << ".");
332 : }
333 : detail::recompute_completion_state(make_not_null(&box), cache);
334 : }
335 : };
336 :
337 : /// Initializes the element-local completion-request tag.
338 1 : struct InitializeElementCompletionRequested
339 : : tt::ConformsTo<db::protocols::Mutator> {
340 0 : using return_tags = tmpl::list<gh::bbh::Tags::ElementCompletionRequested>;
341 0 : using argument_tags = tmpl::list<>;
342 :
343 0 : static void apply(const gsl::not_null<bool*> element_completion_requested) {
344 : *element_completion_requested = false;
345 : }
346 : };
347 : } // namespace Actions
348 :
349 : template <class Metavariables>
350 0 : struct CompletionSingleton {
351 0 : using chare_type = Parallel::Algorithms::Singleton;
352 0 : static constexpr bool checkpoint_data = true;
353 0 : using metavariables = Metavariables;
354 0 : using const_global_cache_tags =
355 : tmpl::list<gh::bbh::Tags::MinCommonHorizonSuccessesBeforeChecks,
356 : gh::bbh::Tags::MaxCommonHorizonSuccesses,
357 : gh::bbh::Tags::GaugeConstraintLinfThreshold,
358 : gh::bbh::Tags::ThreeIndexConstraintLinfThreshold,
359 : gh::bbh::Tags::CommonHorizonLMaxThreshold,
360 : gh::bbh::Tags::ConstraintCheckVerbose>;
361 0 : using phase_dependent_action_list = tmpl::list<Parallel::PhaseActions<
362 : Parallel::Phase::Initialization,
363 : tmpl::list<Initialization::Actions::AddSimpleTags<
364 : gh::bbh::detail::InitializeCompletionState>,
365 : Parallel::Actions::TerminatePhase>>>;
366 0 : using simple_tags_from_options = tmpl::list<>;
367 :
368 0 : static void execute_next_phase(
369 : const Parallel::Phase next_phase,
370 : Parallel::CProxy_GlobalCache<Metavariables>& global_cache) {
371 : auto& local_cache = *Parallel::local_branch(global_cache);
372 : Parallel::get_parallel_component<CompletionSingleton<Metavariables>>(
373 : local_cache)
374 : .start_phase(next_phase);
375 : }
376 : };
377 : } // namespace gh::bbh
|
