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