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 <limits>
8 : #include <string>
9 :
10 : #include "DataStructures/Tensor/TypeAliases.hpp"
11 : #include "Options/String.hpp"
12 : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
13 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Minkowski.hpp"
14 : #include "PointwiseFunctions/AnalyticSolutions/GrMhd/Solutions.hpp"
15 : #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
16 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
17 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
18 : #include "Utilities/Serialization/CharmPupable.hpp"
19 : #include "Utilities/TMPL.hpp"
20 : #include "Utilities/TaggedTuple.hpp"
21 :
22 : /// \cond
23 : namespace PUP {
24 : class er;
25 : } // namespace PUP
26 : /// \endcond
27 :
28 : namespace grmhd::Solutions {
29 :
30 : /*!
31 : * \brief A one-dimensional shock solution for an ideal fluid in Minkowski
32 : * spacetime
33 : *
34 : * This solution consists of a left state for \f$x<0\f$ and a right state for
35 : * \f$x\ge 0\f$, each with constant fluid variables. The interface between these
36 : * states moves with the shock speed \f$\mu\f$ as described in
37 : * \cite Komissarov1999.
38 : *
39 : * \note We do not currently support 1D RMHD, so this class provides a 3D
40 : * solution with \f$x\f$-dependence only. Therefore the computational domain can
41 : * be represented by a single element with periodic boundary conditions in the
42 : * \f$y\f$ and \f$z\f$ directions.
43 : */
44 1 : class KomissarovShock
45 : : public evolution::initial_data::InitialData,
46 : public AnalyticSolution,
47 : public hydro::TemperatureInitialization<KomissarovShock>,
48 : public MarkAsAnalyticSolution {
49 : public:
50 0 : using equation_of_state_type = EquationsOfState::IdealFluid<true>;
51 :
52 0 : struct AdiabaticIndex {
53 0 : using type = double;
54 0 : static constexpr Options::String help = {
55 : "The adiabatic index of the ideal fluid"};
56 0 : static type lower_bound() { return 1.0; }
57 : };
58 0 : struct LeftRestMassDensity {
59 0 : using type = double;
60 0 : static std::string name() { return "LeftDensity"; };
61 0 : static constexpr Options::String help = {
62 : "Fluid rest mass density in the left half-domain"};
63 0 : static type lower_bound() { return 0.0; }
64 : };
65 0 : struct RightRestMassDensity {
66 0 : using type = double;
67 0 : static std::string name() { return "RightDensity"; };
68 0 : static constexpr Options::String help = {
69 : "Fluid rest mass density in the right half-domain"};
70 0 : static type lower_bound() { return 0.0; }
71 : };
72 :
73 0 : struct LeftElectronFraction {
74 0 : using type = double;
75 0 : static std::string name() { return "LeftElectronFraction"; };
76 0 : static constexpr Options::String help = {
77 : "Fluid electron fraction in the left half-domain"};
78 0 : static type lower_bound() { return 0.0; }
79 0 : static type upper_bound() { return 1.0; }
80 : };
81 0 : struct RightElectronFraction {
82 0 : using type = double;
83 0 : static std::string name() { return "RightElectronFraction"; };
84 0 : static constexpr Options::String help = {
85 : "Fluid electron fraction in the right half-domain"};
86 0 : static type lower_bound() { return 0.0; }
87 0 : static type upper_bound() { return 1.0; }
88 : };
89 0 : struct LeftPressure {
90 0 : using type = double;
91 0 : static constexpr Options::String help = {
92 : "Fluid pressure in the left half-domain"};
93 0 : static type lower_bound() { return 0.0; }
94 : };
95 0 : struct RightPressure {
96 0 : using type = double;
97 0 : static constexpr Options::String help = {
98 : "Fluid pressure in the right half-domain"};
99 0 : static type lower_bound() { return 0.0; }
100 : };
101 0 : struct LeftSpatialVelocity {
102 0 : using type = std::array<double, 3>;
103 0 : static std::string name() { return "LeftVelocity"; };
104 0 : static constexpr Options::String help = {
105 : "Fluid spatial velocity in the left half-domain"};
106 : };
107 0 : struct RightSpatialVelocity {
108 0 : using type = std::array<double, 3>;
109 0 : static std::string name() { return "RightVelocity"; };
110 0 : static constexpr Options::String help = {
111 : "Fluid spatial velocity in the right half-domain"};
112 : };
113 0 : struct LeftMagneticField {
114 0 : using type = std::array<double, 3>;
115 0 : static constexpr Options::String help = {
116 : "Magnetic field in the left half-domain"};
117 : };
118 0 : struct RightMagneticField {
119 0 : using type = std::array<double, 3>;
120 0 : static constexpr Options::String help = {
121 : "Magnetic field in the right half-domain"};
122 : };
123 0 : struct ShockSpeed {
124 0 : using type = double;
125 0 : static constexpr Options::String help = {"Propagation speed of the shock"};
126 : };
127 :
128 0 : using options =
129 : tmpl::list<AdiabaticIndex, LeftRestMassDensity, RightRestMassDensity,
130 : LeftElectronFraction, RightElectronFraction, LeftPressure,
131 : RightPressure, LeftSpatialVelocity, RightSpatialVelocity,
132 : LeftMagneticField, RightMagneticField, ShockSpeed>;
133 :
134 0 : static constexpr Options::String help = {
135 : "Analytic initial data for a Komissarov shock test. The fluid variables "
136 : "are set homogeneously on either half of the domain left and right of "
137 : "x=0."};
138 :
139 0 : KomissarovShock() = default;
140 0 : KomissarovShock(const KomissarovShock& /*rhs*/) = default;
141 0 : KomissarovShock& operator=(const KomissarovShock& /*rhs*/) = default;
142 0 : KomissarovShock(KomissarovShock&& /*rhs*/) = default;
143 0 : KomissarovShock& operator=(KomissarovShock&& /*rhs*/) = default;
144 0 : ~KomissarovShock() override = default;
145 :
146 0 : KomissarovShock(double adiabatic_index, double left_rest_mass_density,
147 : double right_rest_mass_density, double left_electron_fraction,
148 : double right_electron_fraction, double left_pressure,
149 : double right_pressure,
150 : const std::array<double, 3>& left_spatial_velocity,
151 : const std::array<double, 3>& right_spatial_velocity,
152 : const std::array<double, 3>& left_magnetic_field,
153 : const std::array<double, 3>& right_magnetic_field,
154 : double shock_speed);
155 :
156 0 : auto get_clone() const
157 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
158 :
159 : /// \cond
160 : explicit KomissarovShock(CkMigrateMessage* msg);
161 : using PUP::able::register_constructor;
162 : WRAPPED_PUPable_decl_template(KomissarovShock);
163 : /// \endcond
164 :
165 : /// @{
166 : /// Retrieve the GRMHD variables at a given position.
167 : template <typename DataType>
168 1 : auto variables(const tnsr::I<DataType, 3>& x, double t,
169 : tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/)
170 : const -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
171 :
172 : template <typename DataType>
173 1 : auto variables(const tnsr::I<DataType, 3>& x, double t,
174 : tmpl::list<hydro::Tags::ElectronFraction<DataType>> /*meta*/)
175 : const -> tuples::TaggedTuple<hydro::Tags::ElectronFraction<DataType>>;
176 :
177 : template <typename DataType>
178 1 : auto variables(
179 : const tnsr::I<DataType, 3>& x, double t,
180 : tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/) const
181 : -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
182 :
183 : template <typename DataType>
184 1 : auto variables(const tnsr::I<DataType, 3>& x, double t,
185 : tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/) const
186 : -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
187 :
188 : template <typename DataType>
189 1 : auto variables(const tnsr::I<DataType, 3>& x, double t,
190 : tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/)
191 : const -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>>;
192 :
193 : template <typename DataType>
194 1 : auto variables(const tnsr::I<DataType, 3>& x, double t,
195 : tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/)
196 : const -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
197 :
198 : template <typename DataType>
199 1 : auto variables(
200 : const tnsr::I<DataType, 3>& x, double t,
201 : tmpl::list<hydro::Tags::DivergenceCleaningField<DataType>> /*meta*/) const
202 : -> tuples::TaggedTuple<hydro::Tags::DivergenceCleaningField<DataType>>;
203 :
204 : template <typename DataType>
205 1 : auto variables(const tnsr::I<DataType, 3>& x, double t,
206 : tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/)
207 : const -> tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataType>>;
208 :
209 : template <typename DataType>
210 1 : auto variables(const tnsr::I<DataType, 3>& x, double t,
211 : tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/)
212 : const -> tuples::TaggedTuple<hydro::Tags::SpecificEnthalpy<DataType>>;
213 :
214 : template <typename DataType>
215 1 : auto variables(const tnsr::I<DataType, 3>& x, double t,
216 : tmpl::list<hydro::Tags::Temperature<DataType>> /*meta*/) const
217 : -> tuples::TaggedTuple<hydro::Tags::Temperature<DataType>> {
218 : return TemperatureInitialization::variables(
219 : x, t, tmpl::list<hydro::Tags::Temperature<DataType>>{});
220 : }
221 : /// @}
222 :
223 : /// Retrieve a collection of hydro variables at `(x, t)`
224 : template <typename DataType, typename Tag1, typename Tag2, typename... Tags>
225 1 : tuples::TaggedTuple<Tag1, Tag2, Tags...> variables(
226 : const tnsr::I<DataType, 3>& x, double t,
227 : tmpl::list<Tag1, Tag2, Tags...> /*meta*/) const {
228 : return {tuples::get<Tag1>(variables(x, t, tmpl::list<Tag1>{})),
229 : tuples::get<Tag2>(variables(x, t, tmpl::list<Tag2>{})),
230 : tuples::get<Tags>(variables(x, t, tmpl::list<Tags>{}))...};
231 : }
232 :
233 : /// Retrieve the metric variables
234 : template <typename DataType, typename Tag,
235 : Requires<tmpl::list_contains_v<
236 : gr::analytic_solution_tags<3, DataType>, Tag>> = nullptr>
237 1 : tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x, double t,
238 : tmpl::list<Tag> /*meta*/) const {
239 : return background_spacetime_.variables(x, t, tmpl::list<Tag>{});
240 : }
241 :
242 0 : const EquationsOfState::IdealFluid<true>& equation_of_state() const {
243 : return equation_of_state_;
244 : }
245 :
246 : // NOLINTNEXTLINE(google-runtime-references)
247 0 : void pup(PUP::er& /*p*/) override;
248 :
249 : protected:
250 0 : EquationsOfState::IdealFluid<true> equation_of_state_{};
251 0 : gr::Solutions::Minkowski<3> background_spacetime_{};
252 :
253 0 : double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
254 0 : double left_rest_mass_density_ = std::numeric_limits<double>::signaling_NaN();
255 0 : double right_rest_mass_density_ =
256 : std::numeric_limits<double>::signaling_NaN();
257 0 : double left_electron_fraction_ = std::numeric_limits<double>::signaling_NaN();
258 0 : double right_electron_fraction_ =
259 : std::numeric_limits<double>::signaling_NaN();
260 0 : double left_pressure_ = std::numeric_limits<double>::signaling_NaN();
261 0 : double right_pressure_ = std::numeric_limits<double>::signaling_NaN();
262 0 : std::array<double, 3> left_spatial_velocity_{
263 : {std::numeric_limits<double>::signaling_NaN(),
264 : std::numeric_limits<double>::signaling_NaN(),
265 : std::numeric_limits<double>::signaling_NaN()}};
266 0 : std::array<double, 3> right_spatial_velocity_{
267 : {std::numeric_limits<double>::signaling_NaN(),
268 : std::numeric_limits<double>::signaling_NaN(),
269 : std::numeric_limits<double>::signaling_NaN()}};
270 0 : std::array<double, 3> left_magnetic_field_{
271 : {std::numeric_limits<double>::signaling_NaN(),
272 : std::numeric_limits<double>::signaling_NaN(),
273 : std::numeric_limits<double>::signaling_NaN()}};
274 0 : std::array<double, 3> right_magnetic_field_{
275 : {std::numeric_limits<double>::signaling_NaN(),
276 : std::numeric_limits<double>::signaling_NaN(),
277 : std::numeric_limits<double>::signaling_NaN()}};
278 0 : double shock_speed_ = std::numeric_limits<double>::signaling_NaN();
279 :
280 0 : friend bool operator==(const KomissarovShock& lhs,
281 : const KomissarovShock& rhs);
282 :
283 0 : friend bool operator!=(const KomissarovShock& lhs,
284 : const KomissarovShock& rhs);
285 : };
286 :
287 : } // namespace grmhd::Solutions
|
