ryujin 2.1.1 revision 3913ebe5e1d74c8fef21bcd286906a244031ce08
initial_state_isentropic_vortex.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#include <simd.h>
12
13namespace ryujin
14{
15 namespace EulerInitialStates
16 {
27 template <typename Description, int dim, typename Number>
28 class IsentropicVortex : public InitialState<Description, dim, Number>
29 {
30 public:
32 using View =
33 typename Description::template HyperbolicSystemView<dim, Number>;
34 using state_type = typename View::state_type;
35
36 IsentropicVortex(const HyperbolicSystem &hyperbolic_system,
37 const std::string subsection)
38 : InitialState<Description, dim, Number>("isentropic vortex",
39 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 mach_number_ = 2.0;
48 this->add_parameter(
49 "mach number", mach_number_, "Mach number of isentropic vortex");
50
51 beta_ = 5.0;
52 this->add_parameter("beta", beta_, "vortex strength beta");
53 }
54
55 state_type compute(const dealii::Point<dim> &point, Number t) final
56 {
57 const auto view = hyperbolic_system_.template view<dim, Number>();
58
59 if constexpr (View::have_gamma) {
60 gamma_ = view.gamma();
61 }
62
63 /* In 3D we simply project onto the 2d plane: */
64 dealii::Point<2> point_bar;
65 point_bar[0] = point[0] - mach_number_ * t;
66 point_bar[1] = point[1];
67
68 const Number r_square = Number(point_bar.norm_square());
69
70 const Number factor = beta_ / Number(2. * M_PI) *
71 exp(Number(0.5) - Number(0.5) * r_square);
72
73 const Number T = Number(1.) - (gamma_ - Number(1.)) /
74 (Number(2.) * gamma_) * factor *
75 factor;
76
77 const Number u = mach_number_ - factor * Number(point_bar[1]);
78 const Number v = factor * Number(point_bar[0]);
79
80 const Number rho = ryujin::pow(T, Number(1.) / (Number(gamma_ - 1.)));
81 const Number p = ryujin::pow(rho, Number(gamma_));
82 const Number E =
83 p / (gamma_ - Number(1.)) + Number(0.5) * rho * (u * u + v * v);
84
85 if constexpr (dim == 2)
86 return state_type({rho, rho * u, rho * v, E});
87 else if constexpr (dim == 3)
88 return state_type({rho, rho * u, rho * v, Number(0.), E});
89 else {
90 AssertThrow(false, dealii::ExcNotImplemented());
91 __builtin_trap();
92 }
93 }
94
95 private:
96 const HyperbolicSystem &hyperbolic_system_;
97
98 Number gamma_;
99 Number mach_number_;
100 Number beta_;
101 };
102 } // namespace EulerInitialStates
103} // namespace ryujin
IsentropicVortex(const HyperbolicSystem &hyperbolic_system, const std::string subsection)
typename Description::HyperbolicSystem HyperbolicSystem
state_type compute(const dealii::Point< dim > &point, Number t) final
typename Description::template HyperbolicSystemView< dim, Number > View
T pow(const T x, const T b)
Euler::HyperbolicSystem HyperbolicSystem
Definition: description.h:34