10#include <compile_time_options.h>
19#include <deal.II/base/parameter_acceptor.h>
20#include <deal.II/base/tensor.h>
36 template <
typename Number>
37 DEAL_II_ALWAYS_INLINE
inline Number
safe_division(
const Number &numerator,
38 const Number &denominator)
41 constexpr ScalarNumber
min = std::numeric_limits<ScalarNumber>::min();
43 return std::max(numerator, Number(0.)) /
44 std::max(denominator, Number(
min));
48 template <
int dim,
typename Number>
70 "Compressible Euler equations (arbitrary EOS)";
83 template <
int dim,
typename Number>
95 std::string equation_of_state_;
96 double reference_density_;
97 double vacuum_state_relaxation_small_;
98 double vacuum_state_relaxation_large_;
99 bool compute_strict_bounds_;
102 equation_of_state_list_;
105 std::shared_ptr<EquationOfState> selected_equation_of_state_;
107 template <
int dim,
typename Number>
129 template <
int dim,
typename Number>
138 : hyperbolic_system_(hyperbolic_system)
145 template <
int dim2,
typename Number2>
163 return hyperbolic_system_.equation_of_state_;
168 return hyperbolic_system_.reference_density_;
174 return hyperbolic_system_.vacuum_state_relaxation_small_;
180 return hyperbolic_system_.vacuum_state_relaxation_large_;
185 return hyperbolic_system_.compute_strict_bounds_;
199 const Number &e)
const
201 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
203 if constexpr (std::is_same_v<ScalarNumber, Number>) {
207 for (
unsigned int k = 0; k < Number::size(); ++k) {
218 DEAL_II_ALWAYS_INLINE
inline Number
221 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
223 if constexpr (std::is_same_v<ScalarNumber, Number>) {
224 return ScalarNumber(eos->specific_internal_energy(rho, p));
227 for (
unsigned int k = 0; k < Number::size(); ++k) {
228 e[k] =
ScalarNumber(eos->specific_internal_energy(rho[k], p[k]));
239 const Number &e)
const
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) {
258 DEAL_II_ALWAYS_INLINE
inline Number
261 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
263 if constexpr (std::is_same_v<ScalarNumber, Number>) {
267 for (
unsigned int k = 0; k < Number::size(); ++k) {
279 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
288 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
298 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
338 using state_type = dealii::Tensor<1, problem_dimension, Number>;
344 dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim, Number>>;
356 []() -> std::array<std::string, problem_dimension> {
357 if constexpr (dim == 1)
358 return {
"rho",
"m",
"E"};
359 else if constexpr (dim == 2)
360 return {
"rho",
"m_1",
"m_2",
"E"};
361 else if constexpr (dim == 3)
362 return {
"rho",
"m_1",
"m_2",
"m_3",
"E"};
371 []() -> std::array<std::string, problem_dimension> {
372 if constexpr (dim == 1)
373 return {
"rho",
"v",
"e"};
374 else if constexpr (dim == 2)
375 return {
"rho",
"v_1",
"v_2",
"e"};
376 else if constexpr (dim == 3)
377 return {
"rho",
"v_1",
"v_2",
"v_3",
"e"};
395 std::array<std::string, n_precomputed_values>{
397 "surrogate_gamma_min",
398 "surrogate_specific_entropy",
399 "surrogate_harten_entropy"}};
410 std::array<Number, n_initial_precomputed_values>;
416 std::array<std::string, n_initial_precomputed_values>{};
422 StateVector<ScalarNumber, problem_dimension, n_precomputed_values>;
459 template <
typename DISPATCH,
typename SPARSITY>
461 const DISPATCH &dispatch_check,
462 const SPARSITY &sparsity_simd,
466 const bool skip_constrained_dofs =
true)
const;
531 const Number &gamma_min)
const;
544 const Number &gamma_min)
const;
559 const Number &gamma_min)
const;
600 const Number &gamma)
const;
621 template <
int component>
627 const dealii::Tensor<1, dim, Number> &normal)
const;
647 template <
typename Lambda>
651 const dealii::Tensor<1, dim, Number> &normal,
652 const Lambda &get_dirichlet_data)
const;
694 const unsigned int i,
700 const unsigned int *js,
710 const dealii::Tensor<1, dim, Number> &c_ij)
const;
720 const dealii::Tensor<1, dim, Number> &c_ij)
const =
delete;
731 const unsigned int i,
736 const unsigned int *js,
756 template <
typename ST>
771 template <
typename ST>
791 template <
typename Lambda>
793 const Lambda &lambda)
const;
806 const std::string &subsection )
807 : ParameterAcceptor(subsection)
809 equation_of_state_ =
"polytropic gas";
813 "The equation of state. Valid names are given by any of the "
814 "subsections defined below");
816 compute_strict_bounds_ =
true;
818 "compute strict bounds",
819 compute_strict_bounds_,
820 "Compute strict, but significantly more expensive bounds at various "
821 "places: (a) an expensive, but better upper wavespeed estimate in "
822 "the approximate RiemannSolver; (b) entropy viscosity-commutator "
823 "with correct gamma_min over the stencil; (c) mathematically correct "
824 "surrogate specific entropy minimum with gamma_min over the "
827 reference_density_ = 1.;
828 add_parameter(
"reference density",
830 "Problem specific density reference");
832 vacuum_state_relaxation_small_ = 1.e2;
833 add_parameter(
"vacuum state relaxation small",
834 vacuum_state_relaxation_small_,
835 "Problem specific vacuum relaxation parameter");
837 vacuum_state_relaxation_large_ = 1.e4;
838 add_parameter(
"vacuum state relaxation large",
839 vacuum_state_relaxation_large_,
840 "Problem specific vacuum relaxation parameter");
847 equation_of_state_list_, subsection);
849 const auto populate_functions = [
this]() {
850 bool initialized =
false;
851 for (
auto &it : equation_of_state_list_)
854 if (it->name() == equation_of_state_) {
855 selected_equation_of_state_ = it;
857 "Compressible Euler equations (" + it->name() +
" EOS)";
865 "Could not find an equation of state description with name \"" +
866 equation_of_state_ +
"\""));
869 ParameterAcceptor::parse_parameters_call_back.connect(populate_functions);
870 populate_functions();
874 template <
int dim,
typename Number>
875 template <
typename DISPATCH,
typename SPARSITY>
876 DEAL_II_ALWAYS_INLINE
inline void
878 unsigned int cycle [[maybe_unused]],
879 const DISPATCH &dispatch_check,
880 const SPARSITY &sparsity_simd,
884 const bool skip_constrained_dofs )
const
886 Assert(cycle == 0 || cycle == 1, dealii::ExcInternalError());
888 const auto &U = std::get<0>(state_vector);
889 auto &precomputed = std::get<1>(state_vector);
893 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
894 unsigned int stride_size = get_stride_size<Number>;
897 if (eos->prefer_vector_interface()) {
902 const auto offset = left;
903 const auto size = right - left;
905 static std::vector<double> p;
906 static std::vector<double> rho;
907 static std::vector<double> e;
916 for (
unsigned int i = 0; i < size; i += stride_size) {
917 const auto U_i = U.template get_tensor<Number>(offset + i);
918 const auto rho_i = density(U_i);
919 const auto e_i = internal_energy(U_i) / rho_i;
925 write_entry<Number>(rho, rho_i, i);
926 write_entry<Number>(e, e_i, i);
932 eos->pressure(p, rho, e);
936 for (
unsigned int i = 0; i < size; i += stride_size) {
938 const unsigned int row_length = sparsity_simd.row_length(i);
939 if (skip_constrained_dofs && row_length == 1)
945 const auto U_i = U.template get_tensor<Number>(offset + i);
946 const auto p_i = get_entry<Number>(p, i);
947 const auto gamma_i = surrogate_gamma(U_i, p_i);
948 const PT prec_i{p_i, gamma_i, Number(0.), Number(0.)};
949 precomputed.template write_tensor<Number>(prec_i, offset + i);
958 for (
unsigned int i = left; i < right; i += stride_size) {
960 const unsigned int row_length = sparsity_simd.row_length(i);
961 if (skip_constrained_dofs && row_length == 1)
966 const auto U_i = U.template get_tensor<Number>(i);
967 const auto rho_i = density(U_i);
968 const auto e_i = internal_energy(U_i) / rho_i;
969 const auto p_i = eos_pressure(rho_i, e_i);
971 const auto gamma_i = surrogate_gamma(U_i, p_i);
973 const PT prec_i{p_i, gamma_i, Number(0.), Number(0.)};
974 precomputed.template write_tensor<Number>(prec_i, i);
981 for (
unsigned int i = left; i < right; i += stride_size) {
985 const unsigned int row_length = sparsity_simd.row_length(i);
986 if (skip_constrained_dofs && row_length == 1)
991 const auto U_i = U.template get_tensor<Number>(i);
992 auto prec_i = precomputed.template get_tensor<Number, PT>(i);
993 auto &[p_i, gamma_min_i, s_i, eta_i] = prec_i;
995 const unsigned int *js = sparsity_simd.columns(i) + stride_size;
996 for (
unsigned int col_idx = 1; col_idx < row_length;
997 ++col_idx, js += stride_size) {
999 const auto U_j = U.template get_tensor<Number>(js);
1000 const auto prec_j = precomputed.template get_tensor<Number, PT>(js);
1001 auto &[p_j, gamma_min_j, s_j, eta_j] = prec_j;
1002 const auto gamma_j = surrogate_gamma(U_j, p_j);
1003 gamma_min_i = std::min(gamma_min_i, gamma_j);
1006 s_i = surrogate_specific_entropy(U_i, gamma_min_i);
1007 eta_i = surrogate_harten_entropy(U_i, gamma_min_i);
1008 precomputed.template write_tensor<Number>(prec_i, i);
1014 template <
int dim,
typename Number>
1015 DEAL_II_ALWAYS_INLINE
inline Number
1022 template <
int dim,
typename Number>
1023 DEAL_II_ALWAYS_INLINE
inline Number
1025 const Number &rho)
const
1027 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
1028 const Number rho_cutoff_large =
1029 reference_density() * vacuum_state_relaxation_large() * eps;
1031 return dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
1032 std::abs(rho), rho_cutoff_large, Number(0.), rho);
1036 template <
int dim,
typename Number>
1037 DEAL_II_ALWAYS_INLINE
inline dealii::Tensor<1, dim, Number>
1040 dealii::Tensor<1, dim, Number> result;
1041 for (
unsigned int i = 0; i < dim; ++i)
1042 result[i] = U[1 + i];
1047 template <
int dim,
typename Number>
1048 DEAL_II_ALWAYS_INLINE
inline Number
1055 template <
int dim,
typename Number>
1056 DEAL_II_ALWAYS_INLINE
inline Number
1062 const Number rho_inverse =
ScalarNumber(1.) / density(U);
1063 const auto m = momentum(U);
1064 const Number E = total_energy(U);
1065 return E -
ScalarNumber(0.5) * m.norm_square() * rho_inverse;
1069 template <
int dim,
typename Number>
1070 DEAL_II_ALWAYS_INLINE
inline auto
1081 const Number rho_inverse =
ScalarNumber(1.) / density(U);
1082 const auto u = momentum(U) * rho_inverse;
1087 for (
unsigned int i = 0; i < dim; ++i) {
1088 result[1 + i] = -u[i];
1096 template <
int dim,
typename Number>
1097 DEAL_II_ALWAYS_INLINE
inline Number
1099 const state_type &U,
const Number &gamma_min)
const
1101 const auto b = Number(eos_interpolation_b());
1102 const auto pinf = Number(eos_interpolation_pinfty());
1103 const auto q = Number(eos_interpolation_q());
1105 const auto rho = density(U);
1108 const auto covolume = Number(1.) - b * rho;
1110 const auto shift = internal_energy(U) - rho * q - pinf * covolume;
1112 return shift *
ryujin::pow(rho_inverse - b, gamma_min) / covolume;
1116 template <
int dim,
typename Number>
1117 DEAL_II_ALWAYS_INLINE
inline Number
1119 const state_type &U,
const Number &gamma_min)
const
1121 const auto b = Number(eos_interpolation_b());
1122 const auto pinf = Number(eos_interpolation_pinfty());
1123 const auto q = Number(eos_interpolation_q());
1125 const auto rho = density(U);
1126 const auto m = momentum(U);
1127 const auto E = total_energy(U);
1128 const auto rho_rho_e_q =
1129 rho * E -
ScalarNumber(0.5) * m.norm_square() - rho * rho * q;
1131 const auto exponent =
ScalarNumber(1.) / (gamma_min + Number(1.));
1133 const auto covolume = Number(1.) - b * rho;
1134 const auto covolume_term =
ryujin::pow(covolume, gamma_min - Number(1.));
1136 const auto rho_pinfcov = rho * pinf * covolume;
1139 positive_part(rho_rho_e_q - rho_pinfcov) * covolume_term, exponent);
1143 template <
int dim,
typename Number>
1144 DEAL_II_ALWAYS_INLINE
inline auto
1146 const state_type &U,
const Number &eta,
const Number &gamma_min)
const
1163 const auto b = Number(eos_interpolation_b());
1164 const auto pinf = Number(eos_interpolation_pinfty());
1165 const auto q = Number(eos_interpolation_q());
1167 const auto rho = density(U);
1168 const auto m = momentum(U);
1169 const auto E = total_energy(U);
1171 const auto covolume = Number(1.) - b * rho;
1172 const auto covolume_inverse =
ScalarNumber(1.) / covolume;
1174 const auto shift = rho * E -
ScalarNumber(0.5) * m.norm_square() -
1175 rho * rho * q - rho * pinf * covolume;
1177 constexpr auto eps = std::numeric_limits<ScalarNumber>::epsilon();
1178 const auto regularization = m.norm() * eps;
1180 std::max(regularization, eta * covolume_inverse), gamma_min);
1182 factor *= fixed_power<2>(covolume_inverse) / (gamma_min + Number(1.));
1186 const auto first_term = E -
ScalarNumber(2.) * rho * q -
1188 const auto second_term = -(gamma_min - Number(1.)) * shift * b;
1190 result[0] = factor * (covolume * first_term + second_term);
1191 for (
unsigned int i = 0; i < dim; ++i)
1192 result[1 + i] = -factor * covolume * m[i];
1193 result[dim + 1] = factor * covolume * rho;
1199 template <
int dim,
typename Number>
1200 DEAL_II_ALWAYS_INLINE
inline Number
1202 const Number &p)
const
1204 const auto b = Number(eos_interpolation_b());
1205 const auto pinf = Number(eos_interpolation_pinfty());
1206 const auto q = Number(eos_interpolation_q());
1208 const auto rho = density(U);
1209 const auto rho_e = internal_energy(U);
1210 const auto covolume = Number(1.) - b * rho;
1212 const auto numerator = (p + pinf) * covolume;
1213 const auto denominator = rho_e - rho * q - covolume * pinf;
1218 template <
int dim,
typename Number>
1219 DEAL_II_ALWAYS_INLINE
inline Number
1221 const state_type &U,
const Number &gamma)
const
1223 const auto b = Number(eos_interpolation_b());
1224 const auto pinf = Number(eos_interpolation_pinfty());
1225 const auto q = Number(eos_interpolation_q());
1227 const auto rho = density(U);
1228 const auto rho_e = internal_energy(U);
1229 const auto covolume = Number(1.) - b * rho;
1237 template <
int dim,
typename Number>
1238 DEAL_II_ALWAYS_INLINE
inline Number
1240 const state_type &U,
const Number &gamma)
const
1242 const auto b = Number(eos_interpolation_b());
1243 const auto pinf = Number(eos_interpolation_pinfty());
1244 const auto q = Number(eos_interpolation_q());
1246 const auto rho = density(U);
1247 const auto rho_e = internal_energy(U);
1248 const auto covolume = Number(1.) - b * rho;
1251 (rho_e - rho * q - pinf * covolume) / (covolume * covolume * rho);
1252 radicand *= gamma * (gamma - 1.);
1257 template <
int dim,
typename Number>
1258 DEAL_II_ALWAYS_INLINE
inline bool
1261 const auto b = Number(eos_interpolation_b());
1262 const auto pinf = Number(eos_interpolation_pinfty());
1263 const auto q = Number(eos_interpolation_q());
1265 const auto rho = density(U);
1266 const auto rho_e = internal_energy(U);
1267 const auto covolume = Number(1.) - b * rho;
1269 const auto shift = rho_e - rho * q - pinf * covolume;
1271 constexpr auto gt = dealii::SIMDComparison::greater_than;
1274 dealii::compare_and_apply_mask<gt>(rho, T(0.), T(0.), T(-1.)) +
1275 dealii::compare_and_apply_mask<gt>(shift, T(0.), T(0.), T(-1.));
1278 if (!(test == Number(0.))) {
1279 std::cout << std::fixed << std::setprecision(16);
1280 std::cout <<
"Bounds violation: Negative state [rho, e] detected!\n";
1281 std::cout <<
"\t\trho: " << rho <<
"\n";
1282 std::cout <<
"\t\tint (shifted): " << shift <<
"\n";
1286 return (test == Number(0.));
1290 template <
int dim,
typename Number>
1291 template <
int component>
1292 DEAL_II_ALWAYS_INLINE
inline auto
1297 const Number &p_bar,
1298 const dealii::Tensor<1, dim, Number> &normal)
const ->
state_type
1300 static_assert(component == 1 || component == 2,
1301 "component has to be 1 or 2");
1303 const auto b = Number(eos_interpolation_b());
1304 const auto pinf = Number(eos_interpolation_pinfty());
1305 const auto q = Number(eos_interpolation_q());
1317 const auto m = momentum(U);
1318 const auto rho = density(U);
1319 const auto vn = m * normal / rho;
1321 const auto gamma = surrogate_gamma(U, p);
1322 const auto a = surrogate_speed_of_sound(U, gamma);
1323 const auto covolume = 1. - b * rho;
1325 const auto m_bar = momentum(U_bar);
1326 const auto rho_bar = density(U_bar);
1327 const auto vn_bar = m_bar * normal / rho_bar;
1329 const auto gamma_bar = surrogate_gamma(U_bar, p_bar);
1330 const auto a_bar = surrogate_speed_of_sound(U_bar, gamma_bar);
1331 const auto covolume_bar = 1. - b * rho_bar;
1344 component == 1 ? vn_bar - 2. * a_bar / (gamma_bar - 1.) * covolume_bar
1345 : vn - 2. * a / (gamma - 1.) * covolume;
1348 component == 2 ? vn_bar + 2. * a_bar / (gamma_bar - 1.) * covolume_bar
1349 : vn + 2. * a / (gamma - 1.) * covolume;
1358 dealii::ExcMessage(
"Encountered R_2 < R_1 in dynamic boundary value "
1359 "enforcement. This implies that the interpolation "
1360 "with Riemann characteristics failed."));
1362 const auto vperp = m / rho - vn * normal;
1364 const auto S = (p + pinf) *
ryujin::pow(Number(1.) / rho - b, gamma);
1381 const auto vn_new = Number(0.5) * (R_1 + R_2);
1411 const auto a_new_square =
1412 ryujin::fixed_power<2>((gamma - 1.) * (R_2 - R_1) / (4. * covolume));
1414 auto term =
ryujin::pow(a_new_square / (gamma * S), 1. / (gamma - 1.));
1416 term *=
std::pow(covolume, 2. / (gamma - 1.));
1419 const auto rho_new = term / (1. + b * term);
1421 const auto covolume_new = (1. - b * rho_new);
1422 const auto p_new = a_new_square / gamma * rho_new * covolume_new - pinf;
1428 const auto rho_e_new =
1429 rho_new * q + (p_new + gamma * pinf) * covolume_new / (gamma - 1.);
1433 for (
unsigned int d = 0; d < dim; ++d) {
1434 U_new[1 + d] = rho_new * (vn_new * normal + vperp)[d];
1437 rho_e_new + 0.5 * rho_new * (vn_new * vn_new + vperp.norm_square());
1443 template <
int dim,
typename Number>
1444 template <
typename Lambda>
1445 DEAL_II_ALWAYS_INLINE
inline auto
1447 dealii::types::boundary_id
id,
1449 const dealii::Tensor<1, dim, Number> &normal,
1450 const Lambda &get_dirichlet_data)
const ->
state_type
1455 result = get_dirichlet_data();
1459 auto m_dirichlet = momentum(get_dirichlet_data());
1460 for (
unsigned int k = 0; k < dim; ++k)
1461 result[k + 1] = m_dirichlet[k];
1464 auto m = momentum(U);
1465 m -= 1. * (m * normal) * normal;
1466 for (
unsigned int k = 0; k < dim; ++k)
1467 result[k + 1] = m[k];
1470 for (
unsigned int k = 0; k < dim; ++k)
1471 result[k + 1] = Number(0.);
1486 const auto m = momentum(U);
1487 const auto rho = density(U);
1488 const auto rho_e = internal_energy(U);
1495 const auto p = eos_pressure(rho, rho_e / rho);
1496 const auto gamma = surrogate_gamma(U, p);
1497 const auto a = surrogate_speed_of_sound(U, gamma);
1498 const auto vn = m * normal / rho;
1502 result = get_dirichlet_data();
1506 if (vn >= -a && vn <= 0.) {
1507 const auto U_dirichlet = get_dirichlet_data();
1508 const auto rho_dirichlet = density(U_dirichlet);
1509 const auto rho_e_dirichlet = internal_energy(U_dirichlet);
1510 const auto p_dirichlet =
1511 eos_pressure(rho_dirichlet, rho_e_dirichlet / rho_dirichlet);
1513 result = prescribe_riemann_characteristic<2>(
1514 U_dirichlet, p_dirichlet, U, p, normal);
1518 if (vn > 0. && vn <= a) {
1519 const auto U_dirichlet = get_dirichlet_data();
1520 const auto rho_dirichlet = density(U_dirichlet);
1521 const auto rho_e_dirichlet = internal_energy(U_dirichlet);
1522 const auto p_dirichlet =
1523 eos_pressure(rho_dirichlet, rho_e_dirichlet / rho_dirichlet);
1525 result = prescribe_riemann_characteristic<1>(
1526 U, p, U_dirichlet, p_dirichlet, normal);
1531 AssertThrow(
false, dealii::ExcNotImplemented());
1538 template <
int dim,
typename Number>
1539 DEAL_II_ALWAYS_INLINE
inline auto
1543 const auto rho_inverse =
ScalarNumber(1.) / density(U);
1544 const auto m = momentum(U);
1545 const auto E = total_energy(U);
1550 for (
unsigned int i = 0; i < dim; ++i) {
1551 result[1 + i] = m * (m[i] * rho_inverse);
1552 result[1 + i][i] += p;
1554 result[dim + 1] = m * (rho_inverse * (E + p));
1560 template <
int dim,
typename Number>
1561 DEAL_II_ALWAYS_INLINE
inline auto
1565 const unsigned int i,
1568 const auto &[p_i, gamma_min_i, s_i, eta_i] =
1569 pv.template get_tensor<Number, precomputed_type>(i);
1574 template <
int dim,
typename Number>
1575 DEAL_II_ALWAYS_INLINE
inline auto
1579 const unsigned int *js,
1582 const auto &[p_j, gamma_min_j, s_j, eta_j] =
1583 pv.template get_tensor<Number, precomputed_type>(js);
1588 template <
int dim,
typename Number>
1589 DEAL_II_ALWAYS_INLINE
inline auto
1593 const dealii::Tensor<1, dim, Number> &c_ij)
const ->
state_type
1599 template <
int dim,
typename Number>
1600 template <
typename ST>
1604 using T =
typename ST::value_type;
1605 static_assert(std::is_same_v<Number, T>,
"template mismatch");
1607 constexpr auto dim2 = ST::dimension - 2;
1608 static_assert(dim >= dim2,
1609 "the space dimension of the argument state must not be "
1610 "larger than the one of the target state");
1613 result[0] = state[0];
1614 result[dim + 1] = state[dim2 + 1];
1615 for (
unsigned int i = 1; i < dim2 + 1; ++i)
1616 result[i] = state[i];
1622 template <
int dim,
typename Number>
1623 template <
typename ST>
1624 DEAL_II_ALWAYS_INLINE
inline auto
1628 auto primitive_state = expand_state(initial_state);
1631 const auto rho = density(primitive_state);
1632 const auto p = total_energy(primitive_state);
1633 const auto e = eos_specific_internal_energy(rho, p);
1634 primitive_state[dim + 1] = e;
1636 return from_primitive_state(primitive_state);
1640 template <
int dim,
typename Number>
1641 DEAL_II_ALWAYS_INLINE
inline auto
1645 const auto rho = density(primitive_state);
1647 const auto u = momentum(primitive_state);
1649 const auto &e = total_energy(primitive_state);
1651 auto state = primitive_state;
1653 for (
unsigned int i = 1; i < dim + 1; ++i)
1657 state[dim + 1] = rho * e + Number(0.5) * rho * u * u;
1663 template <
int dim,
typename Number>
1664 DEAL_II_ALWAYS_INLINE
inline auto
1668 const auto rho = density(state);
1669 const auto rho_inverse = Number(1.) / rho;
1670 const auto rho_e = internal_energy(state);
1672 auto primitive_state = state;
1674 for (
unsigned int i = 1; i < dim + 1; ++i)
1675 primitive_state[i] *= rho_inverse;
1677 primitive_state[dim + 1] = rho_e * rho_inverse;
1679 return primitive_state;
1683 template <
int dim,
typename Number>
1684 template <
typename Lambda>
1688 auto result = state;
1689 const auto M = lambda(momentum(state));
1690 for (
unsigned int d = 0; d < dim; ++d)
1691 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)