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);