8#include <compile_time_options.h>
20#include <deal.II/base/parameter_acceptor.h>
21#include <deal.II/base/tensor.h>
37 template <
typename Number>
38 DEAL_II_ALWAYS_INLINE
inline Number
safe_division(
const Number &numerator,
39 const Number &denominator)
42 constexpr ScalarNumber
min = std::numeric_limits<ScalarNumber>::min();
44 return std::max(numerator, Number(0.)) /
45 std::max(denominator, Number(
min));
49 template <
int dim,
typename Number>
71 "Compressible Euler equations (arbitrary EOS)";
84 template <
int dim,
typename Number>
96 std::string equation_of_state_;
97 double reference_density_;
98 double vacuum_state_relaxation_small_;
99 double vacuum_state_relaxation_large_;
100 bool compute_strict_bounds_;
103 equation_of_state_list_;
106 std::shared_ptr<EquationOfState> selected_equation_of_state_;
108 template <
int dim,
typename Number>
130 template <
int dim,
typename Number>
139 : hyperbolic_system_(hyperbolic_system)
146 template <
int dim2,
typename Number2>
164 return hyperbolic_system_.equation_of_state_;
169 return hyperbolic_system_.reference_density_;
175 return hyperbolic_system_.vacuum_state_relaxation_small_;
181 return hyperbolic_system_.vacuum_state_relaxation_large_;
186 return hyperbolic_system_.compute_strict_bounds_;
200 const Number &e)
const
202 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
204 if constexpr (std::is_same_v<ScalarNumber, Number>) {
208 for (
unsigned int k = 0; k < Number::size(); ++k) {
219 DEAL_II_ALWAYS_INLINE
inline Number
222 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
224 if constexpr (std::is_same_v<ScalarNumber, Number>) {
225 return ScalarNumber(eos->specific_internal_energy(rho, p));
228 for (
unsigned int k = 0; k < Number::size(); ++k) {
229 e[k] =
ScalarNumber(eos->specific_internal_energy(rho[k], p[k]));
240 const Number &e)
const
242 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
244 if constexpr (std::is_same_v<ScalarNumber, Number>) {
248 for (
unsigned int k = 0; k < Number::size(); ++k) {
259 DEAL_II_ALWAYS_INLINE
inline Number
262 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
264 if constexpr (std::is_same_v<ScalarNumber, Number>) {
268 for (
unsigned int k = 0; k < Number::size(); ++k) {
280 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
289 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
299 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
339 using state_type = dealii::Tensor<1, problem_dimension, Number>;
345 dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim, Number>>;
357 []() -> std::array<std::string, problem_dimension> {
358 if constexpr (dim == 1)
359 return {
"rho",
"m",
"E"};
360 else if constexpr (dim == 2)
361 return {
"rho",
"m_1",
"m_2",
"E"};
362 else if constexpr (dim == 3)
363 return {
"rho",
"m_1",
"m_2",
"m_3",
"E"};
372 []() -> std::array<std::string, problem_dimension> {
373 if constexpr (dim == 1)
374 return {
"rho",
"v",
"e"};
375 else if constexpr (dim == 2)
376 return {
"rho",
"v_1",
"v_2",
"e"};
377 else if constexpr (dim == 3)
378 return {
"rho",
"v_1",
"v_2",
"v_3",
"e"};
396 std::array<std::string, n_precomputed_values>{
398 "surrogate_gamma_min",
399 "surrogate_specific_entropy",
400 "surrogate_harten_entropy"}};
411 std::array<Number, n_initial_precomputed_values>;
417 std::array<std::string, n_initial_precomputed_values>{};
423 StateVector<ScalarNumber, problem_dimension, n_precomputed_values>;
460 template <
typename DISPATCH,
typename SPARSITY>
462 const DISPATCH &dispatch_check,
463 const SPARSITY &sparsity_simd,
467 const bool skip_constrained_dofs =
true)
const;
532 const Number &gamma_min)
const;
545 const Number &gamma_min)
const;
560 const Number &gamma_min)
const;
601 const Number &gamma)
const;
622 template <
int component>
628 const dealii::Tensor<1, dim, Number> &normal)
const;
648 template <
typename Lambda>
652 const dealii::Tensor<1, dim, Number> &normal,
653 const Lambda &get_dirichlet_data)
const;
695 const unsigned int i,
701 const unsigned int *js,
711 const dealii::Tensor<1, dim, Number> &c_ij)
const;
721 const dealii::Tensor<1, dim, Number> &c_ij)
const =
delete;
732 const unsigned int i,
737 const unsigned int *js,
757 template <
typename ST>
772 template <
typename ST>
792 template <
typename Lambda>
794 const Lambda &lambda)
const;
807 const std::string &subsection )
808 : ParameterAcceptor(subsection)
810 equation_of_state_ =
"polytropic gas";
814 "The equation of state. Valid names are given by any of the "
815 "subsections defined below");
817 compute_strict_bounds_ =
true;
819 "compute strict bounds",
820 compute_strict_bounds_,
821 "Compute strict, but significantly more expensive bounds at various "
822 "places: (a) an expensive, but better upper wavespeed estimate in "
823 "the approximate RiemannSolver; (b) entropy viscosity-commutator "
824 "with correct gamma_min over the stencil; (c) mathematically correct "
825 "surrogate specific entropy minimum with gamma_min over the "
828 reference_density_ = 1.;
829 add_parameter(
"reference density",
831 "Problem specific density reference");
833 vacuum_state_relaxation_small_ = 1.e2;
834 add_parameter(
"vacuum state relaxation small",
835 vacuum_state_relaxation_small_,
836 "Problem specific vacuum relaxation parameter");
838 vacuum_state_relaxation_large_ = 1.e4;
839 add_parameter(
"vacuum state relaxation large",
840 vacuum_state_relaxation_large_,
841 "Problem specific vacuum relaxation parameter");
848 equation_of_state_list_, subsection);
850 const auto populate_functions = [
this]() {
851 bool initialized =
false;
852 for (
auto &it : equation_of_state_list_)
855 if (it->name() == equation_of_state_) {
856 selected_equation_of_state_ = it;
858 "Compressible Euler equations (" + it->name() +
" EOS)";
866 "Could not find an equation of state description with name \"" +
867 equation_of_state_ +
"\""));
870 ParameterAcceptor::parse_parameters_call_back.connect(populate_functions);
871 populate_functions();
875 template <
int dim,
typename Number>
876 template <
typename DISPATCH,
typename SPARSITY>
877 DEAL_II_ALWAYS_INLINE
inline void
879 unsigned int cycle [[maybe_unused]],
880 const DISPATCH &dispatch_check,
881 const SPARSITY &sparsity_simd,
885 const bool skip_constrained_dofs )
const
887 Assert(cycle == 0 || cycle == 1, dealii::ExcInternalError());
889 const auto &U = std::get<0>(state_vector);
890 auto &precomputed = std::get<1>(state_vector);
894 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
895 unsigned int stride_size = get_stride_size<Number>;
898 if (eos->prefer_vector_interface()) {
903 const auto offset = left;
904 const auto size = right - left;
906 static std::vector<double> p;
907 static std::vector<double> rho;
908 static std::vector<double> e;
917 for (
unsigned int i = 0; i < size; i += stride_size) {
918 const auto U_i = U.template get_tensor<Number>(offset + i);
919 const auto rho_i = density(U_i);
920 const auto e_i = internal_energy(U_i) / rho_i;
926 write_entry<Number>(rho, rho_i, i);
927 write_entry<Number>(e, e_i, i);
933 eos->pressure(p, rho, e);
937 for (
unsigned int i = 0; i < size; i += stride_size) {
939 const unsigned int row_length = sparsity_simd.row_length(i);
940 if (skip_constrained_dofs && row_length == 1)
946 const auto U_i = U.template get_tensor<Number>(offset + i);
947 const auto p_i = get_entry<Number>(p, i);
948 const auto gamma_i = surrogate_gamma(U_i, p_i);
949 const PT prec_i{p_i, gamma_i, Number(0.), Number(0.)};
950 precomputed.template write_tensor<Number>(prec_i, offset + i);
959 for (
unsigned int i = left; i < right; i += stride_size) {
961 const unsigned int row_length = sparsity_simd.row_length(i);
962 if (skip_constrained_dofs && row_length == 1)
967 const auto U_i = U.template get_tensor<Number>(i);
968 const auto rho_i = density(U_i);
969 const auto e_i = internal_energy(U_i) / rho_i;
970 const auto p_i = eos_pressure(rho_i, e_i);
972 const auto gamma_i = surrogate_gamma(U_i, p_i);
974 const PT prec_i{p_i, gamma_i, Number(0.), Number(0.)};
975 precomputed.template write_tensor<Number>(prec_i, i);
982 for (
unsigned int i = left; i < right; i += stride_size) {
986 const unsigned int row_length = sparsity_simd.row_length(i);
987 if (skip_constrained_dofs && row_length == 1)
992 const auto U_i = U.template get_tensor<Number>(i);
993 auto prec_i = precomputed.template get_tensor<Number, PT>(i);
994 auto &[p_i, gamma_min_i, s_i, eta_i] = prec_i;
996 const unsigned int *js = sparsity_simd.columns(i) + stride_size;
997 for (
unsigned int col_idx = 1; col_idx < row_length;
998 ++col_idx, js += stride_size) {
1000 const auto U_j = U.template get_tensor<Number>(js);
1001 const auto prec_j = precomputed.template get_tensor<Number, PT>(js);
1002 auto &[p_j, gamma_min_j, s_j, eta_j] = prec_j;
1003 const auto gamma_j = surrogate_gamma(U_j, p_j);
1004 gamma_min_i = std::min(gamma_min_i, gamma_j);
1007 s_i = surrogate_specific_entropy(U_i, gamma_min_i);
1008 eta_i = surrogate_harten_entropy(U_i, gamma_min_i);
1009 precomputed.template write_tensor<Number>(prec_i, i);
1015 template <
int dim,
typename Number>
1016 DEAL_II_ALWAYS_INLINE
inline Number
1023 template <
int dim,
typename Number>
1024 DEAL_II_ALWAYS_INLINE
inline Number
1026 const Number &rho)
const
1028 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
1029 const Number rho_cutoff_large =
1030 reference_density() * vacuum_state_relaxation_large() * eps;
1032 return dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
1033 std::abs(rho), rho_cutoff_large, Number(0.), rho);
1037 template <
int dim,
typename Number>
1038 DEAL_II_ALWAYS_INLINE
inline dealii::Tensor<1, dim, Number>
1041 dealii::Tensor<1, dim, Number> result;
1042 for (
unsigned int i = 0; i < dim; ++i)
1043 result[i] = U[1 + i];
1048 template <
int dim,
typename Number>
1049 DEAL_II_ALWAYS_INLINE
inline Number
1056 template <
int dim,
typename Number>
1057 DEAL_II_ALWAYS_INLINE
inline Number
1063 const Number rho_inverse =
ScalarNumber(1.) / density(U);
1064 const auto m = momentum(U);
1065 const Number E = total_energy(U);
1066 return E -
ScalarNumber(0.5) * m.norm_square() * rho_inverse;
1070 template <
int dim,
typename Number>
1071 DEAL_II_ALWAYS_INLINE
inline auto
1082 const Number rho_inverse =
ScalarNumber(1.) / density(U);
1083 const auto u = momentum(U) * rho_inverse;
1088 for (
unsigned int i = 0; i < dim; ++i) {
1089 result[1 + i] = -u[i];
1097 template <
int dim,
typename Number>
1098 DEAL_II_ALWAYS_INLINE
inline Number
1100 const state_type &U,
const Number &gamma_min)
const
1102 const auto b = Number(eos_interpolation_b());
1103 const auto pinf = Number(eos_interpolation_pinfty());
1104 const auto q = Number(eos_interpolation_q());
1106 const auto rho = density(U);
1109 const auto covolume = Number(1.) - b * rho;
1111 const auto shift = internal_energy(U) - rho * q - pinf * covolume;
1113 return shift *
ryujin::pow(rho_inverse - b, gamma_min) / covolume;
1117 template <
int dim,
typename Number>
1118 DEAL_II_ALWAYS_INLINE
inline Number
1120 const state_type &U,
const Number &gamma_min)
const
1122 const auto b = Number(eos_interpolation_b());
1123 const auto pinf = Number(eos_interpolation_pinfty());
1124 const auto q = Number(eos_interpolation_q());
1126 const auto rho = density(U);
1127 const auto m = momentum(U);
1128 const auto E = total_energy(U);
1129 const auto rho_rho_e_q =
1130 rho * E -
ScalarNumber(0.5) * m.norm_square() - rho * rho * q;
1132 const auto exponent =
ScalarNumber(1.) / (gamma_min + Number(1.));
1134 const auto covolume = Number(1.) - b * rho;
1135 const auto covolume_term =
ryujin::pow(covolume, gamma_min - Number(1.));
1137 const auto rho_pinfcov = rho * pinf * covolume;
1140 positive_part(rho_rho_e_q - rho_pinfcov) * covolume_term, exponent);
1144 template <
int dim,
typename Number>
1145 DEAL_II_ALWAYS_INLINE
inline auto
1147 const state_type &U,
const Number &eta,
const Number &gamma_min)
const
1164 const auto b = Number(eos_interpolation_b());
1165 const auto pinf = Number(eos_interpolation_pinfty());
1166 const auto q = Number(eos_interpolation_q());
1168 const auto rho = density(U);
1169 const auto m = momentum(U);
1170 const auto E = total_energy(U);
1172 const auto covolume = Number(1.) - b * rho;
1173 const auto covolume_inverse =
ScalarNumber(1.) / covolume;
1175 const auto shift = rho * E -
ScalarNumber(0.5) * m.norm_square() -
1176 rho * rho * q - rho * pinf * covolume;
1178 constexpr auto eps = std::numeric_limits<ScalarNumber>::epsilon();
1179 const auto regularization = m.norm() * eps;
1181 std::max(regularization, eta * covolume_inverse), gamma_min);
1183 factor *= fixed_power<2>(covolume_inverse) / (gamma_min + Number(1.));
1187 const auto first_term = E -
ScalarNumber(2.) * rho * q -
1189 const auto second_term = -(gamma_min - Number(1.)) * shift * b;
1191 result[0] = factor * (covolume * first_term + second_term);
1192 for (
unsigned int i = 0; i < dim; ++i)
1193 result[1 + i] = -factor * covolume * m[i];
1194 result[dim + 1] = factor * covolume * rho;
1200 template <
int dim,
typename Number>
1201 DEAL_II_ALWAYS_INLINE
inline Number
1203 const Number &p)
const
1205 const auto b = Number(eos_interpolation_b());
1206 const auto pinf = Number(eos_interpolation_pinfty());
1207 const auto q = Number(eos_interpolation_q());
1209 const auto rho = density(U);
1210 const auto rho_e = internal_energy(U);
1211 const auto covolume = Number(1.) - b * rho;
1213 const auto numerator = (p + pinf) * covolume;
1214 const auto denominator = rho_e - rho * q - covolume * pinf;
1219 template <
int dim,
typename Number>
1220 DEAL_II_ALWAYS_INLINE
inline Number
1222 const state_type &U,
const Number &gamma)
const
1224 const auto b = Number(eos_interpolation_b());
1225 const auto pinf = Number(eos_interpolation_pinfty());
1226 const auto q = Number(eos_interpolation_q());
1228 const auto rho = density(U);
1229 const auto rho_e = internal_energy(U);
1230 const auto covolume = Number(1.) - b * rho;
1238 template <
int dim,
typename Number>
1239 DEAL_II_ALWAYS_INLINE
inline Number
1241 const state_type &U,
const Number &gamma)
const
1243 const auto b = Number(eos_interpolation_b());
1244 const auto pinf = Number(eos_interpolation_pinfty());
1245 const auto q = Number(eos_interpolation_q());
1247 const auto rho = density(U);
1248 const auto rho_e = internal_energy(U);
1249 const auto covolume = Number(1.) - b * rho;
1252 (rho_e - rho * q - pinf * covolume) / (covolume * covolume * rho);
1253 radicand *= gamma * (gamma - 1.);
1258 template <
int dim,
typename Number>
1259 DEAL_II_ALWAYS_INLINE
inline bool
1262 const auto b = Number(eos_interpolation_b());
1263 const auto pinf = Number(eos_interpolation_pinfty());
1264 const auto q = Number(eos_interpolation_q());
1266 const auto rho = density(U);
1267 const auto rho_e = internal_energy(U);
1268 const auto covolume = Number(1.) - b * rho;
1270 const auto shift = rho_e - rho * q - pinf * covolume;
1272 constexpr auto gt = dealii::SIMDComparison::greater_than;
1275 dealii::compare_and_apply_mask<gt>(rho, T(0.), T(0.), T(-1.)) +
1276 dealii::compare_and_apply_mask<gt>(shift, T(0.), T(0.), T(-1.));
1279 if (!(test == Number(0.))) {
1280 std::cout << std::fixed << std::setprecision(16);
1281 std::cout <<
"Bounds violation: Negative state [rho, e] detected!\n";
1282 std::cout <<
"\t\trho: " << rho <<
"\n";
1283 std::cout <<
"\t\tint (shifted): " << shift <<
"\n";
1287 return (test == Number(0.));
1291 template <
int dim,
typename Number>
1292 template <
int component>
1293 DEAL_II_ALWAYS_INLINE
inline auto
1298 const Number &p_bar,
1299 const dealii::Tensor<1, dim, Number> &normal)
const ->
state_type
1301 static_assert(component == 1 || component == 2,
1302 "component has to be 1 or 2");
1304 const auto b = Number(eos_interpolation_b());
1305 const auto pinf = Number(eos_interpolation_pinfty());
1306 const auto q = Number(eos_interpolation_q());
1318 const auto m = momentum(U);
1319 const auto rho = density(U);
1320 const auto vn = m * normal / rho;
1322 const auto gamma = surrogate_gamma(U, p);
1323 const auto a = surrogate_speed_of_sound(U, gamma);
1324 const auto covolume = 1. - b * rho;
1326 const auto m_bar = momentum(U_bar);
1327 const auto rho_bar = density(U_bar);
1328 const auto vn_bar = m_bar * normal / rho_bar;
1330 const auto gamma_bar = surrogate_gamma(U_bar, p_bar);
1331 const auto a_bar = surrogate_speed_of_sound(U_bar, gamma_bar);
1332 const auto covolume_bar = 1. - b * rho_bar;
1345 component == 1 ? vn_bar - 2. * a_bar / (gamma_bar - 1.) * covolume_bar
1346 : vn - 2. * a / (gamma - 1.) * covolume;
1349 component == 2 ? vn_bar + 2. * a_bar / (gamma_bar - 1.) * covolume_bar
1350 : vn + 2. * a / (gamma - 1.) * covolume;
1359 dealii::ExcMessage(
"Encountered R_2 < R_1 in dynamic boundary value "
1360 "enforcement. This implies that the interpolation "
1361 "with Riemann characteristics failed."));
1363 const auto vperp = m / rho - vn * normal;
1365 const auto S = (p + pinf) *
ryujin::pow(Number(1.) / rho - b, gamma);
1382 const auto vn_new = Number(0.5) * (R_1 + R_2);
1412 const auto a_new_square =
1413 ryujin::fixed_power<2>((gamma - 1.) * (R_2 - R_1) / (4. * covolume));
1415 auto term =
ryujin::pow(a_new_square / (gamma * S), 1. / (gamma - 1.));
1417 term *=
std::pow(covolume, 2. / (gamma - 1.));
1420 const auto rho_new = term / (1. + b * term);
1422 const auto covolume_new = (1. - b * rho_new);
1423 const auto p_new = a_new_square / gamma * rho_new * covolume_new - pinf;
1429 const auto rho_e_new =
1430 rho_new * q + (p_new + gamma * pinf) * covolume_new / (gamma - 1.);
1434 for (
unsigned int d = 0; d < dim; ++d) {
1435 U_new[1 + d] = rho_new * (vn_new * normal + vperp)[d];
1438 rho_e_new + 0.5 * rho_new * (vn_new * vn_new + vperp.norm_square());
1444 template <
int dim,
typename Number>
1445 template <
typename Lambda>
1446 DEAL_II_ALWAYS_INLINE
inline auto
1448 dealii::types::boundary_id
id,
1450 const dealii::Tensor<1, dim, Number> &normal,
1451 const Lambda &get_dirichlet_data)
const ->
state_type
1456 result = get_dirichlet_data();
1460 auto m_dirichlet = momentum(get_dirichlet_data());
1461 for (
unsigned int k = 0; k < dim; ++k)
1462 result[k + 1] = m_dirichlet[k];
1465 auto m = momentum(U);
1466 m -= 1. * (m * normal) * normal;
1467 for (
unsigned int k = 0; k < dim; ++k)
1468 result[k + 1] = m[k];
1471 for (
unsigned int k = 0; k < dim; ++k)
1472 result[k + 1] = Number(0.);
1487 const auto m = momentum(U);
1488 const auto rho = density(U);
1489 const auto rho_e = internal_energy(U);
1496 const auto p = eos_pressure(rho, rho_e / rho);
1497 const auto gamma = surrogate_gamma(U, p);
1498 const auto a = surrogate_speed_of_sound(U, gamma);
1499 const auto vn = m * normal / rho;
1503 result = get_dirichlet_data();
1507 if (vn >= -a && vn <= 0.) {
1508 const auto U_dirichlet = get_dirichlet_data();
1509 const auto rho_dirichlet = density(U_dirichlet);
1510 const auto rho_e_dirichlet = internal_energy(U_dirichlet);
1511 const auto p_dirichlet =
1512 eos_pressure(rho_dirichlet, rho_e_dirichlet / rho_dirichlet);
1514 result = prescribe_riemann_characteristic<2>(
1515 U_dirichlet, p_dirichlet, U, p, normal);
1519 if (vn > 0. && vn <= a) {
1520 const auto U_dirichlet = get_dirichlet_data();
1521 const auto rho_dirichlet = density(U_dirichlet);
1522 const auto rho_e_dirichlet = internal_energy(U_dirichlet);
1523 const auto p_dirichlet =
1524 eos_pressure(rho_dirichlet, rho_e_dirichlet / rho_dirichlet);
1526 result = prescribe_riemann_characteristic<1>(
1527 U, p, U_dirichlet, p_dirichlet, normal);
1532 AssertThrow(
false, dealii::ExcNotImplemented());
1539 template <
int dim,
typename Number>
1540 DEAL_II_ALWAYS_INLINE
inline auto
1544 const auto rho_inverse =
ScalarNumber(1.) / density(U);
1545 const auto m = momentum(U);
1546 const auto E = total_energy(U);
1551 for (
unsigned int i = 0; i < dim; ++i) {
1552 result[1 + i] = m * (m[i] * rho_inverse);
1553 result[1 + i][i] += p;
1555 result[dim + 1] = m * (rho_inverse * (E + p));
1561 template <
int dim,
typename Number>
1562 DEAL_II_ALWAYS_INLINE
inline auto
1566 const unsigned int i,
1569 const auto &[p_i, gamma_min_i, s_i, eta_i] =
1570 pv.template get_tensor<Number, precomputed_type>(i);
1575 template <
int dim,
typename Number>
1576 DEAL_II_ALWAYS_INLINE
inline auto
1580 const unsigned int *js,
1583 const auto &[p_j, gamma_min_j, s_j, eta_j] =
1584 pv.template get_tensor<Number, precomputed_type>(js);
1589 template <
int dim,
typename Number>
1590 DEAL_II_ALWAYS_INLINE
inline auto
1594 const dealii::Tensor<1, dim, Number> &c_ij)
const ->
state_type
1600 template <
int dim,
typename Number>
1601 template <
typename ST>
1605 using T =
typename ST::value_type;
1606 static_assert(std::is_same_v<Number, T>,
"template mismatch");
1608 constexpr auto dim2 = ST::dimension - 2;
1609 static_assert(dim >= dim2,
1610 "the space dimension of the argument state must not be "
1611 "larger than the one of the target state");
1614 result[0] = state[0];
1615 result[dim + 1] = state[dim2 + 1];
1616 for (
unsigned int i = 1; i < dim2 + 1; ++i)
1617 result[i] = state[i];
1623 template <
int dim,
typename Number>
1624 template <
typename ST>
1625 DEAL_II_ALWAYS_INLINE
inline auto
1629 auto primitive_state = expand_state(initial_state);
1632 const auto rho = density(primitive_state);
1633 const auto p = total_energy(primitive_state);
1634 const auto e = eos_specific_internal_energy(rho, p);
1635 primitive_state[dim + 1] = e;
1637 return from_primitive_state(primitive_state);
1641 template <
int dim,
typename Number>
1642 DEAL_II_ALWAYS_INLINE
inline auto
1646 const auto rho = density(primitive_state);
1648 const auto u = momentum(primitive_state);
1650 const auto &e = total_energy(primitive_state);
1652 auto state = primitive_state;
1654 for (
unsigned int i = 1; i < dim + 1; ++i)
1658 state[dim + 1] = rho * e + Number(0.5) * rho * u * u;
1664 template <
int dim,
typename Number>
1665 DEAL_II_ALWAYS_INLINE
inline auto
1669 const auto rho = density(state);
1670 const auto rho_inverse = Number(1.) / rho;
1671 const auto rho_e = internal_energy(state);
1673 auto primitive_state = state;
1675 for (
unsigned int i = 1; i < dim + 1; ++i)
1676 primitive_state[i] *= rho_inverse;
1678 primitive_state[dim + 1] = rho_e * rho_inverse;
1680 return primitive_state;
1684 template <
int dim,
typename Number>
1685 template <
typename Lambda>
1689 auto result = state;
1690 const auto M = lambda(momentum(state));
1691 for (
unsigned int d = 0; d < dim; ++d)
1692 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
Number surrogate_speed_of_sound(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
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
void precomputation_loop(unsigned int cycle, const DISPATCH &dispatch_check, const SPARSITY &sparsity_simd, StateVector &state_vector, unsigned int left, unsigned int right, const bool skip_constrained_dofs=true) 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 prescribe_riemann_characteristic(const state_type &U, const Number &p, const state_type &U_bar, const Number &p_bar, const dealii::Tensor< 1, dim, Number > &normal) 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
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)
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)
std::set< std::shared_ptr< EquationOfState > > equation_of_state_list_type
DEAL_II_ALWAYS_INLINE Number safe_division(const Number &numerator, const Number &denominator)
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)