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/Strahlkorper/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