ryujin 2.1.1 revision 350e54cc11f3d21282bddcf3e6517944dbc508bf
hyperbolic_system.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2020 - 2024 by the ryujin authors
4//
5
6#pragma once
7
9
10#include <compile_time_options.h>
11#include <convenience_macros.h>
12#include <discretization.h>
14#include <openmp.h>
15#include <patterns_conversion.h>
16#include <simd.h>
17#include <state_vector.h>
18
19#include <deal.II/base/parameter_acceptor.h>
20#include <deal.II/base/tensor.h>
21
22#include <array>
23
24namespace ryujin
25{
26 namespace EulerAEOS
27 {
28 /*
29 * For various divisions in the arbirtray equation of state module we
30 * have a mathematical guarantee that the numerator and denominator are
31 * nonnegative and the limit (of zero numerator and denominator) must
32 * converge to zero. The following function takes care of rounding
33 * issues when computing such quotients by (a) avoiding division by
34 * zero and (b) ensuring non-negativity of the result.
35 */
36 template <typename Number>
37 DEAL_II_ALWAYS_INLINE inline Number safe_division(const Number &numerator,
38 const Number &denominator)
39 {
40 using ScalarNumber = typename get_value_type<Number>::type;
41 constexpr ScalarNumber min = std::numeric_limits<ScalarNumber>::min();
42
43 return std::max(numerator, Number(0.)) /
44 std::max(denominator, Number(min));
45 }
46
47
48 template <int dim, typename Number>
50
63 class HyperbolicSystem final : public dealii::ParameterAcceptor
64 {
65 public:
69 static inline std::string problem_name =
70 "Compressible Euler equations (arbitrary EOS)";
71
75 HyperbolicSystem(const std::string &subsection = "/HyperbolicSystem");
76
83 template <int dim, typename Number>
84 auto view() const
85 {
87 }
88
89 private:
94
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_;
100
102 equation_of_state_list_;
103
104 using EquationOfState = EquationOfStateLibrary::EquationOfState;
105 std::shared_ptr<EquationOfState> selected_equation_of_state_;
106
107 template <int dim, typename Number>
110 }; /* HyperbolicSystem */
111
112
129 template <int dim, typename Number>
131 {
132 public:
137 HyperbolicSystemView(const HyperbolicSystem &hyperbolic_system)
138 : hyperbolic_system_(hyperbolic_system)
139 {
140 }
141
145 template <int dim2, typename Number2>
146 auto view() const
147 {
148 return HyperbolicSystemView<dim2, Number2>{hyperbolic_system_};
149 }
150
155
160
161 DEAL_II_ALWAYS_INLINE inline const std::string &equation_of_state() const
162 {
163 return hyperbolic_system_.equation_of_state_;
164 }
165
166 DEAL_II_ALWAYS_INLINE inline ScalarNumber reference_density() const
167 {
168 return hyperbolic_system_.reference_density_;
169 }
170
171 DEAL_II_ALWAYS_INLINE inline ScalarNumber
173 {
174 return hyperbolic_system_.vacuum_state_relaxation_small_;
175 }
176
177 DEAL_II_ALWAYS_INLINE inline ScalarNumber
179 {
180 return hyperbolic_system_.vacuum_state_relaxation_large_;
181 }
182
183 DEAL_II_ALWAYS_INLINE inline bool compute_strict_bounds() const
184 {
185 return hyperbolic_system_.compute_strict_bounds_;
186 }
187
189
193
198 DEAL_II_ALWAYS_INLINE inline Number eos_pressure(const Number &rho,
199 const Number &e) const
200 {
201 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
202
203 if constexpr (std::is_same_v<ScalarNumber, Number>) {
204 return ScalarNumber(eos->pressure(rho, e));
205 } else {
206 Number p;
207 for (unsigned int k = 0; k < Number::size(); ++k) {
208 p[k] = ScalarNumber(eos->pressure(rho[k], e[k]));
209 }
210 return p;
211 }
212 }
213
218 DEAL_II_ALWAYS_INLINE inline Number
219 eos_specific_internal_energy(const Number &rho, const Number &p) const
220 {
221 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
222
223 if constexpr (std::is_same_v<ScalarNumber, Number>) {
224 return ScalarNumber(eos->specific_internal_energy(rho, p));
225 } else {
226 Number e;
227 for (unsigned int k = 0; k < Number::size(); ++k) {
228 e[k] = ScalarNumber(eos->specific_internal_energy(rho[k], p[k]));
229 }
230 return e;
231 }
232 }
233
238 DEAL_II_ALWAYS_INLINE inline Number eos_temperature(const Number &rho,
239 const Number &e) const
240 {
241 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
242
243 if constexpr (std::is_same_v<ScalarNumber, Number>) {
244 return ScalarNumber(eos->temperature(rho, e));
245 } else {
246 Number temp;
247 for (unsigned int k = 0; k < Number::size(); ++k) {
248 temp[k] = ScalarNumber(eos->temperature(rho[k], e[k]));
249 }
250 return temp;
251 }
252 }
253
258 DEAL_II_ALWAYS_INLINE inline Number
259 eos_speed_of_sound(const Number &rho, const Number &e) const
260 {
261 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
262
263 if constexpr (std::is_same_v<ScalarNumber, Number>) {
264 return ScalarNumber(eos->speed_of_sound(rho, e));
265 } else {
266 Number c;
267 for (unsigned int k = 0; k < Number::size(); ++k) {
268 c[k] = ScalarNumber(eos->speed_of_sound(rho[k], e[k]));
269 }
270 return c;
271 }
272 }
273
277 DEAL_II_ALWAYS_INLINE inline ScalarNumber eos_interpolation_b() const
278 {
279 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
280 return ScalarNumber(eos->interpolation_b());
281 }
282
286 DEAL_II_ALWAYS_INLINE inline ScalarNumber eos_interpolation_pinfty() const
287 {
288 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
289 return ScalarNumber(eos->interpolation_pinfty());
290 }
291
296 DEAL_II_ALWAYS_INLINE inline ScalarNumber eos_interpolation_q() const
297 {
298 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
299 return ScalarNumber(eos->interpolation_q());
300 }
301
302
306 static constexpr bool have_gamma = false;
307
311 static constexpr bool have_eos_interpolation_b = true;
312
313
315
319
320 private:
321 const HyperbolicSystem &hyperbolic_system_;
322
323 public:
325
329
333 static constexpr unsigned int problem_dimension = 2 + dim;
334
338 using state_type = dealii::Tensor<1, problem_dimension, Number>;
339
343 using flux_type =
344 dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim, Number>>;
345
350
355 static inline const auto component_names =
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"};
363 __builtin_trap();
364 }();
365
370 static inline const auto primitive_component_names =
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"};
378 __builtin_trap();
379 }();
380
384 static constexpr unsigned int n_precomputed_values = 4;
385
389 using precomputed_type = std::array<Number, n_precomputed_values>;
390
394 static inline const auto precomputed_names =
395 std::array<std::string, n_precomputed_values>{
396 {"p",
397 "surrogate_gamma_min",
398 "surrogate_specific_entropy",
399 "surrogate_harten_entropy"}};
400
404 static constexpr unsigned int n_initial_precomputed_values = 0;
405
410 std::array<Number, n_initial_precomputed_values>;
411
415 static inline const auto initial_precomputed_names =
416 std::array<std::string, n_initial_precomputed_values>{};
417
421 using StateVector = Vectors::
422 StateVector<ScalarNumber, problem_dimension, n_precomputed_values>;
423
429
435
443
445
449
453 static constexpr unsigned int n_precomputation_cycles = 2;
454
459 template <typename DISPATCH, typename SPARSITY>
460 void precomputation_loop(unsigned int cycle,
461 const DISPATCH &dispatch_check,
462 const SPARSITY &sparsity_simd,
463 StateVector &state_vector,
464 unsigned int left,
465 unsigned int right,
466 const bool skip_constrained_dofs = true) const;
467
469
473
478 static Number density(const state_type &U);
479
486 Number filter_vacuum_density(const Number &rho) const;
487
492 static dealii::Tensor<1, dim, Number> momentum(const state_type &U);
493
498 static Number total_energy(const state_type &U);
499
504 static Number internal_energy(const state_type &U);
505
512
514
520
531 const Number &gamma_min) const;
532
543 Number surrogate_harten_entropy(const state_type &U,
544 const Number &gamma_min) const;
545
558 const Number &eta,
559 const Number &gamma_min) const;
560
572 Number surrogate_gamma(const state_type &U, const Number &p) const;
573
589 Number surrogate_pressure(const state_type &U, const Number &gamma) const;
590
599 Number surrogate_speed_of_sound(const state_type &U,
600 const Number &gamma) const;
601
607 bool is_admissible(const state_type &U) const;
608
610
614
621 template <int component>
623 const state_type &U,
624 const Number &p,
625 const state_type &U_bar,
626 const Number &p_bar,
627 const dealii::Tensor<1, dim, Number> &normal) const;
628
647 template <typename Lambda>
649 apply_boundary_conditions(const dealii::types::boundary_id id,
650 const state_type &U,
651 const dealii::Tensor<1, dim, Number> &normal,
652 const Lambda &get_dirichlet_data) const;
653
655
659
670 flux_type f(const state_type &U, const Number &p) const;
671
693 const InitialPrecomputedVector &piv,
694 const unsigned int i,
695 const state_type &U_i) const;
696
699 const InitialPrecomputedVector &piv,
700 const unsigned int *js,
701 const state_type &U_j) const;
702
709 const flux_contribution_type &flux_j,
710 const dealii::Tensor<1, dim, Number> &c_ij) const;
711
715 static constexpr bool have_high_order_flux = false;
716
718 const flux_contribution_type &flux_i,
719 const flux_contribution_type &flux_j,
720 const dealii::Tensor<1, dim, Number> &c_ij) const = delete;
721
726
728 static constexpr bool have_source_terms = false;
729
731 const unsigned int i,
732 const state_type &U_i,
733 const ScalarNumber tau) const = delete;
734
736 const unsigned int *js,
737 const state_type &U_j,
738 const ScalarNumber tau) const = delete;
739
741
745
756 template <typename ST>
757 state_type expand_state(const ST &state) const;
758
771 template <typename ST>
772 state_type from_initial_state(const ST &initial_state) const;
773
778 state_type from_primitive_state(const state_type &primitive_state) const;
779
784 state_type to_primitive_state(const state_type &state) const;
785
791 template <typename Lambda>
793 const Lambda &lambda) const;
795 }; /* HyperbolicSystemView */
796
797
798 /*
799 * -------------------------------------------------------------------------
800 * Inline definitions
801 * -------------------------------------------------------------------------
802 */
803
804
806 const std::string &subsection /*= "HyperbolicSystem"*/)
807 : ParameterAcceptor(subsection)
808 {
809 equation_of_state_ = "polytropic gas";
810 add_parameter(
811 "equation of state",
812 equation_of_state_,
813 "The equation of state. Valid names are given by any of the "
814 "subsections defined below");
815
816 compute_strict_bounds_ = true;
817 add_parameter(
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 "
825 "stencil.");
826
827 reference_density_ = 1.;
828 add_parameter("reference density",
829 reference_density_,
830 "Problem specific density reference");
831
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");
836
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");
841
842 /*
843 * And finally populate the equation of state list with all equation of
844 * state configurations defined in the EquationOfState namespace:
845 */
847 equation_of_state_list_, subsection);
848
849 const auto populate_functions = [this]() {
850 bool initialized = false;
851 for (auto &it : equation_of_state_list_)
852
853 /* Populate EOS-specific quantities and functions */
854 if (it->name() == equation_of_state_) {
855 selected_equation_of_state_ = it;
857 "Compressible Euler equations (" + it->name() + " EOS)";
858 initialized = true;
859 break;
860 }
861
862 AssertThrow(
863 initialized,
864 dealii::ExcMessage(
865 "Could not find an equation of state description with name \"" +
866 equation_of_state_ + "\""));
867 };
868
869 ParameterAcceptor::parse_parameters_call_back.connect(populate_functions);
870 populate_functions();
871 }
872
873
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,
881 StateVector &state_vector,
882 unsigned int left,
883 unsigned int right,
884 const bool skip_constrained_dofs /*= true*/) const
885 {
886 Assert(cycle == 0 || cycle == 1, dealii::ExcInternalError());
887
888 const auto &U = std::get<0>(state_vector);
889 auto &precomputed = std::get<1>(state_vector);
890
891 /* We are inside a thread parallel context */
892
893 const auto &eos = hyperbolic_system_.selected_equation_of_state_;
894 unsigned int stride_size = get_stride_size<Number>;
895
896 if (cycle == 0) {
897 if (eos->prefer_vector_interface()) {
898 /*
899 * Set up temporary storage for p, rho, e and make two calls into
900 * the eos library.
901 */
902 const auto offset = left;
903 const auto size = right - left;
904
905 static /* shared */ std::vector<double> p;
906 static /* shared */ std::vector<double> rho;
907 static /* shared */ std::vector<double> e;
909 {
910 p.resize(size);
911 rho.resize(size);
912 e.resize(size);
913 }
914
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;
920 /*
921 * Populate rho and e also for interpolated values from
922 * constrainted degrees of freedom so that the vectors contain
923 * physically admissible entries throughout.
924 */
925 write_entry<Number>(rho, rho_i, i);
926 write_entry<Number>(e, e_i, i);
927 }
928
929 /* Make sure the call into eospac (and others) is single threaded. */
931 {
932 eos->pressure(p, rho, e);
933 }
934
936 for (unsigned int i = 0; i < size; i += stride_size) {
937 /* Skip constrained degrees of freedom: */
938 const unsigned int row_length = sparsity_simd.row_length(i);
939 if (skip_constrained_dofs && row_length == 1)
940 continue;
941
942 dispatch_check(i);
943
944 using PT = precomputed_type;
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);
950 }
951 } else {
952 /*
953 * This is the variant with slightly better performance provided
954 * that a call to the eos is not too expensive. This variant
955 * calls into the eos library for every single degree of freedom.
956 */
958 for (unsigned int i = left; i < right; i += stride_size) {
959 /* Skip constrained degrees of freedom: */
960 const unsigned int row_length = sparsity_simd.row_length(i);
961 if (skip_constrained_dofs && row_length == 1)
962 continue;
963
964 dispatch_check(i);
965
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);
970
971 const auto gamma_i = surrogate_gamma(U_i, p_i);
972 using PT = precomputed_type;
973 const PT prec_i{p_i, gamma_i, Number(0.), Number(0.)};
974 precomputed.template write_tensor<Number>(prec_i, i);
975 }
976 } /* prefer_vector_interface */
977 } /* cycle == 0 */
978
979 if (cycle == 1) {
981 for (unsigned int i = left; i < right; i += stride_size) {
982 using PT = precomputed_type;
983
984 /* Skip constrained degrees of freedom: */
985 const unsigned int row_length = sparsity_simd.row_length(i);
986 if (skip_constrained_dofs && row_length == 1)
987 continue;
988
989 dispatch_check(i);
990
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;
994
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) {
998
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);
1004 }
1005
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);
1009 }
1010 }
1011 }
1012
1013
1014 template <int dim, typename Number>
1015 DEAL_II_ALWAYS_INLINE inline Number
1017 {
1018 return U[0];
1019 }
1020
1021
1022 template <int dim, typename Number>
1023 DEAL_II_ALWAYS_INLINE inline Number
1025 const Number &rho) const
1026 {
1027 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
1028 const Number rho_cutoff_large =
1029 reference_density() * vacuum_state_relaxation_large() * eps;
1030
1031 return dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
1032 std::abs(rho), rho_cutoff_large, Number(0.), rho);
1033 }
1034
1035
1036 template <int dim, typename Number>
1037 DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim, Number>
1039 {
1040 dealii::Tensor<1, dim, Number> result;
1041 for (unsigned int i = 0; i < dim; ++i)
1042 result[i] = U[1 + i];
1043 return result;
1044 }
1045
1046
1047 template <int dim, typename Number>
1048 DEAL_II_ALWAYS_INLINE inline Number
1050 {
1051 return U[1 + dim];
1052 }
1053
1054
1055 template <int dim, typename Number>
1056 DEAL_II_ALWAYS_INLINE inline Number
1058 {
1059 /*
1060 * rho e = (E - 1/2*m^2/rho)
1061 */
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;
1066 }
1067
1068
1069 template <int dim, typename Number>
1070 DEAL_II_ALWAYS_INLINE inline auto
1072 const state_type &U) -> state_type
1073 {
1074 /*
1075 * With
1076 * rho e = E - 1/2 |m|^2 / rho
1077 * we get
1078 * (rho e)' = (1/2m^2/rho^2, -m/rho , 1 )^T
1079 */
1080
1081 const Number rho_inverse = ScalarNumber(1.) / density(U);
1082 const auto u = momentum(U) * rho_inverse;
1083
1084 state_type result;
1085
1086 result[0] = ScalarNumber(0.5) * u.norm_square();
1087 for (unsigned int i = 0; i < dim; ++i) {
1088 result[1 + i] = -u[i];
1089 }
1090 result[dim + 1] = ScalarNumber(1.);
1091
1092 return result;
1093 }
1094
1095
1096 template <int dim, typename Number>
1097 DEAL_II_ALWAYS_INLINE inline Number
1099 const state_type &U, const Number &gamma_min) const
1100 {
1101 const auto b = Number(eos_interpolation_b());
1102 const auto pinf = Number(eos_interpolation_pinfty());
1103 const auto q = Number(eos_interpolation_q());
1104
1105 const auto rho = density(U);
1106 const auto rho_inverse = ScalarNumber(1.) / rho;
1107
1108 const auto covolume = Number(1.) - b * rho;
1109
1110 const auto shift = internal_energy(U) - rho * q - pinf * covolume;
1111
1112 return shift * ryujin::pow(rho_inverse - b, gamma_min) / covolume;
1113 }
1114
1115
1116 template <int dim, typename Number>
1117 DEAL_II_ALWAYS_INLINE inline Number
1119 const state_type &U, const Number &gamma_min) const
1120 {
1121 const auto b = Number(eos_interpolation_b());
1122 const auto pinf = Number(eos_interpolation_pinfty());
1123 const auto q = Number(eos_interpolation_q());
1124
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;
1130
1131 const auto exponent = ScalarNumber(1.) / (gamma_min + Number(1.));
1132
1133 const auto covolume = Number(1.) - b * rho;
1134 const auto covolume_term = ryujin::pow(covolume, gamma_min - Number(1.));
1135
1136 const auto rho_pinfcov = rho * pinf * covolume;
1137
1138 return ryujin::pow(
1139 positive_part(rho_rho_e_q - rho_pinfcov) * covolume_term, exponent);
1140 }
1141
1142
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
1147 -> state_type
1148 {
1149 /*
1150 * With
1151 * eta = (shift * (1-b*rho)^{gamma-1}) ^ {1/(gamma+1)},
1152 * shift = rho * E - 1/2 |m|^2 - rho^2 * q - p_infty * rho * (1 - b rho)
1153 *
1154 * shift' = [E - 2 * rho * q - p_infty * (1 - 2 b rho), -m, rho]^T
1155 * factor = 1/(gamma+1) * (eta/(1-b rho))^-gamma / (1-b rho)^2
1156 *
1157 * we get
1158 *
1159 * eta' = factor * (1-b*rho) * shift' -
1160 * factor * shift * (gamma - 1) * b * [1, 0, 0]^T
1161 *
1162 */
1163 const auto b = Number(eos_interpolation_b());
1164 const auto pinf = Number(eos_interpolation_pinfty());
1165 const auto q = Number(eos_interpolation_q());
1166
1167 const auto rho = density(U);
1168 const auto m = momentum(U);
1169 const auto E = total_energy(U);
1170
1171 const auto covolume = Number(1.) - b * rho;
1172 const auto covolume_inverse = ScalarNumber(1.) / covolume;
1173
1174 const auto shift = rho * E - ScalarNumber(0.5) * m.norm_square() -
1175 rho * rho * q - rho * pinf * covolume;
1176
1177 constexpr auto eps = std::numeric_limits<ScalarNumber>::epsilon();
1178 const auto regularization = m.norm() * eps;
1179 const auto max_val = ryujin::pow(
1180 std::max(regularization, eta * covolume_inverse), gamma_min);
1181 auto factor = safe_division(Number(1.0), max_val);
1182 factor *= fixed_power<2>(covolume_inverse) / (gamma_min + Number(1.));
1183
1184 state_type result;
1185
1186 const auto first_term = E - ScalarNumber(2.) * rho * q -
1187 pinf * (Number(1.) - ScalarNumber(2.) * b * rho);
1188 const auto second_term = -(gamma_min - Number(1.)) * shift * b;
1189
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;
1194
1195 return result;
1196 }
1197
1198
1199 template <int dim, typename Number>
1200 DEAL_II_ALWAYS_INLINE inline Number
1202 const Number &p) const
1203 {
1204 const auto b = Number(eos_interpolation_b());
1205 const auto pinf = Number(eos_interpolation_pinfty());
1206 const auto q = Number(eos_interpolation_q());
1207
1208 const auto rho = density(U);
1209 const auto rho_e = internal_energy(U);
1210 const auto covolume = Number(1.) - b * rho;
1211
1212 const auto numerator = (p + pinf) * covolume;
1213 const auto denominator = rho_e - rho * q - covolume * pinf;
1214 return Number(1.) + safe_division(numerator, denominator);
1215 }
1216
1217
1218 template <int dim, typename Number>
1219 DEAL_II_ALWAYS_INLINE inline Number
1221 const state_type &U, const Number &gamma) const
1222 {
1223 const auto b = Number(eos_interpolation_b());
1224 const auto pinf = Number(eos_interpolation_pinfty());
1225 const auto q = Number(eos_interpolation_q());
1226
1227 const auto rho = density(U);
1228 const auto rho_e = internal_energy(U);
1229 const auto covolume = Number(1.) - b * rho;
1230
1231 return positive_part(gamma - Number(1.)) *
1232 safe_division(rho_e - rho * q, covolume) -
1233 gamma * pinf;
1234 }
1235
1236
1237 template <int dim, typename Number>
1238 DEAL_II_ALWAYS_INLINE inline Number
1240 const state_type &U, const Number &gamma) const
1241 {
1242 const auto b = Number(eos_interpolation_b());
1243 const auto pinf = Number(eos_interpolation_pinfty());
1244 const auto q = Number(eos_interpolation_q());
1245
1246 const auto rho = density(U);
1247 const auto rho_e = internal_energy(U);
1248 const auto covolume = Number(1.) - b * rho;
1249
1250 auto radicand =
1251 (rho_e - rho * q - pinf * covolume) / (covolume * covolume * rho);
1252 radicand *= gamma * (gamma - 1.);
1253 return std::sqrt(positive_part(radicand));
1254 }
1255
1256
1257 template <int dim, typename Number>
1258 DEAL_II_ALWAYS_INLINE inline bool
1260 {
1261 const auto b = Number(eos_interpolation_b());
1262 const auto pinf = Number(eos_interpolation_pinfty());
1263 const auto q = Number(eos_interpolation_q());
1264
1265 const auto rho = density(U);
1266 const auto rho_e = internal_energy(U);
1267 const auto covolume = Number(1.) - b * rho;
1268
1269 const auto shift = rho_e - rho * q - pinf * covolume;
1270
1271 constexpr auto gt = dealii::SIMDComparison::greater_than;
1272 using T = Number;
1273 const auto test =
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.));
1276
1277#ifdef DEBUG_OUTPUT
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";
1283 }
1284#endif
1285
1286 return (test == Number(0.));
1287 }
1288
1289
1290 template <int dim, typename Number>
1291 template <int component>
1292 DEAL_II_ALWAYS_INLINE inline auto
1294 const state_type &U,
1295 const Number &p,
1296 const state_type &U_bar,
1297 const Number &p_bar,
1298 const dealii::Tensor<1, dim, Number> &normal) const -> state_type
1299 {
1300 static_assert(component == 1 || component == 2,
1301 "component has to be 1 or 2");
1302
1303 const auto b = Number(eos_interpolation_b());
1304 const auto pinf = Number(eos_interpolation_pinfty());
1305 const auto q = Number(eos_interpolation_q());
1306
1307 /*
1308 * The "four" Riemann characteristics are formed under the assumption
1309 * of a locally isentropic flow. For this, we first transform both
1310 * states into {rho, vn, vperp, gamma, a}, where we use the NASG EOS
1311 * interpolation to derive a surrogate gamma and speed of sound a.
1312 *
1313 * See, e.g., https://arxiv.org/pdf/2004.08750, "Compressible flow in
1314 * a NOble-Abel Stiffened-Gas fluid", M. I. Radulescu.
1315 */
1316
1317 const auto m = momentum(U);
1318 const auto rho = density(U);
1319 const auto vn = m * normal / rho;
1320
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;
1324
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;
1328
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;
1332
1333 /*
1334 * Now compute the Riemann characteristics {R_1, R_2, vperp, s}:
1335 * R_1 = v * n - 2 / (gamma - 1) * a * (1 - b * rho)
1336 * R_2 = v * n + 2 / (gamma - 1) * a * (1 - b * rho)
1337 * vperp
1338 * S = (p + p_infty) / rho^gamma * (1 - b * rho)^gamma
1339 *
1340 * Here, we replace either R_1, or R_2 with values coming from U_bar:
1341 */
1342
1343 const auto R_1 =
1344 component == 1 ? vn_bar - 2. * a_bar / (gamma_bar - 1.) * covolume_bar
1345 : vn - 2. * a / (gamma - 1.) * covolume;
1346
1347 const auto R_2 =
1348 component == 2 ? vn_bar + 2. * a_bar / (gamma_bar - 1.) * covolume_bar
1349 : vn + 2. * a / (gamma - 1.) * covolume;
1350
1351 /*
1352 * Note that we are really hoping for the best here... We require
1353 * that R_2 >= R_1 so that we can extract a valid sound speed...
1354 */
1355
1356 Assert(
1357 R_2 >= R_1,
1358 dealii::ExcMessage("Encountered R_2 < R_1 in dynamic boundary value "
1359 "enforcement. This implies that the interpolation "
1360 "with Riemann characteristics failed."));
1361
1362 const auto vperp = m / rho - vn * normal;
1363
1364 const auto S = (p + pinf) * ryujin::pow(Number(1.) / rho - b, gamma);
1365
1366 /*
1367 * Now, we have to reconstruct the actual conserved state U from the
1368 * Riemann characteristics R_1, R_2, vperp, and s. We first set up
1369 * {vn_new, vperp_new, a_new, S} and then solve for {rho_new, p_new}
1370 * with the help of the NASG EOS surrogate formulas:
1371 *
1372 * S = (p + p_infty) / rho^gamma * (1 - b * rho)^gamma
1373 *
1374 * a^2 = gamma * (p + p_infty) / (rho * cov)
1375 *
1376 * This implies:
1377 *
1378 * a^2 / (gamma * S) = rho^{gamma - 1} / (1 - b * rho)^{1 + gamma}
1379 */
1380
1381 const auto vn_new = Number(0.5) * (R_1 + R_2);
1382
1383 /*
1384 * Technically, we would need to solve for rho subject to a number of
1385 * nonlinear relationships:
1386 *
1387 * a = (gamma - 1) * (R_2 - R_1) / (4. * (1 - b * rho))
1388 *
1389 * a^2 / (gamma * S) = rho^{gamma - 1} / (1 - b * rho)^{gamma + 1}
1390 *
1391 * This seems to be a bit expensive for the fact that our dynamic
1392 * boundary conditions are already terribly heuristic...
1393 *
1394 * So instead, we rewrite this system as:
1395 *
1396 * a * (1 - b * rho) = (gamma - 1) * (R_2 - R_1) / 4.
1397 *
1398 * a^2 / (gamma * S) (1 - b * rho)^2
1399 * = (rho / (1 - b * rho))^{gamma - 1}
1400 *
1401 * And compute the terms on the left simply with the old covolume and
1402 * solving an easier easier nonlinear equation for the density. The
1403 * resulting system reads:
1404 *
1405 * a = (gamma - 1) * (R_2 - R_1) / (4. * (1 - b * rho_old))
1406 * A = {a^2 / (gamma * S) (1 - b * rho_old)^{2 gamma}}^{1/(gamma - 1)}
1407 *
1408 * rho = A / (1 + b * A)
1409 */
1410
1411 const auto a_new_square =
1412 ryujin::fixed_power<2>((gamma - 1.) * (R_2 - R_1) / (4. * covolume));
1413
1414 auto term = ryujin::pow(a_new_square / (gamma * S), 1. / (gamma - 1.));
1415 if (b != ScalarNumber(0.)) {
1416 term *= std::pow(covolume, 2. / (gamma - 1.));
1417 }
1418
1419 const auto rho_new = term / (1. + b * term);
1420
1421 const auto covolume_new = (1. - b * rho_new);
1422 const auto p_new = a_new_square / gamma * rho_new * covolume_new - pinf;
1423
1424 /*
1425 * And translate back into conserved quantities:
1426 */
1427
1428 const auto rho_e_new =
1429 rho_new * q + (p_new + gamma * pinf) * covolume_new / (gamma - 1.);
1430
1431 state_type U_new;
1432 U_new[0] = rho_new;
1433 for (unsigned int d = 0; d < dim; ++d) {
1434 U_new[1 + d] = rho_new * (vn_new * normal + vperp)[d];
1435 }
1436 U_new[1 + dim] =
1437 rho_e_new + 0.5 * rho_new * (vn_new * vn_new + vperp.norm_square());
1438
1439 return U_new;
1440 }
1441
1442
1443 template <int dim, typename Number>
1444 template <typename Lambda>
1445 DEAL_II_ALWAYS_INLINE inline auto
1447 dealii::types::boundary_id id,
1448 const state_type &U,
1449 const dealii::Tensor<1, dim, Number> &normal,
1450 const Lambda &get_dirichlet_data) const -> state_type
1451 {
1452 state_type result = U;
1453
1454 if (id == Boundary::dirichlet) {
1455 result = get_dirichlet_data();
1456
1457 } else if (id == Boundary::dirichlet_momentum) {
1458 /* Only enforce Dirichlet conditions on the momentum: */
1459 auto m_dirichlet = momentum(get_dirichlet_data());
1460 for (unsigned int k = 0; k < dim; ++k)
1461 result[k + 1] = m_dirichlet[k];
1462
1463 } else if (id == Boundary::slip) {
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];
1468
1469 } else if (id == Boundary::no_slip) {
1470 for (unsigned int k = 0; k < dim; ++k)
1471 result[k + 1] = Number(0.);
1472
1473 } else if (id == Boundary::dynamic) {
1474 /*
1475 * On dynamic boundary conditions, we distinguish four cases:
1476 *
1477 * - supersonic inflow: prescribe full state
1478 * - subsonic inflow:
1479 * decompose into Riemann invariants and leave R_2
1480 * characteristic untouched.
1481 * - supersonic outflow: do nothing
1482 * - subsonic outflow:
1483 * decompose into Riemann invariants and prescribe incoming
1484 * R_1 characteristic.
1485 */
1486 const auto m = momentum(U);
1487 const auto rho = density(U);
1488 const auto rho_e = internal_energy(U);
1489
1490 /*
1491 * We do not have precomputed values available. Thus, simply query
1492 * the pressure oracle and compute a surrogate speed of sound from
1493 * there:
1494 */
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;
1499
1500 /* Supersonic inflow: */
1501 if (vn < -a) {
1502 result = get_dirichlet_data();
1503 }
1504
1505 /* Subsonic inflow: */
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);
1512
1513 result = prescribe_riemann_characteristic<2>(
1514 U_dirichlet, p_dirichlet, U, p, normal);
1515 }
1516
1517 /* Subsonic outflow: */
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);
1524
1525 result = prescribe_riemann_characteristic<1>(
1526 U, p, U_dirichlet, p_dirichlet, normal);
1527 }
1528 /* Supersonic outflow: do nothing, i.e., keep U as is */
1529
1530 } else {
1531 AssertThrow(false, dealii::ExcNotImplemented());
1532 }
1533
1534 return result;
1535 }
1536
1537
1538 template <int dim, typename Number>
1539 DEAL_II_ALWAYS_INLINE inline auto
1541 const Number &p) const -> flux_type
1542 {
1543 const auto rho_inverse = ScalarNumber(1.) / density(U);
1544 const auto m = momentum(U);
1545 const auto E = total_energy(U);
1546
1547 flux_type result;
1548
1549 result[0] = m;
1550 for (unsigned int i = 0; i < dim; ++i) {
1551 result[1 + i] = m * (m[i] * rho_inverse);
1552 result[1 + i][i] += p;
1553 }
1554 result[dim + 1] = m * (rho_inverse * (E + p));
1555
1556 return result;
1557 }
1558
1559
1560 template <int dim, typename Number>
1561 DEAL_II_ALWAYS_INLINE inline auto
1563 const PrecomputedVector &pv,
1564 const InitialPrecomputedVector & /*piv*/,
1565 const unsigned int i,
1566 const state_type &U_i) const -> flux_contribution_type
1567 {
1568 const auto &[p_i, gamma_min_i, s_i, eta_i] =
1569 pv.template get_tensor<Number, precomputed_type>(i);
1570 return f(U_i, p_i);
1571 }
1572
1573
1574 template <int dim, typename Number>
1575 DEAL_II_ALWAYS_INLINE inline auto
1577 const PrecomputedVector &pv,
1578 const InitialPrecomputedVector & /*piv*/,
1579 const unsigned int *js,
1580 const state_type &U_j) const -> flux_contribution_type
1581 {
1582 const auto &[p_j, gamma_min_j, s_j, eta_j] =
1583 pv.template get_tensor<Number, precomputed_type>(js);
1584 return f(U_j, p_j);
1585 }
1586
1587
1588 template <int dim, typename Number>
1589 DEAL_II_ALWAYS_INLINE inline auto
1591 const flux_contribution_type &flux_i,
1592 const flux_contribution_type &flux_j,
1593 const dealii::Tensor<1, dim, Number> &c_ij) const -> state_type
1594 {
1595 return -contract(add(flux_i, flux_j), c_ij);
1596 }
1597
1598
1599 template <int dim, typename Number>
1600 template <typename ST>
1602 -> state_type
1603 {
1604 using T = typename ST::value_type;
1605 static_assert(std::is_same_v<Number, T>, "template mismatch");
1606
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");
1611
1612 state_type result;
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];
1617
1618 return result;
1619 }
1620
1621
1622 template <int dim, typename Number>
1623 template <typename ST>
1624 DEAL_II_ALWAYS_INLINE inline auto
1626 const ST &initial_state) const -> state_type
1627 {
1628 auto primitive_state = expand_state(initial_state);
1629
1630 /* pressure into specific internal energy: */
1631 const auto rho = density(primitive_state);
1632 const auto p = /*SIC!*/ total_energy(primitive_state);
1633 const auto e = eos_specific_internal_energy(rho, p);
1634 primitive_state[dim + 1] = e;
1635
1636 return from_primitive_state(primitive_state);
1637 }
1638
1639
1640 template <int dim, typename Number>
1641 DEAL_II_ALWAYS_INLINE inline auto
1643 const state_type &primitive_state) const -> state_type
1644 {
1645 const auto rho = density(primitive_state);
1646 /* extract velocity: */
1647 const auto u = /*SIC!*/ momentum(primitive_state);
1648 /* extract specific internal energy: */
1649 const auto &e = /*SIC!*/ total_energy(primitive_state);
1650
1651 auto state = primitive_state;
1652 /* Fix up momentum: */
1653 for (unsigned int i = 1; i < dim + 1; ++i)
1654 state[i] *= rho;
1655
1656 /* Compute total energy: */
1657 state[dim + 1] = rho * e + Number(0.5) * rho * u * u;
1658
1659 return state;
1660 }
1661
1662
1663 template <int dim, typename Number>
1664 DEAL_II_ALWAYS_INLINE inline auto
1666 const state_type &state) const -> state_type
1667 {
1668 const auto rho = density(state);
1669 const auto rho_inverse = Number(1.) / rho;
1670 const auto rho_e = internal_energy(state);
1671
1672 auto primitive_state = state;
1673 /* Fix up velocity: */
1674 for (unsigned int i = 1; i < dim + 1; ++i)
1675 primitive_state[i] *= rho_inverse;
1676 /* Set specific internal energy: */
1677 primitive_state[dim + 1] = rho_e * rho_inverse;
1678
1679 return primitive_state;
1680 }
1681
1682
1683 template <int dim, typename Number>
1684 template <typename Lambda>
1686 const state_type &state, const Lambda &lambda) const -> state_type
1687 {
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];
1692 return result;
1693 }
1694 } // namespace EulerAEOS
1695} // namespace ryujin
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
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
state_type from_primitive_state(const state_type &primitive_state) const
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
DEAL_II_ALWAYS_INLINE Number eos_pressure(const Number &rho, const Number &e) const
static constexpr unsigned int problem_dimension
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
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
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
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")
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)
@ dirichlet_momentum
#define RYUJIN_OMP_FOR
Definition: openmp.h:70
#define RYUJIN_OMP_SINGLE
Definition: openmp.h:100
T pow(const T x, const T b)
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)
Definition: simd.h:112
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)