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 <memory>
       8             : #include <string>
       9             : 
      10             : #include "DataStructures/SpinWeighted.hpp"
      11             : #include "DataStructures/Tensor/TypeAliases.hpp"
      12             : #include "Evolution/Systems/Cce/Initialize/InitializeJ.hpp"
      13             : #include "Utilities/Gsl.hpp"
      14             : #include "Utilities/Serialization/CharmPupable.hpp"
      15             : #include "Utilities/TMPL.hpp"
      16             : 
      17             : /// \cond
      18             : class ComplexDataVector;
      19             : /// \endcond
      20             : 
      21             : namespace Cce {
      22             : namespace InitializeJ {
      23             : 
      24             : /*!
      25             :  * \brief Initialize \f$J\f$ on the first hypersurface by constraining
      26             :  * \f$\Psi_0 = 0\f$.
      27             :  *
      28             :  * \details This algorithm first radially evolves the \f$\Psi_0 = 0\f$
      29             :  * condition, which can be converted to a second-order radial ODE for J. Then,
      30             :  * the initial data generator performs an iterative solve for the angular
      31             :  * coordinates necessary to ensure asymptotic flatness. The parameters for the
      32             :  * iterative procedure are determined by options
      33             :  * `AngularCoordinateTolerance` and `MaxIterations`.
      34             :  */
      35           1 : struct NoIncomingRadiation : InitializeJ<false> {
      36           0 :   struct AngularCoordinateTolerance {
      37           0 :     using type = double;
      38           0 :     static std::string name() { return "AngularCoordTolerance"; }
      39           0 :     static constexpr Options::String help = {
      40             :         "Tolerance of initial angular coordinates for CCE"};
      41           0 :     static type lower_bound() { return 1.0e-14; }
      42           0 :     static type upper_bound() { return 1.0e-3; }
      43           0 :     static type suggested_value() { return 1.0e-10; }
      44             :   };
      45             : 
      46           0 :   struct MaxIterations {
      47           0 :     using type = size_t;
      48           0 :     static constexpr Options::String help = {
      49             :         "Number of linearized inversion iterations."};
      50           0 :     static type lower_bound() { return 10; }
      51           0 :     static type upper_bound() { return 1000; }
      52           0 :     static type suggested_value() { return 300; }
      53             :   };
      54             : 
      55           0 :   struct RequireConvergence {
      56           0 :     using type = bool;
      57           0 :     static constexpr Options::String help = {
      58             :       "If true, initialization will error if it hits MaxIterations"};
      59           0 :     static type suggested_value() { return true; }
      60             :   };
      61             : 
      62           0 :   using options =
      63             :       tmpl::list<AngularCoordinateTolerance, MaxIterations, RequireConvergence>;
      64           0 :   static constexpr Options::String help = {
      65             :       "Initialization process where J is set so Psi0 is vanishing\n"
      66             :       "(roughly a no incoming radiation condition)"};
      67             : 
      68           0 :   WRAPPED_PUPable_decl_template(NoIncomingRadiation);  // NOLINT
      69           0 :   explicit NoIncomingRadiation(CkMigrateMessage* /*unused*/) {}
      70             : 
      71           0 :   NoIncomingRadiation(double angular_coordinate_tolerance,
      72             :                       size_t max_iterations, bool require_convergence = false);
      73             : 
      74           0 :   NoIncomingRadiation() = default;
      75             : 
      76           0 :   std::unique_ptr<InitializeJ> get_clone() const override;
      77             : 
      78           0 :   void operator()(
      79             :       gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> j,
      80             :       gsl::not_null<tnsr::i<DataVector, 3>*> cartesian_cauchy_coordinates,
      81             :       gsl::not_null<
      82             :           tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>*>
      83             :           angular_cauchy_coordinates,
      84             :       const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_j,
      85             :       const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_dr_j,
      86             :       const Scalar<SpinWeighted<ComplexDataVector, 0>>& r,
      87             :       const Scalar<SpinWeighted<ComplexDataVector, 0>>& beta, size_t l_max,
      88             :       size_t number_of_radial_points,
      89             :       gsl::not_null<Parallel::NodeLock*> hdf5_lock) const override;
      90             : 
      91           0 :   void pup(PUP::er& p) override;
      92             : 
      93             :  private:
      94           0 :   bool require_convergence_ = false;
      95           0 :   double angular_coordinate_tolerance_ = 1.0e-10;
      96           0 :   size_t max_iterations_ = 300;
      97             : };
      98             : }  // namespace InitializeJ
      99             : }  // namespace Cce