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