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 <limits>
9 :
10 : #include "DataStructures/DataVector.hpp"
11 : #include "DataStructures/Tensor/Tensor.hpp"
12 : #include "DataStructures/Tensor/TypeAliases.hpp"
13 : #include "Evolution/Systems/NewtonianEuler/Sources/NoSource.hpp"
14 : #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
15 : #include "Options/String.hpp"
16 : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
17 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
18 : #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
19 : #include "Utilities/MakeArray.hpp"
20 : #include "Utilities/TMPL.hpp"
21 : #include "Utilities/TaggedTuple.hpp"
22 :
23 : /// \cond
24 : namespace PUP {
25 : class er; // IWYU pragma: keep
26 : } // namespace PUP
27 : /// \endcond
28 :
29 : namespace NewtonianEuler {
30 : namespace Solutions {
31 :
32 : /*!
33 : * \brief Analytic solution to the Riemann Problem
34 : *
35 : * This class implements the exact Riemann solver described in detail in
36 : * Chapter 4 of \cite Toro2009. We follow the notation there.
37 : * The algorithm implemented here allows for 1, 2 and 3D wave propagation
38 : * along any coordinate axis. Typical initial data for test cases
39 : * (see \cite Toro2009) include:
40 : *
41 : * - Sod's Shock Tube (shock on the right, rarefaction on the left):
42 : * - \f$(\rho_L, u_L, p_L) = (1.0, 0.0, 1.0)\f$
43 : * - \f$(\rho_R, u_R, p_R) = (0.125, 0.0, 0.1)\f$
44 : * - Recommended setup for sample run:
45 : * - InitialTimeStep: 0.0001
46 : * - Final time: 0.2
47 : * - DomainCreator along wave propagation (no AMR):
48 : * - Interval of length 1
49 : * - InitialRefinement: 8
50 : * - InitialGridPoints: 2
51 : *
52 : * - "123" problem (two symmetric rarefaction waves):
53 : * - \f$(\rho_L, u_L, p_L) = (1.0, -2.0, 0.4)\f$
54 : * - \f$(\rho_R, u_R, p_R) = (1.0, 2.0, 0.4)\f$
55 : * - Recommended setup for sample run:
56 : * - InitialTimeStep: 0.0001
57 : * - Final time: 0.15
58 : * - DomainCreator along wave propagation (no AMR):
59 : * - Interval of length 1
60 : * - InitialRefinement: 8
61 : * - InitialGridPoints: 2
62 : *
63 : * - Collision of two blast waves (this test is challenging):
64 : * - \f$(\rho_L, u_L, p_L) = (5.99924, 19.5975, 460.894)\f$
65 : * - \f$(\rho_R, u_R, p_R) = (5.99242, -6.19633, 46.0950)\f$
66 : * - Recommended setup for sample run:
67 : * - InitialTimeStep: 0.00001
68 : * - Final time: 0.012
69 : * - DomainCreator along wave propagation (no AMR):
70 : * - Interval of length 1
71 : * - InitialRefinement: 8
72 : * - InitialGridPoints: 2
73 : *
74 : * - Lax problem:
75 : * - \f$(\rho_L, u_L, p_L) = (0.445, 0.698, 3.528)\f$
76 : * - \f$(\rho_R, u_R, p_R) = (0.5, 0.0, 0.571)\f$
77 : * - Recommended setup for sample run:
78 : * - InitialTimeStep: 0.00001
79 : * - Final time: 0.1
80 : * - DomainCreator along wave propagation (no AMR):
81 : * - Interval of length 1
82 : * - InitialRefinement: 8
83 : * - InitialGridPoints: 2
84 : *
85 : * where \f$\rho\f$ is the mass density, \f$p\f$ is the pressure, and
86 : * \f$u\f$ denotes the normal velocity.
87 : *
88 : * \note Currently the propagation axis must be hard-coded as a `size_t`
89 : * private member variable `propagation_axis_`, which can take one of
90 : * the three values `PropagationAxis::X`, `PropagationAxis::Y`, and
91 : * `PropagationAxis::Z`.
92 : *
93 : * \details The algorithm makes use of the following recipe:
94 : *
95 : * - Given the initial data on both sides of the initial interface of the
96 : * discontinuity (here called "left" and "right" sides, where a coordinate
97 : * axis points from left to right), we compute the pressure,
98 : * \f$p_*\f$, and the normal velocity,
99 : * \f$u_*\f$, in the so-called star region. This is done in the constructor.
100 : * Here "normal" refers to the normal direction to the initial interface.
101 : *
102 : * - Given the pressure and the normal velocity in the star region, two
103 : * `Wave` `struct`s are created, which represent the waves propagating
104 : * at later times on each side of the contact discontinuity. Each `Wave`
105 : * is equipped with two `struct`s named `Shock` and `Rarefaction` which
106 : * contain functions that compute the primitive variables depending on whether
107 : * the wave is a shock or a rarefaction.
108 : *
109 : * - If \f$p_* > p_K\f$, the wave is a shock, otherwise the wave is a
110 : * rarefaction. Here \f$K\f$ stands for \f$L\f$ or \f$R\f$: the left
111 : * and right initial pressure, respectively. Since this comparison can't be
112 : * performed at compile time, each `Wave` holds a `bool` member `is_shock_`
113 : * which is `true` if it is a shock, and `false` if it is a
114 : * rarefaction wave. This variable is used to evaluate the correct functions
115 : * at run time.
116 : *
117 : * - In order to obtain the primitives at a certain time and spatial location,
118 : * we evaluate whether the spatial location is on the left of the propagating
119 : * contact discontinuity \f$(x < u_* t)\f$ or on the right \f$(x > u_* t)\f$,
120 : * and we use the corresponding functions for left or right `Wave`s,
121 : * respectively.
122 : *
123 : * \note The characterization of each propagating wave will only
124 : * depend on the normal velocity, while the initial jump in the components of
125 : * the velocity transverse to the wave propagation will be advected at the
126 : * speed of the contact discontinuity (\f$u_*\f$).
127 : */
128 : template <size_t Dim>
129 1 : class RiemannProblem : public MarkAsAnalyticSolution {
130 0 : enum class Side { Left, Right };
131 0 : enum PropagationAxis { X = 0, Y = 1, Z = 2 };
132 :
133 : struct Wave;
134 : struct Shock;
135 : struct Rarefaction;
136 :
137 : public:
138 0 : using equation_of_state_type = EquationsOfState::IdealFluid<false>;
139 0 : using source_term_type = Sources::NoSource;
140 :
141 : /// The adiabatic index of the fluid.
142 1 : struct AdiabaticIndex {
143 0 : using type = double;
144 0 : static constexpr Options::String help = {
145 : "The adiabatic index of the fluid."};
146 : };
147 :
148 : /// Initial position of the discontinuity
149 1 : struct InitialPosition {
150 0 : using type = double;
151 0 : static constexpr Options::String help = {
152 : "The initial position of the discontinuity."};
153 : };
154 :
155 : /// The mass density on the left of the initial discontinuity
156 1 : struct LeftMassDensity {
157 0 : using type = double;
158 0 : static constexpr Options::String help = {"The left mass density."};
159 : };
160 :
161 : /// The velocity on the left of the initial discontinuity
162 1 : struct LeftVelocity {
163 0 : using type = std::array<double, Dim>;
164 0 : static constexpr Options::String help = {"The left velocity."};
165 : };
166 :
167 : /// The pressure on the left of the initial discontinuity
168 1 : struct LeftPressure {
169 0 : using type = double;
170 0 : static constexpr Options::String help = {"The left pressure."};
171 : };
172 :
173 : /// The mass density on the right of the initial discontinuity
174 1 : struct RightMassDensity {
175 0 : using type = double;
176 0 : static constexpr Options::String help = {"The right mass density."};
177 : };
178 :
179 : /// The velocity on the right of the initial discontinuity
180 1 : struct RightVelocity {
181 0 : using type = std::array<double, Dim>;
182 0 : static constexpr Options::String help = {"The right velocity."};
183 : };
184 :
185 : /// The pressure on the right of the initial discontinuity
186 1 : struct RightPressure {
187 0 : using type = double;
188 0 : static constexpr Options::String help = {"The right pressure."};
189 : };
190 :
191 : /// The tolerance for solving for \f$p_*\f$.
192 1 : struct PressureStarTol {
193 0 : using type = double;
194 0 : static constexpr Options::String help = {
195 : "The tolerance for the numerical solution for p star"};
196 0 : static type suggested_value() { return 1.e-9; }
197 : };
198 :
199 : // Any of the two states that constitute the initial data, including
200 : // some derived quantities that are used repeatedly at each evaluation of the
201 : // different waves of the solution.
202 : /// Holds initial data on a side of the discontinuity and related quantities
203 1 : struct InitialData {
204 0 : InitialData() = default;
205 0 : InitialData(const InitialData& /*rhs*/) = default;
206 0 : InitialData& operator=(const InitialData& /*rhs*/) = default;
207 0 : InitialData(InitialData&& /*rhs*/) = default;
208 0 : InitialData& operator=(InitialData&& /*rhs*/) = default;
209 0 : ~InitialData() = default;
210 :
211 0 : InitialData(double mass_density, const std::array<double, Dim>& velocity,
212 : double pressure, double adiabatic_index,
213 : size_t propagation_axis);
214 :
215 : // clang-tidy: no runtime references
216 0 : void pup(PUP::er& /*p*/); // NOLINT
217 :
218 0 : double mass_density_ = std::numeric_limits<double>::signaling_NaN();
219 0 : std::array<double, Dim> velocity_ =
220 : make_array<Dim>(std::numeric_limits<double>::signaling_NaN());
221 0 : double pressure_ = std::numeric_limits<double>::signaling_NaN();
222 0 : double sound_speed_ = std::numeric_limits<double>::signaling_NaN();
223 0 : double normal_velocity_ = std::numeric_limits<double>::signaling_NaN();
224 :
225 : // Data-dependent constants A and B in Eqns. (4.8) of Toro.
226 0 : double constant_a_ = std::numeric_limits<double>::signaling_NaN();
227 0 : double constant_b_ = std::numeric_limits<double>::signaling_NaN();
228 :
229 0 : friend bool operator==(const InitialData& lhs, const InitialData& rhs) {
230 : return lhs.mass_density_ == rhs.mass_density_ and
231 : lhs.velocity_ == rhs.velocity_ and lhs.pressure_ == rhs.pressure_;
232 : }
233 : };
234 :
235 0 : using options = tmpl::list<AdiabaticIndex, InitialPosition, LeftMassDensity,
236 : LeftVelocity, LeftPressure, RightMassDensity,
237 : RightVelocity, RightPressure, PressureStarTol>;
238 :
239 0 : static constexpr Options::String help = {
240 : "Riemann Problem in 1, 2 or 3D along any coordinate axis."};
241 :
242 0 : RiemannProblem() = default;
243 0 : RiemannProblem(const RiemannProblem& /*rhs*/) = default;
244 0 : RiemannProblem& operator=(const RiemannProblem& /*rhs*/) = default;
245 0 : RiemannProblem(RiemannProblem&& /*rhs*/) = default;
246 0 : RiemannProblem& operator=(RiemannProblem&& /*rhs*/) = default;
247 0 : ~RiemannProblem() = default;
248 :
249 0 : RiemannProblem(double adiabatic_index, double initial_position,
250 : double left_mass_density,
251 : const std::array<double, Dim>& left_velocity,
252 : double left_pressure, double right_mass_density,
253 : const std::array<double, Dim>& right_velocity,
254 : double right_pressure,
255 : double pressure_star_tol = PressureStarTol::suggested_value());
256 :
257 : /// Retrieve a collection of hydrodynamic variables at position `x`
258 : /// and time `t`
259 : template <typename DataType, typename... Tags>
260 1 : tuples::TaggedTuple<Tags...> variables(
261 : const tnsr::I<DataType, Dim, Frame::Inertial>& x, double t,
262 : tmpl::list<Tags...> /*meta*/) const {
263 : const Wave left(left_initial_data_, pressure_star_, velocity_star_,
264 : adiabatic_index_, Side::Left);
265 : const Wave right(right_initial_data_, pressure_star_, velocity_star_,
266 : adiabatic_index_, Side::Right);
267 :
268 : tnsr::I<DataType, Dim, Frame::Inertial> x_shifted(x);
269 : x_shifted.get(propagation_axis_) -= initial_position_;
270 :
271 : return {tuples::get<Tags>(
272 : variables(x_shifted, t, tmpl::list<Tags>{}, left, right))...};
273 : }
274 :
275 0 : const EquationsOfState::IdealFluid<false>& equation_of_state() const {
276 : return equation_of_state_;
277 : }
278 :
279 : // clang-tidy: no runtime references
280 0 : void pup(PUP::er& /*p*/); // NOLINT
281 :
282 : // Retrieve these member variables for testing purposes.
283 0 : constexpr std::array<double, 2> diagnostic_star_region_values() const {
284 : return make_array(pressure_star_, velocity_star_);
285 : }
286 :
287 : private:
288 : /// @{
289 : /// Retrieve hydro variable at `(x, t)`
290 : template <typename DataType>
291 1 : auto variables(const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted,
292 : double t, tmpl::list<Tags::MassDensity<DataType>> /*meta*/,
293 : const Wave& left, const Wave& right) const
294 : -> tuples::TaggedTuple<Tags::MassDensity<DataType>>;
295 :
296 : template <typename DataType>
297 1 : auto variables(const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted,
298 : double t, tmpl::list<Tags::Velocity<DataType, Dim>> /*meta*/,
299 : const Wave& left, const Wave& right) const
300 : -> tuples::TaggedTuple<Tags::Velocity<DataType, Dim>>;
301 :
302 : template <typename DataType>
303 1 : auto variables(const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted,
304 : double t, tmpl::list<Tags::Pressure<DataType>> /*meta*/,
305 : const Wave& left, const Wave& right) const
306 : -> tuples::TaggedTuple<Tags::Pressure<DataType>>;
307 :
308 : template <typename DataType>
309 1 : auto variables(const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted,
310 : double t,
311 : tmpl::list<Tags::SpecificInternalEnergy<DataType>> /*meta*/,
312 : const Wave& left, const Wave& right) const
313 : -> tuples::TaggedTuple<Tags::SpecificInternalEnergy<DataType>>;
314 : /// @}
315 :
316 : // Any of the two waves propagating on each side of the contact discontinuity.
317 : // Depending on whether p_* is larger or smaller
318 : // than the initial pressure on the corresponding side, the wave is a
319 : // shock (larger) or a rarefaction (smaller) wave. Since p_*
320 : // is computed at run time, the characterization of the wave
321 : // must also be done at run time. Shock and rarefaction waves will provide
322 : // different functions to retrieve the primitive variables at a given (x, t).
323 : // Here normal velocity means velocity along the wave propagation.
324 0 : struct Wave {
325 0 : Wave(const InitialData& data, double pressure_star, double velocity_star,
326 : double adiabatic_index, const Side& side);
327 :
328 0 : double mass_density(double x_shifted, double t) const;
329 :
330 0 : double normal_velocity(double x_shifted, double t,
331 : double velocity_star) const;
332 :
333 0 : double pressure(double x_shifted, double t, double pressure_star) const;
334 :
335 : private:
336 : // p_* over initial pressure on the corresponding side
337 0 : double pressure_ratio_ = std::numeric_limits<double>::signaling_NaN();
338 : // false if rarefaction wave
339 0 : bool is_shock_ = true;
340 :
341 0 : InitialData data_{};
342 0 : Shock shock_{};
343 0 : Rarefaction rarefaction_{};
344 : };
345 :
346 0 : struct Shock {
347 0 : Shock(const InitialData& data, double pressure_ratio,
348 : double adiabatic_index, const Side& side);
349 :
350 0 : double mass_density(double x_shifted, double t,
351 : const InitialData& data) const;
352 0 : double normal_velocity(double x_shifted, double t, const InitialData& data,
353 : double velocity_star) const;
354 :
355 0 : double pressure(double x_shifted, double t, const InitialData& data,
356 : double pressure_star) const;
357 :
358 : private:
359 0 : double direction_ = std::numeric_limits<double>::signaling_NaN();
360 0 : double mass_density_star_ = std::numeric_limits<double>::signaling_NaN();
361 0 : double shock_speed_ = std::numeric_limits<double>::signaling_NaN();
362 : };
363 :
364 0 : struct Rarefaction {
365 0 : Rarefaction(const InitialData& data, double pressure_ratio,
366 : double velocity_star, double adiabatic_index, const Side& side);
367 :
368 0 : double mass_density(double x_shifted, double t,
369 : const InitialData& data) const;
370 0 : double normal_velocity(double x_shifted, double t, const InitialData& data,
371 : double velocity_star) const;
372 :
373 0 : double pressure(double x_shifted, double t, const InitialData& data,
374 : double pressure_star) const;
375 :
376 : private:
377 0 : double direction_ = std::numeric_limits<double>::signaling_NaN();
378 0 : double gamma_mm_ = std::numeric_limits<double>::signaling_NaN();
379 0 : double gamma_pp_ = std::numeric_limits<double>::signaling_NaN();
380 0 : double mass_density_star_ = std::numeric_limits<double>::signaling_NaN();
381 0 : double sound_speed_star_ = std::numeric_limits<double>::signaling_NaN();
382 0 : double head_speed_ = std::numeric_limits<double>::signaling_NaN();
383 0 : double tail_speed_ = std::numeric_limits<double>::signaling_NaN();
384 : };
385 :
386 : template <size_t SpatialDim>
387 : friend bool
388 0 : operator==( // NOLINT (clang-tidy: readability-redundant-declaration)
389 : const RiemannProblem<SpatialDim>& lhs,
390 : const RiemannProblem<SpatialDim>& rhs);
391 :
392 0 : double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
393 0 : double initial_position_ = std::numeric_limits<double>::signaling_NaN();
394 0 : size_t propagation_axis_ = PropagationAxis::X;
395 0 : InitialData left_initial_data_{};
396 0 : InitialData right_initial_data_{};
397 :
398 0 : double pressure_star_tol_ = std::numeric_limits<double>::signaling_NaN();
399 : // the pressure in the star region, p_*
400 0 : double pressure_star_ = std::numeric_limits<double>::signaling_NaN();
401 : // the velocity in the star region, u_*
402 0 : double velocity_star_ = std::numeric_limits<double>::signaling_NaN();
403 :
404 0 : EquationsOfState::IdealFluid<false> equation_of_state_{};
405 : };
406 :
407 : template <size_t Dim>
408 0 : bool operator!=(const RiemannProblem<Dim>& lhs, const RiemannProblem<Dim>& rhs);
409 :
410 : } // namespace Solutions
411 : } // namespace NewtonianEuler
|
