SpECTRE Documentation Coverage Report
Current view: top level - PointwiseFunctions/AnalyticSolutions/RelativisticEuler - RotatingStar.hpp Coverage Total Hit
Commit: 0d1d717af84a81fb1284a20c4013708fe85cc2f0 Lines: 4.3 % 94 4
Test Date: 2026-10-22 17:29:42
Legend: Lines:     hit not hit

            Line data    Source code
       1            0 : // Distributed under the MIT License.
       2              : // See LICENSE.txt for details.
       3              : 
       4              : #pragma once
       5              : 
       6              : #include <array>
       7              : #include <cstddef>
       8              : #include <optional>
       9              : #include <string>
      10              : 
      11              : #include "DataStructures/TaggedTuple.hpp"
      12              : #include "DataStructures/Tensor/TypeAliases.hpp"
      13              : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
      14              : #include "Options/Options.hpp"
      15              : #include "Options/String.hpp"
      16              : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
      17              : #include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Solutions.hpp"
      18              : #include "PointwiseFunctions/Hydro/EquationsOfState/Factory.hpp"
      19              : #include "PointwiseFunctions/Hydro/Tags.hpp"
      20              : #include "PointwiseFunctions/Hydro/Temperature.hpp"
      21              : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
      22              : #include "Utilities/TMPL.hpp"
      23              : 
      24              : namespace RelativisticEuler::Solutions {
      25              : namespace detail {
      26              : /*!
      27              :  * \brief Read the CST rotating neutron star solution from file, rescaling the
      28              :  * solution.  If the EoS in RotNS is set to PT (polytrope), then the user
      29              :  * should use `is_polytrope=true`, so that the normalization is done using the
      30              :  * method in \cite Cook1992. In this case, RotNS already uses geometric units,
      31              :  * letting \f$K\f$ (i.e.  \f$P=K\rho^\gamma\f$) be a free parameter. RotNS
      32              :  * assumes \f$K=1\f$ for its calculations. If the user doesn't want to rescale,
      33              :  * they can just set `polytropic_constant=1.` and nothing will change.
      34              :  *
      35              :  * If any other EoS is used, then it's done according to the method in
      36              :  * \cite Cook1994. In this case, RotNS instead normalizes with CGS, so
      37              :  * `is_polytrope=false` converts from CGS into geometric units. The user would
      38              :  * not want to use this for a PT EoS, even if they don't want to rescale.
      39              :  */
      40              : class CstSolution {
      41              :  public:
      42              :   CstSolution() = default;
      43              :   CstSolution(const std::string& filename, bool is_polytrope,
      44              :               double polytropic_constant);
      45              : 
      46              :   std::array<double, 6> interpolate(double target_radius,
      47              :                                     double target_cos_theta,
      48              :                                     bool interpolate_hydro_vars) const;
      49              : 
      50              :   double polytropic_index() const { return polytropic_index_; }
      51              : 
      52              :   double equatorial_radius() const { return equatorial_radius_; }
      53              : 
      54              :   void pup(PUP::er& p);
      55              : 
      56              :  private:
      57              :   double maximum_radius_{std::numeric_limits<double>::signaling_NaN()};
      58              :   double max_density_ratio_for_linear_interpolation_ = 1.0e2;
      59              : 
      60              :   double equatorial_radius_{std::numeric_limits<double>::signaling_NaN()};
      61              :   double polytropic_index_{std::numeric_limits<double>::signaling_NaN()};
      62              :   double central_angular_speed_{std::numeric_limits<double>::signaling_NaN()};
      63              :   double rotation_profile_{std::numeric_limits<double>::signaling_NaN()};
      64              :   size_t num_radial_points_{std::numeric_limits<size_t>::max()};
      65              :   size_t num_angular_points_{std::numeric_limits<size_t>::max()};
      66              :   size_t num_grid_points_{std::numeric_limits<size_t>::max()};
      67              : 
      68              :   DataVector radius_;
      69              :   DataVector cos_theta_;  // Note that cos(theta) is between 0 and 1.
      70              :   DataVector rest_mass_density_;
      71              :   DataVector fluid_velocity_;
      72              :   DataVector alpha_;
      73              :   DataVector rho_;
      74              :   DataVector gamma_;
      75              :   DataVector omega_;
      76              : };
      77              : }  // namespace detail
      78              : 
      79              : /*!
      80              :  * \brief A solution obtained by reading in rotating neutron star initial data
      81              :  * from the RotNS code based on \cite Cook1992 and \cite Cook1994.
      82              :  *
      83              :  * The code that generates the initial data is part of a private SXS repository
      84              :  * called `RotNS`.
      85              :  *
      86              :  * The metric in spherical coordinates is given by \cite Cook1992
      87              :  *
      88              :  * \f{align}{
      89              :  *   ds^2=-e^{\gamma+\rho}dt^2+e^{2\alpha}(dr^2+r^2d\theta^2)
      90              :  *   +e^{\gamma-\rho}r^2\sin^2(\theta)(d\phi-\omega dt)^2.
      91              :  * \f}
      92              :  *
      93              :  * We use rotation about the \f$z\f$-axis. That is,
      94              :  *
      95              :  * \f{align}{
      96              :  *   g_{tt}
      97              :  *   &=-e^{\gamma+\rho} + e^{\gamma-\rho}r^2\sin^2(\theta)\omega^2 \\
      98              :  *   g_{rr}
      99              :  *   &=e^{2\alpha} \\
     100              :  *   g_{\theta\theta}
     101              :  *   &=e^{2\alpha}r^2 \\
     102              :  *   g_{\phi\phi}
     103              :  *   &=e^{\gamma-\rho}r^2\sin^2(\theta) \\
     104              :  *   g_{t\phi}
     105              :  *   &=-e^{\gamma-\rho}r^2\sin^2(\theta)\omega.
     106              :  * \f}
     107              :  *
     108              :  * We can transform from spherical to Cartesian coordinates using
     109              :  *
     110              :  * \f{align}{
     111              :  *   \label{eq:Jacobian}
     112              :  *   \frac{\partial (r,\theta,\phi)}{\partial (x,y,z)}=
     113              :  *   \begin{pmatrix}
     114              :  *     \cos(\phi) \sin(\theta) & \sin(\theta)\sin(\phi) & \cos(\theta) \\
     115              :  *     \tfrac{\cos(\phi)\cos(\theta)}{r} & \tfrac{\sin(\phi)\cos(\theta)}{r}
     116              :  *     & -\tfrac{\sin(\theta)}{r} \\
     117              :  *     -\tfrac{\sin(\phi)}{r\sin(\theta)} & \tfrac{\cos(\phi)}{r\sin(\theta)}
     118              :  *     & 0
     119              :  *   \end{pmatrix}
     120              :  * \f}
     121              :  *
     122              :  * and
     123              :  *
     124              :  * \f{align}{
     125              :  *   \frac{\partial (x,y,z)}{\partial (r,\theta,\phi)}=
     126              :  *   \begin{pmatrix}
     127              :  *     \cos(\phi) \sin(\theta) & r\cos(\phi)\cos(\theta) &
     128              :  * -r\sin(\theta)\sin(\phi)
     129              :  *     \\
     130              :  *     \sin(\theta)\sin(\phi) & r\sin(\phi)\cos(\theta) &
     131              :  * r\cos(\phi)\sin(\theta)
     132              :  *     \\
     133              :  *     \cos(\theta) & -r\sin(\theta) & 0
     134              :  *   \end{pmatrix}
     135              :  * \f}
     136              :  *
     137              :  *
     138              :  * We denote the lapse as \f$N\f$ since \f$\alpha\f$ is already being used,
     139              :  *
     140              :  * \f{align}{
     141              :  *   N = e^{(\gamma+\rho)/2}.
     142              :  * \f}
     143              :  *
     144              :  * The shift is
     145              :  *
     146              :  * \f{align}{
     147              :  *   \beta^\phi =& -\omega
     148              :  * \f}
     149              :  *
     150              :  * so
     151              :  *
     152              :  * \f{align}{
     153              :  *   \beta^x &= \partial_\phi x \beta^\phi = \sin(\phi)\omega r \sin(\theta), \\
     154              :  *   \beta^y &= \partial_\phi y \beta^\phi = -\cos(\phi)\omega r \sin(\theta),
     155              :  * \\ \beta^z &= 0. \f}
     156              :  *
     157              :  * The spatial metric is
     158              :  *
     159              :  * \f{align}{
     160              :  *   \gamma_{xx}
     161              :  *   &= \sin^2(\theta)\cos^2(\phi) e^{(2 \alpha)}
     162              :  *      + \cos^2(\theta)\cos^2(\phi) e^{(2 \alpha)}
     163              :  *      + \sin^2(\phi) e^{(\gamma-\rho)} \notag \\
     164              :  *   &=\cos^2(\phi)e^{2\alpha} + \sin^2(\phi) e^{(\gamma-\rho)} \\
     165              :  *   \gamma_{xy}
     166              :  *   &= \sin^2(\theta)\cos(\phi)\sin(\phi) e^{(2 \alpha)}
     167              :  *      + \cos^2(\theta)\cos(\phi)\sin(\phi) e^{(2 \alpha)}
     168              :  *      - \sin(\phi)\cos(\phi) e^{(\gamma-\rho)} \notag \\
     169              :  *   &=\cos(\phi)\sin(\phi)e^{2\alpha} - \sin(\phi)\cos(\phi) e^{(\gamma-\rho)}
     170              :  * \\ \gamma_{yy}
     171              :  *   &= \sin^2(\theta)\sin^2(\phi) e^{(2 \alpha)}
     172              :  *      + \cos^2(\theta)\sin^2(\phi) e^{(2 \alpha)}
     173              :  *      + \cos^2(\phi) e^{(\gamma-\rho)} \notag \\
     174              :  *   &=\sin^2(\phi)e^{2\alpha} + \cos^2(\phi) e^{(\gamma-\rho)} \\
     175              :  *   \gamma_{xz}
     176              :  *   &= \sin(\theta)\cos(\phi)\cos(\theta) e^{2\alpha}
     177              :  *      -\cos(\theta)\cos(\phi)\sin(\theta)e^{2\alpha} = 0 \\
     178              :  *   \gamma_{yz}
     179              :  *   &= \sin(\theta)\sin(\phi)\cos(\theta) e^{2\alpha}
     180              :  *      -\cos(\theta)\sin(\phi)\sin(\theta)e^{2\alpha} = 0 \\
     181              :  *   \gamma_{zz}
     182              :  *   &=\cos^2(\theta)e^{2\alpha} + \sin^2(\theta) e^{2\alpha}
     183              :  *      = e^{2\alpha}
     184              :  * \f}
     185              :  *
     186              :  * and its determinant is
     187              :  *
     188              :  * \f{align}{
     189              :  *   \gamma = e^{4\alpha + (\gamma-\rho)} = e^{4\alpha}e^{(\gamma-\rho)}.
     190              :  * \f}
     191              :  *
     192              :  * At \f$r=0\f$ we have \f$2\alpha=\gamma-\rho\f$ and so the
     193              :  * \f$\gamma_{xx}=\gamma_{yy}=\gamma_{zz} = e^{2\alpha}\f$ and all other
     194              :  * components are zero. The inverse spatial metric is given by
     195              :  *
     196              :  * \f{align}{
     197              :  *   \gamma^{xx}
     198              :  *   &= \frac{\gamma_{yy}}{e^{2\alpha}e^{(\gamma-\rho)}} =
     199              :  *      \left[\sin^2(\phi)e^{2\alpha} + \cos^2(\phi) e^{(\gamma-\rho)}\right]
     200              :  *      e^{-2\alpha} e^{-(\gamma-\rho)} \notag \\
     201              :  *   &=\sin^2(\phi)e^{-(\gamma-\rho)} + \cos^2(\phi) e^{-2\alpha} \\
     202              :  *   \gamma^{yy}
     203              :  *   &= \frac{\gamma_{xx}}{e^{2\alpha}e^{(\gamma-\rho)}} =
     204              :  *      \left[\cos^2(\phi)e^{2\alpha} + \sin^2(\phi) e^{(\gamma-\rho)}\right]
     205              :  *      e^{-2\alpha} e^{-(\gamma-\rho)} \notag \\
     206              :  *   &=\cos^2(\phi) e^{-(\gamma-\rho)} + \sin^2(\phi) e^{-2\alpha} \\
     207              :  *   \gamma^{xy}
     208              :  *   &=\frac{-\gamma_{xy}}{e^{2\alpha}e^{(\gamma-\rho)}} =
     209              :  *     -\left[\cos(\phi)\sin(\phi)e^{2\alpha} - \sin(\phi)\cos(\phi)
     210              :  *      e^{(\gamma-\rho)}\right] e^{-2\alpha} e^{-(\gamma-\rho)} \notag \\
     211              :  *   &=-\cos(\phi)\sin(\phi)e^{-(\gamma-\rho)} -
     212              :  *     \sin(\phi)\cos(\phi)e^{-2\alpha} \notag \\
     213              :  *   &=\cos(\phi)\sin(\phi) \left[e^{-2\alpha} -
     214              :  *     e^{-(\gamma-\rho)}\right] \\
     215              :  *   \gamma^{xz}
     216              :  *   &= 0 \\
     217              :  *   \gamma^{yz}
     218              :  *   &= 0 \\
     219              :  *   \gamma^{zz} &= e^{-2\alpha}.
     220              :  * \f}
     221              :  *
     222              :  * The 4-velocity in spherical coordinates is given by
     223              :  *
     224              :  * \f{align}{
     225              :  *   u^{\bar{a}}=\frac{e^{-(\rho+\gamma)/2}}{\sqrt{1-v^2}}
     226              :  *   \left[1,0,0,\Omega\right],
     227              :  * \f}
     228              :  *
     229              :  * where
     230              :  *
     231              :  * \f{align}{
     232              :  *   v=(\Omega-\omega)r\sin(\theta)e^{-\rho}.
     233              :  * \f}
     234              :  *
     235              :  * Transforming to Cartesian coordinates we have
     236              :  *
     237              :  * \f{align}{
     238              :  *   u^t
     239              :  *   &=\frac{e^{-(\rho+\gamma)/2}}{\sqrt{1-v^2}} \\
     240              :  *   u^x
     241              :  *   &=\partial_\phi x u^\phi = -r\sin(\theta)\sin(\phi) u^t\Omega \\
     242              :  *   u^y
     243              :  *   &=\partial_\phi y u^\phi = r\sin(\theta)\cos(\phi) u^t\Omega \\
     244              :  *   u^z &= 0.
     245              :  * \f}
     246              :  *
     247              :  * The Lorentz factor is given by
     248              :  *
     249              :  * \f{align}{
     250              :  *   W
     251              :  *   &=Nu^t=e^{(\gamma+\rho)/2}\frac{e^{-(\rho+\gamma)/2}}{\sqrt{1-v^2}} \notag
     252              :  * \\
     253              :  *   &=\frac{1}{\sqrt{1-v^2}}.
     254              :  * \f}
     255              :  *
     256              :  * Using
     257              :  *
     258              :  * \f{align}{
     259              :  *   v^i = \frac{1}{N}\left(\frac{u^i}{u^t} + \beta^i\right)
     260              :  * \f}
     261              :  *
     262              :  * we get
     263              :  *
     264              :  * \f{align}{
     265              :  *   v^x
     266              :  *   &= -e^{-(\gamma+\rho)/2}r\sin(\theta)\sin(\phi)(\Omega-\omega)
     267              :  *   =-e^{-(\gamma-\rho)/2}\sin(\phi) v\\
     268              :  *   v^y
     269              :  *   &= e^{-(\gamma+\rho)/2}r\sin(\theta)\cos(\phi)(\Omega-\omega)
     270              :  *   = e^{-(\gamma-\rho)/2}\cos(\phi)v \\
     271              :  *   v^z&=0.
     272              :  * \f}
     273              :  *
     274              :  * Lowering with the spatial metric we get
     275              :  *
     276              :  * \f{align}{
     277              :  *   v_x
     278              :  *   &=\gamma_{xx} v^x + \gamma_{xy} v^y \notag \\
     279              :  *   &=-\left[\cos^2(\phi)e^{2\alpha}+\sin^2(\phi)e^{\gamma-\rho}\right]
     280              :  *     e^{-(\gamma-\rho)/2}\sin(\phi) v \notag \\
     281              :  *   &+\left[\cos(\phi)\sin(\phi)e^{2\alpha} -
     282              :  *     \sin(\phi)\cos(\phi)e^{\gamma-\rho}\right]
     283              :  *     e^{-(\gamma-\rho)/2}\cos(\phi)v \notag \\
     284              :  *   &=-e^{-(\gamma-\rho)/2}v\sin(\phi)e^{\gamma-\rho}
     285              :  *     \left[\sin^2(\phi)+\cos^2(\phi)\right] \notag \\
     286              :  *   &=-e^{(\gamma-\rho)/2}v\sin(\phi) \\
     287              :  *   v_y
     288              :  *   &=\gamma_{yx} v^x + \gamma_{yy} v^y \notag \\
     289              :  *   &=-\left[\cos(\phi)\sin(\phi)e^{2\alpha} - \sin(\phi)\cos(\phi)
     290              :  *     e^{(\gamma-\rho)}\right] e^{-(\gamma-\rho)/2}\sin(\phi) v \notag \\
     291              :  *   &+\left[\sin^2(\phi)e^{2\alpha} + \cos^2(\phi) e^{(\gamma-\rho)}\right]
     292              :  *     e^{-(\gamma-\rho)/2}\cos(\phi)v \notag \\
     293              :  *   &=e^{(\gamma-\rho)/2}v\cos(\phi) \\
     294              :  *   v_z &= 0.
     295              :  * \f}
     296              :  *
     297              :  * This is consistent with the Lorentz factor read off from \f$u^t\f$ since
     298              :  * \f$v^iv_i=v^2\f$. For completeness, \f$u_i=Wv_i\f$ so
     299              :  *
     300              :  * \f{align}{
     301              :  *   u_x
     302              :  *   &=-\frac{e^{(\gamma-\rho)/2}v\sin(\phi)}{\sqrt{1-v^2}} \\
     303              :  *   u_y
     304              :  *   &=\frac{e^{(\gamma-\rho)/2}v\cos(\phi)}{\sqrt{1-v^2}} \\
     305              :  *   u_z&=0.
     306              :  * \f}
     307              :  *
     308              :  * \warning Near (within `1e-2`) \f$r=0\f$ the numerical errors from
     309              :  * interpolation and computing the metric derivatives by finite difference no
     310              :  * longer cancel out and so the `tilde_s` time derivative only vanishes to
     311              :  * roughly `1e-8` rather than machine precision. Computing the Cartesian
     312              :  * derivatives from analytic differentiation of the radial and angular
     313              :  * polynomial fits might improve the situation but is a decent about of work to
     314              :  * implement.
     315              :  */
     316            1 : class RotatingStar : public virtual evolution::initial_data::InitialData,
     317              :                      public MarkAsAnalyticSolution,
     318              :                      public AnalyticSolution<3>,
     319              :                      public hydro::TemperatureInitialization<RotatingStar> {
     320              :   template <typename DataType>
     321            0 :   struct IntermediateVariables {
     322            0 :     IntermediateVariables(
     323              :         const tnsr::I<DataType, 3, Frame::Inertial>& in_coords,
     324              :         double in_delta_r);
     325              : 
     326            0 :     struct MetricData {
     327            0 :       MetricData() = default;
     328            0 :       MetricData(size_t num_points);
     329            0 :       DataType alpha;
     330            0 :       DataType rho;
     331            0 :       DataType gamma;
     332            0 :       DataType omega;
     333              :     };
     334              : 
     335            0 :     const tnsr::I<DataType, 3, Frame::Inertial>& coords;
     336            0 :     DataType radius;
     337            0 :     DataType phi;
     338            0 :     DataType cos_theta;
     339            0 :     DataType sin_theta;
     340              :     // Data at coords
     341            0 :     std::optional<DataType> rest_mass_density;
     342            0 :     std::optional<DataType> fluid_velocity;
     343            0 :     std::optional<MetricData> metric_data;
     344              :     // Data for 2nd-order FD derivatives.
     345              :     //
     346              :     // Technically we could do full-order accurate derivatives by using
     347              :     // Jacobians to transform quantities, but since we really want the initial
     348              :     // data to be solved natively in SpECTRE and satisfy the Einstein equations
     349              :     // with neutrinos and magnetic fields, doing the 2nd order FD like SpEC
     350              :     // should be fine.
     351              :     //
     352              :     // Note: We guard all the upper/lower values by checking if
     353              :     // metric_data_upper is computed. While not perfect, this simplifies the
     354              :     // code a lot.
     355            0 :     double delta_r;
     356            0 :     std::optional<std::array<MetricData, 3>> metric_data_upper{};
     357            0 :     std::optional<std::array<MetricData, 3>> metric_data_lower{};
     358            0 :     std::array<DataType, 3> radius_upper{};
     359            0 :     std::array<DataType, 3> radius_lower{};
     360            0 :     std::array<DataType, 3> cos_theta_upper{};
     361            0 :     std::array<DataType, 3> cos_theta_lower{};
     362            0 :     std::array<DataType, 3> sin_theta_upper{};
     363            0 :     std::array<DataType, 3> sin_theta_lower{};
     364            0 :     std::array<DataType, 3> phi_upper{};
     365            0 :     std::array<DataType, 3> phi_lower{};
     366              :   };
     367              : 
     368              :  public:
     369              :   /// The path to the RotNS data file.
     370            1 :   struct RotNsFilename {
     371            0 :     using type = std::string;
     372            0 :     static constexpr Options::String help = {
     373              :         "The path to the RotNS data file."};
     374              :   };
     375              : 
     376              :   /// The polytropic constant of the fluid.
     377              :   ///
     378              :   /// The data in the RotNS file will be rescaled.
     379            1 :   struct PolytropicConstant {
     380            0 :     using type = double;
     381            0 :     static constexpr Options::String help = {
     382              :         "The polytropic constant of the fluid."};
     383            0 :     static type lower_bound() { return 0.; }
     384              :   };
     385              : 
     386            0 :   using options = tmpl::list<Options::Alternatives<
     387              :       tmpl::list<RotNsFilename, PolytropicConstant>,
     388              :       tmpl::list<RotNsFilename,
     389              :                  hydro::OptionTags::InitialDataEquationOfState<true, 1>>>>;
     390              : 
     391            0 :   static constexpr Options::String help = {
     392              :       "Rotating neutron star initial data solved by the RotNS solver. The data "
     393              :       "is read in from disk. If RotNS used a polytropic equation of state, use "
     394              :       "the PolytropicConstant option, so that the polytropic index may be read "
     395              :       "directly from the RotNS output. Otherwise, set the equation of state to "
     396              :       "whichever one you would like to use to load the initial data. Note that "
     397              :       "if a polytrope is put into the EquationOfState option rather than using "
     398              :       "the PolytropicConstant option, the generic unit conversion will be used "
     399              :       "instead of the specific polytrope one."};
     400              : 
     401            0 :   RotatingStar() = default;
     402            0 :   RotatingStar(const RotatingStar& /*rhs*/);
     403            0 :   RotatingStar& operator=(const RotatingStar& /*rhs*/);
     404            0 :   RotatingStar(RotatingStar&& /*rhs*/) = default;
     405            0 :   RotatingStar& operator=(RotatingStar&& /*rhs*/) = default;
     406            0 :   ~RotatingStar() override = default;
     407              : 
     408            0 :   RotatingStar(std::string rot_ns_filename,
     409              :                std::unique_ptr<EquationsOfState::EquationOfState<true, 1>>
     410              :                    equation_of_state);
     411            0 :   RotatingStar(std::string rot_ns_filename, double polytropic_constant);
     412              : 
     413            0 :   auto get_clone() const
     414              :       -> std::unique_ptr<evolution::initial_data::InitialData> override;
     415              : 
     416              :   /// \cond
     417              :   explicit RotatingStar(CkMigrateMessage* msg);
     418              :   using PUP::able::register_constructor;
     419              :   WRAPPED_PUPable_decl_template(RotatingStar);
     420              :   /// \endcond
     421              : 
     422              :   /// Retrieve a collection of variables at `(x, t)`
     423              :   template <typename DataType, typename... Tags>
     424            1 :   tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
     425              :                                          const double /*t*/,
     426              :                                          tmpl::list<Tags...> /*meta*/) const {
     427              :     IntermediateVariables<DataType> intermediate_vars{
     428              :         x, 1.0e-4 * cst_solution_.equatorial_radius()};
     429              :     return {get<Tags>(variables(make_not_null(&intermediate_vars), x,
     430              :                                 tmpl::list<Tags>{}))...};
     431              :   }
     432              : 
     433              :   // NOLINTNEXTLINE(google-runtime-references)
     434            0 :   void pup(PUP::er& p) override;
     435              : 
     436            0 :   const EquationsOfState::EquationOfState<true, 1>& equation_of_state() const {
     437              :     return *equation_of_state_;
     438              :   }
     439              : 
     440            0 :   double equatorial_radius() const {
     441              :     return cst_solution_.equatorial_radius();
     442              :   }
     443              : 
     444              :  protected:
     445              :   template <typename DataType>
     446            0 :   using DerivLapse = ::Tags::deriv<gr::Tags::Lapse<DataType>, tmpl::size_t<3>,
     447              :                                    Frame::Inertial>;
     448              : 
     449              :   template <typename DataType>
     450            0 :   using DerivShift = ::Tags::deriv<gr::Tags::Shift<DataType, 3>,
     451              :                                    tmpl::size_t<3>, Frame::Inertial>;
     452              : 
     453              :   template <typename DataType>
     454            0 :   using DerivSpatialMetric = ::Tags::deriv<gr::Tags::SpatialMetric<DataType, 3>,
     455              :                                            tmpl::size_t<3>, Frame::Inertial>;
     456              : 
     457              :   template <typename DataType>
     458            0 :   void interpolate_vars_if_necessary(
     459              :       gsl::not_null<IntermediateVariables<DataType>*> vars) const;
     460              : 
     461              :   template <typename DataType>
     462            0 :   void interpolate_deriv_vars_if_necessary(
     463              :       gsl::not_null<IntermediateVariables<DataType>*> vars) const;
     464              : 
     465              :   template <typename DataType>
     466            0 :   Scalar<DataType> lapse(const DataType& gamma, const DataType& rho) const;
     467              : 
     468              :   template <typename DataType>
     469            0 :   tnsr::I<DataType, 3, Frame::Inertial> shift(const DataType& omega,
     470              :                                               const DataType& phi,
     471              :                                               const DataType& radius,
     472              :                                               const DataType& sin_theta) const;
     473              : 
     474              :   template <typename DataType>
     475            0 :   tnsr::ii<DataType, 3, Frame::Inertial> spatial_metric(
     476              :       const DataType& gamma, const DataType& rho, const DataType& alpha,
     477              :       const DataType& phi) const;
     478              : 
     479              :   template <typename DataType>
     480            0 :   tnsr::II<DataType, 3, Frame::Inertial> inverse_spatial_metric(
     481              :       const DataType& gamma, const DataType& rho, const DataType& alpha,
     482              :       const DataType& phi) const;
     483              : 
     484              :   template <typename DataType>
     485            0 :   Scalar<DataType> sqrt_det_spatial_metric(const DataType& gamma,
     486              :                                            const DataType& rho,
     487              :                                            const DataType& alpha) const;
     488              : 
     489              :   template <typename DataType>
     490            0 :   auto make_metric_data(size_t num_points) const ->
     491              :       typename IntermediateVariables<DataType>::MetricData;
     492              : 
     493              :   template <typename DataType>
     494            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     495              :                  const tnsr::I<DataType, 3>& x,
     496              :                  tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/)
     497              :       const -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
     498              : 
     499              :   template <typename DataType>
     500            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     501              :                  const tnsr::I<DataType, 3>& x,
     502              :                  tmpl::list<hydro::Tags::ElectronFraction<DataType>> /*meta*/)
     503              :       const -> tuples::TaggedTuple<hydro::Tags::ElectronFraction<DataType>>;
     504              : 
     505              :   template <typename DataType>
     506            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     507              :                  const tnsr::I<DataType, 3>& x,
     508              :                  tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/)
     509              :       const -> tuples::TaggedTuple<hydro::Tags::SpecificEnthalpy<DataType>>;
     510              : 
     511              :   template <typename DataType>
     512            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     513              :                  const tnsr::I<DataType, 3>& x,
     514              :                  tmpl::list<hydro::Tags::Temperature<DataType>> /*meta*/) const
     515              :       -> tuples::TaggedTuple<hydro::Tags::Temperature<DataType>>;
     516              : 
     517              :   template <typename DataType>
     518            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     519              :                  const tnsr::I<DataType, 3>& x,
     520              :                  tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/) const
     521              :       -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
     522              : 
     523              :   template <typename DataType>
     524            0 :   auto variables(
     525              :       gsl::not_null<IntermediateVariables<DataType>*> vars,
     526              :       const tnsr::I<DataType, 3>& x,
     527              :       tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/) const
     528              :       -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
     529              : 
     530              :   template <typename DataType>
     531            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     532              :                  const tnsr::I<DataType, 3>& x,
     533              :                  tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/)
     534              :       const -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>>;
     535              : 
     536              :   template <typename DataType>
     537            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     538              :                  const tnsr::I<DataType, 3>& x,
     539              :                  tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/)
     540              :       const -> tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataType>>;
     541              : 
     542              :   template <typename DataType>
     543            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     544              :                  const tnsr::I<DataType, 3>& x,
     545              :                  tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/)
     546              :       const -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
     547              : 
     548              :   template <typename DataType>
     549            0 :   auto variables(
     550              :       gsl::not_null<IntermediateVariables<DataType>*> vars,
     551              :       const tnsr::I<DataType, 3>& x,
     552              :       tmpl::list<hydro::Tags::DivergenceCleaningField<DataType>> /*meta*/) const
     553              :       -> tuples::TaggedTuple<hydro::Tags::DivergenceCleaningField<DataType>>;
     554              : 
     555              :   template <typename DataType>
     556            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     557              :                  const tnsr::I<DataType, 3>& x,
     558              :                  tmpl::list<gr::Tags::Lapse<DataType>> /*meta*/) const
     559              :       -> tuples::TaggedTuple<gr::Tags::Lapse<DataType>>;
     560              : 
     561              :   template <typename DataType>
     562            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     563              :                  const tnsr::I<DataType, 3>& x,
     564              :                  tmpl::list<gr::Tags::Shift<DataType, 3>> /*meta*/) const
     565              :       -> tuples::TaggedTuple<gr::Tags::Shift<DataType, 3>>;
     566              : 
     567              :   template <typename DataType>
     568            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     569              :                  const tnsr::I<DataType, 3>& x,
     570              :                  tmpl::list<gr::Tags::SpatialMetric<DataType, 3>> /*meta*/)
     571              :       const -> tuples::TaggedTuple<gr::Tags::SpatialMetric<DataType, 3>>;
     572              : 
     573              :   template <typename DataType>
     574            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     575              :                  const tnsr::I<DataType, 3>& x,
     576              :                  tmpl::list<gr::Tags::SqrtDetSpatialMetric<DataType>> /*meta*/)
     577              :       const -> tuples::TaggedTuple<gr::Tags::SqrtDetSpatialMetric<DataType>>;
     578              : 
     579              :   template <typename DataType>
     580            0 :   auto variables(
     581              :       gsl::not_null<IntermediateVariables<DataType>*> vars,
     582              :       const tnsr::I<DataType, 3>& x,
     583              :       tmpl::list<gr::Tags::InverseSpatialMetric<DataType, 3>> /*meta*/) const
     584              :       -> tuples::TaggedTuple<gr::Tags::InverseSpatialMetric<DataType, 3>>;
     585              : 
     586              :   template <typename DataType>
     587            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     588              :                  const tnsr::I<DataType, 3>& x,
     589              :                  tmpl::list<::Tags::dt<gr::Tags::Lapse<DataType>>> /*meta*/)
     590              :       const -> tuples::TaggedTuple<::Tags::dt<gr::Tags::Lapse<DataType>>>;
     591              : 
     592              :   template <typename DataType>
     593            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     594              :                  const tnsr::I<DataType, 3>& x,
     595              :                  tmpl::list<DerivLapse<DataType>> /*meta*/) const
     596              :       -> tuples::TaggedTuple<DerivLapse<DataType>>;
     597              : 
     598              :   template <typename DataType>
     599            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     600              :                  const tnsr::I<DataType, 3>& x,
     601              :                  tmpl::list<::Tags::dt<gr::Tags::Shift<DataType, 3>>> /*meta*/)
     602              :       const -> tuples::TaggedTuple<::Tags::dt<gr::Tags::Shift<DataType, 3>>>;
     603              : 
     604              :   template <typename DataType>
     605            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     606              :                  const tnsr::I<DataType, 3>& x,
     607              :                  tmpl::list<DerivShift<DataType>> /*meta*/) const
     608              :       -> tuples::TaggedTuple<DerivShift<DataType>>;
     609              : 
     610              :   template <typename DataType>
     611            0 :   auto variables(
     612              :       gsl::not_null<IntermediateVariables<DataType>*> vars,
     613              :       const tnsr::I<DataType, 3>& x,
     614              :       tmpl::list<::Tags::dt<gr::Tags::SpatialMetric<DataType, 3>>> /*meta*/)
     615              :       const
     616              :       -> tuples::TaggedTuple<::Tags::dt<gr::Tags::SpatialMetric<DataType, 3>>>;
     617              : 
     618              :   template <typename DataType>
     619            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     620              :                  const tnsr::I<DataType, 3>& x,
     621              :                  tmpl::list<DerivSpatialMetric<DataType>> /*meta*/) const
     622              :       -> tuples::TaggedTuple<DerivSpatialMetric<DataType>>;
     623              : 
     624              :   template <typename DataType>
     625            0 :   auto variables(gsl::not_null<IntermediateVariables<DataType>*> vars,
     626              :                  const tnsr::I<DataType, 3>& x,
     627              :                  tmpl::list<gr::Tags::ExtrinsicCurvature<DataType, 3>> /*meta*/)
     628              :       const -> tuples::TaggedTuple<gr::Tags::ExtrinsicCurvature<DataType, 3>>;
     629              : 
     630            0 :   friend bool operator==(const RotatingStar& lhs, const RotatingStar& rhs);
     631              : 
     632            0 :   std::string rot_ns_filename_{};
     633            0 :   detail::CstSolution cst_solution_{};
     634            0 :   double polytropic_constant_ = std::numeric_limits<double>::signaling_NaN();
     635            0 :   double polytropic_exponent_ = std::numeric_limits<double>::signaling_NaN();
     636            0 :   bool is_polytrope_{};
     637              :   std::unique_ptr<EquationsOfState::EquationOfState<true, 1>>
     638            0 :       equation_of_state_;
     639              :   // Floor value to protect EoS from encountering FPEs when computing state
     640              :   // variables in the atmosphere
     641            0 :   static constexpr double atmosphere_floor_ = 1.e-50;
     642              : };
     643              : 
     644            0 : bool operator!=(const RotatingStar& lhs, const RotatingStar& rhs);
     645              : }  // namespace RelativisticEuler::Solutions
        

Generated by: LCOV version 2.0-1