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/SphericalHarmonics/Strahlkorper.hpp"
      12             : #include "Options/Context.hpp"
      13             : #include "Options/String.hpp"
      14             : #include "Parallel/GlobalCache.hpp"
      15             : #include "ParallelAlgorithms/ApparentHorizonFinder/Criteria/Criterion.hpp"
      16             : #include "ParallelAlgorithms/ApparentHorizonFinder/FastFlow.hpp"
      17             : #include "ParallelAlgorithms/ApparentHorizonFinder/Tags.hpp"
      18             : #include "Utilities/ErrorHandling/Assert.hpp"
      19             : #include "Utilities/Gsl.hpp"
      20             : #include "Utilities/TMPL.hpp"
      21             : 
      22             : namespace ah::Criteria {
      23             : /*!
      24             :  * \brief Recommend an updated resolution $L_{\rm new}$ based on the
      25             :  * resolution $L$ of the Strahlkorper and its residual $\delta$.
      26             :  *
      27             :  * \details The returned recommended resolution $L_{\rm new}$ depends on
      28             :  * the following options:
      29             :  * - MinResidual: a minimum residual $\delta_{\rm min}$
      30             :  * - MaxResidual: a maximum residual $\delta_{\rm max}$
      31             :  * - MinResolutionL: a minimum resolution $L_{\rm min}$
      32             :  *
      33             :  * The maximum resolution $L_{\rm max}$ is provided by the global cache entry
      34             :  * `ah::Tags::LMax`. The value returned for $L_{\rm new}$ is
      35             :  * then determined as follows:
      36             :  * - If $L > L_{\rm min}$ and $\delta < \delta_{\rm min}$, then
      37             :  * $L_{\rm new} = L - 1$
      38             :  * - If $L < L_{\rm max}$ and $\delta > \delta_{\rm max}$, then
      39             :  * $L_{\rm new} = L + 1$
      40             :  * - Otherwise, $L_{\rm new} = L$
      41             :  * Note that the residual is obtained from the provided FastFlow::IterInfo. The
      42             :  * residual is the gr::surfaces::expansion weighted by the FastFlow weights.
      43             :  */
      44           1 : class Residual : public Criterion {
      45             :  public:
      46           0 :   struct MinResidual {
      47           0 :     using type = double;
      48           0 :     static constexpr Options::String help = {"The minimum residual."};
      49           0 :     static type lower_bound() { return 0.0; }
      50             :   };
      51           0 :   struct MaxResidual {
      52           0 :     using type = double;
      53           0 :     static constexpr Options::String help = {"The maximum residual."};
      54           0 :     static type lower_bound() { return 0.0; }
      55             :   };
      56           0 :   struct MinResolutionL {
      57           0 :     using type = size_t;
      58           0 :     static constexpr Options::String help = {"The minimum resolution."};
      59             :     // Strahlkorper default constructor sets L=2, so don't allow L below that
      60           0 :     static type lower_bound() { return 2; }
      61             :   };
      62           0 :   using options = tmpl::list<MinResidual, MaxResidual, MinResolutionL>;
      63           0 :   static constexpr Options::String help = {
      64             :       "Use Strahlkorper residual and resolution to suggest a new "
      65             :       "resolution."};
      66             : 
      67           0 :   Residual() = default;
      68           0 :   Residual(double min_residual, double max_residual, size_t min_resolution_l,
      69             :            const Options::Context& context = {});
      70             : 
      71             :   /// \cond
      72             :   explicit Residual(CkMigrateMessage* msg);
      73             :   using PUP::able::register_constructor;
      74             :   WRAPPED_PUPable_decl_template(Residual);  // NOLINT
      75             :   /// \endcond
      76             : 
      77           0 :   std::string observation_name() override { return "Residual"; }
      78             : 
      79           0 :   using argument_tags = tmpl::list<>;
      80           0 :   using compute_tags_for_observation_box = tmpl::list<>;
      81             : 
      82             :   template <typename Metavariables, typename Frame>
      83           0 :   size_t operator()(const Parallel::GlobalCache<Metavariables>& /*cache*/,
      84             :                     const ylm::Strahlkorper<Frame>& strahlkorper,
      85             :                     const FastFlow::IterInfo& info) const;
      86             : 
      87           0 :   void pup(PUP::er& p) override;
      88             : 
      89           0 :   bool is_equal(const Criterion& other) const override;
      90             : 
      91             :  private:
      92           0 :   double min_residual_{std::numeric_limits<double>::signaling_NaN()};
      93           0 :   double max_residual_{std::numeric_limits<double>::signaling_NaN()};
      94           0 :   size_t min_resolution_l_{};
      95             : };
      96             : 
      97             : // Out-of-line definition
      98             : /// \cond
      99             : template <typename Metavariables, typename Frame>
     100             : size_t Residual::operator()(const Parallel::GlobalCache<Metavariables>& cache,
     101             :                             const ylm::Strahlkorper<Frame>& strahlkorper,
     102             :                             const FastFlow::IterInfo& info) const {
     103             :   const auto& max_resolution_l = Parallel::get<ah::Tags::LMax>(cache);
     104             :   ASSERT(min_resolution_l_ <= max_resolution_l,
     105             :          "MinResolutionL (" << min_resolution_l_ << ") must not exceed LMax ("
     106             :                             << max_resolution_l << ").");
     107             :   if (UNLIKELY(min_resolution_l_ == max_resolution_l)) {
     108             :     ASSERT(min_resolution_l_ == strahlkorper.l_max(),
     109             :            "If MinResolutionL == LMax, strahlkorper "
     110             :            "resolution must also equal MinResolutionL, but here MinResolutionL "
     111             :            "is "
     112             :                << min_resolution_l_ << ", LMax is " << max_resolution_l
     113             :                << ", and the current resolution is " << strahlkorper.l_max());
     114             :     return strahlkorper.l_max();
     115             :   }
     116             :   if (strahlkorper.l_max() > min_resolution_l_ and
     117             :       info.residual_ylm < min_residual_ and strahlkorper.l_max() > 0) {
     118             :     return strahlkorper.l_max() - 1;
     119             :   }
     120             :   if (strahlkorper.l_max() < max_resolution_l and
     121             :       info.residual_ylm > max_residual_) {
     122             :     return strahlkorper.l_max() + 1;
     123             :   }
     124             :   return strahlkorper.l_max();
     125             : }
     126             : /// \endcond
     127             : }  // namespace ah::Criteria