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