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 <pup.h>
9 : #include <string>
10 :
11 : #include "NumericalAlgorithms/LinearOperators/PowerMonitors.hpp"
12 : #include "NumericalAlgorithms/Strahlkorper/Strahlkorper.hpp"
13 : #include "NumericalAlgorithms/Strahlkorper/StrahlkorperFunctions.hpp"
14 : #include "Options/Context.hpp"
15 : #include "Options/String.hpp"
16 : #include "Parallel/GlobalCache.hpp"
17 : #include "ParallelAlgorithms/ApparentHorizonFinder/Criteria/Criterion.hpp"
18 : #include "ParallelAlgorithms/ApparentHorizonFinder/FastFlow.hpp"
19 : #include "ParallelAlgorithms/ApparentHorizonFinder/Tags.hpp"
20 : #include "Utilities/ErrorHandling/Assert.hpp"
21 : #include "Utilities/Gsl.hpp"
22 : #include "Utilities/TMPL.hpp"
23 :
24 : namespace ah::Criteria {
25 : /*!
26 : * \brief Recommend an updated resolution $L_{\rm new}$ based on how well
27 : * the shape of the apparent horizon is resolved in terms of its
28 : * truncation error and pile up modes.
29 : *
30 : * \details The returned recommended resolution $L_{\rm new}$ depends on
31 : * the following options:
32 : * - MinTruncationError: a minimum truncation error $\mathcal{T}_{\rm min}$
33 : * - MaxTruncationError: a maximum truncation error $\mathcal{T}_{\rm max}$
34 : * - MaxPileUpModes: a maximum number of pile up modes $\mathcal{P}_{\rm max}$
35 : * - MinResolutionL: a minimum resolution $L_{\rm min}$
36 : *
37 : * The maximum resolution $L_{\rm max}$ is supplied by the global cache entry
38 : * `ah::Tags::LMax`. The value returned for $L_{\rm new}$ is
39 : * then determined as follows:
40 : * The value returned for $L_{\rm new}$ is then determined as follows:
41 : * - If $L > L_{\rm min}$ and $\mathcal{T} < \mathcal{T}_{\rm min}$, then
42 : * $L_{\rm new} = L - 1$
43 : * - If $L < L_{\rm max}$ and $\mathcal{T} > \mathcal{T}_{\rm max}$ and
44 : * $\mathcal{P} \le \mathcal{P}_{\rm max}$, then
45 : * $L_{\rm new} = L + 1$
46 : * - Otherwise, $L_{\rm new} = L$
47 : * Note that the truncation error is obtained via
48 : * PowerMonitors::relative_truncation_error, and the number of pileup modes
49 : * is given by PowerMonitors::convergence_rate_and_number_of_pile_up_modes
50 : */
51 1 : class Shape : public Criterion {
52 : public:
53 0 : struct MinTruncationError {
54 0 : using type = double;
55 0 : static constexpr Options::String help = {"The minimum truncation error."};
56 0 : static type lower_bound() { return 0.0; }
57 : };
58 0 : struct MaxTruncationError {
59 0 : using type = double;
60 0 : static constexpr Options::String help = {"The maximum truncation error."};
61 0 : static type lower_bound() { return 0.0; }
62 : };
63 0 : struct MaxPileUpModes {
64 0 : using type = size_t;
65 0 : static constexpr Options::String help = {
66 : "The maximum number of pile up modes."};
67 0 : static type lower_bound() { return 0; }
68 : };
69 0 : struct MinResolutionL {
70 0 : using type = size_t;
71 0 : static constexpr Options::String help = {"The minimum resolution."};
72 : // Power montior requires at least L=4, so don't allow L below that
73 0 : static type lower_bound() { return 4; }
74 : };
75 0 : using options = tmpl::list<MinTruncationError, MaxTruncationError,
76 : MaxPileUpModes, MinResolutionL>;
77 0 : static constexpr Options::String help = {
78 : "Use Strahlkorper truncation error, number of pile up modes, and "
79 : "resolution to suggest a new resolution."};
80 :
81 0 : Shape() = default;
82 0 : Shape(double min_truncation_error, double max_truncation_error,
83 : size_t max_pile_up_modes, size_t min_resolution_l,
84 : const Options::Context& context = {});
85 :
86 : /// \cond
87 : explicit Shape(CkMigrateMessage* msg);
88 : using PUP::able::register_constructor;
89 : WRAPPED_PUPable_decl_template(Shape); // NOLINT
90 : /// \endcond
91 :
92 0 : std::string observation_name() override { return "Shape"; }
93 :
94 0 : using argument_tags = tmpl::list<>;
95 0 : using compute_tags_for_observation_box = tmpl::list<>;
96 :
97 : template <typename Metavariables, typename Frame>
98 0 : size_t operator()(const Parallel::GlobalCache<Metavariables>& cache,
99 : const ylm::Strahlkorper<Frame>& strahlkorper,
100 : const FastFlow::IterInfo& /*info*/) const;
101 :
102 0 : void pup(PUP::er& p) override;
103 :
104 0 : bool is_equal(const Criterion& other) const override;
105 :
106 : private:
107 0 : double min_truncation_error_{std::numeric_limits<double>::signaling_NaN()};
108 0 : double max_truncation_error_{std::numeric_limits<double>::signaling_NaN()};
109 0 : size_t max_pile_up_modes_{};
110 0 : size_t min_resolution_l_{};
111 : };
112 :
113 : // Out-of-line definition
114 : /// \cond
115 : template <typename Metavariables, typename Frame>
116 : size_t Shape::operator()(const Parallel::GlobalCache<Metavariables>& cache,
117 : const ylm::Strahlkorper<Frame>& strahlkorper,
118 : const FastFlow::IterInfo& /*info*/) const {
119 : const auto& max_resolution_l = Parallel::get<ah::Tags::LMax>(cache);
120 : ASSERT(min_resolution_l_ <= max_resolution_l,
121 : "MinResolutionL (" << min_resolution_l_ << ") must not exceed LMax ("
122 : << max_resolution_l << ").");
123 : if (UNLIKELY(min_resolution_l_ == max_resolution_l)) {
124 : ASSERT(min_resolution_l_ == strahlkorper.l_max(),
125 : "If MinResolutionL == LMax, strahlkorper "
126 : "resolution must also equal MinResolutionL, but here MinResolutionL "
127 : "is "
128 : << min_resolution_l_ << ", LMax is " << max_resolution_l
129 : << ", and the current resolution is " << strahlkorper.l_max());
130 : return strahlkorper.l_max();
131 : }
132 :
133 : const DataVector& power_monitor = ylm::power_monitor(strahlkorper);
134 : const double truncation_error = PowerMonitors::relative_truncation_error(
135 : power_monitor, power_monitor.size());
136 : const size_t pile_up_modes =
137 : PowerMonitors::convergence_rate_and_number_of_pile_up_modes(power_monitor)
138 : .number_of_pile_up_modes;
139 : if (strahlkorper.l_max() > min_resolution_l_ and
140 : truncation_error < min_truncation_error_ and strahlkorper.l_max() > 0) {
141 : return strahlkorper.l_max() - 1;
142 : }
143 : if (strahlkorper.l_max() < max_resolution_l and
144 : truncation_error > max_truncation_error_ and
145 : pile_up_modes <= max_pile_up_modes_) {
146 : return strahlkorper.l_max() + 1;
147 : }
148 : return strahlkorper.l_max();
149 : }
150 : /// \endcond
151 : } // namespace ah::Criteria
|
