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