FastFlow.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <limits>
8 #include <ostream>
9 #include <utility>
10 
12 #include "Options/Options.hpp"
14 #include "Utilities/TMPL.hpp"
15 
16 /// \cond
17 class DataVector;
18 namespace PUP {
19 class er;
20 } // namespace PUP
21 namespace gsl {
22 template <class T>
23 class not_null;
24 } // namespace gsl
25 template <typename>
26 class Strahlkorper;
27 /// \endcond
28 
29 /// \ingroup SurfacesGroup
30 /// \brief Fast flow method for finding apparent horizons.
31 ///
32 /// \details Based on \cite Gundlach1997us.
33 // The method is iterative.
34 class FastFlow {
35  public:
36  enum class FlowType { Jacobi, Curvature, Fast };
37 
38  // Error conditions are negative, successes are nonnegative
39  enum class Status {
40  // SuccessfulIteration means we should iterate again.
41  SuccessfulIteration = 0,
42  // The following indicate convergence, can stop iterating.
43  AbsTol = 1,
44  TruncationTol = 2,
45  // The following indicate errors.
46  MaxIts = -1,
47  NegativeRadius = -2,
48  DivergenceError = -3
49  };
50 
51  /// Holds information about an iteration of the algorithm.
52  struct IterInfo {
53  size_t iteration{std::numeric_limits<size_t>::max()};
60  };
61 
62  struct Flow {
63  using type = FlowType;
64  static constexpr OptionString help = {
65  "Flow method: Jacobi, Curvature, or Fast"};
66  static type default_value() noexcept { return FlowType::Fast; }
67  };
68 
69  struct Alpha {
70  using type = double;
71  static constexpr OptionString help = {
72  "Alpha parameter in PRD 57, 863 (1998)"};
73  static type default_value() noexcept { return 1.0; }
74  };
75 
76  struct Beta {
77  using type = double;
78  static constexpr OptionString help = {
79  "Beta parameter in PRD 57, 863 (1998)"};
80  static type default_value() noexcept { return 0.5; }
81  };
82 
83  struct AbsTol {
84  using type = double;
85  static constexpr OptionString help = {
86  "Convergence found if R_{Y_lm} < AbsTol"};
87  static type default_value() noexcept { return 1.e-12; }
88  };
89 
90  struct TruncationTol {
91  using type = double;
92  static constexpr OptionString help = {
93  "Convergence found if R_{Y_lm} < TruncationTol*R_{mesh}"};
94  static type default_value() noexcept { return 1.e-2; }
95  };
96 
97  struct DivergenceTol {
98  using type = double;
99  static constexpr OptionString help = {
100  "Fraction that residual can increase before dying"};
101  static type default_value() noexcept { return 1.2; }
102  static type lower_bound() noexcept { return 1.0; }
103  };
104 
105  struct DivergenceIter {
106  using type = size_t;
107  static constexpr OptionString help = {
108  "Num iterations residual can increase before dying"};
109  static type default_value() noexcept { return 5; }
110  };
111 
112  struct MaxIts {
113  using type = size_t;
114  static constexpr OptionString help = {"Maximum number of iterations."};
115  static type default_value() noexcept { return 100; }
116  };
117 
118  using options = tmpl::list<Flow, Alpha, Beta, AbsTol, TruncationTol,
120 
121  static constexpr OptionString help{
122  "Find a Strahlkorper using a 'fast flow' method.\n"
123  "Based on Gundlach, PRD 57, 863 (1998).\n"
124  "Expands the surface in terms of spherical harmonics Y_lm up to a given\n"
125  "l_surface, and varies the coefficients S_lm where 0<=l<=l_surface to\n"
126  "minimize the residual of the apparent horizon equation. Also keeps\n"
127  "another representation of the surface that is expanded up to\n"
128  "l_mesh > l_surface. Let R_{Y_lm} be the residual computed using the\n"
129  "surface represented up to l_surface; this residual can in principle be\n"
130  "lowered to machine roundoff by enough iterations. Let R_{mesh} be the\n"
131  "residual computed using the surface represented up to l_mesh; this\n"
132  "residual represents the truncation error, since l_mesh>l_surface and\n"
133  "since coefficients S_lm with l>l_surface are not modified in the\n"
134  "iteration.\n\n"
135  "Convergence is achieved if R_{Y_lm}< TruncationTol*R_{mesh}, or if\n"
136  "R_{Y_lm}<AbsTol, where TruncationTol and AbsTol are input parameters.\n"
137  "If instead |R_{mesh}|_i > DivergenceTol * min_{j}(|R_{mesh}|_j) where\n"
138  "i is the iteration index and j runs from 0 to i-DivergenceIter, then\n"
139  "FastFlow exits with Status::DivergenceError. Here DivergenceIter and\n"
140  "DivergenceTol are input parameters."};
141 
142  FastFlow(Flow::type flow, Alpha::type alpha, Beta::type beta,
143  AbsTol::type abs_tol, TruncationTol::type trunc_tol,
144  DivergenceTol::type divergence_tol,
145  DivergenceIter::type divergence_iter, MaxIts::type max_its) noexcept;
146 
147  FastFlow() noexcept
148  : FastFlow(FlowType::Fast, 1.0, 0.5, 1.e-12, 1.e-2, 1.2, 5, 100) {}
149 
150  FastFlow(const FastFlow& /*rhs*/) = default;
151  FastFlow& operator=(const FastFlow& /*rhs*/) = default;
152  FastFlow(FastFlow&& /*rhs*/) noexcept = default;
153  FastFlow& operator=(FastFlow&& /*rhs*/) noexcept = default;
154  ~FastFlow() = default;
155 
156  // clang-tidy: no runtime references
157  void pup(PUP::er& p) noexcept; // NOLINT
158 
159  /// Evaluate residuals and compute the next iteration. If
160  /// Status==SuccessfulIteration, then `current_strahlkorper` is
161  /// modified and `current_iteration()` is incremented. Otherwise, we
162  /// end with success or failure, and neither `current_strahlkorper`
163  /// nor `current_iteration()` is changed.
164  template <typename Frame>
165  std::pair<Status, IterInfo> iterate_horizon_finder(
166  gsl::not_null<Strahlkorper<Frame>*> current_strahlkorper,
167  const tnsr::II<DataVector, 3, Frame>& upper_spatial_metric,
168  const tnsr::ii<DataVector, 3, Frame>& extrinsic_curvature,
169  const tnsr::Ijj<DataVector, 3, Frame>& christoffel_2nd_kind) noexcept;
170 
171  size_t current_iteration() const noexcept { return current_iter_; }
172 
173  /// Given a Strahlkorper defined up to some maximum Y_lm l called
174  /// l_surface, returns a larger value of l, l_mesh, that is used for
175  /// evaluating convergence.
176  template <typename Frame>
177  size_t current_l_mesh(const Strahlkorper<Frame>& strahlkorper) const noexcept;
178 
179  /// Resets the finder.
181  current_iter_ = 0;
182  }
183 
184  private:
185  friend bool operator==(const FastFlow& /*lhs*/,
186  const FastFlow& /*rhs*/) noexcept;
187  double alpha_, beta_, abs_tol_, trunc_tol_, divergence_tol_;
188  size_t divergence_iter_, max_its_;
189  FlowType flow_;
190  size_t current_iter_;
191  double previous_residual_mesh_norm_, min_residual_mesh_norm_;
192  size_t iter_at_min_residual_mesh_norm_;
193 };
194 
195 SPECTRE_ALWAYS_INLINE bool converged(const FastFlow::Status& status) noexcept {
196  return static_cast<int>(status) > 0;
197 }
198 
199 std::ostream& operator<<(std::ostream& os,
200  const FastFlow::FlowType& flow_type) noexcept;
201 
202 std::ostream& operator<<(std::ostream& os,
203  const FastFlow::Status& status) noexcept;
204 
205 SPECTRE_ALWAYS_INLINE bool operator!=(const FastFlow& lhs,
206  const FastFlow& rhs) noexcept {
207  return not(lhs == rhs);
208 }
209 
210 template <>
211 struct create_from_yaml<FastFlow::FlowType> {
212  static FastFlow::FlowType create(const Option& options);
213 };
Holds information about an iteration of the algorithm.
Definition: FastFlow.hpp:52
Definition: Strahlkorper.hpp:14
The type that options are passed around as. Contains YAML node data and an OptionContext.
Definition: Options.hpp:103
Implementations from the Guideline Support Library.
Definition: ConservativeFromPrimitive.hpp:10
T signaling_NaN(T... args)
Used by the parser to create an object. The default action is to parse options using T::options...
Definition: Options.hpp:143
Definition: FastFlow.hpp:62
Defines classes and functions for making classes creatable from input files.
Definition: FastFlow.hpp:90
constexpr auto create(Args &&... args)
Create a new DataBox.
Definition: DataBox.hpp:1259
void reset_for_next_find() noexcept
Resets the finder.
Definition: FastFlow.hpp:180
Definition: FastFlow.hpp:97
Definition: FastFlow.hpp:83
tnsr::ii< DataVector, 3, Frame > extrinsic_curvature(const tnsr::ii< DataVector, 3, Frame > &grad_normal, const tnsr::i< DataVector, 3, Frame > &unit_normal_one_form, const tnsr::I< DataVector, 3, Frame > &unit_normal_vector) noexcept
Extrinsic curvature of a 2D Strahlkorper embedded in a 3D space.
Definition: FastFlow.hpp:112
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:26
#define SPECTRE_ALWAYS_INLINE
Always inline a function. Only use this if you benchmarked the code.
Definition: ForceInline.hpp:20
Definition: FastFlow.hpp:76
T max(T... args)
Definition: FastFlow.hpp:105
Defines macro to always inline a function.
Defines a list of useful type aliases for tensors.
A star-shaped surface expanded in spherical harmonics.
Definition: Strahlkorper.hpp:21
Stores a collection of function values.
Definition: DataVector.hpp:46
Definition: FastFlow.hpp:69
Wraps the template metaprogramming library used (brigand)
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12
Fast flow method for finding apparent horizons.
Definition: FastFlow.hpp:34