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 <memory>
       9             : #include <string>
      10             : 
      11             : #include "DataStructures/SpinWeighted.hpp"
      12             : #include "DataStructures/Tensor/TypeAliases.hpp"
      13             : #include "Evolution/Systems/Cce/Initialize/InitializeJ.hpp"
      14             : #include "Utilities/Gsl.hpp"
      15             : #include "Utilities/Serialization/CharmPupable.hpp"
      16             : #include "Utilities/TMPL.hpp"
      17             : 
      18             : /// \cond
      19             : class ComplexDataVector;
      20             : /// \endcond
      21             : 
      22           1 : namespace Cce::InitializeJ {
      23             : 
      24             : /*!
      25             :  * \brief Initialize \f$J\f$ on the first hypersurface using a second-order
      26             :  * matching at the worldtube.
      27             :  *
      28             :  * \details The volume \f$J\f$ is built from the worldtube values of
      29             :  * \f$J\f$, \f$\partial_r J\f$, and \f$\partial_y^2 J\f$ computed from the
      30             :  * H hypersurface equation. The remaining angular coordinates are determined
      31             :  * iteratively to ensure asymptotic flatness. As a safeguard, the
      32             :  * initialization aborts if the second radial derivative of \f$J\f$ at scri+
      33             :  * of the final solution exceeds `MaxScriSecondDerivative`.
      34             :  */
      35           1 : struct CauchySecondOrder : 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-12; }
      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 :   struct MaxScriSecondDerivative {
      63           0 :     using type = double;
      64           0 :     static constexpr Options::String help = {
      65             :         "Abort initialization if the largest second radial derivative of J at "
      66             :         "scri+ of the final initial data exceeds this threshold. The "
      67             :         "second-order construction drives this derivative to (near) zero, so a "
      68             :         "large value indicates a poorly matched solution. Set to a large value "
      69             :         "to effectively disable the check."};
      70           0 :     static type lower_bound() { return 1.0e-14; }
      71           0 :     static type upper_bound() { return 1.0e2; }
      72           0 :     static type suggested_value() { return 1.0e-8; }
      73             :   };
      74             : 
      75           0 :   using options = tmpl::list<AngularCoordinateTolerance, MaxIterations,
      76             :                              RequireConvergence, MaxScriSecondDerivative>;
      77           0 :   static constexpr Options::String help = {
      78             :       "Second-order initial data generator for the Cauchy CCE evolution."};
      79             : 
      80           0 :   WRAPPED_PUPable_decl_template(CauchySecondOrder);  // NOLINT
      81           0 :   explicit CauchySecondOrder(CkMigrateMessage* /*unused*/) {}
      82             : 
      83           0 :   CauchySecondOrder(double angular_coordinate_tolerance, size_t max_iterations,
      84             :                     bool require_convergence,
      85             :                     double max_scri_second_derivative);
      86             : 
      87           0 :   CauchySecondOrder() = default;
      88             : 
      89           0 :   std::unique_ptr<InitializeJ> get_clone() const override;
      90             : 
      91             :   // Per-class tag lists. The flexible dispatch in `InitializeJ<false>` reads
      92             :   // these via `call_with_dynamic_type` so this generator can request more
      93             :   // worldtube boundary values than the simpler sibling classes do.
      94           0 :   using return_tags = tmpl::list<Tags::BondiJ, Tags::CauchyCartesianCoords,
      95             :                                  Tags::CauchyAngularCoords>;
      96           0 :   using argument_tags = tmpl::list<
      97             :       Tags::BoundaryValue<Tags::BondiJ>, Tags::BoundaryValue<Tags::BondiU>,
      98             :       Tags::BoundaryValue<Tags::BondiW>, Tags::BoundaryValue<Tags::BondiBeta>,
      99             :       Tags::BoundaryValue<Tags::BondiQ>,
     100             :       Tags::BoundaryValue<Tags::Du<Tags::BondiJ>>,
     101             :       Tags::BoundaryValue<Tags::Dr<Tags::BondiJ>>,
     102             :       Tags::BoundaryValue<Tags::Du<Tags::Dr<Tags::BondiJ>>>,
     103             :       Tags::BoundaryValue<Tags::Du<Tags::BondiR>>,
     104             :       Tags::BoundaryValue<Tags::BondiR>, Tags::LMax,
     105             :       Tags::NumberOfRadialPoints>;
     106             : 
     107           0 :   void operator()(
     108             :       gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> j,
     109             :       gsl::not_null<tnsr::i<DataVector, 3>*> cartesian_cauchy_coordinates,
     110             :       gsl::not_null<
     111             :           tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>*>
     112             :           angular_cauchy_coordinates,
     113             :       const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_j,
     114             :       const Scalar<SpinWeighted<ComplexDataVector, 1>>& boundary_u,
     115             :       const Scalar<SpinWeighted<ComplexDataVector, 0>>& boundary_w,
     116             :       const Scalar<SpinWeighted<ComplexDataVector, 0>>& boundary_beta,
     117             :       const Scalar<SpinWeighted<ComplexDataVector, 1>>& boundary_q,
     118             :       const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_du_j,
     119             :       const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_dr_j,
     120             :       const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_du_dr_j,
     121             :       const Scalar<SpinWeighted<ComplexDataVector, 0>>& boundary_du_r,
     122             :       const Scalar<SpinWeighted<ComplexDataVector, 0>>& r, size_t l_max,
     123             :       size_t number_of_radial_points,
     124             :       gsl::not_null<Parallel::NodeLock*> hdf5_lock) const;
     125             : 
     126           0 :   void pup(PUP::er& p) override;
     127             : 
     128             :  private:
     129           0 :   bool require_convergence_ = true;
     130           0 :   double angular_coordinate_tolerance_ =
     131             :       std::numeric_limits<double>::signaling_NaN();
     132           0 :   size_t max_iterations_ = 0;
     133           0 :   double max_scri_second_derivative_ =
     134             :       std::numeric_limits<double>::signaling_NaN();
     135             : };
     136             : }  // namespace Cce::InitializeJ