ryujin 2.1.1 revision ef0fcd4010d109b860652ece4a7b8963fb7d46b1
hyperbolic_system.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0
3// [LANL Copyright Statement]
4// Copyright (C) 2020 - 2024 by the ryujin authors
5// Copyright (C) 2023 - 2024 by Triad National Security, LLC
6//
7
8#pragma once
9
10#include <compile_time_options.h>
11
12#include <convenience_macros.h>
13#include <deal.II/base/config.h>
14#include <discretization.h>
16#include <openmp.h>
17#include <patterns_conversion.h>
18#include <simd.h>
19#include <state_vector.h>
20
21#include <deal.II/base/parameter_acceptor.h>
22#include <deal.II/base/tensor.h>
23
24#include <array>
25
26namespace ryujin
27{
28 namespace ShallowWater
29 {
30 template <int dim, typename Number>
31 class HyperbolicSystemView;
32
42 class HyperbolicSystem final : public dealii::ParameterAcceptor
43 {
44 public:
48 static inline const std::string problem_name = "Shallow water equations";
49
53 HyperbolicSystem(const std::string &subsection = "/HyperbolicSystem");
54
61 template <int dim, typename Number>
62 auto view() const
63 {
65 }
66
67 private:
72 double gravity_;
73 double manning_friction_coefficient_;
74
75 double reference_water_depth_;
76 double dry_state_relaxation_small_;
77 double dry_state_relaxation_large_;
78
79 template <int dim, typename Number>
82 }; /* HyperbolicSystem */
83
84
101 template <int dim, typename Number>
103 {
104 public:
109 HyperbolicSystemView(const HyperbolicSystem &hyperbolic_system)
110 : hyperbolic_system_(hyperbolic_system)
111 {
112 }
113
117 template <int dim2, typename Number2>
118 auto view() const
119 {
120 return HyperbolicSystemView<dim2, Number2>{hyperbolic_system_};
121 }
122
127
132
133 DEAL_II_ALWAYS_INLINE inline ScalarNumber gravity() const
134 {
135 return hyperbolic_system_.gravity_;
136 }
137
138 DEAL_II_ALWAYS_INLINE inline ScalarNumber
140 {
141 return hyperbolic_system_.manning_friction_coefficient_;
142 }
143
144 DEAL_II_ALWAYS_INLINE inline ScalarNumber reference_water_depth() const
145 {
146 return hyperbolic_system_.reference_water_depth_;
147 }
148
149 DEAL_II_ALWAYS_INLINE inline ScalarNumber
151 {
152 return hyperbolic_system_.dry_state_relaxation_small_;
153 }
154
155 DEAL_II_ALWAYS_INLINE inline ScalarNumber
157 {
158 return hyperbolic_system_.dry_state_relaxation_large_;
159 }
160
162
166
167 private:
168 const HyperbolicSystem &hyperbolic_system_;
169
170 public:
172
176
180 static constexpr unsigned int problem_dimension = 1 + dim;
181
185 using state_type = dealii::Tensor<1, problem_dimension, Number>;
186
190 using flux_type =
191 dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim, Number>>;
192
196 using flux_contribution_type = std::tuple<state_type, Number>;
197
202 static inline const auto component_names =
203 []() -> std::array<std::string, problem_dimension> {
204 if constexpr (dim == 1)
205 return {"h", "m"};
206 else if constexpr (dim == 2)
207 return {"h", "m_1", "m_2"};
208 else if constexpr (dim == 3)
209 return {"h", "m_1", "m_2", "m_3"};
210 __builtin_trap();
211 }();
212
217 static inline const auto primitive_component_names =
218 []() -> std::array<std::string, problem_dimension> {
219 if constexpr (dim == 1)
220 return {"h", "v"};
221 else if constexpr (dim == 2)
222 return {"h", "v_1", "v_2"};
223 else if constexpr (dim == 3)
224 return {"h", "v_1", "v_2", "v_3"};
225 __builtin_trap();
226 }();
227
231 static constexpr unsigned int n_precomputed_values = 2;
232
236 using precomputed_type = std::array<Number, n_precomputed_values>;
237
241 static inline const auto precomputed_names =
242 std::array<std::string, n_precomputed_values>{"eta_m", "h_star"};
243
247 static constexpr unsigned int n_initial_precomputed_values = 1;
248
253 std::array<Number, n_initial_precomputed_values>;
254
258 static inline const auto initial_precomputed_names =
259 std::array<std::string, n_initial_precomputed_values>{"bathymetry"};
260
264 using StateVector = Vectors::
265 StateVector<ScalarNumber, problem_dimension, n_precomputed_values>;
266
272
278
286
288
292
296 static constexpr unsigned int n_precomputation_cycles = 1;
297
302 template <typename DISPATCH, typename SPARSITY>
303 void precomputation_loop(unsigned int cycle,
304 const DISPATCH &dispatch_check,
305 const SPARSITY &sparsity_simd,
306 StateVector &state_vector,
307 unsigned int left,
308 unsigned int right,
309 const bool skip_constrained_dofs = true) const;
310
312
316
321 static Number water_depth(const state_type &U);
322
329 Number inverse_water_depth_mollified(const state_type &U) const;
330
337 Number water_depth_sharp(const state_type &U) const;
338
345 Number inverse_water_depth_sharp(const state_type &U) const;
346
353 Number filter_dry_water_depth(const Number &h) const;
354
359 static dealii::Tensor<1, dim, Number> momentum(const state_type &U);
360
368 Number kinetic_energy(const state_type &U) const;
369
377 Number pressure(const state_type &U) const;
378
386 Number speed_of_sound(const state_type &U) const;
387
392 Number mathematical_entropy(const state_type &U) const;
393
399
405 bool is_admissible(const state_type &U) const;
406
408
412
418 template <int component>
420 const state_type &U,
421 const state_type &U_bar,
422 const dealii::Tensor<1, dim, Number> &normal) const;
423
427 template <typename Lambda>
429 apply_boundary_conditions(const dealii::types::boundary_id id,
430 const state_type &U,
431 const dealii::Tensor<1, dim, Number> &normal,
432 const Lambda &get_dirichlet_data) const;
433
435
439
449 flux_type f(const state_type &U) const;
450
460 flux_type g(const state_type &U) const;
461
468 const Number &Z_left,
469 const Number &Z_right) const;
470
476 std::array<state_type, 2>
478 const flux_contribution_type &) const;
479
502 const InitialPrecomputedVector &piv,
503 const unsigned int i,
504 const state_type &U_i) const;
505
508 const InitialPrecomputedVector &piv,
509 const unsigned int *js,
510 const state_type &U_j) const;
511
519 const flux_contribution_type &flux_j,
520 const dealii::Tensor<1, dim, Number> &c_ij) const;
521
525 static constexpr bool have_high_order_flux = true;
526
533 const flux_contribution_type &flux_i,
534 const flux_contribution_type &flux_j,
535 const dealii::Tensor<1, dim, Number> &c_ij) const;
536
543 const flux_contribution_type &flux_j,
544 const dealii::Tensor<1, dim, Number> &c_ij,
545 const Number &d_ij) const;
546
548
552
554 static constexpr bool have_source_terms = true;
555
562 const Number &h_star,
563 const ScalarNumber tau) const;
564
566 const unsigned int i,
567 const state_type &U_i,
568 const ScalarNumber tau) const;
569
571 const unsigned int *js,
572 const state_type &U_j,
573 const ScalarNumber tau) const;
574
576
581
592 template <typename ST>
593 state_type expand_state(const ST &state) const;
594
606 template <typename ST>
607 state_type from_initial_state(const ST &initial_state) const;
608
613 state_type from_primitive_state(const state_type &primitive_state) const;
614
618 state_type to_primitive_state(const state_type &state) const;
619
625 template <typename Lambda>
627 const Lambda &lambda) const;
628
629 }; /* HyperbolicSystemView */
630
631
632 /*
633 * -------------------------------------------------------------------------
634 * Inline definitions
635 * -------------------------------------------------------------------------
636 */
637
638
639 inline HyperbolicSystem::HyperbolicSystem(const std::string &subsection)
640 : ParameterAcceptor(subsection)
641 {
642 gravity_ = 9.81;
643 add_parameter("gravity", gravity_, "Gravitational constant [m/s^2]");
644
645 manning_friction_coefficient_ = 0.;
646 add_parameter("manning friction coefficient",
647 manning_friction_coefficient_,
648 "Roughness coefficient for friction source");
649
650 reference_water_depth_ = 1.;
651 add_parameter("reference water depth",
652 reference_water_depth_,
653 "Problem specific water depth reference");
654
655 dry_state_relaxation_small_ = 1.e2;
656 add_parameter("dry state relaxation small",
657 dry_state_relaxation_small_,
658 "Problem specific dry-state relaxation parameter");
659
660 dry_state_relaxation_large_ = 1.e4;
661 add_parameter("dry state relaxation large",
662 dry_state_relaxation_large_,
663 "Problem specific dry-state relaxation parameter");
664 }
665
666
667 template <int dim, typename Number>
668 template <typename DISPATCH, typename SPARSITY>
669 DEAL_II_ALWAYS_INLINE inline void
671 unsigned int cycle [[maybe_unused]],
672 const DISPATCH &dispatch_check,
673 const SPARSITY &sparsity_simd,
674 StateVector &state_vector,
675 unsigned int left,
676 unsigned int right,
677 const bool skip_constrained_dofs /*= true*/) const
678 {
679 Assert(cycle == 0, dealii::ExcInternalError());
680
681 const auto &U = std::get<0>(state_vector);
682 auto &precomputed = std::get<1>(state_vector);
683
684 /* We are inside a thread parallel context */
685
686 unsigned int stride_size = get_stride_size<Number>;
687
689 for (unsigned int i = left; i < right; i += stride_size) {
690
691 /* Skip constrained degrees of freedom: */
692 const unsigned int row_length = sparsity_simd.row_length(i);
693 if (skip_constrained_dofs && row_length == 1)
694 continue;
695
696 dispatch_check(i);
697
698 const auto U_i = U.template get_tensor<Number>(i);
699
700 const auto eta_m = mathematical_entropy(U_i);
701
702 const auto h_sharp = water_depth_sharp(U_i);
703 const auto h_star = ryujin::pow(h_sharp, ScalarNumber(4. / 3.));
704
705 const precomputed_type prec_i{eta_m, h_star};
706 precomputed.template write_tensor<Number>(prec_i, i);
707 }
708 }
709
710
711 template <int dim, typename Number>
712 DEAL_II_ALWAYS_INLINE inline Number
714 {
715 return U[0];
716 }
717
718
719 template <int dim, typename Number>
720 DEAL_II_ALWAYS_INLINE inline Number
722 const state_type &U) const
723 {
724 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
725
726 const Number h_cutoff_mollified =
727 reference_water_depth() * dry_state_relaxation_large() * Number(eps);
728
729 const Number h = water_depth(U);
730 const Number h_pos = positive_part(water_depth(U));
731 const Number h_max = std::max(h, h_cutoff_mollified);
732 const Number denom = h * h + h_max * h_max;
733 return ScalarNumber(2.) * h_pos / denom;
734 }
735
736
737 template <int dim, typename Number>
738 DEAL_II_ALWAYS_INLINE inline Number
740 const state_type &U) const
741 {
742 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
743
744 const Number h_cutoff_small =
745 reference_water_depth() * dry_state_relaxation_small() * Number(eps);
746
747 const Number h = water_depth(U);
748 const Number h_max = std::max(h, h_cutoff_small);
749 return h_max;
750 }
751
752
753 template <int dim, typename Number>
754 DEAL_II_ALWAYS_INLINE inline Number
756 const state_type &U) const
757 {
758 return ScalarNumber(1.) / water_depth_sharp(U);
759 }
760
761
762 template <int dim, typename Number>
763 DEAL_II_ALWAYS_INLINE inline Number
765 const Number &h) const
766 {
768 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
769
770 const Number h_cutoff_large =
771 reference_water_depth() * dry_state_relaxation_large() * Number(eps);
772
773 return dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
774 std::abs(h), h_cutoff_large, Number(0.), h);
775 }
776
777
778 template <int dim, typename Number>
779 DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim, Number>
781 {
782 dealii::Tensor<1, dim, Number> result;
783
784 for (unsigned int i = 0; i < dim; ++i)
785 result[i] = U[1 + i];
786 return result;
787 }
788
789
790 template <int dim, typename Number>
791 DEAL_II_ALWAYS_INLINE inline Number
793 {
794 const auto h = water_depth(U);
795 const auto vel = momentum(U) * inverse_water_depth_sharp(U);
796
797 /* KE = 1/2 h |v|^2 */
798 return ScalarNumber(0.5) * h * vel.norm_square();
799 }
800
801
802 template <int dim, typename Number>
803 DEAL_II_ALWAYS_INLINE inline Number
805 {
806 const Number h_sqd = U[0] * U[0];
807
808 /* p = 1/2 g h^2 */
809 return ScalarNumber(0.5) * gravity() * h_sqd;
810 }
811
812
813 template <int dim, typename Number>
814 DEAL_II_ALWAYS_INLINE inline Number
816 {
817 /* c^2 = g * h */
818 return std::sqrt(gravity() * U[0]);
819 }
820
821
822 template <int dim, typename Number>
823 DEAL_II_ALWAYS_INLINE inline Number
825 const state_type &U) const
826 {
827 const auto p = pressure(U);
828 const auto k_e = kinetic_energy(U);
829 return p + k_e;
830 }
831
832
833 template <int dim, typename Number>
834 DEAL_II_ALWAYS_INLINE inline auto
836 const state_type &U) const -> state_type
837 {
838 /*
839 * With
840 * eta = 1/2 g h^2 + 1/2 |m|^2 / h
841 *
842 * we get
843 *
844 * eta' = (g h - 1/2 |vel|^2, vel)
845 *
846 * where vel = m / h
847 */
848
849 state_type result;
850
851 const Number &h = U[0];
852 const auto vel = momentum(U) * inverse_water_depth_sharp(U);
853
854 // water depth component
855 result[0] = gravity() * h - ScalarNumber(0.5) * vel.norm_square();
856
857 // momentum components
858 for (unsigned int i = 0; i < dim; ++i) {
859 result[1 + i] = vel[i];
860 }
861
862 return result;
863 }
864
865
866 template <int dim, typename Number>
867 DEAL_II_ALWAYS_INLINE inline bool
869 {
870 const auto h = filter_dry_water_depth(water_depth(U));
871
872 constexpr auto gte = dealii::SIMDComparison::greater_than_or_equal;
873 const auto test = dealii::compare_and_apply_mask<gte>(
874 h, Number(0.), Number(0.), Number(-1.));
875
876#ifdef DEBUG_OUTPUT
877 if (!(test == Number(0.))) {
878 std::cout << std::fixed << std::setprecision(16);
879 std::cout << "Bounds violation: Negative state [h] detected!\n";
880 std::cout << "\t\th: " << h << "\n" << std::endl;
881 __builtin_trap();
882 }
883#endif
884
885 return (test == Number(0.));
886 }
887
888
889 template <int dim, typename Number>
890 template <int component>
891 DEAL_II_ALWAYS_INLINE inline auto
893 const state_type &U,
894 const state_type &U_bar,
895 const dealii::Tensor<1, dim, Number> &normal) const -> state_type
896 {
897 /* Note that U_bar are the dirichlet values that are prescribed */
898 static_assert(component == 1 || component == 2,
899 "component has to be 1 or 2");
900
902
903 const auto m = momentum(U);
904 const auto a = speed_of_sound(U);
905 const auto vn = m * normal * inverse_water_depth_sharp(U);
906
907 const auto m_bar = momentum(U_bar);
908 const auto a_bar = speed_of_sound(U_bar);
909 const auto vn_bar = m_bar * normal * inverse_water_depth_sharp(U_bar);
910
911 /* First Riemann characteristic: v * n - 2 * a */
912
913 const auto R_1 = component == 1 ? vn_bar - ScalarNumber(2.) * a_bar
914 : vn - ScalarNumber(2.) * a;
915
916 /* Second Riemann characteristic: v * n + 2 * a */
917
918 const auto R_2 = component == 2 ? vn_bar + ScalarNumber(2.) * a_bar
919 : vn + ScalarNumber(2.) * a;
920
921 const auto vperp = m * inverse_water_depth_sharp(U) - vn * normal;
922
923 const auto vn_new = ScalarNumber(0.5) * (R_1 + R_2);
924
925 const auto h_new =
926 ryujin::fixed_power<2>((R_2 - R_1) / ScalarNumber(4.)) / gravity();
927
928 state_type U_new;
929 U_new[0] = h_new;
930 for (unsigned int d = 0; d < dim; ++d) {
931 U_new[1 + d] = h_new * (vn_new * normal + vperp)[d];
932 }
933
934 return U_new;
935 }
936
937
938 template <int dim, typename Number>
939 template <typename Lambda>
940 DEAL_II_ALWAYS_INLINE inline auto
942 const dealii::types::boundary_id id,
943 const state_type &U,
944 const dealii::Tensor<1, dim, Number> &normal,
945 const Lambda &get_dirichlet_data) const -> state_type
946 {
947 state_type result = U;
948
949 if (id == Boundary::dirichlet) {
950 result = get_dirichlet_data();
951
952 } else if (id == Boundary::dirichlet_momentum) {
953 /* Only enforce Dirichlet conditions on the momentum: */
954 auto m_dirichlet = momentum(get_dirichlet_data());
955 for (unsigned int k = 0; k < dim; ++k)
956 result[k + 1] = m_dirichlet[k];
957
958 } else if (id == Boundary::slip) {
959 auto m = momentum(U);
960 m -= 1. * (m * normal) * normal;
961 for (unsigned int k = 0; k < dim; ++k)
962 result[k + 1] = m[k];
963
964 } else if (id == Boundary::no_slip) {
965 for (unsigned int k = 0; k < dim; ++k)
966 result[k + 1] = Number(0.);
967
968 } else if (id == Boundary::dynamic) {
969 /*
970 * On dynamic boundary conditions, we distinguish four cases:
971 *
972 * - supersonic inflow: prescribe full state
973 * - subsonic inflow:
974 * decompose into Riemann invariants and leave R_2
975 * characteristic untouched.
976 * - supersonic outflow: do nothing
977 * - subsonic outflow:
978 * decompose into Riemann invariants and prescribe incoming
979 * R_1 characteristic.
980 */
981 const auto m = momentum(U);
982 const auto h_inverse = inverse_water_depth_sharp(U);
983 const auto a = speed_of_sound(U);
984 const auto vn = m * normal * h_inverse;
985
986 /* Supersonic inflow: */
987 if (vn < -a) {
988 result = get_dirichlet_data();
989 }
990
991 /* Subsonic inflow: */
992 if (vn >= -a && vn <= 0.) {
993 const auto U_dirichlet = get_dirichlet_data();
994 result = prescribe_riemann_characteristic<2>(U_dirichlet, U, normal);
995 }
996
997 /* Subsonic outflow: */
998 if (vn > 0. && vn <= a) {
999 const auto U_dirichlet = get_dirichlet_data();
1000 result = prescribe_riemann_characteristic<1>(U, U_dirichlet, normal);
1001 }
1002
1003 /* Supersonic outflow: do nothing, i.e., keep U as is */
1004
1005 } else {
1006 AssertThrow(false, dealii::ExcNotImplemented());
1007 }
1008
1009 return result;
1010 }
1011
1012
1013 template <int dim, typename Number>
1014 DEAL_II_ALWAYS_INLINE inline auto
1016 {
1017 const auto h_inverse = inverse_water_depth_sharp(U);
1018 const auto m = momentum(U);
1019 const auto p = pressure(U);
1020
1021 flux_type result;
1022
1023 result[0] = (m * h_inverse) * U[0];
1024 for (unsigned int i = 0; i < dim; ++i) {
1025 result[1 + i] = (m * h_inverse) * m[i];
1026 result[1 + i][i] += p;
1027 }
1028 return result;
1029 }
1030
1031
1032 template <int dim, typename Number>
1033 DEAL_II_ALWAYS_INLINE inline auto
1035 {
1036 const auto h_inverse = inverse_water_depth_sharp(U);
1037 const auto m = momentum(U);
1038
1039 flux_type result;
1040
1041 result[0] = (m * h_inverse) * U[0];
1042 for (unsigned int i = 0; i < dim; ++i) {
1043 result[1 + i] = (m * h_inverse) * m[i];
1044 }
1045 return result;
1046 }
1047
1048
1049 template <int dim, typename Number>
1050 DEAL_II_ALWAYS_INLINE inline auto
1052 const Number &Z_left,
1053 const Number &Z_right) const
1054 -> state_type
1055 {
1056 const Number Z_max = std::max(Z_left, Z_right);
1057 const Number h = water_depth(U);
1058 const Number H_star = std::max(Number(0.), h + Z_left - Z_max);
1059
1060 return U * H_star * inverse_water_depth_mollified(U);
1061 }
1062
1063
1064 template <int dim, typename Number>
1065 DEAL_II_ALWAYS_INLINE inline auto
1067 const flux_contribution_type &flux_i,
1068 const flux_contribution_type &flux_j) const -> std::array<state_type, 2>
1069 {
1070 const auto &[U_i, Z_i] = flux_i;
1071 const auto &[U_j, Z_j] = flux_j;
1072
1073 const auto U_star_ij = star_state(U_i, Z_i, Z_j);
1074 const auto U_star_ji = star_state(U_j, Z_j, Z_i);
1075
1076 return {U_star_ij, U_star_ji};
1077 }
1078
1079
1080 template <int dim, typename Number>
1081 DEAL_II_ALWAYS_INLINE inline auto
1083 const PrecomputedVector & /*pv*/,
1084 const InitialPrecomputedVector &piv,
1085 const unsigned int i,
1086 const state_type &U_i) const -> flux_contribution_type
1087 {
1088 const auto Z_i = piv.template get_tensor<Number>(i)[0];
1089 return {U_i, Z_i};
1090 }
1091
1092
1093 template <int dim, typename Number>
1094 DEAL_II_ALWAYS_INLINE inline auto
1096 const PrecomputedVector & /*pv*/,
1097 const InitialPrecomputedVector &piv,
1098 const unsigned int *js,
1099 const state_type &U_j) const -> flux_contribution_type
1100 {
1101 const auto Z_j = piv.template get_tensor<Number>(js)[0];
1102 return {U_j, Z_j};
1103 }
1104
1105
1106 template <int dim, typename Number>
1107 DEAL_II_ALWAYS_INLINE inline auto
1109 const flux_contribution_type &flux_i,
1110 const flux_contribution_type &flux_j,
1111 const dealii::Tensor<1, dim, Number> &c_ij) const -> state_type
1112 {
1113 const auto &[U_i, Z_i] = flux_i;
1114 const auto &[U_star_ij, U_star_ji] = equilibrated_states(flux_i, flux_j);
1115
1116 const auto H_i = water_depth(U_i);
1117 const auto H_star_ij = water_depth(U_star_ij);
1118 const auto H_star_ji = water_depth(U_star_ji);
1119
1120 const auto g_i = g(U_star_ij);
1121 const auto g_j = g(U_star_ji);
1122
1123 auto result = -add(g_i, g_j);
1124
1125 const auto factor =
1126 (ScalarNumber(0.5) * (H_star_ji * H_star_ji - H_star_ij * H_star_ij) +
1127 H_i * H_i) *
1128 gravity();
1129
1130 for (unsigned int i = 0; i < dim; ++i) {
1131 result[1 + i][i] -= factor;
1132 }
1133
1134 return contract(result, c_ij);
1135 }
1136
1137
1138 template <int dim, typename Number>
1139 DEAL_II_ALWAYS_INLINE inline auto
1141 const flux_contribution_type &flux_i,
1142 const flux_contribution_type &flux_j,
1143 const dealii::Tensor<1, dim, Number> &c_ij) const -> state_type
1144 {
1145 const auto &[U_i, Z_i] = flux_i;
1146 const auto &[U_j, Z_j] = flux_j;
1147
1148 const auto H_i = water_depth(U_i);
1149 const auto H_j = water_depth(U_j);
1150
1151 const auto g_i = g(U_i);
1152 const auto g_j = g(U_j);
1153
1154 auto result = -add(g_i, g_j);
1155
1156 const auto factor = gravity() * H_i * (H_j + Z_j - Z_i);
1157 for (unsigned int i = 0; i < dim; ++i) {
1158 result[1 + i][i] -= factor;
1159 }
1160
1161 return contract(result, c_ij);
1162 }
1163
1164
1165 template <int dim, typename Number>
1166 DEAL_II_ALWAYS_INLINE inline auto
1168 const flux_contribution_type &flux_i,
1169 const flux_contribution_type &flux_j,
1170 const dealii::Tensor<1, dim, Number> &c_ij,
1171 const Number &d_ij) const -> state_type
1172 {
1173 const auto &[U_i, Z_i] = flux_i;
1174 const auto &[U_j, Z_j] = flux_j;
1175 const auto U_star_ij = star_state(U_i, Z_i, Z_j);
1176
1177 const auto h_inverse = inverse_water_depth_sharp(U_i);
1178 const auto m = momentum(U_i);
1179 const auto factor = ScalarNumber(2.) * (d_ij + h_inverse * (m * c_ij));
1180
1181 return -factor * (U_star_ij - U_i);
1182 }
1183
1184
1185 template <int dim, typename Number>
1186 DEAL_II_ALWAYS_INLINE inline auto
1188 const state_type &U, const Number &h_star, const ScalarNumber tau) const
1189 -> state_type
1190 {
1191 state_type result;
1192
1193 const auto g = gravity();
1194 const auto n = manning_friction_coefficient();
1195
1196 const auto h_inverse = inverse_water_depth_mollified(U);
1197
1198 const auto m = momentum(U);
1199 const auto v_norm = (m * h_inverse).norm();
1200 const auto factor = ScalarNumber(2.) * g * n * n * v_norm;
1201
1202 const auto denominator = h_star + std::max(h_star, tau * factor);
1203 const auto denominator_inverse = ScalarNumber(1.) / denominator;
1204
1205 for (unsigned int d = 0; d < dim; ++d)
1206 result[d + 1] = -factor * denominator_inverse * m[d];
1207
1208 return result;
1209 }
1210
1211
1212 template <int dim, typename Number>
1213 DEAL_II_ALWAYS_INLINE inline auto
1215 const PrecomputedVector &pv,
1216 const unsigned int i,
1217 const state_type &U_i,
1218 const ScalarNumber tau) const -> state_type
1219 {
1220 const auto &[eta_m, h_star] =
1221 pv.template get_tensor<Number, precomputed_type>(i);
1222
1223 return manning_friction(U_i, h_star, tau);
1224 }
1225
1226
1227 template <int dim, typename Number>
1228 DEAL_II_ALWAYS_INLINE inline auto
1230 const PrecomputedVector &pv,
1231 const unsigned int *js,
1232 const state_type &U_j,
1233 const ScalarNumber tau) const -> state_type
1234 {
1235 const auto &[eta_m, h_star] =
1236 pv.template get_tensor<Number, precomputed_type>(js);
1237
1238 return manning_friction(U_j, h_star, tau);
1239 }
1240
1241
1242 template <int dim, typename Number>
1243 template <typename ST>
1244 DEAL_II_ALWAYS_INLINE inline auto
1246 -> state_type
1247 {
1248 using T = typename ST::value_type;
1249 static_assert(std::is_same_v<Number, T>, "template mismatch");
1250
1251 constexpr auto dim2 = ST::dimension - 1;
1252 static_assert(dim >= dim2,
1253 "the space dimension of the argument state must not be "
1254 "larger than the one of the target state");
1255
1256 state_type result;
1257 result[0] = state[0];
1258 for (unsigned int i = 1; i < dim2 + 1; ++i)
1259 result[i] = state[i];
1260
1261 return result;
1262 }
1263
1264 template <int dim, typename Number>
1265 template <typename ST>
1266 DEAL_II_ALWAYS_INLINE inline auto
1268 const ST &initial_state) const -> state_type
1269 {
1270 const auto primitive_state = expand_state(initial_state);
1271 return from_primitive_state(primitive_state);
1272 }
1273
1274
1275 template <int dim, typename Number>
1276 DEAL_II_ALWAYS_INLINE inline auto
1278 const state_type &primitive_state) const -> state_type
1279 {
1280 const auto &h = primitive_state[0];
1281
1282 auto state = primitive_state;
1283 /* Fix up momentum: */
1284 for (unsigned int i = 1; i < dim + 1; ++i)
1285 state[i] *= h;
1286
1287 return state;
1288 }
1289
1290
1291 template <int dim, typename Number>
1292 DEAL_II_ALWAYS_INLINE inline auto
1294 const state_type &state) const -> state_type
1295 {
1296 const auto h_inverse = inverse_water_depth_sharp(state);
1297
1298 auto primitive_state = state;
1299 /* Fix up velocity: */
1300 for (unsigned int i = 1; i < dim + 1; ++i)
1301 primitive_state[i] *= h_inverse;
1302
1303 return primitive_state;
1304 }
1305
1306
1307 template <int dim, typename Number>
1308 template <typename Lambda>
1309 DEAL_II_ALWAYS_INLINE inline auto
1311 const state_type &state, const Lambda &lambda) const -> state_type
1312 {
1313 auto result = state;
1314 auto M = lambda(momentum(state));
1315 for (unsigned int d = 0; d < dim; ++d)
1316 result[1 + d] = M[d];
1317 return result;
1318 }
1319
1320 } // namespace ShallowWater
1321} // namespace ryujin
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
dealii::Tensor< 1, problem_dimension, dealii::Tensor< 1, dim, Number > > flux_type
Number pressure(const state_type &U) 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 ScalarNumber gravity() const
bool is_admissible(const state_type &U) const
state_type mathematical_entropy_derivative(const state_type &U) const
typename get_value_type< Number >::type ScalarNumber
static Number water_depth(const state_type &U)
DEAL_II_ALWAYS_INLINE ScalarNumber dry_state_relaxation_large() const
static constexpr unsigned int n_initial_precomputed_values
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
DEAL_II_ALWAYS_INLINE ScalarNumber dry_state_relaxation_small() const
Number inverse_water_depth_mollified(const state_type &U) const
Number water_depth_sharp(const state_type &U) const
Vectors::StateVector< ScalarNumber, problem_dimension, n_precomputed_values > StateVector
state_type nodal_source(const PrecomputedVector &pv, const unsigned int i, const state_type &U_i, const ScalarNumber tau) const
static constexpr unsigned int n_precomputed_values
Number inverse_water_depth_sharp(const state_type &U) const
Number kinetic_energy(const state_type &U) const
Number filter_dry_water_depth(const Number &h) const
state_type from_primitive_state(const state_type &primitive_state) const
state_type manning_friction(const state_type &U, const Number &h_star, const ScalarNumber tau) const
state_type expand_state(const ST &state) const
std::array< Number, n_initial_precomputed_values > initial_precomputed_type
DEAL_II_ALWAYS_INLINE ScalarNumber manning_friction_coefficient() const
Number speed_of_sound(const state_type &U) const
state_type star_state(const state_type &U, const Number &Z_left, const Number &Z_right) const
dealii::Tensor< 1, problem_dimension, Number > state_type
static dealii::Tensor< 1, dim, Number > momentum(const state_type &U)
state_type affine_shift(const flux_contribution_type &flux_i, const flux_contribution_type &flux_j, const dealii::Tensor< 1, dim, Number > &c_ij, const Number &d_ij) const
std::array< Number, n_precomputed_values > precomputed_type
state_type to_primitive_state(const state_type &state) const
state_type from_initial_state(const ST &initial_state) const
static constexpr unsigned int problem_dimension
flux_contribution_type flux_contribution(const PrecomputedVector &pv, const InitialPrecomputedVector &piv, const unsigned int i, const state_type &U_i) const
static constexpr unsigned int n_precomputation_cycles
flux_type g(const state_type &U) const
Number mathematical_entropy(const state_type &U) const
HyperbolicSystemView(const HyperbolicSystem &hyperbolic_system)
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
std::array< state_type, 2 > equilibrated_states(const flux_contribution_type &, const flux_contribution_type &) 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
std::tuple< state_type, Number > flux_contribution_type
flux_type f(const state_type &U) const
state_type prescribe_riemann_characteristic(const state_type &U, const state_type &U_bar, const dealii::Tensor< 1, dim, Number > &normal) const
state_type apply_galilei_transform(const state_type &state, const Lambda &lambda) const
DEAL_II_ALWAYS_INLINE ScalarNumber reference_water_depth() const
HyperbolicSystem(const std::string &subsection="/HyperbolicSystem")
@ dirichlet_momentum
#define RYUJIN_OMP_FOR
Definition: openmp.h:70
T pow(const T x, const T b)
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)
Definition: simd.h:112
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)