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 <ostream>
9 : #include <utility>
10 :
11 : #include "DataStructures/Tensor/TypeAliases.hpp"
12 : #include "Options/Options.hpp"
13 : #include "Options/String.hpp"
14 : #include "Utilities/ForceInline.hpp"
15 : #include "Utilities/TMPL.hpp"
16 :
17 : /// \cond
18 : class DataVector;
19 : namespace PUP {
20 : class er;
21 : } // namespace PUP
22 : namespace gsl {
23 : template <class T>
24 : class not_null;
25 : } // namespace gsl
26 : namespace ylm {
27 : template <typename>
28 : class Strahlkorper;
29 : } // namespace ylm
30 : /// \endcond
31 :
32 : /// \ingroup SurfacesGroup
33 : /// \brief Fast flow method for finding apparent horizons.
34 : ///
35 : /// \details Based on \cite Gundlach1997us.
36 : // The method is iterative.
37 1 : class FastFlow {
38 : public:
39 0 : enum class FlowType { Jacobi, Curvature, Fast };
40 :
41 : // Error conditions are negative, successes are nonnegative
42 0 : enum class Status {
43 : // SuccessfulIteration means we should iterate again.
44 : SuccessfulIteration = 0,
45 : // The following indicate convergence, can stop iterating.
46 : AbsTol = 1,
47 : TruncationTol = 2,
48 : // The following indicate errors.
49 : MaxIts = -1,
50 : NegativeRadius = -2,
51 : DivergenceError = -3,
52 : InterpolationFailure = -4
53 : };
54 :
55 : /// Holds information about an iteration of the algorithm.
56 1 : struct IterInfo {
57 0 : size_t iteration{std::numeric_limits<size_t>::max()};
58 0 : double r_min{std::numeric_limits<double>::signaling_NaN()},
59 0 : r_max{std::numeric_limits<double>::signaling_NaN()},
60 0 : min_residual{std::numeric_limits<double>::signaling_NaN()},
61 0 : max_residual{std::numeric_limits<double>::signaling_NaN()},
62 0 : residual_ylm{std::numeric_limits<double>::signaling_NaN()},
63 0 : residual_mesh{std::numeric_limits<double>::signaling_NaN()};
64 : };
65 :
66 0 : struct Flow {
67 0 : using type = FlowType;
68 0 : static constexpr Options::String help = {
69 : "Flow method: Jacobi, Curvature, or Fast"};
70 0 : static type suggested_value() { return FlowType::Fast; }
71 : };
72 :
73 0 : struct Alpha {
74 0 : using type = double;
75 0 : static constexpr Options::String help = {
76 : "Alpha parameter in PRD 57, 863 (1998)"};
77 0 : static type suggested_value() { return 1.0; }
78 : };
79 :
80 0 : struct Beta {
81 0 : using type = double;
82 0 : static constexpr Options::String help = {
83 : "Beta parameter in PRD 57, 863 (1998)"};
84 0 : static type suggested_value() { return 0.5; }
85 : };
86 :
87 0 : struct AbsTol {
88 0 : using type = double;
89 0 : static constexpr Options::String help = {
90 : "Convergence found if R_{Y_lm} < AbsTol"};
91 0 : static type suggested_value() { return 1.e-12; }
92 : };
93 :
94 0 : struct TruncationTol {
95 0 : using type = double;
96 0 : static constexpr Options::String help = {
97 : "Convergence found if R_{Y_lm} < TruncationTol*R_{mesh}"};
98 0 : static type suggested_value() { return 1.e-2; }
99 : };
100 :
101 0 : struct DivergenceTol {
102 0 : using type = double;
103 0 : static constexpr Options::String help = {
104 : "Fraction that residual can increase before dying"};
105 0 : static type suggested_value() { return 1.2; }
106 0 : static type lower_bound() { return 1.0; }
107 : };
108 :
109 0 : struct DivergenceIter {
110 0 : using type = size_t;
111 0 : static constexpr Options::String help = {
112 : "Num iterations residual can increase before dying"};
113 0 : static type suggested_value() { return 5; }
114 : };
115 :
116 0 : struct MaxIts {
117 0 : using type = size_t;
118 0 : static constexpr Options::String help = {"Maximum number of iterations."};
119 0 : static type suggested_value() { return 100; }
120 : };
121 :
122 0 : using options = tmpl::list<Flow, Alpha, Beta, AbsTol, TruncationTol,
123 : DivergenceTol, DivergenceIter, MaxIts>;
124 :
125 0 : static constexpr Options::String help{
126 : "Find a Strahlkorper using a 'fast flow' method.\n"
127 : "Based on Gundlach, PRD 57, 863 (1998).\n"
128 : "Expands the surface in terms of spherical harmonics Y_lm up to a given\n"
129 : "l_surface, and varies the coefficients S_lm where 0<=l<=l_surface to\n"
130 : "minimize the residual of the apparent horizon equation. Also keeps\n"
131 : "another representation of the surface that is expanded up to\n"
132 : "l_mesh > l_surface. Let R_{Y_lm} be the residual computed using the\n"
133 : "surface represented up to l_surface; this residual can in principle be\n"
134 : "lowered to machine roundoff by enough iterations. Let R_{mesh} be the\n"
135 : "residual computed using the surface represented up to l_mesh; this\n"
136 : "residual represents the truncation error, since l_mesh>l_surface and\n"
137 : "since coefficients S_lm with l>l_surface are not modified in the\n"
138 : "iteration.\n\n"
139 : "Convergence is achieved if R_{Y_lm}< TruncationTol*R_{mesh}, or if\n"
140 : "R_{Y_lm}<AbsTol, where TruncationTol and AbsTol are input parameters.\n"
141 : "If instead |R_{mesh}|_i > DivergenceTol * min_{j}(|R_{mesh}|_j) where\n"
142 : "i is the iteration index and j runs from 0 to i-DivergenceIter, then\n"
143 : "FastFlow exits with Status::DivergenceError. Here DivergenceIter and\n"
144 : "DivergenceTol are input parameters."};
145 :
146 0 : FastFlow(Flow::type flow, Alpha::type alpha, Beta::type beta,
147 : AbsTol::type abs_tol, TruncationTol::type trunc_tol,
148 : DivergenceTol::type divergence_tol,
149 : DivergenceIter::type divergence_iter, MaxIts::type max_its);
150 :
151 0 : FastFlow() : FastFlow(FlowType::Fast, 1.0, 0.5, 1.e-12, 1.e-2, 1.2, 5, 100) {}
152 :
153 0 : FastFlow(const FastFlow& /*rhs*/) = default;
154 0 : FastFlow& operator=(const FastFlow& /*rhs*/) = default;
155 0 : FastFlow(FastFlow&& /*rhs*/) = default;
156 0 : FastFlow& operator=(FastFlow&& /*rhs*/) = default;
157 0 : ~FastFlow() = default;
158 :
159 : // NOLINTNEXTLINE(google-runtime-references)
160 0 : void pup(PUP::er& p);
161 :
162 : /// Evaluate residuals and compute the next iteration. If
163 : /// Status==SuccessfulIteration, then `current_strahlkorper` is
164 : /// modified and `current_iteration()` is incremented. Otherwise, we
165 : /// end with success or failure, and neither `current_strahlkorper`
166 : /// nor `current_iteration()` is changed.
167 : template <typename Frame>
168 1 : std::pair<Status, IterInfo> iterate_horizon_finder(
169 : gsl::not_null<ylm::Strahlkorper<Frame>*> current_strahlkorper,
170 : const tnsr::II<DataVector, 3, Frame>& upper_spatial_metric,
171 : const tnsr::ii<DataVector, 3, Frame>& extrinsic_curvature,
172 : const tnsr::Ijj<DataVector, 3, Frame>& christoffel_2nd_kind);
173 :
174 0 : size_t current_iteration() const { return current_iter_; }
175 :
176 : /// Given a Strahlkorper defined up to some maximum Y_lm l called
177 : /// l_surface, returns a larger value of l, l_mesh, that is used for
178 : /// evaluating convergence.
179 : template <typename Frame>
180 1 : size_t current_l_mesh(const ylm::Strahlkorper<Frame>& strahlkorper) const;
181 :
182 : /// Resets the finder.
183 1 : SPECTRE_ALWAYS_INLINE void reset_for_next_find() {
184 : current_iter_ = 0;
185 : previous_residual_mesh_norm_ = 0.0;
186 : min_residual_mesh_norm_ = std::numeric_limits<double>::max();
187 : iter_at_min_residual_mesh_norm_ = 0;
188 : }
189 :
190 : private:
191 0 : friend bool operator==(const FastFlow& /*lhs*/, const FastFlow& /*rhs*/);
192 0 : double alpha_, beta_, abs_tol_, trunc_tol_, divergence_tol_;
193 0 : size_t divergence_iter_, max_its_;
194 0 : FlowType flow_;
195 0 : size_t current_iter_;
196 0 : double previous_residual_mesh_norm_, min_residual_mesh_norm_;
197 0 : size_t iter_at_min_residual_mesh_norm_;
198 : };
199 :
200 0 : SPECTRE_ALWAYS_INLINE bool converged(const FastFlow::Status& status) {
201 : return static_cast<int>(status) > 0;
202 : }
203 :
204 0 : std::ostream& operator<<(std::ostream& os, const FastFlow::FlowType& flow_type);
205 :
206 0 : std::ostream& operator<<(std::ostream& os, const FastFlow::Status& status);
207 :
208 0 : SPECTRE_ALWAYS_INLINE bool operator!=(const FastFlow& lhs,
209 : const FastFlow& rhs) {
210 : return not(lhs == rhs);
211 : }
212 :
213 : template <>
214 0 : struct Options::create_from_yaml<FastFlow::FlowType> {
215 : template <typename Metavariables>
216 0 : static FastFlow::FlowType create(const Options::Option& options) {
217 : return create<void>(options);
218 : }
219 : };
220 : template <>
221 0 : FastFlow::FlowType Options::create_from_yaml<FastFlow::FlowType>::create<void>(
222 : const Options::Option& options);
|
