ryujin 2.1.1 revision ef0fcd4010d109b860652ece4a7b8963fb7d46b1
initial_state_rarefaction.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2023 - 2024 by the ryujin authors
4//
5
6#pragma once
7
8#include <compile_time_options.h>
9
11
12namespace ryujin
13{
14 namespace EulerInitialStates
15 {
28 template <typename Description, int dim, typename Number>
29 class Rarefaction : public InitialState<Description, dim, Number>
30 {
31 public:
33 using View =
34 typename Description::template HyperbolicSystemView<dim, Number>;
35 using state_type = typename View::state_type;
36
37 using state_type_1d = std::array<Number, 4>;
38
39 Rarefaction(const HyperbolicSystem &hyperbolic_system,
40 const std::string subsection)
41 : InitialState<Description, dim, Number>("rarefaction", subsection)
42 , hyperbolic_system_(hyperbolic_system)
43 {
44 gamma_ = 1.4;
45 if constexpr (!View::have_gamma) {
46 this->add_parameter("gamma", gamma_, "The ratio of specific heats");
47 }
48
49 /*
50 * Compute the speed of sound:
51 */
52 const auto speed_of_sound = [&](const Number rho, const Number p) {
53 return std::sqrt(gamma_ * p / rho);
54 };
55
56 /*
57 * Compute the rarefaction right side:
58 */
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;
63 state_type_1d primitive_right{{rho_right, 0., 0.}};
64
65 /* Isentropic condition: pR = (rhoR/rhoL)^{gamma} * pL */
66 primitive_right[2] = std::pow(rho_right / rho_left, gamma_) * p_left;
67
68 const auto c_right = speed_of_sound(rho_right, primitive_right[2]);
69 primitive_right[3] = c_right;
70
71 /* 1-Riemann invariant: uR + 2 cR/(gamma -1) = uL + 2 cL/(gamma -1) */
72 primitive_right[1] =
73 u_left + 2.0 * (c_left - c_right) / (gamma_ - 1.0);
74
75 return primitive_right;
76 };
77
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();
83 }
84
85 /*
86 * Initial left and right states (rho, u, p, c):
87 */
88
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; /* verify */
93 const Number rho_right = 0.5;
94
95 primitive_left_ = {rho_left, c_left, p_left, c_left};
96 primitive_right_ =
97 rarefaction_right_state(primitive_left_, rho_right);
98
99 /*
100 * Populate constants:
101 */
102
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);
108 };
109
110 this->parse_parameters_call_back.connect(compute_constants);
111 compute_constants();
112 } /* Constructor */
113
114 state_type compute(const dealii::Point<dim> &point, Number delta_t) final
115 {
116 /*
117 * Compute rarefaction solution:
118 */
119
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_;
122
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;
126
127 state_type_1d primitive;
128
129 if (x <= t * (u_left - c_left)) {
130 primitive = primitive_left_;
131
132 } else if (x <= t * (u_right - c_right)) {
133
134 /* Self-similar variable: */
135 const double chi = x / t;
136
137 primitive[0] =
138 rho_left * std::pow(k1 + k2 * (u_left - chi), density_exponent);
139 primitive[1] = k1 * (k3 + chi);
140 primitive[2] =
141 p_left * std::pow(k1 + k2 * (u_left - chi), pressure_exponent);
142
143 } else {
144 primitive = primitive_right_;
145 }
146
147 state_type conserved_state;
148 {
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;
154 }
155 return conserved_state;
156 }
157
158 private:
159 const HyperbolicSystem &hyperbolic_system_;
160 Number gamma_;
161
162 state_type_1d primitive_left_;
163 state_type_1d primitive_right_;
164 Number k1;
165 Number k2;
166 Number density_exponent;
167 Number k3;
168 Number pressure_exponent;
169 };
170 } // namespace EulerInitialStates
171} // namespace ryujin
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
T pow(const T x, const T b)
Euler::HyperbolicSystem HyperbolicSystem
Definition: description.h:34