ryujin 2.1.1 revision 46bf70e400e423a8ffffe8300887eeb35b8dfb2c
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
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 ShockFront : 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 ShockFront(const HyperbolicSystem &hyperbolic_system,
38 const std::string subsection)
39 : InitialState<Description, dim, Number>("shock front", subsection)
40 , hyperbolic_system_(hyperbolic_system)
41 {
42 gamma_ = 1.4;
43 if constexpr (!View::have_gamma) {
44 this->add_parameter("gamma", gamma_, "The ratio of specific heats");
45 }
46
47 primitive_right_[0] = 1.4;
48 primitive_right_[1] = 0.0;
49 primitive_right_[2] = 1.;
50 this->add_parameter("primitive state",
51 primitive_right_,
52 "Initial 1d primitive state (rho, u, p) before the "
53 "shock (to the right)");
54
55 mach_number_ = 2.0;
56 this->add_parameter(
57 "mach number",
58 mach_number_,
59 "Mach number of shock front (S1, S3 = mach * a_L/R)");
60
61 const auto compute_and_convert_states = [&]() {
62 const auto view = hyperbolic_system_.template view<dim, Number>();
63
64 if constexpr (View::have_gamma) {
65 gamma_ = view.gamma();
66 }
67
68 /* Compute post-shock state and S3: */
69
70 auto b = Number(0.);
71 if constexpr (View::have_eos_interpolation_b)
72 b = view.eos_interpolation_b();
73
74 const auto &rho_R = primitive_right_[0];
75 const auto &u_R = primitive_right_[1];
76 const auto &p_R = primitive_right_[2];
77 /* a_R^2 = gamma * p / rho / (1 - b * rho) */
78 const Number a_R = std::sqrt(gamma_ * p_R / rho_R / (1 - b * rho_R));
79 const Number mach_R = u_R / a_R;
80
81 S3_ = mach_number_ * a_R;
82 const Number delta_mach = mach_R - mach_number_;
83
84 const Number rho_L =
85 rho_R * (gamma_ + Number(1.)) * delta_mach * delta_mach /
86 ((gamma_ - Number(1.)) * delta_mach * delta_mach + Number(2.));
87 const Number u_L =
88 (Number(1.) - rho_R / rho_L) * S3_ + rho_R / rho_L * u_R;
89 const Number p_L = p_R *
90 (Number(2.) * gamma_ * delta_mach * delta_mach -
91 (gamma_ - Number(1.))) /
92 (gamma_ + Number(1.));
93
94 primitive_left_[0] = rho_L;
95 primitive_left_[1] = u_L;
96 primitive_left_[2] = p_L;
97
98 state_left_ = view.from_initial_state(primitive_left_);
99 state_right_ = view.from_initial_state(primitive_right_);
100 };
101
102 this->parse_parameters_call_back.connect(compute_and_convert_states);
103 compute_and_convert_states();
104 }
105
106
107 state_type compute(const dealii::Point<dim> &point, Number t) final
108 {
109 const Number position_1d = Number(point[0] - S3_ * t);
110 return (position_1d > 0. ? state_right_ : state_left_);
111 }
112
113 private:
114 const HyperbolicSystem &hyperbolic_system_;
115
116 Number gamma_;
117
118 dealii::Tensor<1, 3, Number> primitive_left_;
119 dealii::Tensor<1, 3, Number> primitive_right_;
120 Number mach_number_;
121 Number S3_;
122
123 state_type state_left_;
124 state_type state_right_;
125 };
126 } // namespace EulerInitialStates
127} // 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:34