10#include <compile_time_options.h>
19#include <deal.II/base/parameter_acceptor.h>
20#include <deal.II/base/tensor.h>
28 template <
int dim,
typename Number>
29 class HyperbolicSystemView;
50 "Compressible Euler equations (arbitrary EOS)";
63 template <
int dim,
typename Number>
75 std::string equation_of_state_;
76 double reference_density_;
77 double vacuum_state_relaxation_small_;
78 double vacuum_state_relaxation_large_;
79 bool compute_strict_bounds_;
82 equation_of_state_list_;
85 std::shared_ptr<EquationOfState> selected_equation_of_state_;
87 template <
int dim,
typename Number>
109 template <
int dim,
typename Number>
118 : hyperbolic_system_(hyperbolic_system)
125 template <
int dim2,
typename Number2>
143 return hyperbolic_system_.equation_of_state_;
148 return hyperbolic_system_.reference_density_;
154 return hyperbolic_system_.vacuum_state_relaxation_small_;
160 return hyperbolic_system_.vacuum_state_relaxation_large_;
165 return hyperbolic_system_.compute_strict_bounds_;
179 const Number &e)
const
181 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
183 if constexpr (std::is_same_v<ScalarNumber, Number>) {
187 for (
unsigned int k = 0; k < Number::size(); ++k) {
198 DEAL_II_ALWAYS_INLINE
inline Number
201 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
203 if constexpr (std::is_same_v<ScalarNumber, Number>) {
204 return ScalarNumber(eos->specific_internal_energy(rho, p));
207 for (
unsigned int k = 0; k < Number::size(); ++k) {
208 e[k] =
ScalarNumber(eos->specific_internal_energy(rho[k], p[k]));
219 const Number &e)
const
221 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
223 if constexpr (std::is_same_v<ScalarNumber, Number>) {
227 for (
unsigned int k = 0; k < Number::size(); ++k) {
238 DEAL_II_ALWAYS_INLINE
inline Number
241 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
243 if constexpr (std::is_same_v<ScalarNumber, Number>) {
247 for (
unsigned int k = 0; k < Number::size(); ++k) {
259 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
268 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
278 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
318 using state_type = dealii::Tensor<1, problem_dimension, Number>;
324 dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim, Number>>;
336 []() -> std::array<std::string, problem_dimension> {
337 if constexpr (dim == 1)
338 return {
"rho",
"m",
"E"};
339 else if constexpr (dim == 2)
340 return {
"rho",
"m_1",
"m_2",
"E"};
341 else if constexpr (dim == 3)
342 return {
"rho",
"m_1",
"m_2",
"m_3",
"E"};
351 []() -> std::array<std::string, problem_dimension> {
352 if constexpr (dim == 1)
353 return {
"rho",
"v",
"e"};
354 else if constexpr (dim == 2)
355 return {
"rho",
"v_1",
"v_2",
"e"};
356 else if constexpr (dim == 3)
357 return {
"rho",
"v_1",
"v_2",
"v_3",
"e"};
375 std::array<std::string, n_precomputed_values>{
378 "surrogate_specific_entropy",
379 "surrogate_harten_entropy"}};
390 std::array<Number, n_initial_precomputed_values>;
396 std::array<std::string, n_initial_precomputed_values>{};
402 StateVector<ScalarNumber, problem_dimension, n_precomputed_values>;
439 template <
typename DISPATCH,
typename SPARSITY>
441 const DISPATCH &dispatch_check,
442 const SPARSITY &sparsity_simd,
445 unsigned int right)
const;
510 const Number &gamma_min)
const;
523 const Number &gamma_min)
const;
538 const Number &gamma_min)
const;
589 template <
int component>
593 const dealii::Tensor<1, dim, Number> &normal)
const;
613 template <
typename Lambda>
617 const dealii::Tensor<1, dim, Number> &normal,
618 const Lambda &get_dirichlet_data)
const;
660 const unsigned int i,
666 const unsigned int *js,
676 const dealii::Tensor<1, dim, Number> &c_ij)
const;
686 const dealii::Tensor<1, dim, Number> &c_ij)
const =
delete;
697 const unsigned int i,
702 const unsigned int *js,
722 template <
typename ST>
737 template <
typename ST>
757 template <
typename Lambda>
759 const Lambda &lambda)
const;
772 const std::string &subsection )
773 : ParameterAcceptor(subsection)
775 equation_of_state_ =
"polytropic gas";
779 "The equation of state. Valid names are given by any of the "
780 "subsections defined below");
782 compute_strict_bounds_ =
true;
784 "compute strict bounds",
785 compute_strict_bounds_,
786 "Compute strict, but significantly more expensive bounds at various "
787 "places: (a) an expensive, but better upper wavespeed estimate in "
788 "the approximate RiemannSolver; (b) entropy viscosity-commutator "
789 "with correct gamma_min over the stencil; (c) mathematically correct "
790 "surrogate specific entropy minimum with gamma_min over the "
793 reference_density_ = 1.;
794 add_parameter(
"reference density",
796 "Problem specific density reference");
798 vacuum_state_relaxation_small_ = 1.e2;
799 add_parameter(
"vacuum state relaxation small",
800 vacuum_state_relaxation_small_,
801 "Problem specific vacuum relaxation parameter");
803 vacuum_state_relaxation_large_ = 1.e4;
804 add_parameter(
"vacuum state relaxation large",
805 vacuum_state_relaxation_large_,
806 "Problem specific vacuum relaxation parameter");
813 equation_of_state_list_, subsection);
815 const auto populate_functions = [
this]() {
816 bool initialized =
false;
817 for (
auto &it : equation_of_state_list_)
820 if (it->name() == equation_of_state_) {
821 selected_equation_of_state_ = it;
823 "Compressible Euler equations (" + it->name() +
" EOS)";
831 "Could not find an equation of state description with name \"" +
832 equation_of_state_ +
"\""));
835 ParameterAcceptor::parse_parameters_call_back.connect(populate_functions);
836 populate_functions();
840 template <
int dim,
typename Number>
841 template <
typename DISPATCH,
typename SPARSITY>
842 DEAL_II_ALWAYS_INLINE
inline void
844 unsigned int cycle [[maybe_unused]],
845 const DISPATCH &dispatch_check,
846 const SPARSITY &sparsity_simd,
849 unsigned int right)
const
851 Assert(cycle == 0 || cycle == 1, dealii::ExcInternalError());
853 const auto &U = std::get<0>(state_vector);
854 auto &precomputed = std::get<1>(state_vector);
858 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
859 unsigned int stride_size = get_stride_size<Number>;
862 if (eos->prefer_vector_interface()) {
867 const auto offset = left;
868 const auto size = right - left;
870 static std::vector<double> p;
871 static std::vector<double> rho;
872 static std::vector<double> e;
881 for (
unsigned int i = 0; i < size; i += stride_size) {
882 const auto U_i = U.template get_tensor<Number>(offset + i);
883 const auto rho_i = density(U_i);
884 const auto e_i = internal_energy(U_i) / rho_i;
890 write_entry<Number>(rho, rho_i, i);
891 write_entry<Number>(e, e_i, i);
897 eos->pressure(p, rho, e);
901 for (
unsigned int i = 0; i < size; i += stride_size) {
903 const unsigned int row_length = sparsity_simd.row_length(i);
910 const auto U_i = U.template get_tensor<Number>(offset + i);
911 const auto p_i = get_entry<Number>(p, i);
912 const auto gamma_i = surrogate_gamma(U_i, p_i);
913 const PT prec_i{p_i, gamma_i, Number(0.), Number(0.)};
914 precomputed.template write_tensor<Number>(prec_i, offset + i);
923 for (
unsigned int i = left; i < right; i += stride_size) {
925 const unsigned int row_length = sparsity_simd.row_length(i);
931 const auto U_i = U.template get_tensor<Number>(i);
932 const auto rho_i = density(U_i);
933 const auto e_i = internal_energy(U_i) / rho_i;
934 const auto p_i = eos_pressure(rho_i, e_i);
936 const auto gamma_i = surrogate_gamma(U_i, p_i);
938 const PT prec_i{p_i, gamma_i, Number(0.), Number(0.)};
939 precomputed.template write_tensor<Number>(prec_i, i);
946 for (
unsigned int i = left; i < right; i += stride_size) {
950 const unsigned int row_length = sparsity_simd.row_length(i);
956 const auto U_i = U.template get_tensor<Number>(i);
957 auto prec_i = precomputed.template get_tensor<Number, PT>(i);
958 auto &[p_i, gamma_min_i, s_i, eta_i] = prec_i;
960 const unsigned int *js = sparsity_simd.columns(i) + stride_size;
961 for (
unsigned int col_idx = 1; col_idx < row_length;
962 ++col_idx, js += stride_size) {
964 const auto U_j = U.template get_tensor<Number>(js);
965 const auto prec_j = precomputed.template get_tensor<Number, PT>(js);
966 auto &[p_j, gamma_min_j, s_j, eta_j] = prec_j;
967 const auto gamma_j = surrogate_gamma(U_j, p_j);
968 gamma_min_i = std::min(gamma_min_i, gamma_j);
971 s_i = surrogate_specific_entropy(U_i, gamma_min_i);
972 eta_i = surrogate_harten_entropy(U_i, gamma_min_i);
973 precomputed.template write_tensor<Number>(prec_i, i);
979 template <
int dim,
typename Number>
980 DEAL_II_ALWAYS_INLINE
inline Number
987 template <
int dim,
typename Number>
988 DEAL_II_ALWAYS_INLINE
inline Number
990 const Number &rho)
const
992 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
993 const Number rho_cutoff_large =
994 reference_density() * vacuum_state_relaxation_large() * eps;
996 return dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
997 std::abs(rho), rho_cutoff_large, Number(0.), rho);
1001 template <
int dim,
typename Number>
1002 DEAL_II_ALWAYS_INLINE
inline dealii::Tensor<1, dim, Number>
1005 dealii::Tensor<1, dim, Number> result;
1006 for (
unsigned int i = 0; i < dim; ++i)
1007 result[i] = U[1 + i];
1012 template <
int dim,
typename Number>
1013 DEAL_II_ALWAYS_INLINE
inline Number
1020 template <
int dim,
typename Number>
1021 DEAL_II_ALWAYS_INLINE
inline Number
1027 const Number rho_inverse =
ScalarNumber(1.) / density(U);
1028 const auto m = momentum(U);
1029 const Number E = total_energy(U);
1030 return E -
ScalarNumber(0.5) * m.norm_square() * rho_inverse;
1034 template <
int dim,
typename Number>
1035 DEAL_II_ALWAYS_INLINE
inline auto
1046 const Number rho_inverse =
ScalarNumber(1.) / density(U);
1047 const auto u = momentum(U) * rho_inverse;
1052 for (
unsigned int i = 0; i < dim; ++i) {
1053 result[1 + i] = -u[i];
1061 template <
int dim,
typename Number>
1062 DEAL_II_ALWAYS_INLINE
inline Number
1064 const state_type &U,
const Number &gamma_min)
const
1066 const auto b = Number(eos_interpolation_b());
1067 const auto pinf = Number(eos_interpolation_pinfty());
1068 const auto q = Number(eos_interpolation_q());
1070 const auto rho = density(U);
1073 const auto covolume = Number(1.) - b * rho;
1075 const auto shift = internal_energy(U) - rho * q - pinf * covolume;
1077 return shift *
ryujin::pow(rho_inverse - b, gamma_min) / covolume;
1081 template <
int dim,
typename Number>
1082 DEAL_II_ALWAYS_INLINE
inline Number
1084 const state_type &U,
const Number &gamma_min)
const
1086 const auto b = Number(eos_interpolation_b());
1087 const auto pinf = Number(eos_interpolation_pinfty());
1088 const auto q = Number(eos_interpolation_q());
1090 const auto rho = density(U);
1091 const auto m = momentum(U);
1092 const auto E = total_energy(U);
1093 const auto rho_rho_e_q =
1094 rho * E -
ScalarNumber(0.5) * m.norm_square() - rho * rho * q;
1096 const auto exponent =
ScalarNumber(1.) / (gamma_min + Number(1.));
1098 const auto covolume = Number(1.) - b * rho;
1099 const auto covolume_term =
ryujin::pow(covolume, gamma_min - Number(1.));
1101 const auto rho_pinfcov = rho * pinf * covolume;
1103 return ryujin::pow((rho_rho_e_q - rho_pinfcov) * covolume_term, exponent);
1107 template <
int dim,
typename Number>
1108 DEAL_II_ALWAYS_INLINE
inline auto
1110 const state_type &U,
const Number &eta,
const Number &gamma_min)
const
1127 const auto b = Number(eos_interpolation_b());
1128 const auto pinf = Number(eos_interpolation_pinfty());
1129 const auto q = Number(eos_interpolation_q());
1131 const auto rho = density(U);
1132 const auto m = momentum(U);
1133 const auto E = total_energy(U);
1135 const auto covolume = Number(1.) - b * rho;
1136 const auto covolume_inverse =
ScalarNumber(1.) / covolume;
1138 const auto shift = rho * E -
ScalarNumber(0.5) * m.norm_square() -
1139 rho * rho * q - rho * pinf * covolume;
1141 const auto factor =
ryujin::pow(eta * covolume_inverse, -gamma_min) *
1142 fixed_power<2>(covolume_inverse) /
1143 (gamma_min + Number(1.));
1147 const auto first_term = E -
ScalarNumber(2.) * rho * q -
1149 const auto second_term = -(gamma_min - Number(1.)) * shift * b;
1151 result[0] = factor * (covolume * first_term + second_term);
1152 for (
unsigned int i = 0; i < dim; ++i)
1153 result[1 + i] = -factor * covolume * m[i];
1154 result[dim + 1] = factor * covolume * rho;
1160 template <
int dim,
typename Number>
1161 DEAL_II_ALWAYS_INLINE
inline Number
1163 const Number &p)
const
1165 const auto b = Number(eos_interpolation_b());
1166 const auto pinf = Number(eos_interpolation_pinfty());
1167 const auto q = Number(eos_interpolation_q());
1169 const auto rho = density(U);
1170 const auto rho_e = internal_energy(U);
1171 const auto covolume = Number(1.) - b * rho;
1173 const auto numerator = (p + pinf) * covolume;
1174 const auto denominator = rho_e - rho * q - covolume * pinf;
1175 return Number(1.) + numerator / denominator;
1179 template <
int dim,
typename Number>
1180 DEAL_II_ALWAYS_INLINE
inline Number
1182 const state_type &U,
const Number &gamma)
const
1184 const auto b = Number(eos_interpolation_b());
1185 const auto pinf = Number(eos_interpolation_pinfty());
1186 const auto q = Number(eos_interpolation_q());
1188 const auto rho = density(U);
1189 const auto rho_e = internal_energy(U);
1190 const auto covolume = Number(1.) - b * rho;
1193 return (gamma - Number(1.)) * (rho_e - rho * q) / covolume - gamma * pinf;
1197 template <
int dim,
typename Number>
1198 DEAL_II_ALWAYS_INLINE
inline bool
1201 const auto b = Number(eos_interpolation_b());
1202 const auto pinf = Number(eos_interpolation_pinfty());
1203 const auto q = Number(eos_interpolation_q());
1205 const auto rho = density(U);
1206 const auto rho_e = internal_energy(U);
1207 const auto covolume = Number(1.) - b * rho;
1209 const auto shift = rho_e - rho * q - pinf * covolume;
1211 constexpr auto gt = dealii::SIMDComparison::greater_than;
1214 dealii::compare_and_apply_mask<gt>(rho, T(0.), T(0.), T(-1.)) +
1215 dealii::compare_and_apply_mask<gt>(shift, T(0.), T(0.), T(-1.));
1218 if (!(test == Number(0.))) {
1219 std::cout << std::fixed << std::setprecision(16);
1220 std::cout <<
"Bounds violation: Negative state [rho, e] detected!\n";
1221 std::cout <<
"\t\trho: " << rho <<
"\n";
1222 std::cout <<
"\t\tint (shifted): " << shift <<
"\n";
1226 return (test == Number(0.));
1230 template <
int dim,
typename Number>
1231 template <
int component>
1232 DEAL_II_ALWAYS_INLINE
inline auto
1236 const dealii::Tensor<1, dim, Number> &normal)
const ->
state_type
1240 static_assert(component == 1 || component == 2,
1241 "component has to be 1 or 2");
1243 const auto m = momentum(U);
1244 const auto rho = density(U);
1245 const auto vn = m * normal / rho;
1247 const auto p = surrogate_pressure(U);
1248 const auto gamma = surrogate_gamma(U, p);
1249 const auto interpolation_b =
ScalarNumber(eos_interpolation_b());
1250 const auto x = Number(1.) - interpolation_b * rho;
1251 const auto a = std::sqrt(gamma * p / (rho * x));
1254 const auto m_bar = momentum(U_bar);
1255 const auto rho_bar = density(U_bar);
1256 const auto vn_bar = m_bar * normal / rho_bar;
1258 const auto p_bar = surrogate_pressure(U_bar);
1259 const auto gamma_bar = surrogate_gamma(U_bar, p_bar);
1260 const auto x_bar = Number(1.) - interpolation_b * rho_bar;
1261 const auto a_bar = std::sqrt(gamma_bar * p_bar / (rho_bar * x_bar));
1265 const auto R_1 = component == 1
1266 ? vn_bar - 2. * a_bar / (gamma_bar - 1.) * x_bar
1267 : vn - 2. * a / (gamma - 1.) * x;
1271 const auto R_2 = component == 2
1272 ? vn_bar + 2. * a_bar / (gamma_bar - 1.) * x_bar
1273 : vn + 2. * a / (gamma - 1.) * x;
1277 const auto vperp = m / rho - vn * normal;
1279 const auto vn_new = 0.5 * (R_1 + R_2);
1282 1. / (gamma * s) *
ryujin::pow((gamma - 1.) / 4. * (R_2 - R_1), 2.);
1284 1. / (
ryujin::pow(rho_new, 1. / (1. - gamma)) + interpolation_b);
1286 const auto x_new = 1. - interpolation_b * rho_new;
1293 for (
unsigned int d = 0; d < dim; ++d) {
1294 U_new[1 + d] = rho_new * (vn_new * normal + vperp)[d];
1296 U_new[1 + dim] = p_new / (gamma - 1.) +
1297 0.5 * rho_new * (vn_new * vn_new + vperp.norm_square());
1303 template <
int dim,
typename Number>
1304 template <
typename Lambda>
1305 DEAL_II_ALWAYS_INLINE
inline auto
1307 dealii::types::boundary_id
id,
1309 const dealii::Tensor<1, dim, Number> &normal,
1310 const Lambda &get_dirichlet_data)
const ->
state_type
1315 result = get_dirichlet_data();
1318 auto m = momentum(U);
1319 m -= 1. * (m * normal) * normal;
1320 for (
unsigned int k = 0; k < dim; ++k)
1321 result[k + 1] = m[k];
1324 for (
unsigned int k = 0; k < dim; ++k)
1325 result[k + 1] = Number(0.);
1342 const auto m = momentum(U);
1343 const auto rho = density(U);
1344 const auto a = speed_of_sound(U);
1345 const auto vn = m * normal / rho;
1349 result = get_dirichlet_data();
1353 if (vn >= -a && vn <= 0.) {
1354 const auto U_dirichlet = get_dirichlet_data();
1355 result = prescribe_riemann_characteristic<2>(U_dirichlet, U, normal);
1359 if (vn > 0. && vn <= a) {
1360 const auto U_dirichlet = get_dirichlet_data();
1361 result = prescribe_riemann_characteristic<1>(U, U_dirichlet, normal);
1371 template <
int dim,
typename Number>
1372 DEAL_II_ALWAYS_INLINE
inline auto
1376 const auto rho_inverse =
ScalarNumber(1.) / density(U);
1377 const auto m = momentum(U);
1378 const auto E = total_energy(U);
1383 for (
unsigned int i = 0; i < dim; ++i) {
1384 result[1 + i] = m * (m[i] * rho_inverse);
1385 result[1 + i][i] += p;
1387 result[dim + 1] = m * (rho_inverse * (E + p));
1393 template <
int dim,
typename Number>
1394 DEAL_II_ALWAYS_INLINE
inline auto
1398 const unsigned int i,
1401 const auto &[p_i, gamma_min_i, s_i, eta_i] =
1402 pv.template get_tensor<Number, precomputed_type>(i);
1407 template <
int dim,
typename Number>
1408 DEAL_II_ALWAYS_INLINE
inline auto
1412 const unsigned int *js,
1415 const auto &[p_j, gamma_min_j, s_j, eta_j] =
1416 pv.template get_tensor<Number, precomputed_type>(js);
1421 template <
int dim,
typename Number>
1422 DEAL_II_ALWAYS_INLINE
inline auto
1426 const dealii::Tensor<1, dim, Number> &c_ij)
const ->
state_type
1432 template <
int dim,
typename Number>
1433 template <
typename ST>
1437 using T =
typename ST::value_type;
1438 static_assert(std::is_same_v<Number, T>,
"template mismatch");
1440 constexpr auto dim2 = ST::dimension - 2;
1441 static_assert(dim >= dim2,
1442 "the space dimension of the argument state must not be "
1443 "larger than the one of the target state");
1446 result[0] = state[0];
1447 result[dim + 1] = state[dim2 + 1];
1448 for (
unsigned int i = 1; i < dim2 + 1; ++i)
1449 result[i] = state[i];
1455 template <
int dim,
typename Number>
1456 template <
typename ST>
1457 DEAL_II_ALWAYS_INLINE
inline auto
1461 auto primitive_state = expand_state(initial_state);
1464 const auto rho = density(primitive_state);
1465 const auto p = total_energy(primitive_state);
1466 const auto e = eos_specific_internal_energy(rho, p);
1467 primitive_state[dim + 1] = e;
1469 return from_primitive_state(primitive_state);
1473 template <
int dim,
typename Number>
1474 DEAL_II_ALWAYS_INLINE
inline auto
1478 const auto rho = density(primitive_state);
1480 const auto u = momentum(primitive_state);
1482 const auto &e = total_energy(primitive_state);
1484 auto state = primitive_state;
1486 for (
unsigned int i = 1; i < dim + 1; ++i)
1490 state[dim + 1] = rho * e + Number(0.5) * rho * u * u;
1496 template <
int dim,
typename Number>
1497 DEAL_II_ALWAYS_INLINE
inline auto
1501 const auto rho = density(state);
1502 const auto rho_inverse = Number(1.) / rho;
1503 const auto rho_e = internal_energy(state);
1505 auto primitive_state = state;
1507 for (
unsigned int i = 1; i < dim + 1; ++i)
1508 primitive_state[i] *= rho_inverse;
1510 primitive_state[dim + 1] = rho_e * rho_inverse;
1512 return primitive_state;
1516 template <
int dim,
typename Number>
1517 template <
typename Lambda>
1521 auto result = state;
1522 const auto M = lambda(momentum(state));
1523 for (
unsigned int d = 0; d < dim; ++d)
1524 result[1 + d] = M[d];
dealii::Tensor< 1, problem_dimension, Number > state_type
Vectors::StateVector< ScalarNumber, problem_dimension, n_precomputed_values > StateVector
state_type from_initial_state(const ST &initial_state) const
std::array< Number, n_initial_precomputed_values > initial_precomputed_type
Number surrogate_harten_entropy(const state_type &U, const Number &gamma_min) const
flux_contribution_type flux_contribution(const PrecomputedVector &pv, const InitialPrecomputedVector &piv, const unsigned int i, const state_type &U_i) const
static Number internal_energy(const state_type &U)
state_type to_primitive_state(const state_type &state) const
state_type high_order_flux_divergence(const flux_contribution_type &flux_i, const flux_contribution_type &flux_j, const dealii::Tensor< 1, dim, Number > &c_ij) const =delete
static constexpr bool have_high_order_flux
typename get_value_type< Number >::type ScalarNumber
static dealii::Tensor< 1, dim, Number > momentum(const state_type &U)
Number surrogate_pressure(const state_type &U, const Number &gamma) const
DEAL_II_ALWAYS_INLINE const std::string & equation_of_state() const
DEAL_II_ALWAYS_INLINE ScalarNumber vacuum_state_relaxation_large() const
static constexpr unsigned int n_precomputed_values
DEAL_II_ALWAYS_INLINE ScalarNumber reference_density() const
static const auto component_names
state_type from_primitive_state(const state_type &primitive_state) const
void precomputation_loop(unsigned int cycle, const DISPATCH &dispatch_check, const SPARSITY &sparsity_simd, StateVector &state_vector, unsigned int left, unsigned int right) const
static constexpr bool have_eos_interpolation_b
static state_type internal_energy_derivative(const state_type &U)
DEAL_II_ALWAYS_INLINE Number eos_temperature(const Number &rho, const Number &e) const
DEAL_II_ALWAYS_INLINE ScalarNumber eos_interpolation_b() const
static constexpr bool have_source_terms
DEAL_II_ALWAYS_INLINE Number eos_pressure(const Number &rho, const Number &e) const
static constexpr unsigned int problem_dimension
flux_type flux_contribution_type
state_type nodal_source(const PrecomputedVector &pv, const unsigned int *js, const state_type &U_j, const ScalarNumber tau) const =delete
Number surrogate_specific_entropy(const state_type &U, const Number &gamma_min) const
static const auto precomputed_names
static const auto primitive_component_names
DEAL_II_ALWAYS_INLINE ScalarNumber eos_interpolation_pinfty() const
static Number total_energy(const state_type &U)
bool is_admissible(const state_type &U) const
DEAL_II_ALWAYS_INLINE bool compute_strict_bounds() const
state_type flux_divergence(const flux_contribution_type &flux_i, const flux_contribution_type &flux_j, const dealii::Tensor< 1, dim, Number > &c_ij) const
DEAL_II_ALWAYS_INLINE Number eos_specific_internal_energy(const Number &rho, const Number &p) const
static constexpr unsigned int n_precomputation_cycles
HyperbolicSystemView(const HyperbolicSystem &hyperbolic_system)
DEAL_II_ALWAYS_INLINE ScalarNumber eos_interpolation_q() const
Number filter_vacuum_density(const Number &rho) const
flux_type f(const state_type &U, const Number &p) const
state_type expand_state(const ST &state) const
std::array< Number, n_precomputed_values > precomputed_type
state_type surrogate_harten_entropy_derivative(const state_type &U, const Number &eta, const Number &gamma_min) const
static const auto initial_precomputed_names
DEAL_II_ALWAYS_INLINE ScalarNumber vacuum_state_relaxation_small() const
Number surrogate_gamma(const state_type &U, const Number &p) const
state_type prescribe_riemann_characteristic(const state_type &U, const state_type &U_bar, const dealii::Tensor< 1, dim, Number > &normal) const
static constexpr unsigned int n_initial_precomputed_values
DEAL_II_ALWAYS_INLINE Number eos_speed_of_sound(const Number &rho, const Number &e) const
dealii::Tensor< 1, problem_dimension, dealii::Tensor< 1, dim, Number > > flux_type
static Number density(const state_type &U)
state_type apply_boundary_conditions(const dealii::types::boundary_id id, const state_type &U, const dealii::Tensor< 1, dim, Number > &normal, const Lambda &get_dirichlet_data) const
static constexpr bool have_gamma
state_type nodal_source(const PrecomputedVector &pv, const unsigned int i, const state_type &U_i, const ScalarNumber tau) const =delete
state_type apply_galilei_transform(const state_type &state, const Lambda &lambda) const
HyperbolicSystem(const std::string &subsection="/HyperbolicSystem")
static std::string problem_name
typename get_value_type< Number >::type ScalarNumber
dealii::Tensor< 1, problem_dimension, Number > state_type
std::array< Number, n_precomputed_values > precomputed_type
Vectors::StateVector< ScalarNumber, problem_dimension, n_precomputed_values > StateVector
void populate_equation_of_state_list(equation_of_state_list_type &equation_of_state_list, const std::string &subsection)
#define RYUJIN_OMP_SINGLE
T pow(const T x, const T b)
std::set< std::shared_ptr< EquationOfState > > equation_of_state_list_type
DEAL_II_ALWAYS_INLINE FT add(const FT &flux_left_ij, const FT &flux_right_ij)
DEAL_II_ALWAYS_INLINE dealii::Tensor< 1, problem_dim, T > contract(const FT &flux_ij, const TT &c_ij)