8#include <compile_time_options.h>
14 namespace EulerInitialStates
28 template <
typename Description,
int dim,
typename Number>
40 const std::string subsection)
42 , hyperbolic_system_(hyperbolic_system)
45 if constexpr (!View::have_gamma) {
46 this->add_parameter(
"gamma", gamma_,
"The ratio of specific heats");
52 const auto speed_of_sound = [&](
const Number rho,
const Number p) {
53 return std::sqrt(gamma_ * p / rho);
59 const auto rarefaction_right_state = [
this, speed_of_sound](
60 const auto primitive_left,
61 const Number rho_right) {
62 const auto &[rho_left, u_left, p_left, c_left] = primitive_left;
66 primitive_right[2] =
std::pow(rho_right / rho_left, gamma_) * p_left;
68 const auto c_right = speed_of_sound(rho_right, primitive_right[2]);
69 primitive_right[3] = c_right;
73 u_left + 2.0 * (c_left - c_right) / (gamma_ - 1.0);
75 return primitive_right;
78 const auto compute_constants =
79 [
this, speed_of_sound, rarefaction_right_state]() {
80 const auto view = hyperbolic_system_.template view<dim, Number>();
81 if constexpr (View::have_gamma) {
82 gamma_ = view.gamma();
89 const Number rho_left = 3.0;
90 const Number p_left = 1.0;
91 const Number c_left = speed_of_sound(rho_left, p_left);
92 const Number u_left = c_left;
93 const Number rho_right = 0.5;
95 primitive_left_ = {rho_left, c_left, p_left, c_left};
97 rarefaction_right_state(primitive_left_, rho_right);
103 k1 = 2.0 / (gamma_ + 1.0);
104 k2 = ((gamma_ - 1.0) / ((gamma_ + 1.0) * c_left));
105 density_exponent = 2.0 / (gamma_ - 1.0);
106 k3 = c_left + ((gamma_ - 1.0) / 2.0) * u_left;
107 pressure_exponent = 2.0 * gamma_ / (gamma_ - 1.0);
110 this->parse_parameters_call_back.connect(compute_constants);
120 const auto &[rho_left, u_left, p_left, c_left] = primitive_left_;
121 const auto &[rho_right, u_right, p_right, c_right] = primitive_right_;
123 const double x = point[0];
124 const auto t_0 = 0.2 / (u_right - u_left);
125 const auto t = t_0 + delta_t;
129 if (x <= t * (u_left - c_left)) {
130 primitive = primitive_left_;
132 }
else if (x <= t * (u_right - c_right)) {
135 const double chi = x / t;
138 rho_left *
std::pow(k1 + k2 * (u_left - chi), density_exponent);
139 primitive[1] = k1 * (k3 + chi);
141 p_left *
std::pow(k1 + k2 * (u_left - chi), pressure_exponent);
144 primitive = primitive_right_;
149 const auto &[rho, u, p, c] = primitive;
150 conserved_state[0] = rho;
151 conserved_state[1] = rho * u;
152 conserved_state[dim + 1] =
153 p / Number(gamma_ - 1.) + Number(0.5) * rho * u * u;
155 return conserved_state;
166 Number density_exponent;
168 Number pressure_exponent;
std::array< Number, 4 > state_type_1d
Rarefaction(const HyperbolicSystem &hyperbolic_system, const std::string subsection)
state_type compute(const dealii::Point< dim > &point, Number delta_t) final
typename Description::template HyperbolicSystemView< dim, Number > View
typename Description::HyperbolicSystem HyperbolicSystem
typename View::state_type state_type
T pow(const T x, const T b)
Euler::HyperbolicSystem HyperbolicSystem