ryujin 2.1.1 revision 863a4d36dcc743d4e1a9b41cfabd03d0aea57863
initial_state_shock_front.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2022 - 2024 by the ryujin authors
4//
5
6#pragma once
7
9
10namespace ryujin
11{
12 namespace EulerInitialStates
13 {
26 template <typename Description, int dim, typename Number>
27 class ShockFront : public InitialState<Description, dim, Number>
28 {
29 public:
31 using View =
32 typename Description::template HyperbolicSystemView<dim, Number>;
33 using state_type = typename View::state_type;
34
35 ShockFront(const HyperbolicSystem &hyperbolic_system,
36 const std::string subsection)
37 : InitialState<Description, dim, Number>("shock front", subsection)
38 , hyperbolic_system_(hyperbolic_system)
39 {
40 gamma_ = 1.4;
41 if constexpr (!View::have_gamma) {
42 this->add_parameter("gamma", gamma_, "The ratio of specific heats");
43 }
44
45 primitive_right_[0] = 1.4;
46 primitive_right_[1] = 0.0;
47 primitive_right_[2] = 1.;
48 this->add_parameter("primitive state",
49 primitive_right_,
50 "Initial 1d primitive state (rho, u, p) before the "
51 "shock (to the right)");
52
53 mach_number_ = 2.0;
54 this->add_parameter(
55 "mach number",
56 mach_number_,
57 "Mach number of shock front (S1, S3 = mach * a_L/R)");
58
59 const auto compute_and_convert_states = [&]() {
60 const auto view = hyperbolic_system_.template view<dim, Number>();
61
62 if constexpr (View::have_gamma) {
63 gamma_ = view.gamma();
64 }
65
66 /* Compute post-shock state and S3: */
67
68 auto b = Number(0.);
69 if constexpr (View::have_eos_interpolation_b)
70 b = view.eos_interpolation_b();
71
72 const auto &rho_R = primitive_right_[0];
73 const auto &u_R = primitive_right_[1];
74 const auto &p_R = primitive_right_[2];
75 /* a_R^2 = gamma * p / rho / (1 - b * rho) */
76 const Number a_R = std::sqrt(gamma_ * p_R / rho_R / (1 - b * rho_R));
77 const Number mach_R = u_R / a_R;
78
79 S3_ = mach_number_ * a_R;
80 const Number delta_mach = mach_R - mach_number_;
81
82 const Number rho_L =
83 rho_R * (gamma_ + Number(1.)) * delta_mach * delta_mach /
84 ((gamma_ - Number(1.)) * delta_mach * delta_mach + Number(2.));
85 const Number u_L =
86 (Number(1.) - rho_R / rho_L) * S3_ + rho_R / rho_L * u_R;
87 const Number p_L = p_R *
88 (Number(2.) * gamma_ * delta_mach * delta_mach -
89 (gamma_ - Number(1.))) /
90 (gamma_ + Number(1.));
91
92 primitive_left_[0] = rho_L;
93 primitive_left_[1] = u_L;
94 primitive_left_[2] = p_L;
95
96 state_left_ = view.from_initial_state(primitive_left_);
97 state_right_ = view.from_initial_state(primitive_right_);
98 };
99
100 this->parse_parameters_call_back.connect(compute_and_convert_states);
101 compute_and_convert_states();
102 }
103
104
105 state_type compute(const dealii::Point<dim> &point, Number t) final
106 {
107 const Number position_1d = Number(point[0] - S3_ * t);
108 return (position_1d > 0. ? state_right_ : state_left_);
109 }
110
111 private:
112 const HyperbolicSystem &hyperbolic_system_;
113
114 Number gamma_;
115
116 dealii::Tensor<1, 3, Number> primitive_left_;
117 dealii::Tensor<1, 3, Number> primitive_right_;
118 Number mach_number_;
119 Number S3_;
120
121 state_type state_left_;
122 state_type state_right_;
123 };
124 } // namespace EulerInitialStates
125} // namespace ryujin
typename Description::HyperbolicSystem HyperbolicSystem
typename Description::template HyperbolicSystemView< dim, Number > View
state_type compute(const dealii::Point< dim > &point, Number t) final
ShockFront(const HyperbolicSystem &hyperbolic_system, const std::string subsection)
Euler::HyperbolicSystem HyperbolicSystem
Definition: description.h:32