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             : 
       8             : #include "DataStructures/Tensor/Tensor.hpp"
       9             : #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
      10             : #include "Options/String.hpp"
      11             : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
      12             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      13             : #include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp"
      14             : #include "PointwiseFunctions/Hydro/Tags.hpp"
      15             : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
      16             : #include "Utilities/TMPL.hpp"
      17             : #include "Utilities/TaggedTuple.hpp"
      18             : 
      19             : /// \cond
      20             : namespace PUP {
      21             : class er;  // IWYU pragma: keep
      22             : }  // namespace PUP
      23             : /// \endcond
      24             : 
      25             : // IWYU pragma: no_include <pup.h>
      26             : 
      27             : namespace NewtonianEuler::Solutions {
      28             : 
      29             : /*!
      30             :  * \brief A static spherically symmetric star in Newtonian gravity
      31             :  *
      32             :  * The solution for a static, spherically-symmetric star in 3 dimensions, found
      33             :  * by solving the Lane-Emden equation \cite Chandrasekhar1939
      34             :  * \cite Shapiro1983 .
      35             :  * The Lane-Emden equation has closed-form solutions for certain equations of
      36             :  * state; this class implements the solution for a polytropic fluid with
      37             :  * polytropic exponent \f$\Gamma=2\f$ (i.e., with polytropic index \f$n=1\f$).
      38             :  * The solution is returned in units where \f$G=1\f$, with \f$G\f$ the
      39             :  * gravitational constant.
      40             :  *
      41             :  * The radius and mass of the star are determined by the polytropic constant
      42             :  * \f$\kappa\f$ and central density \f$\rho_c\f$.
      43             :  * The radius is \f$R = \pi \alpha\f$,
      44             :  * and the mass is \f$M = 4 \pi^2 \alpha^3 \rho_c\f$,
      45             :  * where \f$\alpha = \sqrt{\kappa / (2 \pi)}\f$ and \f$G=1\f$.
      46             :  */
      47           1 : class LaneEmdenStar : public evolution::initial_data::InitialData,
      48             :                       public MarkAsAnalyticSolution {
      49             :  public:
      50           0 :   using equation_of_state_type = EquationsOfState::PolytropicFluid<false>;
      51             : 
      52             :   /// The central mass density of the star.
      53           1 :   struct CentralMassDensity {
      54           0 :     using type = double;
      55           0 :     static constexpr Options::String help = {
      56             :         "The central mass density of the star."};
      57           0 :     static type lower_bound() { return 0.; }
      58             :   };
      59             : 
      60             :   /// The polytropic constant of the polytropic fluid.
      61           1 :   struct PolytropicConstant {
      62           0 :     using type = double;
      63           0 :     static constexpr Options::String help = {
      64             :         "The polytropic constant of the fluid."};
      65           0 :     static type lower_bound() { return 0.; }
      66             :   };
      67             : 
      68           0 :   using options = tmpl::list<CentralMassDensity, PolytropicConstant>;
      69             : 
      70           0 :   static constexpr Options::String help = {
      71             :       "A static, spherically-symmetric star in Newtonian gravity, found by\n"
      72             :       "solving the Lane-Emden equations, with a given central density and\n"
      73             :       "polytropic fluid. The fluid has polytropic index 1, but the polytropic\n"
      74             :       "constant is specifiable"};
      75             : 
      76           0 :   LaneEmdenStar() = default;
      77           0 :   LaneEmdenStar(const LaneEmdenStar& /*rhs*/) = default;
      78           0 :   LaneEmdenStar& operator=(const LaneEmdenStar& /*rhs*/) = default;
      79           0 :   LaneEmdenStar(LaneEmdenStar&& /*rhs*/) = default;
      80           0 :   LaneEmdenStar& operator=(LaneEmdenStar&& /*rhs*/) = default;
      81           0 :   ~LaneEmdenStar() override = default;
      82             : 
      83           0 :   auto get_clone() const
      84             :       -> std::unique_ptr<evolution::initial_data::InitialData> override;
      85             : 
      86             :   /// \cond
      87             :   explicit LaneEmdenStar(CkMigrateMessage* msg);
      88             :   using PUP::able::register_constructor;
      89             :   WRAPPED_PUPable_decl_template(LaneEmdenStar);
      90             :   /// \endcond
      91             : 
      92           0 :   LaneEmdenStar(double central_mass_density, double polytropic_constant);
      93             : 
      94             :   /// Retrieve a collection of variables at `(x, t)`
      95             :   template <typename DataType, typename... Tags>
      96           1 :   tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
      97             :                                          const double /*t*/,
      98             :                                          tmpl::list<Tags...> /*meta*/) const {
      99             :     const auto mass_density = precompute_mass_density(x);
     100             :     return {tuples::get<Tags>(variables(tmpl::list<Tags>{}, mass_density))...};
     101             :   }
     102             : 
     103             :   /// \brief Compute the gravitational field for the corresponding source term,
     104             :   /// LaneEmdenGravitationalField.
     105             :   ///
     106             :   /// The result is the vector-field giving the acceleration due to gravity
     107             :   /// that is felt by a test particle.
     108             :   template <typename DataType>
     109           1 :   tnsr::I<DataType, 3> gravitational_field(const tnsr::I<DataType, 3>& x) const;
     110             : 
     111           0 :   const EquationsOfState::PolytropicFluid<false>& equation_of_state() const {
     112             :     return equation_of_state_;
     113             :   }
     114             : 
     115             :   // NOLINTNEXTLINE(google-runtime-references)
     116           0 :   void pup(PUP::er& /*p*/) override;
     117             : 
     118             :  private:
     119             :   template <typename DataType>
     120           0 :   Scalar<DataType> precompute_mass_density(const tnsr::I<DataType, 3>& x) const;
     121             : 
     122             :   template <typename DataType>
     123           0 :   tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>> variables(
     124             :       tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/,
     125             :       const Scalar<DataType>& mass_density) const;
     126             : 
     127             :   template <typename DataType>
     128           0 :   tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>> variables(
     129             :       tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/,
     130             :       const Scalar<DataType>& mass_density) const;
     131             : 
     132             :   template <typename DataType>
     133           0 :   tuples::TaggedTuple<hydro::Tags::Pressure<DataType>> variables(
     134             :       tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/,
     135             :       const Scalar<DataType>& mass_density) const;
     136             : 
     137             :   template <typename DataType>
     138           0 :   tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>> variables(
     139             :       tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/,
     140             :       const Scalar<DataType>& mass_density) const;
     141             : 
     142           0 :   friend bool operator==(const LaneEmdenStar& lhs, const LaneEmdenStar& rhs);
     143             : 
     144           0 :   double central_mass_density_ = std::numeric_limits<double>::signaling_NaN();
     145           0 :   double polytropic_constant_ = std::numeric_limits<double>::signaling_NaN();
     146           0 :   EquationsOfState::PolytropicFluid<false> equation_of_state_{};
     147             : };
     148             : 
     149           0 : bool operator!=(const LaneEmdenStar& lhs, const LaneEmdenStar& rhs);
     150             : 
     151             : }  // namespace NewtonianEuler::Solutions