SpECTRE Documentation Coverage Report
 Current view: top level - PointwiseFunctions/AnalyticSolutions/NewtonianEuler - RiemannProblem.hpp Hit Total Coverage Commit: 9f349d3c09e1c03107f00c2135ca40e209d3b84c Lines: 16 109 14.7 % Date: 2023-06-09 21:05:06 Legend: Lines: hit not hit
  Line data Source code  1 0 : // Distributed under the MIT License. 2 : // See LICENSE.txt for details. 3 : 4 : #pragma once 5 : 6 : #include 7 : #include 8 : #include 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 structs 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 structs 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 Waves, 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 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; 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; 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; 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& 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::signaling_NaN(); 219 0 : std::array velocity_ = 220 : make_array(std::numeric_limits::signaling_NaN()); 221 0 : double pressure_ = std::numeric_limits::signaling_NaN(); 222 0 : double sound_speed_ = std::numeric_limits::signaling_NaN(); 223 0 : double normal_velocity_ = std::numeric_limits::signaling_NaN(); 224 : 225 : // Data-dependent constants A and B in Eqns. (4.8) of Toro. 226 0 : double constant_a_ = std::numeric_limits::signaling_NaN(); 227 0 : double constant_b_ = std::numeric_limits::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; 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& left_velocity, 252 : double left_pressure, double right_mass_density, 253 : const std::array& 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 260 1 : tuples::TaggedTuple variables( 261 : const tnsr::I& x, double t, 262 : tmpl::list /*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 x_shifted(x); 269 : x_shifted.get(propagation_axis_) -= initial_position_; 270 : 271 : return {tuples::get( 272 : variables(x_shifted, t, tmpl::list{}, left, right))...}; 273 : } 274 : 275 0 : const EquationsOfState::IdealFluid& 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 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 291 1 : auto variables(const tnsr::I& x_shifted, 292 : double t, tmpl::list> /*meta*/, 293 : const Wave& left, const Wave& right) const 294 : -> tuples::TaggedTuple>; 295 : 296 : template 297 1 : auto variables(const tnsr::I& x_shifted, 298 : double t, tmpl::list> /*meta*/, 299 : const Wave& left, const Wave& right) const 300 : -> tuples::TaggedTuple>; 301 : 302 : template 303 1 : auto variables(const tnsr::I& x_shifted, 304 : double t, tmpl::list> /*meta*/, 305 : const Wave& left, const Wave& right) const 306 : -> tuples::TaggedTuple>; 307 : 308 : template 309 1 : auto variables(const tnsr::I& x_shifted, 310 : double t, 311 : tmpl::list> /*meta*/, 312 : const Wave& left, const Wave& right) const 313 : -> tuples::TaggedTuple>; 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::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::signaling_NaN(); 360 0 : double mass_density_star_ = std::numeric_limits::signaling_NaN(); 361 0 : double shock_speed_ = std::numeric_limits::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::signaling_NaN(); 378 0 : double gamma_mm_ = std::numeric_limits::signaling_NaN(); 379 0 : double gamma_pp_ = std::numeric_limits::signaling_NaN(); 380 0 : double mass_density_star_ = std::numeric_limits::signaling_NaN(); 381 0 : double sound_speed_star_ = std::numeric_limits::signaling_NaN(); 382 0 : double head_speed_ = std::numeric_limits::signaling_NaN(); 383 0 : double tail_speed_ = std::numeric_limits::signaling_NaN(); 384 : }; 385 : 386 : template 387 : friend bool 388 0 : operator==( // NOLINT (clang-tidy: readability-redundant-declaration) 389 : const RiemannProblem& lhs, 390 : const RiemannProblem& rhs); 391 : 392 0 : double adiabatic_index_ = std::numeric_limits::signaling_NaN(); 393 0 : double initial_position_ = std::numeric_limits::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::signaling_NaN(); 399 : // the pressure in the star region, p_* 400 0 : double pressure_star_ = std::numeric_limits::signaling_NaN(); 401 : // the velocity in the star region, u_* 402 0 : double velocity_star_ = std::numeric_limits::signaling_NaN(); 403 : 404 0 : EquationsOfState::IdealFluid equation_of_state_{}; 405 : }; 406 : 407 : template 408 0 : bool operator!=(const RiemannProblem& lhs, const RiemannProblem& rhs); 409 : 410 : } // namespace Solutions 411 : } // namespace NewtonianEuler 

 Generated by: LCOV version 1.14