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 <limits>
8 : #include <optional>
9 : #include <string>
10 : #include <type_traits>
11 : #include <utility>
12 : #include <vector>
13 :
14 : #include "DataStructures/DataBox/DataBox.hpp"
15 : #include "DataStructures/DataVector.hpp"
16 : #include "DataStructures/LinkedMessageId.hpp"
17 : #include "DataStructures/Tensor/IndexType.hpp"
18 : #include "DataStructures/Variables.hpp"
19 : #include "Domain/Creators/Tags/Domain.hpp"
20 : #include "Domain/FunctionsOfTime/Tags.hpp"
21 : #include "Domain/Structure/ElementId.hpp"
22 : #include "IO/Logging/Verbosity.hpp"
23 : #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp"
24 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
25 : #include "NumericalAlgorithms/Strahlkorper/Strahlkorper.hpp"
26 : #include "Parallel/GlobalCache.hpp"
27 : #include "ParallelAlgorithms/Actions/FunctionsOfTimeAreReady.hpp"
28 : #include "ParallelAlgorithms/ApparentHorizonFinder/Callbacks/InvokeCallbacks.hpp"
29 : #include "ParallelAlgorithms/ApparentHorizonFinder/CleanUp.hpp"
30 : #include "ParallelAlgorithms/ApparentHorizonFinder/ComputeCoords.hpp"
31 : #include "ParallelAlgorithms/ApparentHorizonFinder/CurrentTime.hpp"
32 : #include "ParallelAlgorithms/ApparentHorizonFinder/Destination.hpp"
33 : #include "ParallelAlgorithms/ApparentHorizonFinder/FastFlow.hpp"
34 : #include "ParallelAlgorithms/ApparentHorizonFinder/HorizonAliases.hpp"
35 : #include "ParallelAlgorithms/ApparentHorizonFinder/InterpolateVolumeVars.hpp"
36 : #include "ParallelAlgorithms/ApparentHorizonFinder/OptionTags.hpp"
37 : #include "ParallelAlgorithms/ApparentHorizonFinder/Storage.hpp"
38 : #include "ParallelAlgorithms/ApparentHorizonFinder/Tags.hpp"
39 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
40 : #include "Utilities/Algorithm.hpp"
41 : #include "Utilities/Gsl.hpp"
42 : #include "Utilities/TMPL.hpp"
43 :
44 : namespace ah {
45 : namespace detail {
46 : template <typename Criterion>
47 : struct get_tags {
48 : using type = typename Criterion::compute_tags_for_observation_box;
49 : };
50 : } // namespace detail
51 :
52 : /*!
53 : * \brief Simple action run on the horizon finder by the Elements which receives
54 : * volume data, finds the apparent horizon, and calls the callbacks after the
55 : * horizon is found.
56 : *
57 : * \details First, the volume variables `ah::vars_to_interpolate_to_target`
58 : * and the \p dependency are put into the `ah::Tags::Storage` for the incoming
59 : * time. Then, we try to set the
60 : * `ah::Tags::CurrentTime` with the `ah::set_current_time` function. Once we
61 : * have a current time, we ensure the functions of time are up-to-date at this
62 : * current time with `ah::check_if_current_time_is_ready`. If they are not,
63 : * we re-trigger this action with `vars_have_already_been_received` set to
64 : * true and not sending any new volume data.
65 : *
66 : * Once we ensure the functions of time are ready, we compute the cartesian
67 : * coordinates for the current fast-flow iteration surface with
68 : * `ah::set_current_iteration_coords`. Once we have the coords, we interpolate
69 : * the volume date onto these coordinates. After that, we do one iteration of
70 : * the `FastFlow` algorithm to get a new surface. If we haven't converged yet,
71 : * we compute new cartesian and start another iteration. Once we have converged,
72 : * we call the callbacks with `ah::invoke_callbacks` and then clean up the
73 : * horizon finer with `ah::clean_up_horizon_finder`. Then we try to set a new
74 : * current time and the process starts over again. If `FastFlow` never
75 : * converges or the interpolation fails, we call the
76 : * `horizon_find_failure_callbacks` from the \p HorizonMetavars.
77 : *
78 : * The \p dependency is the name of the volume data that an
79 : * `observers::ObserverWriter` holds for the incoming time until the horizon
80 : * find at that time succeeds or fails (see
81 : * `ah::callbacks::SendDependencyToObserverWriter`), or `std::nullopt` if there
82 : * is no such data. Since this action can receive volume data for a later time
83 : * while it is still working on the current time, the dependency is stored per
84 : * time and copied from the storage of the current time into
85 : * `ah::Tags::Dependency` right before the callbacks for that time are invoked,
86 : * on success (in `ah::invoke_callbacks`) as well as on failure.
87 : */
88 : template <typename HorizonMetavars>
89 1 : struct FindApparentHorizon {
90 0 : using frame = typename HorizonMetavars::frame;
91 :
92 : template <typename ParallelComponent, typename DbTags, typename Metavariables,
93 : typename ArrayIndex>
94 0 : static void apply(db::DataBox<DbTags>& box,
95 : Parallel::GlobalCache<Metavariables>& cache,
96 : const ArrayIndex& /*array_index*/,
97 : const LinkedMessageId<double>& incoming_time,
98 : const ElementId<3>& incoming_element_id,
99 : const ::Mesh<3>& incoming_mesh,
100 : Variables<ah::vars_to_interpolate_to_target<3, frame>>&&
101 : incoming_vars_to_interpolate,
102 : const std::optional<std::string>& dependency,
103 : const bool vars_have_already_been_received = false) {
104 : const auto& verbosity = db::get<ah::Tags::Verbosity>(box);
105 : const bool quiet_print = verbosity >= ::Verbosity::Quiet;
106 : const bool verbose_print = verbosity >= ::Verbosity::Verbose;
107 : const bool debug_print = verbosity >= ::Verbosity::Debug;
108 : const std::string& name = pretty_type::name<HorizonMetavars>();
109 :
110 : // If we've already completed this horizon find, ignore the volume data and
111 : // just return
112 : auto& completed_times = db::get_mutable_reference<ah::Tags::CompletedTimes>(
113 : make_not_null(&box));
114 : if (completed_times.contains(incoming_time)) {
115 : if (debug_print) {
116 : Parallel::printf(
117 : "%s: Received volume data at finished time %s from Element %s. "
118 : "Ignoring.\n",
119 : name, incoming_time, incoming_element_id);
120 : }
121 : return;
122 : }
123 :
124 : auto& current_time_optional =
125 : db::get_mutable_reference<ah::Tags::CurrentTime>(make_not_null(&box));
126 : auto& pending_times =
127 : db::get_mutable_reference<ah::Tags::PendingTimes>(make_not_null(&box));
128 :
129 : // Only add to pending if it isn't the current time.
130 : if (not(current_time_optional.has_value() and
131 : current_time_optional.value() == incoming_time)) {
132 : pending_times.insert(incoming_time);
133 : }
134 :
135 : auto& all_storage = db::get_mutable_reference<ah::Tags::Storage<frame>>(
136 : make_not_null(&box));
137 : auto& block_search_order =
138 : db::get_mutable_reference<ah::Tags::BlockSearchOrder>(
139 : make_not_null(&box));
140 :
141 : // Add volume variables, destination, and dependency to the storage if they
142 : // haven't already been received
143 : if (not vars_have_already_been_received) {
144 : auto& current_time_storage = all_storage[incoming_time];
145 : current_time_storage.all_volume_variables.emplace(
146 : incoming_element_id,
147 : ah::Storage::VolumeVariables<frame>{
148 : incoming_mesh, std::move(incoming_vars_to_interpolate)});
149 : current_time_storage.destination = HorizonMetavars::destination;
150 : current_time_storage.dependency = dependency;
151 : }
152 :
153 : // ========================================================================
154 : // The incoming vars are moved and used in the above code. Below, we
155 : // only use quantities at the current time, except for forwarding the
156 : // incoming quantities to a callback for checking if the functions of time
157 : // are ready.
158 : // ========================================================================
159 :
160 : // Keep trying to find horizons for as long as we can.
161 : // In the first pass, interpolate only from the incoming element because
162 : // this was missing data. On subsequent passes, interpolate from all
163 : // elements. Notes:
164 : // - The `vars_have_already_been_received` flag is used only in
165 : // `check_if_current_time_is_ready` to re-trigger the AH find in a
166 : // callback but not sending any new data.
167 : // - If the incoming time is not the current time (but later), then we have
168 : // to do the full interpolation for the current time first before moving
169 : // on to the incoming time. In this case we also need to do the full
170 : // interpolation for the incoming time once we get there.
171 : bool interpolate_only_from_incoming_element =
172 : current_time_optional == incoming_time and
173 : not vars_have_already_been_received;
174 : while (true) {
175 : // Set current time using the pending times if it isn't already set
176 : set_current_time(make_not_null(¤t_time_optional),
177 : make_not_null(&pending_times), completed_times,
178 : all_storage, verbosity, name);
179 :
180 : // At this point, if we don't have a current id we can't do anything so
181 : // just exit.
182 : if (not current_time_optional.has_value()) {
183 : if (verbose_print) {
184 : Parallel::printf(
185 : "%s: No current time set. Stopping here. Incoming time: %s\n",
186 : name, incoming_time);
187 : }
188 : return;
189 : }
190 :
191 : const LinkedMessageId<double>& current_time =
192 : current_time_optional.value();
193 : ASSERT(all_storage.contains(current_time),
194 : "Current time doesn't have any data.");
195 :
196 : auto& current_time_storage = all_storage.at(current_time);
197 :
198 : // Makes the dependency of the current time available to the failure
199 : // callbacks, e.g. so `ah::callbacks::SendDependencyToObserverWriter` can
200 : // tell the ObserverWriter to drop the volume data it holds for this
201 : // time. `ah::invoke_callbacks` does the same for the success callbacks.
202 : const auto set_dependency_for_failure_callbacks =
203 : [&box, ¤t_time_storage]() {
204 : db::mutate<ah::Tags::Dependency>(
205 : [¤t_time_storage](
206 : const gsl::not_null<std::optional<std::string>*>
207 : dependency_in_box) {
208 : *dependency_in_box = current_time_storage.dependency;
209 : },
210 : make_not_null(&box));
211 : };
212 :
213 : // Only check if the current time is ready once
214 : if (UNLIKELY(not current_time_storage.time_is_ready)) {
215 : if (check_if_current_time_is_ready<HorizonMetavars>(
216 : current_time, cache, incoming_time, incoming_element_id,
217 : incoming_mesh, dependency)) {
218 : // Set this so we don't check it again
219 : current_time_storage.time_is_ready = true;
220 :
221 : if (verbose_print) {
222 : Parallel::printf("%s: Current time %s is ready.\n", name,
223 : current_time);
224 : }
225 : } else {
226 : if (verbose_print) {
227 : Parallel::printf("%s: Current time %s is NOT ready.\n", name,
228 : current_time);
229 : }
230 : // We need to wait for the functions of time to be updated so we can't
231 : // do anything yet.
232 : return;
233 : }
234 : }
235 :
236 : const auto& domain = Parallel::get<domain::Tags::Domain<3>>(cache);
237 : const auto& functions_of_time =
238 : Parallel::get<domain::Tags::FunctionsOfTime>(cache);
239 : const auto& options =
240 : Parallel::get<ah::Tags::ApparentHorizonOptions<HorizonMetavars>>(
241 : cache);
242 : auto& fast_flow =
243 : db::get_mutable_reference<ah::Tags::FastFlow>(make_not_null(&box));
244 :
245 : auto& all_volume_variables = current_time_storage.all_volume_variables;
246 : auto& current_iteration_storage = current_time_storage.current_iteration;
247 : auto& previous_iteration_surface =
248 : current_time_storage.previous_iteration_surface;
249 : auto& element_order = current_time_storage.element_order;
250 : auto& previous_surfaces =
251 : db::get_mutable_reference<ah::Tags::PreviousSurfaces<frame>>(
252 : make_not_null(&box));
253 :
254 : auto& current_resolution_l =
255 : db::get_mutable_reference<ah::Tags::CurrentResolutionL>(
256 : make_not_null(&box));
257 : // Set current resolution to the initial guess resolution if it is not
258 : // yet specified.
259 : if (UNLIKELY(not current_resolution_l.has_value())) {
260 : ASSERT(fast_flow.current_iteration() == 0,
261 : "Current resolution L is not set, but fast flow iteration > 0. "
262 : "Iteration 0 should set the current resolution L.");
263 : ASSERT(previous_surfaces.empty(),
264 : "Current resolution L is not set, but previous horizons have "
265 : "been found. Iteration 0 of each horizon find should set the "
266 : "current resolution L.");
267 : current_resolution_l = options.initial_guess.l_max();
268 : }
269 :
270 : // Keep finding the horizon until it's found with sufficient resolution
271 : bool rerunning_with_higher_resolution = false;
272 : while (true) {
273 : std::pair<FastFlow::Status, FastFlow::IterInfo> status_and_info;
274 :
275 : // Keep doing fast flow iterations for as long as we can until we find a
276 : // horizon
277 : while (true) {
278 : // Points haven't been set for this iteration so need to do so now
279 : if (not current_iteration_storage.block_coord_holders.has_value() or
280 : rerunning_with_higher_resolution) {
281 : const bool coords_set_successfully = set_current_iteration_coords(
282 : make_not_null(¤t_iteration_storage),
283 : make_not_null(&block_search_order), current_time, fast_flow,
284 : options.initial_guess, previous_iteration_surface,
285 : previous_surfaces, options.max_compute_coords_retries, domain,
286 : functions_of_time, current_resolution_l,
287 : rerunning_with_higher_resolution);
288 :
289 : // Some points are outside the domain and we cannot recover
290 : if (not coords_set_successfully) {
291 : // This will possibly throw an error, depending on which
292 : // callback is used in
293 : // HorizonMetavars::horizon_find_failure_callbacks
294 : // (ErrorOnFailedApparentHorizon or IgnoreFailedApparentHorizon)
295 : set_dependency_for_failure_callbacks();
296 : tmpl::for_each<
297 : typename HorizonMetavars::horizon_find_failure_callbacks>(
298 : [&]<typename Callback>(tmpl::type_<Callback> /*meta*/) {
299 : Callback::apply(box, cache,
300 : FastFlow::Status::InterpolationFailure);
301 : });
302 :
303 : // Still clean up in case we didn't error so we don't keep
304 : // accepting volume data for this time
305 : clean_up_horizon_finder(make_not_null(¤t_time_optional),
306 : make_not_null(&all_storage),
307 : make_not_null(&completed_times),
308 : make_not_null(&fast_flow));
309 :
310 : // If we didn't error, move on to the next horizon find
311 : break;
312 : }
313 :
314 : if (debug_print) {
315 : Parallel::printf(
316 : "%s: For current time %s, at iteration %zu, coords have been "
317 : "set.\n",
318 : name, current_time, fast_flow.current_iteration());
319 : }
320 : }
321 :
322 : // Interpolate to the current iteration points
323 : bool interpolation_is_complete = false;
324 : if (interpolate_only_from_incoming_element) {
325 : // Data is coming in one element at a time, so try to interpolate
326 : // only from the incoming element
327 : ASSERT(all_volume_variables.contains(incoming_element_id),
328 : "Trying to interpolate incoming data from element "
329 : << incoming_element_id
330 : << " but it's not available. This is a bug.");
331 : const bool interpolated_any_points = interpolate_volume_data(
332 : make_not_null(¤t_iteration_storage),
333 : all_volume_variables.at(incoming_element_id),
334 : incoming_element_id);
335 : // Store which elements we found points in for next iteration
336 : if (interpolated_any_points) {
337 : element_order.push_back(incoming_element_id);
338 : }
339 : interpolation_is_complete =
340 : current_iteration_storage.interpolation_is_complete();
341 : } else {
342 : // We already have some element data, so go through the elements in
343 : // the order we found points in before
344 : for (const auto& element_id : element_order) {
345 : ASSERT(all_volume_variables.contains(element_id),
346 : "Trying to interpolate data from element "
347 : << element_id
348 : << " that was found in a previous FastFlow iteration, "
349 : "but it's not available. This is a bug.");
350 : interpolate_volume_data(make_not_null(¤t_iteration_storage),
351 : all_volume_variables.at(element_id),
352 : element_id);
353 : interpolation_is_complete =
354 : current_iteration_storage.interpolation_is_complete();
355 : if (interpolation_is_complete) {
356 : break;
357 : }
358 : }
359 : // If we still need more points, try going through all the elements
360 : // that we have received data from so far
361 : if (not interpolation_is_complete) {
362 : for (const auto& [element_id, volume_vars_storage] :
363 : all_volume_variables) {
364 : if (std::find(element_order.begin(), element_order.end(),
365 : element_id) != element_order.end()) {
366 : // Already tried this element above
367 : continue;
368 : }
369 : const bool interpolated_any_points = interpolate_volume_data(
370 : make_not_null(¤t_iteration_storage),
371 : volume_vars_storage, element_id);
372 : if (interpolated_any_points) {
373 : element_order.push_back(element_id);
374 : }
375 : interpolation_is_complete =
376 : current_iteration_storage.interpolation_is_complete();
377 : if (interpolation_is_complete) {
378 : break;
379 : }
380 : }
381 : }
382 : }
383 :
384 : if (interpolation_is_complete) {
385 : // Next iteration try to interpolate from all elements
386 : interpolate_only_from_incoming_element = false;
387 : } else {
388 : // We need more points, return and wait for more volume data.
389 : // When this action is called again, we'll try to interpolate only
390 : // from the incoming element until we have enough data.
391 : if (debug_print) {
392 : Parallel::printf(
393 : "%s: For current time %s, at iteration %zu, need more volume "
394 : "data.\n",
395 : name, current_time, fast_flow.current_iteration());
396 : }
397 : return;
398 : }
399 :
400 : // Set this before we iterate the fast flow.
401 : previous_iteration_surface = current_iteration_storage.strahlkorper;
402 :
403 : // Do a single fast flow iteration.
404 : {
405 : using InverseSpatialMetric =
406 : gr::Tags::InverseSpatialMetric<DataVector, 3, frame>;
407 : using ExtrinsicCurvature =
408 : gr::Tags::ExtrinsicCurvature<DataVector, 3, frame>;
409 : using SpatialChristoffel =
410 : gr::Tags::SpatialChristoffelSecondKind<DataVector, 3, frame>;
411 :
412 : auto& interpolated_vars =
413 : current_iteration_storage.interpolated_vars;
414 : status_and_info = fast_flow.iterate_horizon_finder(
415 : make_not_null(¤t_iteration_storage.strahlkorper),
416 : get<InverseSpatialMetric>(interpolated_vars),
417 : get<ExtrinsicCurvature>(interpolated_vars),
418 : get<SpatialChristoffel>(interpolated_vars));
419 : }
420 :
421 : const auto& status = status_and_info.first;
422 : const auto& info = status_and_info.second;
423 : const auto has_converged = converged(status);
424 :
425 : if (verbose_print or (quiet_print and has_converged)) {
426 : Parallel::printf(
427 : "%s: t=%.6g: L=%zu: its=%zu: %.1e<R<%.0e, |R|=%.1g, "
428 : "|R_grid|=%.1g, %.4g<r<%.4g\n",
429 : name, current_time.id,
430 : all_storage.at(current_time)
431 : .current_iteration.strahlkorper.l_max(),
432 : info.iteration, info.min_residual, info.max_residual,
433 : info.residual_ylm, info.residual_mesh, info.r_min, info.r_max);
434 : }
435 :
436 : if (status == FastFlow::Status::SuccessfulIteration) {
437 : // Reset so we compute and interpolate to new points
438 : current_iteration_storage.reset_for_next_iteration();
439 :
440 : // Continue in the iteration while loop
441 : continue;
442 : } else if (not has_converged) {
443 : // This will possibly throw an error
444 : set_dependency_for_failure_callbacks();
445 : tmpl::for_each<
446 : typename HorizonMetavars::horizon_find_failure_callbacks>(
447 : [&](auto callback_v) {
448 : using callback = tmpl::type_from<decltype(callback_v)>;
449 : callback::apply(box, cache, status);
450 : });
451 :
452 : // Still clean up in case we didn't error so we don't keep accepting
453 : // volume data for this time
454 : clean_up_horizon_finder(make_not_null(¤t_time_optional),
455 : make_not_null(&all_storage),
456 : make_not_null(&completed_times),
457 : make_not_null(&fast_flow));
458 :
459 : // If we didn't error, move on to the next horizon find
460 : break;
461 : }
462 :
463 : // Break out of the iteration while loop. Go back into the
464 : // while loop for sufficient resolution
465 : break;
466 : } // while (true) [iterations]
467 : if (converged(status_and_info.first)) {
468 : std::vector<size_t> recommended_resolutions{};
469 : const auto& criteria = options.criteria;
470 :
471 : // Create an ObservationBox with the union of all compute tags needed
472 : // by criteria
473 : using compute_tags_for_obs_box =
474 : tmpl::remove_duplicates<tmpl::flatten<tmpl::transform<
475 : tmpl::at<
476 : typename Metavariables::factory_creation::factory_classes,
477 : ah::Criterion>,
478 : detail::get_tags<tmpl::_1>>>>;
479 : const auto observation_box =
480 : make_observation_box<compute_tags_for_obs_box>(
481 : make_not_null(&box));
482 :
483 : for (const auto& criterion : criteria) {
484 : recommended_resolutions.push_back(criterion->evaluate(
485 : observation_box, cache,
486 : current_time_storage.current_iteration.strahlkorper,
487 : status_and_info.second));
488 : }
489 : const size_t next_resolution_l =
490 : recommended_resolutions.empty()
491 : ? all_storage.at(current_time)
492 : .current_iteration.strahlkorper.l_max()
493 : : *std::max_element(recommended_resolutions.begin(),
494 : recommended_resolutions.end());
495 : if (next_resolution_l > *current_resolution_l) {
496 : current_resolution_l = next_resolution_l;
497 : rerunning_with_higher_resolution = true;
498 :
499 : // Prepare for trying the find again with higher resolution:
500 : // Reset fast flow, but don't "clean up," which also erases the
501 : // storage for the current time, resets the current time, and
502 : // updates the completed times. None of those things happen yet,
503 : // because the current time horizon finding is not yet finished.
504 : // Also reset the current iteration storage, since the next find
505 : // is effectively a new iteration.
506 : fast_flow.reset_for_next_find();
507 : current_time_storage.current_iteration.reset_for_next_iteration();
508 :
509 : continue;
510 : } else {
511 : rerunning_with_higher_resolution = false;
512 : }
513 :
514 : // We have converged to the apparent horizon. Invoke the callbacks and
515 : // clean up for the next horizon find
516 : invoke_callbacks<HorizonMetavars>(make_not_null(&box), cache,
517 : status_and_info.first);
518 :
519 : if (debug_print) {
520 : Parallel::printf(
521 : "%s: Horizon find finished at current time %s. Cleaning up and "
522 : "moving to next time\n",
523 : name, current_time);
524 : }
525 :
526 : current_resolution_l = next_resolution_l;
527 : clean_up_horizon_finder(make_not_null(¤t_time_optional),
528 : make_not_null(&all_storage),
529 : make_not_null(&completed_times),
530 : make_not_null(&fast_flow));
531 : }
532 :
533 : // We have either i) converged with sufficient resolution, or
534 : // ii) failed to converge but not thrown an error. In both cases,
535 : // continue to the next time by breaking out of the resolution
536 : // while loop, going back into the overall while loop for time.
537 : break;
538 : } // while (true) [sufficient resolution]
539 : } // while (true) [time]
540 : }
541 : };
542 : } // namespace ah
|
