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 : #include <deque>
8 : #include <utility>
9 :
10 : #include "DataStructures/DataBox/DataBox.hpp"
11 : #include "DataStructures/VariablesTag.hpp"
12 : #include "Domain/Creators/Tags/Domain.hpp"
13 : #include "Domain/FunctionsOfTime/Tags.hpp"
14 : #include "Domain/StrahlkorperTransformations.hpp"
15 : #include "IO/Logging/Tags.hpp"
16 : #include "IO/Logging/Verbosity.hpp"
17 : #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
18 : #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
19 : #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
20 : #include "Parallel/GlobalCache.hpp"
21 : #include "Parallel/Invoke.hpp"
22 : #include "Parallel/ParallelComponentHelpers.hpp"
23 : #include "Parallel/Printf/Printf.hpp"
24 : #include "ParallelAlgorithms/ApparentHorizonFinder/FastFlow.hpp"
25 : #include "ParallelAlgorithms/ApparentHorizonFinder/InterpolationTarget.hpp"
26 : #include "ParallelAlgorithms/ApparentHorizonFinder/Tags.hpp"
27 : #include "ParallelAlgorithms/Interpolation/Actions/SendPointsToInterpolator.hpp"
28 : #include "ParallelAlgorithms/Interpolation/InterpolationTargetDetail.hpp"
29 : #include "ParallelAlgorithms/Interpolation/Protocols/PostInterpolationCallback.hpp"
30 : #include "PointwiseFunctions/GeneralRelativity/Surfaces/Tags.hpp"
31 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
32 : #include "Utilities/ErrorHandling/Error.hpp"
33 : #include "Utilities/Gsl.hpp"
34 : #include "Utilities/PrettyType.hpp"
35 : #include "Utilities/ProtocolHelpers.hpp"
36 : #include "Utilities/TMPL.hpp"
37 :
38 : /// \cond
39 : namespace ylm::Tags {
40 : template <typename Frame>
41 : struct Strahlkorper;
42 : } // namespace ylm::Tags
43 : namespace db {
44 : template <typename TagsList>
45 : class DataBox;
46 : } // namespace db
47 : namespace intrp {
48 : template <class Metavariables, typename InterpolationTargetTag>
49 : struct InterpolationTarget;
50 : namespace Tags {
51 : template <typename Metavariables>
52 : struct TemporalIds;
53 : } // namespace Tags
54 : } // namespace intrp
55 : /// \endcond
56 :
57 : namespace intrp::callbacks {
58 :
59 : /// \brief post interpolation callback (see
60 : /// intrp::protocols::PostInterpolationCallback) that does a FastFlow iteration
61 : /// and triggers another one until convergence.
62 : ///
63 : /// Assumes that InterpolationTargetTag contains an additional
64 : /// type alias called `post_horizon_find_callbacks`, which is a list of
65 : /// structs, each of which has a function
66 : ///
67 : /// \snippet ParallelAlgorithms/ApparentHorizonFinder/Test_ApparentHorizonFinder.cpp post_horizon_find_callback_example
68 : ///
69 : /// that is called if the FastFlow iteration has converged.
70 : /// InterpolationTargetTag also is assumed to contain an additional
71 : /// type alias called `horizon_find_failure_callbacks`, which is a list of
72 : /// structs, each of which has a function
73 : ///
74 : /// \snippet ParallelAlgorithms/ApparentHorizonFinder/Test_ApparentHorizonFinder.cpp horizon_find_failure_callbacks_example
75 : ///
76 : /// that is called if the FastFlow iteration or the interpolation has
77 : /// failed.
78 : ///
79 : /// Uses:
80 : /// - Metavariables:
81 : /// - `temporal_id`
82 : /// - DataBox:
83 : /// - `logging::Tags::Verbosity<InterpolationTargetTag>`
84 : /// - `::gr::Tags::InverseSpatialMetric<DataVector, 3, Frame>`
85 : /// - `::gr::Tags::ExtrinsicCurvature<DataVector, 3, Frame>`
86 : /// - `::gr::Tags::SpatialChristoffelSecondKind<DataVector, 3, Frame>`
87 : /// - `::ah::Tags::FastFlow`
88 : /// - `ylm::Tags::Strahlkorper<Frame>`
89 : ///
90 : /// Modifies:
91 : /// - DataBox:
92 : /// - `::ah::Tags::FastFlow`
93 : /// - `ylm::Tags::Strahlkorper<Frame>`
94 : ///
95 : /// This can be used in InterpolationTargetTag::post_interpolation_callbacks;
96 : /// see intrp::protocols::InterpolationTargetTag for details on
97 : /// InterpolationTargetTag.
98 : ///
99 : /// ### Output
100 : ///
101 : /// Optionally, a single line of output is printed to stdout either on each
102 : /// iteration (if verbosity > Verbosity::Quiet) or on convergence
103 : /// (if verbosity > Verbosity::Silent). The output consists of the
104 : /// following labels, and values associated with each label:
105 : /// - t = time given by value of `temporal_id` argument to `apply`
106 : /// - its = current iteration number.
107 : /// - R = min and max of residual over all prolonged grid points.
108 : /// - |R| = L2 norm of residual, counting only L modes solved for.
109 : /// - |R_mesh| = L2 norm of residual over prolonged grid points.
110 : /// - r = min and max radius of trial horizon surface.
111 : ///
112 : /// #### Difference between |R| and |R_mesh|:
113 : /// The horizon is represented in a \f$Y_{lm}\f$ expansion up
114 : /// to \f$l=l_{\mathrm{surface}}\f$;
115 : /// the residual |R| represents the failure of that surface to satisfy
116 : /// the apparent horizon equation.
117 : ///
118 : /// However, at each iteration we also interpolate the horizon surface
119 : /// to a higher resolution ("prolongation"). The prolonged surface
120 : /// includes \f$Y_{lm}\f$ coefficients up to \f$l=l_{\mathrm{mesh}}\f$,
121 : /// where \f$l_{\mathrm{mesh}} > l_{\mathrm{surface}}\f$.
122 : /// The residual computed on this higher-resolution surface is |R_mesh|.
123 : ///
124 : /// As iterations proceed, |R| should decrease until it reaches
125 : /// numerical roundoff error, because we are varying all the \f$Y_{lm}\f$
126 : /// coefficients up to \f$l=l_{\mathrm{surface}}\f$
127 : /// to minimize the residual. However,
128 : /// |R_mesh| should eventually stop decreasing as iterations proceed,
129 : /// hitting a floor that represents the numerical truncation error of
130 : /// the solution.
131 : ///
132 : /// The convergence criterion looks at both |R| and |R_mesh|: Once
133 : /// |R| is small enough and |R_mesh| has stabilized, it is pointless
134 : /// to keep iterating (even though one could iterate until |R| reaches
135 : /// roundoff).
136 : ///
137 : /// Problems with convergence of the apparent horizon finder can often
138 : /// be diagnosed by looking at the behavior of |R| and |R_mesh| as a
139 : /// function of iteration.
140 : ///
141 : template <typename InterpolationTargetTag, typename Frame>
142 1 : struct FindApparentHorizon
143 : : tt::ConformsTo<intrp::protocols::PostInterpolationCallback> {
144 0 : using const_global_cache_tags =
145 : Parallel::get_const_global_cache_tags_from_actions<
146 : typename InterpolationTargetTag::post_horizon_find_callbacks>;
147 : template <typename DbTags, typename Metavariables, typename TemporalId>
148 0 : static bool apply(
149 : const gsl::not_null<db::DataBox<DbTags>*> box,
150 : const gsl::not_null<Parallel::GlobalCache<Metavariables>*> cache,
151 : const TemporalId& temporal_id) {
152 : bool horizon_finder_failed = false;
153 :
154 : const size_t current_iteration =
155 : get<::ah::Tags::FastFlow>(*box).current_iteration();
156 :
157 : if (current_iteration == 0) {
158 : // If we get here, we are in a new apparent horizon search, as
159 : // opposed to a subsequent iteration of the same horizon search.
160 : //
161 : // So put new initial guess into ylm::Tags::Strahlkorper<Frame>.
162 : // We need to do this now, and not at the end of the previous horizon
163 : // search, because only now do we know the temporal_id of this horizon
164 : // search.
165 : db::mutate<ylm::Tags::Strahlkorper<Frame>,
166 : ylm::Tags::PreviousStrahlkorpers<Frame>,
167 : ah::Tags::PreviousIterationStrahlkorper<Frame>,
168 : ah::Tags::FailedInterpolationIterations>(
169 : [&temporal_id](
170 : const gsl::not_null<ylm::Strahlkorper<Frame>*> strahlkorper,
171 : const gsl::not_null<
172 : std::deque<std::pair<double, ylm::Strahlkorper<Frame>>>*>
173 : previous_strahlkorpers,
174 : const gsl::not_null<ylm::Strahlkorper<Frame>*>
175 : previous_fast_flow_strahlkorper,
176 : const gsl::not_null<size_t*> failed_interpolation_iterations) {
177 : // Reset to zero. Will only be used starting from the second
178 : // iteration
179 : *failed_interpolation_iterations = 0;
180 :
181 : // If we have zero previous_strahlkorpers, then the
182 : // initial guess is already in strahlkorper, so do
183 : // nothing.
184 : //
185 : // If we have one previous_strahlkorper, then we have had
186 : // a successful horizon find, and the initial guess for the
187 : // next horizon find is already in strahlkorper, so
188 : // again we do nothing.
189 : //
190 : // If we have 2 previous_strahlkorpers and the time of the second
191 : // one is a NaN, this means that the corresponding
192 : // previous_strahlkorper is the original initial guess, so
193 : // again we do nothing.
194 : //
195 : // If we have 2 valid previous_strahlkorpers, then
196 : // we set the initial guess by linear extrapolation in time
197 : // using the last 2 previous_strahlkorpers.
198 : //
199 : // If we have 3 valid previous_strahlkorpers, then
200 : // we set the initial guess by quadratic extrapolation in time
201 : // using the last 3 previous_strahlkorpers.
202 : //
203 : // For extrapolation, we assume that
204 : // * Expansion center of all the Strahlkorpers are equal.
205 : // * Maximum L of all the Strahlkorpers are equal.
206 : // It is easy to relax the max L assumption once we start
207 : // adaptively changing the L of the strahlkorpers.
208 : if (previous_strahlkorpers->size() > 1 and
209 : not std::isnan((*previous_strahlkorpers)[1].first)) {
210 : if (previous_strahlkorpers->size() > 2 and
211 : not std::isnan((*previous_strahlkorpers)[2].first)) {
212 : // Quadratic extrapolation
213 : const double new_time =
214 : InterpolationTarget_detail::get_temporal_id_value(
215 : temporal_id);
216 : const double dt_0 =
217 : (*previous_strahlkorpers)[0].first - new_time;
218 : const double dt_1 =
219 : (*previous_strahlkorpers)[1].first - new_time;
220 : const double dt_2 =
221 : (*previous_strahlkorpers)[2].first - new_time;
222 : const double fac_0 =
223 : dt_1 * dt_2 / ((dt_1 - dt_0) * (dt_2 - dt_0));
224 : const double fac_1 =
225 : dt_0 * dt_2 / ((dt_2 - dt_1) * (dt_0 - dt_1));
226 : const double fac_2 = 1.0 - fac_0 - fac_1;
227 : strahlkorper->coefficients() =
228 : fac_0 * (*previous_strahlkorpers)[0].second.coefficients() +
229 : fac_1 * (*previous_strahlkorpers)[1].second.coefficients() +
230 : fac_2 * (*previous_strahlkorpers)[2].second.coefficients();
231 : } else {
232 : // Linear extrapolation
233 : const double new_time =
234 : InterpolationTarget_detail::get_temporal_id_value(
235 : temporal_id);
236 : const double dt_0 =
237 : (*previous_strahlkorpers)[0].first - new_time;
238 : const double dt_1 =
239 : (*previous_strahlkorpers)[1].first - new_time;
240 : const double fac_0 = dt_1 / (dt_1 - dt_0);
241 : const double fac_1 = 1.0 - fac_0;
242 : strahlkorper->coefficients() =
243 : fac_0 * (*previous_strahlkorpers)[0].second.coefficients() +
244 : fac_1 * (*previous_strahlkorpers)[1].second.coefficients();
245 : }
246 : }
247 :
248 : // First iteration the previous is just the initial guess
249 : (*previous_fast_flow_strahlkorper) = *strahlkorper;
250 : },
251 : box);
252 : }
253 :
254 : // Deal with the possibility that some of the points might be
255 : // outside of the Domain.
256 : const auto& indices_of_invalid_pts =
257 : db::get<Tags::IndicesOfInvalidInterpPoints<TemporalId>>(*box);
258 : size_t& failed_interpolation_iterations =
259 : db::get_mutable_reference<ah::Tags::FailedInterpolationIterations>(box);
260 : if (indices_of_invalid_pts.count(temporal_id) > 0 and
261 : not indices_of_invalid_pts.at(temporal_id).empty()) {
262 : ++failed_interpolation_iterations;
263 :
264 : const auto& options = Parallel::get<
265 : intrp::Tags::ApparentHorizon<InterpolationTargetTag, Frame>>(*cache);
266 :
267 : // Can't recover from the first iteration or if we've exceeded our number
268 : // of attempts
269 : if (current_iteration > 0 and failed_interpolation_iterations <=
270 : options.max_interpolation_retries) {
271 : // Move the new trial surface halfway between the current surface and
272 : // the previous surface
273 : db::mutate<ylm::Tags::Strahlkorper<Frame>>(
274 : [](const gsl::not_null<ylm::Strahlkorper<Frame>*> strahlkorper,
275 : const ylm::Strahlkorper<Frame>&
276 : previous_strahlkorper_iteration) {
277 : strahlkorper->coefficients() +=
278 : 0.5 * (previous_strahlkorper_iteration.coefficients() -
279 : strahlkorper->coefficients());
280 : },
281 : box, db::get<ah::Tags::PreviousIterationStrahlkorper<Frame>>(*box));
282 :
283 : // Resend the new coords at the same time and iteration. Use 0 for the
284 : // array index because it's always 0 due to the targets being
285 : // singletons. No need to worry about the invalid points since those get
286 : // reset within SendPointsToInterpolator
287 : Actions::SendPointsToInterpolator<InterpolationTargetTag>::
288 : template apply<::intrp::InterpolationTarget<
289 : Metavariables, InterpolationTargetTag>>(
290 : *box, *cache, 0, temporal_id, current_iteration,
291 : failed_interpolation_iterations);
292 :
293 : // We return false because we don't want this iteration to clean
294 : // up the volume data, since we still need it.
295 : return false;
296 : } else {
297 : tmpl::for_each<
298 : typename InterpolationTargetTag::horizon_find_failure_callbacks>(
299 : [&box, &cache, &temporal_id](auto callback_v) {
300 : using callback = tmpl::type_from<decltype(callback_v)>;
301 : callback::template apply<InterpolationTargetTag>(
302 : *box, *cache, temporal_id,
303 : FastFlow::Status::InterpolationFailure);
304 : });
305 : horizon_finder_failed = true;
306 : }
307 : }
308 :
309 : // Reset now that we were able to interpolate (or if we are abandoning this
310 : // horizon find)
311 : failed_interpolation_iterations = 0;
312 :
313 : if (not horizon_finder_failed) {
314 : const auto& verbosity =
315 : db::get<logging::Tags::Verbosity<InterpolationTargetTag>>(*box);
316 : const auto& inv_g =
317 : db::get<::gr::Tags::InverseSpatialMetric<DataVector, 3, Frame>>(*box);
318 : const auto& ex_curv =
319 : db::get<::gr::Tags::ExtrinsicCurvature<DataVector, 3, Frame>>(*box);
320 : const auto& christoffel = db::get<
321 : ::gr::Tags::SpatialChristoffelSecondKind<DataVector, 3, Frame>>(*box);
322 :
323 : std::pair<FastFlow::Status, FastFlow::IterInfo> status_and_info;
324 :
325 : // Do a FastFlow iteration.
326 : db::mutate<::ah::Tags::FastFlow, ylm::Tags::Strahlkorper<Frame>,
327 : ah::Tags::PreviousIterationStrahlkorper<Frame>>(
328 : [&inv_g, &ex_curv, &christoffel, &status_and_info](
329 : const gsl::not_null<::FastFlow*> fast_flow,
330 : const gsl::not_null<ylm::Strahlkorper<Frame>*> strahlkorper,
331 : const gsl::not_null<ylm::Strahlkorper<Frame>*>
332 : previous_fast_flow_strahlkorper) {
333 : // Set this before we iterate the fast flow. Note that this will
334 : // have to be updated once we have adaptive horizon finding
335 : (*previous_fast_flow_strahlkorper) = *strahlkorper;
336 :
337 : status_and_info = fast_flow->template iterate_horizon_finder<Frame>(
338 : strahlkorper, inv_g, ex_curv, christoffel);
339 : },
340 : box);
341 :
342 : // Determine whether we have converged, whether we need another step,
343 : // or whether we have encountered an error.
344 :
345 : const auto& status = status_and_info.first;
346 : const auto& info = status_and_info.second;
347 : const auto has_converged = converged(status);
348 :
349 : if (verbosity >= ::Verbosity::Verbose or
350 : (verbosity >= ::Verbosity::Quiet and has_converged)) {
351 : Parallel::printf(
352 : "%s: t=%.6g: its=%d: %.1e<R<%.0e, |R|=%.1g, "
353 : "|R_grid|=%.1g, %.4g<r<%.4g\n",
354 : pretty_type::name<InterpolationTargetTag>(),
355 : InterpolationTarget_detail::get_temporal_id_value(temporal_id),
356 : info.iteration, info.min_residual, info.max_residual,
357 : info.residual_ylm, info.residual_mesh, info.r_min, info.r_max);
358 : }
359 :
360 : if (status == FastFlow::Status::SuccessfulIteration) {
361 : // Do another iteration of the same horizon search.
362 : const auto& current_id =
363 : db::get<intrp::Tags::CurrentTemporalId<TemporalId>>(*box);
364 : using ::operator<<;
365 : ASSERT(current_id.has_value(),
366 : "While horizon finding, the current temporal id doesn't have a "
367 : "value, but it should.");
368 : // The iteration of these new coords is the fast flow iteration + 1
369 : // because the zeroth iteration was the initial guess. Use 0 for the
370 : // array index because it's always 0 due to the targets being singletons
371 : Actions::SendPointsToInterpolator<InterpolationTargetTag>::
372 : template apply<::intrp::InterpolationTarget<
373 : Metavariables, InterpolationTargetTag>>(
374 : *box, *cache, 0, current_id.value(), info.iteration + 1);
375 : // We return false because we don't want this iteration to clean
376 : // up the volume data, since we are using it for the next iteration
377 : // (i.e. the simple_action that we just called).
378 : return false;
379 : } else if (not has_converged) {
380 : tmpl::for_each<
381 : typename InterpolationTargetTag::horizon_find_failure_callbacks>(
382 : [&box, &cache, &temporal_id, &status](auto callback_v) {
383 : using callback = tmpl::type_from<decltype(callback_v)>;
384 : callback::template apply<InterpolationTargetTag>(
385 : *box, *cache, temporal_id, status);
386 : });
387 : horizon_finder_failed = true;
388 : }
389 : }
390 :
391 : // If it failed, don't update any variables, just reset the Strahlkorper to
392 : // it's previous value
393 : if (horizon_finder_failed) {
394 : db::mutate<ylm::Tags::Strahlkorper<Frame>>(
395 : [](const gsl::not_null<ylm::Strahlkorper<Frame>*> strahlkorper,
396 : const std::deque<std::pair<double, ylm::Strahlkorper<Frame>>>&
397 : previous_strahlkorpers) {
398 : // Don't keep a partially-converged strahlkorper in the
399 : // DataBox. Reset to either the original initial guess or
400 : // to the last-found Strahlkorper (whichever one happens
401 : // to be in previous_strahlkorpers).
402 : *strahlkorper = previous_strahlkorpers.front().second;
403 : },
404 : box, db::get<ylm::Tags::PreviousStrahlkorpers<Frame>>(*box));
405 : } else {
406 : // The interpolated variables
407 : // Tags::Variables<InterpolationTargetTag::vars_to_interpolate_to_target>
408 : // have been interpolated from the volume to the points on the
409 : // prolonged_strahlkorper, not to the points on the actual
410 : // strahlkorper. So here we do a restriction of these
411 : // quantities onto the actual strahlkorper.
412 :
413 : // Type alias to make code more understandable.
414 : using vars_tags =
415 : typename InterpolationTargetTag::vars_to_interpolate_to_target;
416 : db::mutate_apply<
417 : tmpl::list<::Tags::Variables<vars_tags>>,
418 : tmpl::list<ylm::Tags::Strahlkorper<Frame>, ::ah::Tags::FastFlow>>(
419 : [](const gsl::not_null<Variables<vars_tags>*> vars,
420 : const ylm::Strahlkorper<Frame>& strahlkorper,
421 : const FastFlow& fast_flow) {
422 : const size_t L_mesh = fast_flow.current_l_mesh(strahlkorper);
423 : const auto prolonged_strahlkorper =
424 : ylm::Strahlkorper<Frame>(L_mesh, L_mesh, strahlkorper);
425 : auto new_vars = ::Variables<vars_tags>(
426 : strahlkorper.ylm_spherepack().physical_size());
427 :
428 : tmpl::for_each<vars_tags>([&strahlkorper, &prolonged_strahlkorper,
429 : &vars, &new_vars](auto tag_v) {
430 : using tag = typename decltype(tag_v)::type;
431 : const auto& old_var = get<tag>(*vars);
432 : auto& new_var = get<tag>(new_vars);
433 : auto old_iter = old_var.begin();
434 : auto new_iter = new_var.begin();
435 : for (; old_iter != old_var.end() and new_iter != new_var.end();
436 : ++old_iter, ++new_iter) {
437 : *new_iter = strahlkorper.ylm_spherepack().spec_to_phys(
438 : prolonged_strahlkorper.ylm_spherepack().prolong_or_restrict(
439 : prolonged_strahlkorper.ylm_spherepack().phys_to_spec(
440 : *old_iter),
441 : strahlkorper.ylm_spherepack()));
442 : }
443 : });
444 : *vars = std::move(new_vars);
445 : },
446 : box);
447 :
448 : // Compute Strahlkorper Cartesian coordinates in Inertial frame
449 : // if the current frame is not inertial.
450 : if constexpr (not std::is_same_v<Frame, ::Frame::Inertial>) {
451 : db::mutate_apply<
452 : tmpl::list<ylm::Tags::CartesianCoords<::Frame::Inertial>>,
453 : tmpl::list<ylm::Tags::Strahlkorper<Frame>,
454 : domain::Tags::Domain<Metavariables::volume_dim>>>(
455 : [&cache, &temporal_id](
456 : const gsl::not_null<tnsr::I<DataVector, 3, ::Frame::Inertial>*>
457 : inertial_strahlkorper_coords,
458 : const ylm::Strahlkorper<Frame>& strahlkorper,
459 : const Domain<Metavariables::volume_dim>& domain) {
460 : // Note that functions_of_time must already be up to
461 : // date at temporal_id because they were used in the AH
462 : // search above.
463 : const auto& functions_of_time =
464 : get<domain::Tags::FunctionsOfTime>(*cache);
465 : strahlkorper_coords_in_different_frame(
466 : inertial_strahlkorper_coords, strahlkorper, domain,
467 : functions_of_time,
468 : InterpolationTarget_detail::get_temporal_id_value(
469 : temporal_id));
470 : },
471 : box);
472 : }
473 :
474 : // Update the previous strahlkorpers. We do this before the callbacks
475 : // in case any of the callbacks need the previous strahlkorpers with the
476 : // current strahlkorper already in it.
477 : db::mutate<ylm::Tags::Strahlkorper<Frame>,
478 : ylm::Tags::PreviousStrahlkorpers<Frame>>(
479 : [&temporal_id](
480 : const gsl::not_null<ylm::Strahlkorper<Frame>*> strahlkorper,
481 : const gsl::not_null<
482 : std::deque<std::pair<double, ylm::Strahlkorper<Frame>>>*>
483 : previous_strahlkorpers) {
484 : // This is the number of previous strahlkorpers that we
485 : // keep around.
486 : const size_t num_previous_strahlkorpers = 3;
487 :
488 : // Save a new previous_strahlkorper.
489 : previous_strahlkorpers->emplace_front(
490 : InterpolationTarget_detail::get_temporal_id_value(temporal_id),
491 : *strahlkorper);
492 :
493 : // Remove old previous_strahlkorpers that are no longer relevant.
494 : while (previous_strahlkorpers->size() >
495 : num_previous_strahlkorpers) {
496 : previous_strahlkorpers->pop_back();
497 : }
498 : },
499 : box);
500 :
501 : // Finally call callbacks
502 : tmpl::for_each<
503 : typename InterpolationTargetTag::post_horizon_find_callbacks>(
504 : [&box, &cache, &temporal_id](auto callback_v) {
505 : using callback = tmpl::type_from<decltype(callback_v)>;
506 : callback::apply(*box, *cache, temporal_id);
507 : });
508 : }
509 :
510 : // Prepare for finding horizon at a new time. Regardless of if we failed or
511 : // not, we reset fast flow.
512 : db::mutate<::ah::Tags::FastFlow>(
513 : [](const gsl::not_null<::FastFlow*> fast_flow) {
514 : fast_flow->reset_for_next_find();
515 : },
516 : box);
517 :
518 : // We return true because we are now done with all the volume data
519 : // at this temporal_id, so we want it cleaned up.
520 : return true;
521 : }
522 : };
523 : } // namespace intrp::callbacks
|
