ryujin 2.1.1 revision 3913ebe5e1d74c8fef21bcd286906a244031ce08
initial_state_noh.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
12namespace ryujin
13{
14 namespace EulerInitialStates
15 {
27 template <typename Description, int dim, typename Number>
28 class Noh : 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 Noh(const HyperbolicSystem &hyperbolic_system,
37 const std::string &subsection)
38 : InitialState<Description, dim, Number>("noh", subsection)
39 , hyperbolic_system_(hyperbolic_system)
40 {
41 gamma_ = 1.4;
42 if constexpr (!View::have_gamma) {
43 this->add_parameter("gamma", gamma_, "The ratio of specific heats");
44 }
45
46 rho0_ = 1.0;
47 this->add_parameter(
48 "reference density", rho0_, "The reference density");
49
50 /*
51 * Exact solution assumes this value is negative, but we are just
52 * switching u0 to -u0 by hand in the formulas.
53 */
54 u0_ = 1.0;
55 this->add_parameter("reference velocity magnitude",
56 u0_,
57 "The reference velocity magnitude");
58
59
60 p0_ = 1.e-12;
61 this->add_parameter(
62 "reference pressure", p0_, "The reference pressure");
63
64 this->parse_parameters_call_back.connect([&]() {
65 if constexpr (View::have_gamma) {
66 const auto view = hyperbolic_system_.template view<dim, Number>();
67 gamma_ = view.gamma();
68 }
69 });
70 }
71
72 /* Compute solution */
73 auto compute(const dealii::Point<dim> &point, Number t)
74 -> state_type final
75 {
76 const auto view = hyperbolic_system_.template view<dim, Number>();
77
78 const auto norm = point.norm();
79 const auto min = 10. * std::numeric_limits<Number>::min();
80
81 /* Initialize primitive variables */
82 Number rho = rho0_;
83 auto vel = -u0_ * point / (norm + min);
84 Number p = p0_;
85
86 /* Define exact solutions */
87 const auto D = u0_ * (gamma_ - 1.) / 2.;
88 const bool in_interior = t == Number(0.) ? false : norm / t < D;
89
90 if (in_interior) {
91 rho = rho0_ * std::pow((gamma_ + 1.) / (gamma_ - 1.), dim);
92 vel = 0. * point;
93 p = 0.5 * rho0_ * u0_ * u0_;
94 p *= std::pow(gamma_ + 1., dim) / std::pow(gamma_ - 1., dim - 1);
95 } else {
96 rho = rho0_ * std::pow(1. + t / (norm + min), dim - 1);
97 }
98
99 /* Set final state */
100 if constexpr (dim == 1)
101 return view.from_initial_state(state_type{{rho, Number(vel[0]), p}});
102 else if constexpr (dim == 2)
103 return view.from_initial_state(
104 state_type{{rho, Number(vel[0]), Number(vel[1]), p}});
105 else
106 return view.from_initial_state(state_type{
107 {rho, Number(vel[0]), Number(vel[1]), Number(vel[2]), p}});
108 }
109
110 private:
111 const HyperbolicSystem &hyperbolic_system_;
112 Number gamma_;
113 Number rho0_;
114 Number u0_;
115 Number p0_;
116 };
117 } // namespace EulerInitialStates
118} // namespace ryujin
Noh(const HyperbolicSystem &hyperbolic_system, const std::string &subsection)
typename View::state_type state_type
typename Description::template HyperbolicSystemView< dim, Number > View
auto compute(const dealii::Point< dim > &point, Number t) -> state_type final
typename Description::HyperbolicSystem HyperbolicSystem
T pow(const T x, const T b)
Euler::HyperbolicSystem HyperbolicSystem
Euler::HyperbolicSystem HyperbolicSystem
Definition: description.h:34