FindApparentHorizon.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <utility>
7 
8 #include "ApparentHorizons/FastFlow.hpp"
10 #include "ErrorHandling/Error.hpp"
11 #include "Informer/Verbosity.hpp"
13 #include "Parallel/Invoke.hpp"
14 #include "Parallel/Printf.hpp"
15 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
16 #include "Utilities/Gsl.hpp"
17 #include "Utilities/PrettyType.hpp"
18 
19 /// \cond
20 namespace StrahlkorperTags {
21 template <typename Frame>
22 struct Strahlkorper;
23 } // namespace StrahlkorperTags
24 namespace db {
25 template <typename TagsList>
26 class DataBox;
27 } // namespace db
28 namespace intrp {
29 template <class Metavariables, typename InterpolationTargetTag>
30 struct InterpolationTarget;
31 namespace Tags {
32 template <typename Metavariables>
33 struct TemporalIds;
34 } // namespace Tags
35 } // namespace intrp
36 namespace ah {
37 namespace Tags {
38 struct FastFlow;
39 } // namespace Tags
40 } // namespace ah
41 namespace Tags {
42 struct Verbosity;
43 } // namespace Tags
44 template <typename Frame>
45 class Strahlkorper;
46 /// \endcond
47 
48 namespace intrp {
49 namespace callbacks {
50 
51 /// \brief post interpolation callback (see InterpolationTarget)
52 /// that does a FastFlow iteration and triggers another one until convergence.
53 ///
54 /// Assumes that InterpolationTargetTag contains an additional
55 /// struct called `post_horizon_find_callback`, which has a function
56 ///```
57 /// static void apply(const DataBox<DbTags>&,
58 /// const intrp::ConstGlobalCache<Metavariables>&,
59 /// const Metavariables::temporal_id&) noexcept;
60 ///```
61 /// that is called if the FastFlow iteration has converged.
62 ///
63 /// Uses:
64 /// - Metavariables:
65 /// - `temporal_id`
66 /// - `domain_frame`
67 /// - DataBox:
68 /// - `::Tags::Verbosity`
69 /// - `::gr::Tags::InverseSpatialMetric<3,Frame>`
70 /// - `::gr::Tags::ExtrinsicCurvature<3,Frame>`
71 /// - `::gr::Tags::SpatialChristoffelSecondKind<3,Frame>`
72 /// - `::ah::Tags::FastFlow`
73 /// - `StrahlkorperTags::Strahlkorper<Frame>`
74 ///
75 /// Modifies:
76 /// - DataBox:
77 /// - `::ah::Tags::FastFlow`
78 /// - `StrahlkorperTags::Strahlkorper<Frame>`
79 ///
80 /// This is an InterpolationTargetTag::post_interpolation_callback;
81 /// see InterpolationTarget for a description of InterpolationTargetTag.
82 template <typename InterpolationTargetTag>
84  template <typename DbTags, typename Metavariables>
85  static bool apply(
88  const typename Metavariables::temporal_id::type& temporal_id) noexcept {
89  const auto& verbosity = db::get<::Tags::Verbosity>(*box);
90  const auto& inv_g = db::get<::gr::Tags::InverseSpatialMetric<
91  3, typename Metavariables::domain_frame>>(*box);
92  const auto& ex_curv = db::get<::gr::Tags::ExtrinsicCurvature<
93  3, typename Metavariables::domain_frame>>(*box);
94  const auto& christoffel = db::get<::gr::Tags::SpatialChristoffelSecondKind<
95  3, typename Metavariables::domain_frame>>(*box);
96 
98 
99  // Do a FastFlow iteration.
101  typename Metavariables::domain_frame>>(
102  box, [&inv_g, &ex_curv, &christoffel, &status_and_info ](
103  const gsl::not_null<::FastFlow*> fast_flow,
104  const gsl::not_null<
106  strahlkorper) noexcept {
107  status_and_info = fast_flow->template iterate_horizon_finder<
108  typename Metavariables::domain_frame>(strahlkorper, inv_g,
109  ex_curv, christoffel);
110  });
111 
112  // Determine whether we have converged, whether we need another step,
113  // or whether we have encountered an error.
114 
115  const auto& status = status_and_info.first;
116  const auto& info = status_and_info.second;
117  const auto has_converged = converged(status);
118 
119  if (verbosity > ::Verbosity::Quiet or
120  (verbosity > ::Verbosity::Silent and has_converged)) {
121  // The things printed out here are:
122  // its = current iteration number.
123  // R = min and max of residual over all prolonged grid points.
124  // |R| = L2 norm of residual, counting only L modes solved for.
125  // |R_mesh| = L2 norm of residual over prolonged grid points.
126  // r = min and max radius of trial horizon surface.
127  //
128  // Difference between |R| and |R_mesh|:
129  // The horizon is represented in a Y_lm expansion up to l=l_surface;
130  // the residual |R| represents the failure of that surface to satisfy
131  // the apparent horizon equation.
132  //
133  // However, at each iteration we also interpolate the horizon surface
134  // to a higher resolution ("prolongation"). The prolonged surface
135  // includes Ylm coefficents up to l=l_mesh, where l_mesh > l_surface.
136  // The residual computed on this higher-resolution surface is |R_mesh|.
137  //
138  // As iterations proceed, |R| should decrease until it reaches numerical
139  // roundoff error, because we are varying all the Y_lm coefficients up to
140  // l=l_surface to minimize the residual. However, |R_mesh| should
141  // eventually stop decreasing as iterations proceed, hitting a
142  // floor that represents the numerical truncation error of the solution.
143  //
144  // The convergence criterion looks at both |R| and |R_mesh|: Once
145  // |R| is small enough and |R_mesh| has stabilized, it is pointless
146  // to keep iterating (even though one could iterate until |R| reaches
147  // roundoff).
148  //
149  // Problems with convergence of the apparent horizon finder can often
150  // be diagnosed by looking at the behavior of |R| and |R_mesh| as a
151  // function of iteration.
153  "%s: its=%d: %.1e<R<%.0e, |R|=%.1g, "
154  "|R_grid|=%.1g, %.4g<r<%.4g\n",
155  pretty_type::short_name<InterpolationTargetTag>(), info.iteration,
156  info.min_residual, info.max_residual, info.residual_ylm,
157  info.residual_mesh, info.r_min, info.r_max);
158  }
159 
160  if (status == FastFlow::Status::SuccessfulIteration) {
161  // Do another iteration of the same horizon search.
162  const auto& temporal_ids =
163  db::get<intrp::Tags::TemporalIds<Metavariables>>(*box);
164  auto& interpolation_target = Parallel::get_parallel_component<
166  *cache);
168  typename InterpolationTargetTag::compute_target_points>(
169  interpolation_target, temporal_ids.front());
170  // We return false because we don't want this iteration to clean
171  // up the volume data, since we are using it for the next iteration
172  // (i.e. the simple_action that we just called).
173  return false;
174  } else if (not has_converged) {
175  ERROR("Apparent horizon finder "
176  << pretty_type::short_name<InterpolationTargetTag>()
177  << " failed, reason = " << status);
178  }
179  // If we get here, the horizon finder has converged.
180 
182  temporal_id);
183 
184  // Prepare for finding horizon at a new time.
185  // For now, the initial guess for the new
186  // horizon is the result of the old one, so we don't need
187  // to modify the strahlkorper.
188  // Eventually we will hold more than one previous guess and
189  // will do time-extrapolation to set the next guess.
190  db::mutate<::ah::Tags::FastFlow>(
191  box, [](const gsl::not_null<::FastFlow*> fast_flow) noexcept {
192  fast_flow->reset_for_next_find();
193  });
194  // We return true because we are now done with all the volume data
195  // at this temporal_id, so we want it cleaned up.
196  return true;
197  }
198 };
199 } // namespace callbacks
200 } // namespace intrp
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:35
void mutate(const gsl::not_null< DataBox< TagList > *> box, Invokable &&invokable, Args &&... args) noexcept
Allows changing the state of one or more non-computed elements in the DataBox.
Definition: DataBox.hpp:1099
Definition: AddTemporalIdsToInterpolationTarget.hpp:17
Definition: Tags.hpp:31
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args)
Apply the function f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1609
Defines classes and functions used for manipulating DataBox&#39;s.
Verbosity
Indicates how much informative output a class should output.
Definition: Verbosity.hpp:16
Defines Parallel::printf for writing to stdout.
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
Definition: DataBoxTag.hpp:29
A Charm++ chare that caches constant data once per Charm++ node.
Definition: ConstGlobalCache.hpp:76
Namespace for DataBox related things.
Definition: DataBox.hpp:33
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:1211
Contains a pretty_type library to write types in a "pretty" format.
A star-shaped surface expanded in spherical harmonics.
Definition: Strahlkorper.hpp:21
Definition: ComputeItems.hpp:23
auto get_parallel_component(ConstGlobalCache< Metavariables > &cache) noexcept -> Parallel::proxy_from_parallel_component< ConstGlobalCache_detail::get_component_if_mocked< typename Metavariables::component_list, ParallelComponentTag >> &
Access the Charm++ proxy associated with a ParallelComponent.
Definition: ConstGlobalCache.hpp:163
Defines functions and classes from the GSL.
void simple_action(Proxy &&proxy) noexcept
Invoke a simple action on proxy
Definition: Invoke.hpp:112
post interpolation callback (see InterpolationTarget) that does a FastFlow iteration and triggers ano...
Definition: FindApparentHorizon.hpp:83
Definition: Tags.hpp:26
Defines class template ConstGlobalCache.
Defines macro ERROR.
Holds tags and ComputeItems associated with a Strahlkorper.
Definition: Tags.cpp:15
void printf(const std::string &format, Args &&... args)
Print an atomic message to stdout with C printf usage.
Definition: Printf.hpp:100
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12
ParallelComponent representing a set of points to be interpolated to and a function to call upon inte...
Definition: InterpolationTarget.hpp:99
Fast flow method for finding apparent horizons.
Definition: FastFlow.hpp:34
Definition: Tags.hpp:96
Tag referring to a Strahlkorper
Definition: Tags.hpp:39