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/SphericalHarmonics/Strahlkorper.hpp"
      13             : #include "NumericalAlgorithms/SphericalHarmonics/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