8#include <compile_time_options.h>
17#include <deal.II/base/parameter_acceptor.h>
18#include <deal.II/base/tensor.h>
26 template <
int dim,
typename Number>
27 class HyperbolicSystemView;
46 "Compressible Euler equations (polytropic gas EOS, optimized)";
59 template <
int dim,
typename Number>
72 double reference_density_;
73 double vacuum_state_relaxation_small_;
74 double vacuum_state_relaxation_large_;
76 double gamma_inverse_;
77 double gamma_minus_one_inverse_;
78 double gamma_minus_one_over_gamma_plus_one_;
79 double gamma_plus_one_inverse_;
81 template <
int dim,
typename Number>
103 template <
int dim,
typename Number>
112 : hyperbolic_system_(hyperbolic_system)
119 template <
int dim2,
typename Number2>
137 return hyperbolic_system_.gamma_;
142 return hyperbolic_system_.reference_density_;
148 return hyperbolic_system_.vacuum_state_relaxation_small_;
154 return hyperbolic_system_.vacuum_state_relaxation_large_;
174 return ScalarNumber(hyperbolic_system_.gamma_plus_one_inverse_);
179 return ScalarNumber(hyperbolic_system_.gamma_minus_one_inverse_);
186 hyperbolic_system_.gamma_minus_one_over_gamma_plus_one_);
223 using state_type = dealii::Tensor<1, problem_dimension, Number>;
229 dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim, Number>>;
241 []() -> std::array<std::string, problem_dimension> {
242 if constexpr (dim == 1)
243 return {
"rho",
"m",
"E"};
244 else if constexpr (dim == 2)
245 return {
"rho",
"m_1",
"m_2",
"E"};
246 else if constexpr (dim == 3)
247 return {
"rho",
"m_1",
"m_2",
"m_3",
"E"};
256 []() -> std::array<std::string, problem_dimension> {
257 if constexpr (dim == 1)
258 return {
"rho",
"v",
"p"};
259 else if constexpr (dim == 2)
260 return {
"rho",
"v_1",
"v_2",
"p"};
261 else if constexpr (dim == 3)
262 return {
"rho",
"v_1",
"v_2",
"v_3",
"p"};
280 std::array<std::string, n_precomputed_values>{
"s",
"eta_h"};
291 std::array<Number, n_initial_precomputed_values>;
297 std::array<std::string, n_initial_precomputed_values>{};
303 StateVector<ScalarNumber, problem_dimension, n_precomputed_values>;
340 template <
typename DISPATCH,
typename SPARSITY>
342 const DISPATCH &dispatch_check,
343 const SPARSITY &sparsity_simd,
347 const bool skip_constrained_dofs =
true)
const;
473 template <
int component>
476 const dealii::Tensor<1, dim, Number> &normal)
const;
484 template <
int component>
488 const dealii::Tensor<1, dim, Number> &normal)
const;
508 template <
typename Lambda>
512 const dealii::Tensor<1, dim, Number> &normal,
513 const Lambda &get_dirichlet_data)
const;
555 const unsigned int i,
561 const unsigned int *js,
571 const dealii::Tensor<1, dim, Number> &c_ij)
const;
579 const dealii::Tensor<1, dim, Number> &c_ij)
const =
delete;
591 const unsigned int i,
596 const unsigned int *js,
616 template <
typename ST>
632 template <
typename ST>
652 template <
typename Lambda>
654 const Lambda &lambda)
const;
667 : ParameterAcceptor(subsection)
670 add_parameter(
"gamma", gamma_,
"The ratio of specific heats");
672 reference_density_ = 1.;
673 add_parameter(
"reference density",
675 "Problem specific density reference");
677 vacuum_state_relaxation_small_ = 1.e2;
678 add_parameter(
"vacuum state relaxation small",
679 vacuum_state_relaxation_small_,
680 "Problem specific vacuum relaxation parameter");
682 vacuum_state_relaxation_large_ = 1.e4;
683 add_parameter(
"vacuum state relaxation large",
684 vacuum_state_relaxation_large_,
685 "Problem specific vacuum relaxation parameter");
691 const auto compute_inverses = [
this] {
692 gamma_inverse_ = 1. / gamma_;
693 gamma_plus_one_inverse_ = 1. / (gamma_ + 1.);
694 gamma_minus_one_inverse_ = 1. / (gamma_ - 1.);
695 gamma_minus_one_over_gamma_plus_one_ = (gamma_ - 1.) / (gamma_ + 1.);
699 ParameterAcceptor::parse_parameters_call_back.connect(compute_inverses);
703 template <
int dim,
typename Number>
704 template <
typename DISPATCH,
typename SPARSITY>
705 DEAL_II_ALWAYS_INLINE
inline void
707 unsigned int cycle [[maybe_unused]],
708 const DISPATCH &dispatch_check,
709 const SPARSITY &sparsity_simd,
713 const bool skip_constrained_dofs )
const
715 Assert(cycle == 0, dealii::ExcInternalError());
717 const auto &U = std::get<0>(state_vector);
718 auto &precomputed = std::get<1>(state_vector);
722 unsigned int stride_size = get_stride_size<Number>;
725 for (
unsigned int i = left; i < right; i += stride_size) {
728 const unsigned int row_length = sparsity_simd.row_length(i);
729 if (skip_constrained_dofs && row_length == 1)
734 const auto U_i = U.template get_tensor<Number>(i);
736 harten_entropy(U_i)};
737 precomputed.template write_tensor<Number>(prec_i, i);
742 template <
int dim,
typename Number>
743 DEAL_II_ALWAYS_INLINE
inline Number
750 template <
int dim,
typename Number>
751 DEAL_II_ALWAYS_INLINE
inline Number
753 const Number &rho)
const
755 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
756 const Number rho_cutoff_large =
757 reference_density() * vacuum_state_relaxation_large() * eps;
759 return dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
760 std::abs(rho), rho_cutoff_large, Number(0.), rho);
764 template <
int dim,
typename Number>
765 DEAL_II_ALWAYS_INLINE
inline dealii::Tensor<1, dim, Number>
768 dealii::Tensor<1, dim, Number> result;
769 for (
unsigned int i = 0; i < dim; ++i)
770 result[i] = U[1 + i];
775 template <
int dim,
typename Number>
776 DEAL_II_ALWAYS_INLINE
inline Number
783 template <
int dim,
typename Number>
784 DEAL_II_ALWAYS_INLINE
inline Number
790 const Number rho_inverse =
ScalarNumber(1.) / density(U);
791 const auto m = momentum(U);
792 const Number E = total_energy(U);
793 return E -
ScalarNumber(0.5) * m.norm_square() * rho_inverse;
797 template <
int dim,
typename Number>
798 DEAL_II_ALWAYS_INLINE
inline auto
809 const Number rho_inverse =
ScalarNumber(1.) / density(U);
810 const auto u = momentum(U) * rho_inverse;
815 for (
unsigned int i = 0; i < dim; ++i) {
816 result[1 + i] = -u[i];
824 template <
int dim,
typename Number>
825 DEAL_II_ALWAYS_INLINE
inline Number
829 return (gamma() -
ScalarNumber(1.)) * internal_energy(U);
833 template <
int dim,
typename Number>
834 DEAL_II_ALWAYS_INLINE
inline Number
838 const Number rho_inverse =
ScalarNumber(1.) / density(U);
839 const Number p = pressure(U);
840 return std::sqrt(gamma() * p * rho_inverse);
844 template <
int dim,
typename Number>
845 DEAL_II_ALWAYS_INLINE
inline Number
851 return internal_energy(U) *
ryujin::pow(rho_inverse, gamma());
855 template <
int dim,
typename Number>
856 DEAL_II_ALWAYS_INLINE
inline Number
861 const Number rho = density(U);
862 const auto m = momentum(U);
863 const Number E = total_energy(U);
865 const Number rho_rho_e = rho * E -
ScalarNumber(0.5) * m.norm_square();
866 return ryujin::pow(rho_rho_e, gamma_plus_one_inverse());
870 template <
int dim,
typename Number>
871 DEAL_II_ALWAYS_INLINE
inline auto
886 const Number rho = density(U);
887 const auto m = momentum(U);
888 const Number E = total_energy(U);
890 const Number rho_rho_e = rho * E -
ScalarNumber(0.5) * m.norm_square();
893 gamma_plus_one_inverse() *
894 ryujin::pow(rho_rho_e, -gamma() * gamma_plus_one_inverse());
898 result[0] = factor * E;
899 for (
unsigned int i = 0; i < dim; ++i)
900 result[1 + i] = -factor * m[i];
901 result[dim + 1] = factor * rho;
907 template <
int dim,
typename Number>
908 DEAL_II_ALWAYS_INLINE
inline Number
913 const auto p = pressure(U);
918 template <
int dim,
typename Number>
919 DEAL_II_ALWAYS_INLINE
inline auto
935 const Number rho = density(U);
937 const auto u = momentum(U) * rho_inverse;
938 const auto p = pressure(U);
940 const auto factor = (gamma() -
ScalarNumber(1.0)) * gamma_inverse() *
945 result[0] = factor *
ScalarNumber(0.5) * u.norm_square();
946 result[dim + 1] = factor;
947 for (
unsigned int i = 0; i < dim; ++i) {
948 result[1 + i] = -factor * u[i];
955 template <
int dim,
typename Number>
956 DEAL_II_ALWAYS_INLINE
inline bool
959 const auto rho_new = density(U);
960 const auto e_new = internal_energy(U);
961 const auto s_new = specific_entropy(U);
963 constexpr auto gt = dealii::SIMDComparison::greater_than;
966 dealii::compare_and_apply_mask<gt>(rho_new, T(0.), T(0.), T(-1.)) +
967 dealii::compare_and_apply_mask<gt>(e_new, T(0.), T(0.), T(-1.)) +
968 dealii::compare_and_apply_mask<gt>(s_new, T(0.), T(0.), T(-1.));
971 if (!(test == Number(0.))) {
972 std::cout << std::fixed << std::setprecision(16);
973 std::cout <<
"Bounds violation: Negative state [rho, e, s] detected!\n";
974 std::cout <<
"\t\trho: " << rho_new <<
"\n";
975 std::cout <<
"\t\tint: " << e_new <<
"\n";
976 std::cout <<
"\t\tent: " << s_new <<
"\n" << std::endl;
980 return (test == Number(0.));
984 template <
int dim,
typename Number>
985 template <
int component>
986 DEAL_II_ALWAYS_INLINE
inline auto
988 const state_type &U,
const dealii::Tensor<1, dim, Number> &normal)
const
989 -> std::array<state_type, 2>
991 static_assert(component == 1 || component == problem_dimension,
992 "Only first and last eigenvectors implemented");
994 const auto rho = density(U);
995 const auto m = momentum(U);
996 const auto v = m / rho;
997 const auto a = speed_of_sound(U);
998 const auto gamma = this->gamma();
1003 const auto e_k = 0.5 * v.norm_square();
1005 switch (component) {
1007 b[0] = (gamma - 1.) * e_k + a * v * normal;
1008 for (
unsigned int i = 0; i < dim; ++i)
1009 b[1 + i] = (1. - gamma) * v[i] - a * normal[i];
1010 b[dim + 1] = gamma - 1.;
1014 for (
unsigned int i = 0; i < dim; ++i)
1015 c[1 + i] = v[i] - a * normal[i];
1016 c[dim + 1] = a * a / (gamma - 1) + e_k - a * (v * normal);
1020 case problem_dimension:
1021 b[0] = (gamma - 1.) * e_k - a * v * normal;
1022 for (
unsigned int i = 0; i < dim; ++i)
1023 b[1 + i] = (1. - gamma) * v[i] + a * normal[i];
1024 b[dim + 1] = gamma - 1.;
1028 for (
unsigned int i = 0; i < dim; ++i)
1029 c[1 + i] = v[i] + a * normal[i];
1030 c[dim + 1] = a * a / (gamma - 1) + e_k + a * (v * normal);
1035 __builtin_unreachable();
1039 template <
int dim,
typename Number>
1040 template <
int component>
1041 DEAL_II_ALWAYS_INLINE
inline auto
1045 const dealii::Tensor<1, dim, Number> &normal)
const ->
state_type
1047 static_assert(component == 1 || component == 2,
1048 "component has to be 1 or 2");
1050 const auto m = momentum(U);
1051 const auto rho = density(U);
1052 const auto a = speed_of_sound(U);
1053 const auto vn = m * normal / rho;
1055 const auto m_bar = momentum(U_bar);
1056 const auto rho_bar = density(U_bar);
1057 const auto a_bar = speed_of_sound(U_bar);
1058 const auto vn_bar = m_bar * normal / rho_bar;
1062 const auto R_1 = component == 1
1068 const auto R_2 = component == 2
1072 const auto p = pressure(U);
1075 const auto vperp = m / rho - vn * normal;
1077 const auto vn_new = 0.5 * (R_1 + R_2);
1079 auto rho_new = 1. / (gamma() * s) *
1080 ryujin::fixed_power<2>(
ScalarNumber((gamma() - 1.) / 4.) *
1082 rho_new =
ryujin::pow(rho_new, 1. / (gamma() - 1.));
1084 const auto p_new = s *
std::pow(rho_new, gamma());
1088 for (
unsigned int d = 0; d < dim; ++d) {
1089 U_new[1 + d] = rho_new * (vn_new * normal + vperp)[d];
1092 0.5 * rho_new * (vn_new * vn_new + vperp.norm_square());
1098 template <
int dim,
typename Number>
1099 template <
typename Lambda>
1100 DEAL_II_ALWAYS_INLINE
inline auto
1102 dealii::types::boundary_id
id,
1104 const dealii::Tensor<1, dim, Number> &normal,
1105 const Lambda &get_dirichlet_data)
const ->
state_type
1110 result = get_dirichlet_data();
1114 auto m_dirichlet = momentum(get_dirichlet_data());
1115 for (
unsigned int k = 0; k < dim; ++k)
1116 result[k + 1] = m_dirichlet[k];
1119 auto m = momentum(U);
1120 m -= 1. * (m * normal) * normal;
1121 for (
unsigned int k = 0; k < dim; ++k)
1122 result[k + 1] = m[k];
1125 for (
unsigned int k = 0; k < dim; ++k)
1126 result[k + 1] = Number(0.);
1141 const auto m = momentum(U);
1142 const auto rho = density(U);
1143 const auto a = speed_of_sound(U);
1144 const auto vn = m * normal / rho;
1148 result = get_dirichlet_data();
1152 if (vn >= -a && vn <= 0.) {
1153 const auto U_dirichlet = get_dirichlet_data();
1154 result = prescribe_riemann_characteristic<2>(U_dirichlet, U, normal);
1158 if (vn > 0. && vn <= a) {
1159 const auto U_dirichlet = get_dirichlet_data();
1160 result = prescribe_riemann_characteristic<1>(U, U_dirichlet, normal);
1166 AssertThrow(
false, dealii::ExcNotImplemented());
1173 template <
int dim,
typename Number>
1174 DEAL_II_ALWAYS_INLINE
inline auto
1177 const auto rho_inverse =
ScalarNumber(1.) / density(U);
1178 const auto m = momentum(U);
1179 const auto p = pressure(U);
1180 const auto E = total_energy(U);
1185 for (
unsigned int i = 0; i < dim; ++i) {
1186 result[1 + i] = m * (m[i] * rho_inverse);
1187 result[1 + i][i] += p;
1189 result[dim + 1] = m * (rho_inverse * (E + p));
1195 template <
int dim,
typename Number>
1196 DEAL_II_ALWAYS_INLINE
inline auto
1200 const unsigned int ,
1207 template <
int dim,
typename Number>
1208 DEAL_II_ALWAYS_INLINE
inline auto
1212 const unsigned int * ,
1219 template <
int dim,
typename Number>
1220 DEAL_II_ALWAYS_INLINE
inline auto
1224 const dealii::Tensor<1, dim, Number> &c_ij)
const ->
state_type
1230 template <
int dim,
typename Number>
1231 template <
typename ST>
1235 using T =
typename ST::value_type;
1236 static_assert(std::is_same_v<Number, T>,
"template mismatch");
1238 constexpr auto dim2 = ST::dimension - 2;
1239 static_assert(dim >= dim2,
1240 "the space dimension of the argument state must not be "
1241 "larger than the one of the target state");
1244 result[0] = state[0];
1245 result[dim + 1] = state[dim2 + 1];
1246 for (
unsigned int i = 1; i < dim2 + 1; ++i)
1247 result[i] = state[i];
1253 template <
int dim,
typename Number>
1254 template <
typename ST>
1255 DEAL_II_ALWAYS_INLINE
inline auto
1259 const auto primitive_state = expand_state(initial_state);
1260 return from_primitive_state(primitive_state);
1264 template <
int dim,
typename Number>
1265 DEAL_II_ALWAYS_INLINE
inline auto
1269 const auto &rho = primitive_state[0];
1271 const auto u = momentum(primitive_state);
1272 const auto &p = primitive_state[dim + 1];
1274 auto state = primitive_state;
1276 for (
unsigned int i = 1; i < dim + 1; ++i)
1280 p / (
ScalarNumber(gamma() - 1.)) + Number(0.5) * rho * u * u;
1286 template <
int dim,
typename Number>
1287 DEAL_II_ALWAYS_INLINE
inline auto
1291 const auto &rho = state[0];
1292 const auto rho_inverse = Number(1.) / rho;
1293 const auto p = pressure(state);
1295 auto primitive_state = state;
1297 for (
unsigned int i = 1; i < dim + 1; ++i)
1298 primitive_state[i] *= rho_inverse;
1300 primitive_state[dim + 1] = p;
1302 return primitive_state;
1306 template <
int dim,
typename Number>
1307 template <
typename Lambda>
1311 auto result = state;
1312 const auto M = lambda(momentum(state));
1313 for (
unsigned int d = 0; d < dim; ++d)
1314 result[1 + d] = M[d];
typename get_value_type< Number >::type ScalarNumber
static constexpr unsigned int n_initial_precomputed_values
static constexpr bool have_eos_interpolation_b
state_type from_initial_state(const ST &initial_state) const
static Number internal_energy(const state_type &U)
state_type from_primitive_state(const state_type &primitive_state) const
dealii::Tensor< 1, problem_dimension, dealii::Tensor< 1, dim, Number > > flux_type
Number harten_entropy(const state_type &U) const
flux_contribution_type flux_contribution(const PrecomputedVector &pv, const InitialPrecomputedVector &ipv, const unsigned int i, const state_type &U_i) const
static const auto primitive_component_names
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
DEAL_II_ALWAYS_INLINE ScalarNumber reference_density() const
state_type harten_entropy_derivative(const state_type &U) const
static constexpr unsigned int n_precomputation_cycles
static const auto initial_precomputed_names
static constexpr bool have_high_order_flux
state_type nodal_source(const PrecomputedVector &pv, const unsigned int *js, const state_type &U_j, const ScalarNumber tau) const =delete
static const auto precomputed_names
state_type expand_state(const ST &state) const
HyperbolicSystemView(const HyperbolicSystem &hyperbolic_system)
static constexpr unsigned int problem_dimension
dealii::Tensor< 1, problem_dimension, Number > state_type
flux_type f(const state_type &U) const
static constexpr unsigned int n_precomputed_values
state_type nodal_source(const PrecomputedVector &pv, const unsigned int i, const state_type &U_i, const ScalarNumber tau) const =delete
std::array< Number, n_precomputed_values > precomputed_type
std::array< Number, n_initial_precomputed_values > initial_precomputed_type
state_type prescribe_riemann_characteristic(const state_type &U, const state_type &U_bar, const dealii::Tensor< 1, dim, Number > &normal) const
Number pressure(const state_type &U) const
state_type apply_galilei_transform(const state_type &state, const Lambda &lambda) const
Number mathematical_entropy(const state_type &U) const
DEAL_II_ALWAYS_INLINE ScalarNumber gamma_inverse() const
DEAL_II_ALWAYS_INLINE ScalarNumber gamma() const
static Number total_energy(const state_type &U)
DEAL_II_ALWAYS_INLINE ScalarNumber vacuum_state_relaxation_small() 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 const auto component_names
DEAL_II_ALWAYS_INLINE ScalarNumber gamma_minus_one_inverse() const
DEAL_II_ALWAYS_INLINE ScalarNumber gamma_plus_one_inverse() const
static Number density(const state_type &U)
Number filter_vacuum_density(const Number &rho) const
DEAL_II_ALWAYS_INLINE ScalarNumber vacuum_state_relaxation_large() const
Vectors::StateVector< ScalarNumber, problem_dimension, n_precomputed_values > StateVector
static state_type internal_energy_derivative(const state_type &U)
static constexpr bool have_source_terms
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
std::array< state_type, 2 > linearized_eigenvector(const state_type &U, const dealii::Tensor< 1, dim, Number > &normal) const
DEAL_II_ALWAYS_INLINE ScalarNumber gamma_minus_one_over_gamma_plus_one() const
flux_type flux_contribution_type
state_type mathematical_entropy_derivative(const state_type &U) const
static constexpr bool have_gamma
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
state_type to_primitive_state(const state_type &state) const
static dealii::Tensor< 1, dim, Number > momentum(const state_type &U)
Number speed_of_sound(const state_type &U) const
bool is_admissible(const state_type &U) const
Number specific_entropy(const state_type &U) const
HyperbolicSystem(const std::string &subsection="/HyperbolicSystem")
static const std::string problem_name
T pow(const T x, const T b)
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)