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