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/TaggedTuple.hpp"
12 : #include "DataStructures/Tensor/Tensor.hpp"
13 : #include "DataStructures/Tensor/TypeAliases.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 "PointwiseFunctions/Hydro/Tags.hpp"
20 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
21 : #include "Utilities/MakeArray.hpp"
22 : #include "Utilities/TMPL.hpp"
23 :
24 : /// \cond
25 : namespace PUP {
26 : class er;
27 : } // namespace PUP
28 : /// \endcond
29 :
30 : namespace NewtonianEuler::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 evolution::initial_data::InitialData,
130 : public MarkAsAnalyticSolution {
131 0 : enum class Side { Left, Right };
132 0 : enum PropagationAxis { X = 0, Y = 1, Z = 2 };
133 :
134 : struct Wave;
135 : struct Shock;
136 : struct Rarefaction;
137 :
138 : public:
139 0 : using equation_of_state_type = EquationsOfState::IdealFluid<false>;
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() override = default;
248 :
249 0 : auto get_clone() const
250 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
251 :
252 : /// \cond
253 : explicit RiemannProblem(CkMigrateMessage* msg);
254 : using PUP::able::register_constructor;
255 : WRAPPED_PUPable_decl_template(RiemannProblem);
256 : /// \endcond
257 :
258 0 : RiemannProblem(double adiabatic_index, double initial_position,
259 : double left_mass_density,
260 : const std::array<double, Dim>& left_velocity,
261 : double left_pressure, double right_mass_density,
262 : const std::array<double, Dim>& right_velocity,
263 : double right_pressure,
264 : double pressure_star_tol = PressureStarTol::suggested_value());
265 :
266 : /// Retrieve a collection of hydrodynamic variables at position `x`
267 : /// and time `t`
268 : template <typename DataType, typename... Tags>
269 1 : tuples::TaggedTuple<Tags...> variables(
270 : const tnsr::I<DataType, Dim, Frame::Inertial>& x, double t,
271 : tmpl::list<Tags...> /*meta*/) const {
272 : const Wave left(left_initial_data_, pressure_star_, velocity_star_,
273 : adiabatic_index_, Side::Left);
274 : const Wave right(right_initial_data_, pressure_star_, velocity_star_,
275 : adiabatic_index_, Side::Right);
276 :
277 : tnsr::I<DataType, Dim, Frame::Inertial> x_shifted(x);
278 : x_shifted.get(propagation_axis_) -= initial_position_;
279 :
280 : return {tuples::get<Tags>(
281 : variables(x_shifted, t, tmpl::list<Tags>{}, left, right))...};
282 : }
283 :
284 0 : const EquationsOfState::IdealFluid<false>& equation_of_state() const {
285 : return equation_of_state_;
286 : }
287 :
288 : // NOLINTNEXTLINE(google-runtime-references)
289 0 : void pup(PUP::er& /*p*/) override;
290 :
291 : // Retrieve these member variables for testing purposes.
292 0 : constexpr std::array<double, 2> diagnostic_star_region_values() const {
293 : return make_array(pressure_star_, velocity_star_);
294 : }
295 :
296 : private:
297 : /// @{
298 : /// Retrieve hydro variable at `(x, t)`
299 : template <typename DataType>
300 1 : auto variables(const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted,
301 : double t,
302 : tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/,
303 : const Wave& left, const Wave& right) const
304 : -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
305 :
306 : template <typename DataType>
307 1 : auto variables(
308 : const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted, double t,
309 : tmpl::list<hydro::Tags::SpatialVelocity<DataType, Dim>> /*meta*/,
310 : const Wave& left, const Wave& right) const
311 : -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, Dim>>;
312 :
313 : template <typename DataType>
314 1 : auto variables(const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted,
315 : double t, tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/,
316 : const Wave& left, const Wave& right) const
317 : -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
318 :
319 : template <typename DataType>
320 1 : auto variables(
321 : const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted, double t,
322 : tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/,
323 : const Wave& left, const Wave& right) const
324 : -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
325 : /// @}
326 :
327 : // Any of the two waves propagating on each side of the contact discontinuity.
328 : // Depending on whether p_* is larger or smaller
329 : // than the initial pressure on the corresponding side, the wave is a
330 : // shock (larger) or a rarefaction (smaller) wave. Since p_*
331 : // is computed at run time, the characterization of the wave
332 : // must also be done at run time. Shock and rarefaction waves will provide
333 : // different functions to retrieve the primitive variables at a given (x, t).
334 : // Here normal velocity means velocity along the wave propagation.
335 0 : struct Wave {
336 0 : Wave(const InitialData& data, double pressure_star, double velocity_star,
337 : double adiabatic_index, const Side& side);
338 :
339 0 : double mass_density(double x_shifted, double t) const;
340 :
341 0 : double normal_velocity(double x_shifted, double t,
342 : double velocity_star) const;
343 :
344 0 : double pressure(double x_shifted, double t, double pressure_star) const;
345 :
346 : private:
347 : // p_* over initial pressure on the corresponding side
348 0 : double pressure_ratio_ = std::numeric_limits<double>::signaling_NaN();
349 : // false if rarefaction wave
350 0 : bool is_shock_ = true;
351 :
352 0 : InitialData data_{};
353 0 : Shock shock_{};
354 0 : Rarefaction rarefaction_{};
355 : };
356 :
357 0 : struct Shock {
358 0 : Shock(const InitialData& data, double pressure_ratio,
359 : double adiabatic_index, const Side& side);
360 :
361 0 : double mass_density(double x_shifted, double t,
362 : const InitialData& data) const;
363 0 : double normal_velocity(double x_shifted, double t, const InitialData& data,
364 : double velocity_star) const;
365 :
366 0 : double pressure(double x_shifted, double t, const InitialData& data,
367 : double pressure_star) const;
368 :
369 : private:
370 0 : double direction_ = std::numeric_limits<double>::signaling_NaN();
371 0 : double mass_density_star_ = std::numeric_limits<double>::signaling_NaN();
372 0 : double shock_speed_ = std::numeric_limits<double>::signaling_NaN();
373 : };
374 :
375 0 : struct Rarefaction {
376 0 : Rarefaction(const InitialData& data, double pressure_ratio,
377 : double velocity_star, double adiabatic_index, const Side& side);
378 :
379 0 : double mass_density(double x_shifted, double t,
380 : const InitialData& data) const;
381 0 : double normal_velocity(double x_shifted, double t, const InitialData& data,
382 : double velocity_star) const;
383 :
384 0 : double pressure(double x_shifted, double t, const InitialData& data,
385 : double pressure_star) const;
386 :
387 : private:
388 0 : double direction_ = std::numeric_limits<double>::signaling_NaN();
389 0 : double gamma_mm_ = std::numeric_limits<double>::signaling_NaN();
390 0 : double gamma_pp_ = std::numeric_limits<double>::signaling_NaN();
391 0 : double mass_density_star_ = std::numeric_limits<double>::signaling_NaN();
392 0 : double sound_speed_star_ = std::numeric_limits<double>::signaling_NaN();
393 0 : double head_speed_ = std::numeric_limits<double>::signaling_NaN();
394 0 : double tail_speed_ = std::numeric_limits<double>::signaling_NaN();
395 : };
396 :
397 : template <size_t SpatialDim>
398 : friend bool
399 0 : operator==( // NOLINT (clang-tidy: readability-redundant-declaration)
400 : const RiemannProblem<SpatialDim>& lhs,
401 : const RiemannProblem<SpatialDim>& rhs);
402 :
403 0 : double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
404 0 : double initial_position_ = std::numeric_limits<double>::signaling_NaN();
405 0 : size_t propagation_axis_ = PropagationAxis::X;
406 0 : InitialData left_initial_data_{};
407 0 : InitialData right_initial_data_{};
408 :
409 0 : double pressure_star_tol_ = std::numeric_limits<double>::signaling_NaN();
410 : // the pressure in the star region, p_*
411 0 : double pressure_star_ = std::numeric_limits<double>::signaling_NaN();
412 : // the velocity in the star region, u_*
413 0 : double velocity_star_ = std::numeric_limits<double>::signaling_NaN();
414 :
415 0 : EquationsOfState::IdealFluid<false> equation_of_state_{};
416 : };
417 :
418 : template <size_t Dim>
419 0 : bool operator!=(const RiemannProblem<Dim>& lhs, const RiemannProblem<Dim>& rhs);
420 :
421 : } // namespace NewtonianEuler::Solutions
|
