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