ryujin 2.1.1 revision fa044d3939298da517bf6a5839b34e00b155547d
initial_state_leblanc.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
12#include <cmath>
13
14namespace ryujin
15{
16 namespace EulerInitialStates
17 {
30 template <typename Description, int dim, typename Number>
31 class LeBlanc : public InitialState<Description, dim, Number>
32 {
33 public:
35 using View =
36 typename Description::template HyperbolicSystemView<dim, Number>;
37 using state_type = typename View::state_type;
38
39 using ScalarNumber = typename View::ScalarNumber;
40
41 LeBlanc(const HyperbolicSystem &hyperbolic_system,
42 const std::string subsection)
43 : InitialState<Description, dim, Number>("leblanc", subsection)
44 , hyperbolic_system_(hyperbolic_system)
45 {
46 } /* Constructor */
47
48 state_type compute(const dealii::Point<dim> &point, Number t) final
49 {
50 /*
51 * The Le Blanc shock tube:
52 */
53
54 /* Initial left and right states (rho, u, p): */
55 using state_type_1d = std::array<Number, 3>;
56 constexpr state_type_1d primitive_left{1., 0., Number(2. / 3. * 1.e-1)};
57 constexpr state_type_1d primitive_right{
58 1.e-3, 0., Number(2. / 3. * 1.e-10)};
59
60 /* The intermediate wave-speeds appearing on the Riemann fan: */
61 constexpr Number rarefaction_speed = 0.49578489518897934;
62 constexpr Number contact_velocity = 0.62183867139173454;
63 constexpr Number right_shock_speed = 0.82911836253346982;
64
65 /*
66 * Velocity and pressure are constant across the middle discontinuity,
67 * only the density jumps: it's a contact wave!
68 */
69 constexpr Number pre_contact_density = 5.4079335349316249e-02;
70 constexpr Number post_contact_density = 3.9999980604299963e-03;
71 constexpr Number contact_pressure = 0.51557792765096996e-03;
72
73 state_type_1d primitive_state;
74 const double &x = point[0];
75
76 if (x <= -1.0 / 3.0 * t) {
77 /* Left state: */
78 primitive_state = primitive_left;
79
80 } else if (x < rarefaction_speed * t) {
81 /* Expansion data (with self-similar variable chi): */
82 const double chi = x / t;
83 primitive_state[0] = std::pow(0.75 - 0.75 * chi, 3.0);
84 primitive_state[1] = 0.75 * (1.0 / 3.0 + chi);
85 primitive_state[2] = (1.0 / 15.0) * std::pow(0.75 - 0.75 * chi, 5.0);
86
87 } else if (x < contact_velocity * t) {
88 primitive_state[0] = pre_contact_density;
89 primitive_state[1] = contact_velocity;
90 primitive_state[2] = contact_pressure;
91
92 } else if (x < right_shock_speed * t) {
93 /* Contact-wave data (velocity and pressure are continuous): */
94 primitive_state[0] = post_contact_density;
95 primitive_state[1] = contact_velocity;
96 primitive_state[2] = contact_pressure;
97
98 } else {
99 /* Right state: */
100 primitive_state = primitive_right;
101 }
102
103 state_type conserved_state;
104 {
105 const auto &[rho, u, p] = primitive_state;
106 conserved_state[0] = rho;
107 conserved_state[1] = rho * u;
108 conserved_state[dim + 1] =
109 p / ScalarNumber(5. / 3. - 1.) + ScalarNumber(0.5) * rho * u * u;
110 }
111
112 return conserved_state;
113 }
114
115 private:
116 const HyperbolicSystem &hyperbolic_system_;
117 };
118 } // namespace EulerInitialStates
119} // namespace ryujin
typename Description::template HyperbolicSystemView< dim, Number > View
LeBlanc(const HyperbolicSystem &hyperbolic_system, const std::string subsection)
typename View::ScalarNumber ScalarNumber
typename Description::HyperbolicSystem HyperbolicSystem
state_type compute(const dealii::Point< dim > &point, Number t) final
T pow(const T x, const T b)
Euler::HyperbolicSystem HyperbolicSystem
Definition: description.h:34