Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <limits>
       7             : #include <ostream>
       8             : 
       9             : #include "NumericalAlgorithms/Interpolation/CubicSpline.hpp"
      10             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      11             : 
      12             : /// \cond
      13             : namespace Options {
      14             : struct Option;
      15             : template <typename T>
      16             : struct create_from_yaml;
      17             : }  // namespace Options
      18             : namespace PUP {
      19             : class er;
      20             : }  // namespace PUP
      21             : /// \endcond
      22             : 
      23             : namespace RelativisticEuler::Solutions {
      24             : 
      25             : /// Radial coordinate of a TOV solution.
      26           1 : enum class TovCoordinates {
      27             :   /// Areal radial coordinate, in which the exterior TOV solution coincides with
      28             :   /// the Schwarzschild metric in standard Schwarzschild coordinates.
      29             :   Schwarzschild,
      30             :   /// Isotropic radial coordinate, in which the spatial metric is conformally
      31             :   /// flat.
      32             :   Isotropic
      33             : };
      34             : 
      35           0 : std::ostream& operator<<(std::ostream& os, TovCoordinates coords);
      36             : 
      37             : }  // namespace RelativisticEuler::Solutions
      38             : 
      39             : template <>
      40           0 : struct Options::create_from_yaml<RelativisticEuler::Solutions::TovCoordinates> {
      41             :   template <typename Metavariables>
      42           0 :   static RelativisticEuler::Solutions::TovCoordinates create(
      43             :       const Options::Option& options) {
      44             :     return create<void>(options);
      45             :   }
      46             : };
      47             : 
      48             : template <>
      49           0 : RelativisticEuler::Solutions::TovCoordinates
      50             : Options::create_from_yaml<RelativisticEuler::Solutions::TovCoordinates>::create<
      51             :     void>(const Options::Option& options);
      52             : 
      53             : namespace RelativisticEuler::Solutions {
      54             : 
      55             : /*!
      56             :  * \brief TOV solver based on Lindblom's method
      57             :  *
      58             :  * Uses Lindblom's method of integrating the TOV equations from
      59             :  * \cite Lindblom1998dp .
      60             :  *
      61             :  * Instead of integrating the interior mass \f$m(r)\f$ and pressure \f$p(r)\f$,
      62             :  * Lindblom introduces the variables \f$u=r^2\f$ and \f$v=m/r\f$. Then, the TOV
      63             :  * equations are integrated with the log of the specific enthalpy as the
      64             :  * independent variable, \f$\ln(h)\f$, from the center of the star where
      65             :  * \f$h(r=0) = h_c\f$ to its surface \f$h(r=R) = 1\f$. The ODEs being solved are
      66             :  * Eq. (A2) and (A3) in \cite Lindblom1998dp :
      67             :  *
      68             :  * \f{align}
      69             :  * \frac{\mathrm{d}u}{\mathrm{d}\ln{h}} &= \frac{-2u (1 - 2v)}{4\pi u p + v} \\
      70             :  * \frac{\mathrm{d}v}{\mathrm{d}\ln{h}} &=
      71             :  *   -(1 - 2v) \frac{4\pi u \rho - v}{4\pi u p + v}
      72             :  * \f}
      73             :  *
      74             :  * Note that Lindblom's paper labels the independent variable as \f$h\f$.
      75             :  * However the \f$h\f$ in Lindblom's paper is **not** the specific enthalpy but
      76             :  * its logarithm, \f$\ln(h)\f$.
      77             :  *
      78             :  * The ODEs are solved numerically when this class is constructed, and the
      79             :  * quantities \f$m(r)/r\f$ and \f$\ln(h)\f$ are interpolated and exposed as
      80             :  * member functions. With these quantities the metric can be constructed as:
      81             :  *
      82             :  * \f{equation}
      83             :  * \mathrm{d}s^2 = -\alpha^2 \mathrm{d}t^2 + (1 - 2m/r)^{-1} \mathrm{d}r^2
      84             :  *   + r^2 \mathrm{d}\Omega^2
      85             :  * \f}
      86             :  *
      87             :  * where the lapse is
      88             :  *
      89             :  * \f{equation}
      90             :  * \alpha(r < R) = \mathcal{E} / h = \alpha(r=R) / h
      91             :  * \text{,}
      92             :  * \f}
      93             :  *
      94             :  * with the conserved `injection_energy()` \f$\mathcal{E}\f$, such that the
      95             :  * lapse matches the exterior Schwarzschild solution:
      96             :  *
      97             :  * \f{equation}
      98             :  * \alpha(r \geq R) = \sqrt{1 - \frac{2M}{r}}
      99             :  * \f}
     100             :  *
     101             :  * \par Isotropic radial coordinate
     102             :  * This class also supports transforming to an isotropic radial coordinate. When
     103             :  * you pass `RelativisticEuler::Solutions::TovCoordinates::Isotropic` to the
     104             :  * constructor, an additional ODE is integrated alongside the TOV equations to
     105             :  * determine the conformal factor
     106             :  *
     107             :  * \f{equation}
     108             :  * \psi^2 = \frac{r}{\bar{r}}
     109             :  * \f}
     110             :  *
     111             :  * where \f$r\f$ is the areal (Schwarzschild) radius and \f$\bar{r}\f$ is the
     112             :  * isotropic radius. The additional ODE is:
     113             :  *
     114             :  * \f{equation}
     115             :  * \frac{\mathrm{d}\ln(\psi)}{\mathrm{d}\ln{h}} =
     116             :  *   \frac{\sqrt{1 - 2v}}{1 + \sqrt{1 - 2v}} \frac{v}{4\pi u p + v}
     117             :  * \f}
     118             :  *
     119             :  * In isotropic coordinates, the spatial metric is conformally flat:
     120             :  *
     121             :  * \f{equation}
     122             :  * \mathrm{d}s^2 = -\alpha^2 \mathrm{d}t^2 + \psi^4 (\mathrm{d}\bar{r}^2 +
     123             :  *   \bar{r}^2 \mathrm{d}\Omega^2)
     124             :  * \f}
     125             :  *
     126             :  * When isotropic coordinates are selected, radii returned by member functions
     127             :  * or taken as arguments are isotropic. An exception is `mass_over_radius()`,
     128             :  * which always returns the quantity \f$m / r\f$ because that is the quantity
     129             :  * stored internally and hence most numerically precise. See
     130             :  * `mass_over_radius()` for details.
     131             :  */
     132           1 : class TovSolution {
     133             :  public:
     134           0 :   TovSolution(
     135             :       const EquationsOfState::EquationOfState<true, 1>& equation_of_state,
     136             :       double central_mass_density,
     137             :       const TovCoordinates coordinate_system = TovCoordinates::Schwarzschild,
     138             :       double log_enthalpy_at_outer_radius = 0.0,
     139             :       double absolute_tolerance = 1.e-18, double relative_tolerance = 1.0e-14);
     140             : 
     141           0 :   TovSolution() = default;
     142           0 :   TovSolution(const TovSolution& /*rhs*/) = default;
     143           0 :   TovSolution& operator=(const TovSolution& /*rhs*/) = default;
     144           0 :   TovSolution(TovSolution&& /*rhs*/) = default;
     145           0 :   TovSolution& operator=(TovSolution&& /*rhs*/) = default;
     146           0 :   ~TovSolution() = default;
     147             : 
     148             :   /// The type of radial coordinate.
     149             :   ///
     150             :   /// \see RelativisticEuler::Solutions::TovCoordinates
     151           1 :   TovCoordinates coordinate_system() const { return coordinate_system_; }
     152             : 
     153             :   /// \brief The outer radius of the solution.
     154             :   ///
     155             :   /// This is the outer radius in the specified `coordinate_system()`, i.e.,
     156             :   /// areal or isotropic.
     157             :   ///
     158             :   /// \note This is the radius at which `log_specific_enthalpy` is equal
     159             :   /// to the value of `log_enthalpy_at_outer_radius` that was given when
     160             :   /// constructing this TovSolution
     161           1 :   double outer_radius() const { return outer_radius_; }
     162             : 
     163             :   /// The total mass \f$m(R)\f$, where \f$R\f$ is the outer radius.
     164           1 :   double total_mass() const { return total_mass_; }
     165             : 
     166             :   /*!
     167             :    * \brief The injection energy \f$\mathcal{E}=\alpha(r=R)=\sqrt{1-2M/R}\f$.
     168             :    *
     169             :    * The injection energy of the TOV solution is
     170             :    *
     171             :    * \f{equation}
     172             :    * \mathcal{E} = -h k^a u_a = h \alpha
     173             :    * \text{,}
     174             :    * \f}
     175             :    *
     176             :    * where \f$\boldsymbol{k} = \partial_t\f$ is a Killing vector of the static
     177             :    * solution, \f$h\f$ is the specific enthalpy, \f$u_a\f$ is the fluid
     178             :    * four-velocity, and \f$\alpha\f$ is the lapse (see, e.g., Eqs. (2.19) and
     179             :    * (4.2) in \cite Moldenhauer2014yaa). Since the TOV solution is static, the
     180             :    * injection energy is conserved not only along stream lines but throughout
     181             :    * the star,
     182             :    *
     183             :    * \f{equation}
     184             :    * \nabla_a \mathcal{E} = 0
     185             :    * \text{.}
     186             :    * \f}
     187             :    *
     188             :    * Therefore,
     189             :    *
     190             :    * \f{equation}
     191             :    * \mathcal{E} = \alpha(r=R) = \sqrt{1 - 2M/R}
     192             :    * \f}
     193             :    *
     194             :    * by evaluating the injection energy at the outer (areal) radius \f$R\f$,
     195             :    * where \f$h=1\f$ and where we match the lapse to the outer Schwarzschild
     196             :    * solution. The conservation also implies
     197             :    *
     198             :    * \f{equation}
     199             :    * \alpha = \mathcal{E} / h
     200             :    * \f}
     201             :    *
     202             :    * throughout the star.
     203             :    */
     204           1 :   double injection_energy() const { return injection_energy_; }
     205             : 
     206             :   /// \brief The mass inside the given radius over the areal radius,
     207             :   /// \f$\frac{m(r)}{r}\f$
     208             :   ///
     209             :   /// The argument to this function is the radius in the `coordinate_system()`,
     210             :   /// i.e., areal (Schwarzschild) or isotropic radius. The denominator \f$r\f$
     211             :   /// in the return value is always the areal (Schwarzschild) radius. You can
     212             :   /// use the conformal factor \f$\psi=\sqrt{r / \bar{r}}\f$ returned by the
     213             :   /// `conformal_factor()` function to obtain the mass over the isotropic
     214             :   /// radius, or the mass alone. The reason for this choice is that we represent
     215             :   /// the solution internally as the mass over the areal radius, so this is the
     216             :   /// most numerically precise quantity from which other quantities can be
     217             :   /// derived.
     218             :   ///
     219             :   /// \note `r` should be non-negative and not greater than `outer_radius()`.
     220             :   template <typename DataType>
     221           1 :   DataType mass_over_radius(const DataType& r) const;
     222             : 
     223             :   /// \brief The log of the specific enthalpy at the given radius
     224             :   ///
     225             :   /// \note `r` should be non-negative and not greater than `outer_radius()`
     226             :   template <typename DataType>
     227           1 :   DataType log_specific_enthalpy(const DataType& r) const;
     228             : 
     229             :   /// \brief The conformal factor \f$\psi=\sqrt{r / \bar{r}}\f$.
     230             :   ///
     231             :   /// The conformal factor is computed only when the `coordinate_system()` is
     232             :   /// `RelativisticEuler::Solution::TovCoordinates::Isotropic`. Otherwise, it is
     233             :   /// an error to call this function.
     234             :   ///
     235             :   /// \note `r` should be non-negative and not greater than `outer_radius()`
     236             :   template <typename DataType>
     237           1 :   DataType conformal_factor(const DataType& r) const;
     238             : 
     239           0 :   const intrp::CubicSpline& mass_over_radius_interpolant() const {
     240             :     return mass_over_radius_interpolant_;
     241             :   }
     242           0 :   const intrp::CubicSpline& log_specific_enthalpy_interpolant() const {
     243             :     return log_enthalpy_interpolant_;
     244             :   }
     245           0 :   const intrp::CubicSpline& conformal_factor_interpolant() const {
     246             :     return conformal_factor_interpolant_;
     247             :   }
     248             : 
     249             :   // NOLINTNEXTLINE(google-runtime-references)
     250           0 :   void pup(PUP::er& p);
     251             : 
     252             :  private:
     253             :   template <TovCoordinates CoordSystem>
     254           0 :   void integrate(
     255             :       const EquationsOfState::EquationOfState<true, 1>& equation_of_state,
     256             :       const double central_mass_density,
     257             :       const double log_enthalpy_at_outer_radius,
     258             :       const double absolute_tolerance, const double relative_tolerance);
     259             : 
     260           0 :   TovCoordinates coordinate_system_{};
     261           0 :   double outer_radius_{std::numeric_limits<double>::signaling_NaN()};
     262           0 :   double total_mass_{std::numeric_limits<double>::signaling_NaN()};
     263           0 :   double injection_energy_{std::numeric_limits<double>::signaling_NaN()};
     264           0 :   intrp::CubicSpline mass_over_radius_interpolant_;
     265           0 :   intrp::CubicSpline log_enthalpy_interpolant_;
     266           0 :   intrp::CubicSpline conformal_factor_interpolant_;
     267             : };
     268             : 
     269             : }  // namespace RelativisticEuler::Solutions