8#include <compile_time_options.h>
16 namespace EulerInitialStates
30 template <
typename Description,
int dim,
typename Number>
42 const std::string subsection)
44 , hyperbolic_system_(hyperbolic_system)
55 using state_type_1d = std::array<Number, 3>;
56 constexpr state_type_1d primitive_left{1., 0., Number(2. / 3. * 1.e-1)};
57 constexpr state_type_1d primitive_right{
58 1.e-3, 0., Number(2. / 3. * 1.e-10)};
61 constexpr Number rarefaction_speed = 0.49578489518897934;
62 constexpr Number contact_velocity = 0.62183867139173454;
63 constexpr Number right_shock_speed = 0.82911836253346982;
69 constexpr Number pre_contact_density = 5.4079335349316249e-02;
70 constexpr Number post_contact_density = 3.9999980604299963e-03;
71 constexpr Number contact_pressure = 0.51557792765096996e-03;
73 state_type_1d primitive_state;
74 const double &x = point[0];
76 if (x <= -1.0 / 3.0 * t) {
78 primitive_state = primitive_left;
80 }
else if (x < rarefaction_speed * t) {
82 const double chi = x / t;
83 primitive_state[0] =
std::pow(0.75 - 0.75 * chi, 3.0);
84 primitive_state[1] = 0.75 * (1.0 / 3.0 + chi);
85 primitive_state[2] = (1.0 / 15.0) *
std::pow(0.75 - 0.75 * chi, 5.0);
87 }
else if (x < contact_velocity * t) {
88 primitive_state[0] = pre_contact_density;
89 primitive_state[1] = contact_velocity;
90 primitive_state[2] = contact_pressure;
92 }
else if (x < right_shock_speed * t) {
94 primitive_state[0] = post_contact_density;
95 primitive_state[1] = contact_velocity;
96 primitive_state[2] = contact_pressure;
100 primitive_state = primitive_right;
105 const auto &[rho, u, p] = primitive_state;
106 conserved_state[0] = rho;
107 conserved_state[1] = rho * u;
108 conserved_state[dim + 1] =
112 return conserved_state;
typename Description::template HyperbolicSystemView< dim, Number > View
LeBlanc(const HyperbolicSystem &hyperbolic_system, const std::string subsection)
typename View::state_type state_type
typename View::ScalarNumber ScalarNumber
typename Description::HyperbolicSystem HyperbolicSystem
state_type compute(const dealii::Point< dim > &point, Number t) final
T pow(const T x, const T b)
Euler::HyperbolicSystem HyperbolicSystem