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
8#include <compile_time_options.h>
10#include <discretization.h>
12#include <openmp.h>
13#include <patterns_conversion.h>
14#include <simd.h>
15#include <state_vector.h>
16
17#include <deal.II/base/parameter_acceptor.h>
18#include <deal.II/base/tensor.h>
19
20#include <array>
21
22namespace ryujin
23{
24 namespace Euler
25 {
26 template <int dim, typename Number>
27 class HyperbolicSystemView;
28
39 class HyperbolicSystem final : public dealii::ParameterAcceptor
40 {
41 public:
45 static inline const std::string problem_name =
46 "Compressible Euler equations (polytropic gas EOS, optimized)";
47
51 HyperbolicSystem(const std::string &subsection = "/HyperbolicSystem");
52
59 template <int dim, typename Number>
60 auto view() const
61 {
63 }
64
65 private:
70 double gamma_;
71
72 double reference_density_;
73 double vacuum_state_relaxation_small_;
74 double vacuum_state_relaxation_large_;
75
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_;
80
81 template <int dim, typename Number>
84 }; /* HyperbolicSystem */
85
86
103 template <int dim, typename Number>
105 {
106 public:
111 HyperbolicSystemView(const HyperbolicSystem &hyperbolic_system)
112 : hyperbolic_system_(hyperbolic_system)
113 {
114 }
115
119 template <int dim2, typename Number2>
120 auto view() const
121 {
122 return HyperbolicSystemView<dim2, Number2>{hyperbolic_system_};
123 }
124
129
134
135 DEAL_II_ALWAYS_INLINE inline ScalarNumber gamma() const
136 {
137 return hyperbolic_system_.gamma_;
138 }
139
140 DEAL_II_ALWAYS_INLINE inline ScalarNumber reference_density() const
141 {
142 return hyperbolic_system_.reference_density_;
143 }
144
145 DEAL_II_ALWAYS_INLINE inline ScalarNumber
147 {
148 return hyperbolic_system_.vacuum_state_relaxation_small_;
149 }
150
151 DEAL_II_ALWAYS_INLINE inline ScalarNumber
153 {
154 return hyperbolic_system_.vacuum_state_relaxation_large_;
155 }
156
158
166
167 DEAL_II_ALWAYS_INLINE inline ScalarNumber gamma_inverse() const
168 {
169 return ScalarNumber(hyperbolic_system_.gamma_inverse_);
170 }
171
172 DEAL_II_ALWAYS_INLINE inline ScalarNumber gamma_plus_one_inverse() const
173 {
174 return ScalarNumber(hyperbolic_system_.gamma_plus_one_inverse_);
175 }
176
177 DEAL_II_ALWAYS_INLINE inline ScalarNumber gamma_minus_one_inverse() const
178 {
179 return ScalarNumber(hyperbolic_system_.gamma_minus_one_inverse_);
180 }
181
182 DEAL_II_ALWAYS_INLINE inline ScalarNumber
184 {
185 return ScalarNumber(
186 hyperbolic_system_.gamma_minus_one_over_gamma_plus_one_);
187 }
188
192 static constexpr bool have_gamma = true;
193
197 static constexpr bool have_eos_interpolation_b = false;
198
200
204
205 private:
206 const HyperbolicSystem &hyperbolic_system_;
207
208 public:
210
214
218 static constexpr unsigned int problem_dimension = 2 + dim;
219
223 using state_type = dealii::Tensor<1, problem_dimension, Number>;
224
228 using flux_type =
229 dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim, Number>>;
230
235
240 static inline const auto component_names =
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"};
248 __builtin_trap();
249 }();
250
255 static inline const auto primitive_component_names =
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"};
263 __builtin_trap();
264 }();
265
269 static constexpr unsigned int n_precomputed_values = 2;
270
274 using precomputed_type = std::array<Number, n_precomputed_values>;
275
279 static inline const auto precomputed_names =
280 std::array<std::string, n_precomputed_values>{"s", "eta_h"};
281
285 static constexpr unsigned int n_initial_precomputed_values = 0;
286
291 std::array<Number, n_initial_precomputed_values>;
292
296 static inline const auto initial_precomputed_names =
297 std::array<std::string, n_initial_precomputed_values>{};
298
302 using StateVector = Vectors::
303 StateVector<ScalarNumber, problem_dimension, n_precomputed_values>;
304
310
316
324
326
330
334 static constexpr unsigned int n_precomputation_cycles = 1;
335
340 template <typename DISPATCH, typename SPARSITY>
341 void precomputation_loop(unsigned int cycle,
342 const DISPATCH &dispatch_check,
343 const SPARSITY &sparsity_simd,
344 StateVector &state_vector,
345 unsigned int left,
346 unsigned int right,
347 const bool skip_constrained_dofs = true) const;
348
350
354
359 static Number density(const state_type &U);
360
367 Number filter_vacuum_density(const Number &rho) const;
368
373 static dealii::Tensor<1, dim, Number> momentum(const state_type &U);
374
379 static Number total_energy(const state_type &U);
380
385 static Number internal_energy(const state_type &U);
386
393
404 Number pressure(const state_type &U) const;
405
413 Number speed_of_sound(const state_type &U) const;
414
422 Number specific_entropy(const state_type &U) const;
423
431 Number harten_entropy(const state_type &U) const;
432
441
446 Number mathematical_entropy(const state_type &U) const;
447
454
460 bool is_admissible(const state_type &U) const;
461
463
467
473 template <int component>
474 std::array<state_type, 2> linearized_eigenvector(
475 const state_type &U,
476 const dealii::Tensor<1, dim, Number> &normal) const;
477
484 template <int component>
486 const state_type &U,
487 const state_type &U_bar,
488 const dealii::Tensor<1, dim, Number> &normal) const;
489
508 template <typename Lambda>
510 apply_boundary_conditions(const dealii::types::boundary_id id,
511 const state_type &U,
512 const dealii::Tensor<1, dim, Number> &normal,
513 const Lambda &get_dirichlet_data) const;
514
516
520
531 flux_type f(const state_type &U) const;
532
554 const InitialPrecomputedVector &ipv,
555 const unsigned int i,
556 const state_type &U_i) const;
557
560 const InitialPrecomputedVector &ipv,
561 const unsigned int *js,
562 const state_type &U_j) const;
563
570 const flux_contribution_type &flux_j,
571 const dealii::Tensor<1, dim, Number> &c_ij) const;
572
574 static constexpr bool have_high_order_flux = false;
575
577 const flux_contribution_type &flux_i,
578 const flux_contribution_type &flux_j,
579 const dealii::Tensor<1, dim, Number> &c_ij) const = delete;
580
582
586
588 static constexpr bool have_source_terms = false;
589
591 const unsigned int i,
592 const state_type &U_i,
593 const ScalarNumber tau) const = delete;
594
596 const unsigned int *js,
597 const state_type &U_j,
598 const ScalarNumber tau) const = delete;
599
601
605
616 template <typename ST>
617 state_type expand_state(const ST &state) const;
618
632 template <typename ST>
633 state_type from_initial_state(const ST &initial_state) const;
634
639 state_type from_primitive_state(const state_type &primitive_state) const;
640
645 state_type to_primitive_state(const state_type &state) const;
646
652 template <typename Lambda>
654 const Lambda &lambda) const;
656 }; /* HyperbolicSystemView */
657
658
659 /*
660 * -------------------------------------------------------------------------
661 * Inline definitions
662 * -------------------------------------------------------------------------
663 */
664
665
666 inline HyperbolicSystem::HyperbolicSystem(const std::string &subsection)
667 : ParameterAcceptor(subsection)
668 {
669 gamma_ = 7. / 5.;
670 add_parameter("gamma", gamma_, "The ratio of specific heats");
671
672 reference_density_ = 1.;
673 add_parameter("reference density",
674 reference_density_,
675 "Problem specific density reference");
676
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");
681
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");
686
687 /*
688 * Precompute a number of derived gamma coefficients that contain
689 * divisions:
690 */
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.);
696 };
697
698 compute_inverses();
699 ParameterAcceptor::parse_parameters_call_back.connect(compute_inverses);
700 }
701
702
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,
710 StateVector &state_vector,
711 unsigned int left,
712 unsigned int right,
713 const bool skip_constrained_dofs /*= true*/) const
714 {
715 Assert(cycle == 0, dealii::ExcInternalError());
716
717 const auto &U = std::get<0>(state_vector);
718 auto &precomputed = std::get<1>(state_vector);
719
720 /* We are inside a thread parallel context */
721
722 unsigned int stride_size = get_stride_size<Number>;
723
725 for (unsigned int i = left; i < right; i += stride_size) {
726
727 /* Skip constrained degrees of freedom: */
728 const unsigned int row_length = sparsity_simd.row_length(i);
729 if (skip_constrained_dofs && row_length == 1)
730 continue;
731
732 dispatch_check(i);
733
734 const auto U_i = U.template get_tensor<Number>(i);
735 const precomputed_type prec_i{specific_entropy(U_i),
736 harten_entropy(U_i)};
737 precomputed.template write_tensor<Number>(prec_i, i);
738 }
739 }
740
741
742 template <int dim, typename Number>
743 DEAL_II_ALWAYS_INLINE inline Number
745 {
746 return U[0];
747 }
748
749
750 template <int dim, typename Number>
751 DEAL_II_ALWAYS_INLINE inline Number
753 const Number &rho) const
754 {
755 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
756 const Number rho_cutoff_large =
757 reference_density() * vacuum_state_relaxation_large() * eps;
758
759 return dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
760 std::abs(rho), rho_cutoff_large, Number(0.), rho);
761 }
762
763
764 template <int dim, typename Number>
765 DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim, Number>
767 {
768 dealii::Tensor<1, dim, Number> result;
769 for (unsigned int i = 0; i < dim; ++i)
770 result[i] = U[1 + i];
771 return result;
772 }
773
774
775 template <int dim, typename Number>
776 DEAL_II_ALWAYS_INLINE inline Number
778 {
779 return U[1 + dim];
780 }
781
782
783 template <int dim, typename Number>
784 DEAL_II_ALWAYS_INLINE inline Number
786 {
787 /*
788 * rho e = (E - 1/2*m^2/rho)
789 */
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;
794 }
795
796
797 template <int dim, typename Number>
798 DEAL_II_ALWAYS_INLINE inline auto
800 const state_type &U) -> state_type
801 {
802 /*
803 * With
804 * rho e = E - 1/2 |m|^2 / rho
805 * we get
806 * (rho e)' = (1/2m^2/rho^2, -m/rho , 1 )^T
807 */
808
809 const Number rho_inverse = ScalarNumber(1.) / density(U);
810 const auto u = momentum(U) * rho_inverse;
811
812 state_type result;
813
814 result[0] = ScalarNumber(0.5) * u.norm_square();
815 for (unsigned int i = 0; i < dim; ++i) {
816 result[1 + i] = -u[i];
817 }
818 result[dim + 1] = ScalarNumber(1.);
819
820 return result;
821 }
822
823
824 template <int dim, typename Number>
825 DEAL_II_ALWAYS_INLINE inline Number
827 {
828 /* p = (gamma - 1) * (rho e) */
829 return (gamma() - ScalarNumber(1.)) * internal_energy(U);
830 }
831
832
833 template <int dim, typename Number>
834 DEAL_II_ALWAYS_INLINE inline Number
836 {
837 /* c^2 = gamma * p / rho */
838 const Number rho_inverse = ScalarNumber(1.) / density(U);
839 const Number p = pressure(U);
840 return std::sqrt(gamma() * p * rho_inverse);
841 }
842
843
844 template <int dim, typename Number>
845 DEAL_II_ALWAYS_INLINE inline Number
847 const state_type &U) const
848 {
849 /* exp((gamma - 1)s) = (rho e) / rho ^ gamma */
850 const auto rho_inverse = ScalarNumber(1.) / density(U);
851 return internal_energy(U) * ryujin::pow(rho_inverse, gamma());
852 }
853
854
855 template <int dim, typename Number>
856 DEAL_II_ALWAYS_INLINE inline Number
858 {
859 /* rho^2 e = \rho E - 1/2*m^2 */
860
861 const Number rho = density(U);
862 const auto m = momentum(U);
863 const Number E = total_energy(U);
864
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());
867 }
868
869
870 template <int dim, typename Number>
871 DEAL_II_ALWAYS_INLINE inline auto
873 const state_type &U) const -> state_type
874 {
875 /*
876 * With
877 * eta = (rho^2 e) ^ 1/(gamma+1)
878 * rho^2 e = rho * E - 1/2 |m|^2
879 *
880 * we get
881 *
882 * eta' = 1/(gamma+1) * (rho^2 e) ^ -gamma/(gamma+1) * (E,-m,rho)^T
883 *
884 */
885
886 const Number rho = density(U);
887 const auto m = momentum(U);
888 const Number E = total_energy(U);
889
890 const Number rho_rho_e = rho * E - ScalarNumber(0.5) * m.norm_square();
891
892 const auto factor =
893 gamma_plus_one_inverse() *
894 ryujin::pow(rho_rho_e, -gamma() * gamma_plus_one_inverse());
895
896 state_type result;
897
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;
902
903 return result;
904 }
905
906
907 template <int dim, typename Number>
908 DEAL_II_ALWAYS_INLINE inline Number
910 const state_type &U) const
911 {
913 const auto p = pressure(U);
914 return ryujin::pow(p, gamma_inverse());
915 }
916
917
918 template <int dim, typename Number>
919 DEAL_II_ALWAYS_INLINE inline auto
921 const state_type &U) const -> state_type
922 {
923 /*
924 * With
925 * eta = p ^ (1/gamma)
926 * p = (gamma - 1) * (rho e)
927 * rho e = E - 1/2 |m|^2 / rho
928 *
929 * we get
930 *
931 * eta' = (gamma - 1)/gamma p ^(1/gamma - 1) *
932 *
933 * (1/2m^2/rho^2 , -m/rho , 1 )^T
934 */
935 const Number rho = density(U);
936 const Number rho_inverse = ScalarNumber(1.) / rho;
937 const auto u = momentum(U) * rho_inverse;
938 const auto p = pressure(U);
939
940 const auto factor = (gamma() - ScalarNumber(1.0)) * gamma_inverse() *
941 ryujin::pow(p, gamma_inverse() - ScalarNumber(1.));
942
943 state_type result;
944
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];
949 }
950
951 return result;
952 }
953
954
955 template <int dim, typename Number>
956 DEAL_II_ALWAYS_INLINE inline bool
958 {
959 const auto rho_new = density(U);
960 const auto e_new = internal_energy(U);
961 const auto s_new = specific_entropy(U);
962
963 constexpr auto gt = dealii::SIMDComparison::greater_than;
964 using T = Number;
965 const auto test =
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.));
969
970#ifdef DEBUG_OUTPUT
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;
977 }
978#endif
979
980 return (test == Number(0.));
981 }
982
983
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>
990 {
991 static_assert(component == 1 || component == problem_dimension,
992 "Only first and last eigenvectors implemented");
993
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();
999
1000 state_type b;
1001 state_type c;
1002
1003 const auto e_k = 0.5 * v.norm_square();
1004
1005 switch (component) {
1006 case 1:
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.;
1011 b /= 2. * a * a;
1012
1013 c[0] = 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);
1017
1018 return {b, c};
1019
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.;
1025 b /= 2. * a * a;
1026
1027 c[0] = 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);
1031
1032 return {b, c};
1033 }
1034
1035 __builtin_unreachable();
1036 }
1037
1038
1039 template <int dim, typename Number>
1040 template <int component>
1041 DEAL_II_ALWAYS_INLINE inline auto
1043 const state_type &U,
1044 const state_type &U_bar,
1045 const dealii::Tensor<1, dim, Number> &normal) const -> state_type
1046 {
1047 static_assert(component == 1 || component == 2,
1048 "component has to be 1 or 2");
1049
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;
1054
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;
1059
1060 /* First Riemann characteristic: v* n - 2 / (gamma - 1) * a */
1061
1062 const auto R_1 = component == 1
1063 ? vn_bar - 2. * a_bar / (gamma() - ScalarNumber(1.))
1064 : vn - 2. * a / (gamma() - ScalarNumber(1.));
1065
1066 /* Second Riemann characteristic: v* n + 2 / (gamma() - 1) * a */
1067
1068 const auto R_2 = component == 2
1069 ? vn_bar + 2. * a_bar / (gamma() - ScalarNumber(1.))
1070 : vn + 2. * a / (gamma() - ScalarNumber(1.));
1071
1072 const auto p = pressure(U);
1073 const auto s = p / ryujin::pow(rho, gamma());
1074
1075 const auto vperp = m / rho - vn * normal;
1076
1077 const auto vn_new = 0.5 * (R_1 + R_2);
1078
1079 auto rho_new = 1. / (gamma() * s) *
1080 ryujin::fixed_power<2>(ScalarNumber((gamma() - 1.) / 4.) *
1081 (R_2 - R_1));
1082 rho_new = ryujin::pow(rho_new, 1. / (gamma() - 1.));
1083
1084 const auto p_new = s * std::pow(rho_new, gamma());
1085
1086 state_type U_new;
1087 U_new[0] = rho_new;
1088 for (unsigned int d = 0; d < dim; ++d) {
1089 U_new[1 + d] = rho_new * (vn_new * normal + vperp)[d];
1090 }
1091 U_new[1 + dim] = p_new / ScalarNumber(gamma() - 1.) +
1092 0.5 * rho_new * (vn_new * vn_new + vperp.norm_square());
1093
1094 return U_new;
1095 }
1096
1097
1098 template <int dim, typename Number>
1099 template <typename Lambda>
1100 DEAL_II_ALWAYS_INLINE inline auto
1102 dealii::types::boundary_id id,
1103 const state_type &U,
1104 const dealii::Tensor<1, dim, Number> &normal,
1105 const Lambda &get_dirichlet_data) const -> state_type
1106 {
1107 state_type result = U;
1108
1109 if (id == Boundary::dirichlet) {
1110 result = get_dirichlet_data();
1111
1112 } else if (id == Boundary::dirichlet_momentum) {
1113 /* Only enforce Dirichlet conditions on the momentum: */
1114 auto m_dirichlet = momentum(get_dirichlet_data());
1115 for (unsigned int k = 0; k < dim; ++k)
1116 result[k + 1] = m_dirichlet[k];
1117
1118 } else if (id == Boundary::slip) {
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];
1123
1124 } else if (id == Boundary::no_slip) {
1125 for (unsigned int k = 0; k < dim; ++k)
1126 result[k + 1] = Number(0.);
1127
1128 } else if (id == Boundary::dynamic) {
1129 /*
1130 * On dynamic boundary conditions, we distinguish four cases:
1131 *
1132 * - supersonic inflow: prescribe full state
1133 * - subsonic inflow:
1134 * decompose into Riemann invariants and leave R_2
1135 * characteristic untouched.
1136 * - supersonic outflow: do nothing
1137 * - subsonic outflow:
1138 * decompose into Riemann invariants and prescribe incoming
1139 * R_1 characteristic.
1140 */
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;
1145
1146 /* Supersonic inflow: */
1147 if (vn < -a) {
1148 result = get_dirichlet_data();
1149 }
1150
1151 /* Subsonic inflow: */
1152 if (vn >= -a && vn <= 0.) {
1153 const auto U_dirichlet = get_dirichlet_data();
1154 result = prescribe_riemann_characteristic<2>(U_dirichlet, U, normal);
1155 }
1156
1157 /* Subsonic outflow: */
1158 if (vn > 0. && vn <= a) {
1159 const auto U_dirichlet = get_dirichlet_data();
1160 result = prescribe_riemann_characteristic<1>(U, U_dirichlet, normal);
1161 }
1162
1163 /* Supersonic outflow: do nothing, i.e., keep U as is */
1164
1165 } else {
1166 AssertThrow(false, dealii::ExcNotImplemented());
1167 }
1168
1169 return result;
1170 }
1171
1172
1173 template <int dim, typename Number>
1174 DEAL_II_ALWAYS_INLINE inline auto
1176 {
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);
1181
1182 flux_type result;
1183
1184 result[0] = m;
1185 for (unsigned int i = 0; i < dim; ++i) {
1186 result[1 + i] = m * (m[i] * rho_inverse);
1187 result[1 + i][i] += p;
1188 }
1189 result[dim + 1] = m * (rho_inverse * (E + p));
1190
1191 return result;
1192 }
1193
1194
1195 template <int dim, typename Number>
1196 DEAL_II_ALWAYS_INLINE inline auto
1198 const PrecomputedVector & /*pv*/,
1199 const InitialPrecomputedVector & /*ipv*/,
1200 const unsigned int /*i*/,
1201 const state_type &U_i) const -> flux_contribution_type
1202 {
1203 return f(U_i);
1204 }
1205
1206
1207 template <int dim, typename Number>
1208 DEAL_II_ALWAYS_INLINE inline auto
1210 const PrecomputedVector & /*pv*/,
1211 const InitialPrecomputedVector & /*ipv*/,
1212 const unsigned int * /*js*/,
1213 const state_type &U_j) const -> flux_contribution_type
1214 {
1215 return f(U_j);
1216 }
1217
1218
1219 template <int dim, typename Number>
1220 DEAL_II_ALWAYS_INLINE inline auto
1222 const flux_contribution_type &flux_i,
1223 const flux_contribution_type &flux_j,
1224 const dealii::Tensor<1, dim, Number> &c_ij) const -> state_type
1225 {
1226 return -contract(add(flux_i, flux_j), c_ij);
1227 }
1228
1229
1230 template <int dim, typename Number>
1231 template <typename ST>
1233 -> state_type
1234 {
1235 using T = typename ST::value_type;
1236 static_assert(std::is_same_v<Number, T>, "template mismatch");
1237
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");
1242
1243 state_type result;
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];
1248
1249 return result;
1250 }
1251
1252
1253 template <int dim, typename Number>
1254 template <typename ST>
1255 DEAL_II_ALWAYS_INLINE inline auto
1257 const ST &initial_state) const -> state_type
1258 {
1259 const auto primitive_state = expand_state(initial_state);
1260 return from_primitive_state(primitive_state);
1261 }
1262
1263
1264 template <int dim, typename Number>
1265 DEAL_II_ALWAYS_INLINE inline auto
1267 const state_type &primitive_state) const -> state_type
1268 {
1269 const auto &rho = primitive_state[0];
1270 /* extract velocity: */
1271 const auto u = /*SIC!*/ momentum(primitive_state);
1272 const auto &p = primitive_state[dim + 1];
1273
1274 auto state = primitive_state;
1275 /* Fix up momentum: */
1276 for (unsigned int i = 1; i < dim + 1; ++i)
1277 state[i] *= rho;
1278 /* Compute total energy: */
1279 state[dim + 1] =
1280 p / (ScalarNumber(gamma() - 1.)) + Number(0.5) * rho * u * u;
1281
1282 return state;
1283 }
1284
1285
1286 template <int dim, typename Number>
1287 DEAL_II_ALWAYS_INLINE inline auto
1289 const state_type &state) const -> state_type
1290 {
1291 const auto &rho = state[0];
1292 const auto rho_inverse = Number(1.) / rho;
1293 const auto p = pressure(state);
1294
1295 auto primitive_state = state;
1296 /* Fix up velocity: */
1297 for (unsigned int i = 1; i < dim + 1; ++i)
1298 primitive_state[i] *= rho_inverse;
1299 /* Set pressure: */
1300 primitive_state[dim + 1] = p;
1301
1302 return primitive_state;
1303 }
1304
1305
1306 template <int dim, typename Number>
1307 template <typename Lambda>
1309 const state_type &state, const Lambda &lambda) const -> state_type
1310 {
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];
1315 return result;
1316 }
1317
1318 } // namespace Euler
1319} // namespace ryujin
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
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 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
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
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)
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
state_type mathematical_entropy_derivative(const state_type &U) const
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
@ dirichlet_momentum
#define RYUJIN_OMP_FOR
Definition: openmp.h:70
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)