ryujin 2.1.1 revision 350e54cc11f3d21282bddcf3e6517944dbc508bf
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 <convenience_macros.h>
11#include <deal.II/base/config.h>
12#include <discretization.h>
14#include <openmp.h>
15#include <patterns_conversion.h>
16#include <simd.h>
17#include <state_vector.h>
18
19#include <deal.II/base/parameter_acceptor.h>
20#include <deal.II/base/tensor.h>
21
22#include <array>
23
24namespace ryujin
25{
26 namespace ShallowWater
27 {
28 template <int dim, typename Number>
29 class HyperbolicSystemView;
30
40 class HyperbolicSystem final : public dealii::ParameterAcceptor
41 {
42 public:
46 static inline const std::string problem_name = "Shallow water equations";
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 gravity_;
71 double manning_friction_coefficient_;
72
73 double reference_water_depth_;
74 double dry_state_relaxation_small_;
75 double dry_state_relaxation_large_;
76
77 template <int dim, typename Number>
80 }; /* HyperbolicSystem */
81
82
99 template <int dim, typename Number>
101 {
102 public:
107 HyperbolicSystemView(const HyperbolicSystem &hyperbolic_system)
108 : hyperbolic_system_(hyperbolic_system)
109 {
110 }
111
115 template <int dim2, typename Number2>
116 auto view() const
117 {
118 return HyperbolicSystemView<dim2, Number2>{hyperbolic_system_};
119 }
120
125
130
131 DEAL_II_ALWAYS_INLINE inline ScalarNumber gravity() const
132 {
133 return hyperbolic_system_.gravity_;
134 }
135
136 DEAL_II_ALWAYS_INLINE inline ScalarNumber
138 {
139 return hyperbolic_system_.manning_friction_coefficient_;
140 }
141
142 DEAL_II_ALWAYS_INLINE inline ScalarNumber reference_water_depth() const
143 {
144 return hyperbolic_system_.reference_water_depth_;
145 }
146
147 DEAL_II_ALWAYS_INLINE inline ScalarNumber
149 {
150 return hyperbolic_system_.dry_state_relaxation_small_;
151 }
152
153 DEAL_II_ALWAYS_INLINE inline ScalarNumber
155 {
156 return hyperbolic_system_.dry_state_relaxation_large_;
157 }
158
160
164
165 private:
166 const HyperbolicSystem &hyperbolic_system_;
167
168 public:
170
174
178 static constexpr unsigned int problem_dimension = 1 + dim;
179
183 using state_type = dealii::Tensor<1, problem_dimension, Number>;
184
188 using flux_type =
189 dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim, Number>>;
190
194 using flux_contribution_type = std::tuple<state_type, Number>;
195
200 static inline const auto component_names =
201 []() -> std::array<std::string, problem_dimension> {
202 if constexpr (dim == 1)
203 return {"h", "m"};
204 else if constexpr (dim == 2)
205 return {"h", "m_1", "m_2"};
206 else if constexpr (dim == 3)
207 return {"h", "m_1", "m_2", "m_3"};
208 __builtin_trap();
209 }();
210
215 static inline const auto primitive_component_names =
216 []() -> std::array<std::string, problem_dimension> {
217 if constexpr (dim == 1)
218 return {"h", "v"};
219 else if constexpr (dim == 2)
220 return {"h", "v_1", "v_2"};
221 else if constexpr (dim == 3)
222 return {"h", "v_1", "v_2", "v_3"};
223 __builtin_trap();
224 }();
225
229 static constexpr unsigned int n_precomputed_values = 2;
230
234 using precomputed_type = std::array<Number, n_precomputed_values>;
235
239 static inline const auto precomputed_names =
240 std::array<std::string, n_precomputed_values>{"eta_m", "h_star"};
241
245 static constexpr unsigned int n_initial_precomputed_values = 1;
246
251 std::array<Number, n_initial_precomputed_values>;
252
256 static inline const auto initial_precomputed_names =
257 std::array<std::string, n_initial_precomputed_values>{"bathymetry"};
258
262 using StateVector = Vectors::
263 StateVector<ScalarNumber, problem_dimension, n_precomputed_values>;
264
270
276
284
286
290
294 static constexpr unsigned int n_precomputation_cycles = 1;
295
300 template <typename DISPATCH, typename SPARSITY>
301 void precomputation_loop(unsigned int cycle,
302 const DISPATCH &dispatch_check,
303 const SPARSITY &sparsity_simd,
304 StateVector &state_vector,
305 unsigned int left,
306 unsigned int right,
307 const bool skip_constrained_dofs = true) const;
308
310
314
319 static Number water_depth(const state_type &U);
320
327 Number inverse_water_depth_mollified(const state_type &U) const;
328
335 Number water_depth_sharp(const state_type &U) const;
336
343 Number inverse_water_depth_sharp(const state_type &U) const;
344
351 Number filter_dry_water_depth(const Number &h) const;
352
357 static dealii::Tensor<1, dim, Number> momentum(const state_type &U);
358
366 Number kinetic_energy(const state_type &U) const;
367
375 Number pressure(const state_type &U) const;
376
384 Number speed_of_sound(const state_type &U) const;
385
390 Number mathematical_entropy(const state_type &U) const;
391
397
403 bool is_admissible(const state_type &U) const;
404
406
410
416 template <int component>
418 const state_type &U,
419 const state_type &U_bar,
420 const dealii::Tensor<1, dim, Number> &normal) const;
421
425 template <typename Lambda>
427 apply_boundary_conditions(const dealii::types::boundary_id id,
428 const state_type &U,
429 const dealii::Tensor<1, dim, Number> &normal,
430 const Lambda &get_dirichlet_data) const;
431
433
437
447 flux_type f(const state_type &U) const;
448
458 flux_type g(const state_type &U) const;
459
466 const Number &Z_left,
467 const Number &Z_right) const;
468
474 std::array<state_type, 2>
476 const flux_contribution_type &) const;
477
500 const InitialPrecomputedVector &piv,
501 const unsigned int i,
502 const state_type &U_i) const;
503
506 const InitialPrecomputedVector &piv,
507 const unsigned int *js,
508 const state_type &U_j) const;
509
517 const flux_contribution_type &flux_j,
518 const dealii::Tensor<1, dim, Number> &c_ij) const;
519
523 static constexpr bool have_high_order_flux = true;
524
531 const flux_contribution_type &flux_i,
532 const flux_contribution_type &flux_j,
533 const dealii::Tensor<1, dim, Number> &c_ij) const;
534
541 const flux_contribution_type &flux_j,
542 const dealii::Tensor<1, dim, Number> &c_ij,
543 const Number &d_ij) const;
544
546
550
552 static constexpr bool have_source_terms = true;
553
560 const Number &h_star,
561 const ScalarNumber tau) const;
562
564 const unsigned int i,
565 const state_type &U_i,
566 const ScalarNumber tau) const;
567
569 const unsigned int *js,
570 const state_type &U_j,
571 const ScalarNumber tau) const;
572
574
579
590 template <typename ST>
591 state_type expand_state(const ST &state) const;
592
604 template <typename ST>
605 state_type from_initial_state(const ST &initial_state) const;
606
611 state_type from_primitive_state(const state_type &primitive_state) const;
612
616 state_type to_primitive_state(const state_type &state) const;
617
623 template <typename Lambda>
625 const Lambda &lambda) const;
626
627 }; /* HyperbolicSystemView */
628
629
630 /*
631 * -------------------------------------------------------------------------
632 * Inline definitions
633 * -------------------------------------------------------------------------
634 */
635
636
637 inline HyperbolicSystem::HyperbolicSystem(const std::string &subsection)
638 : ParameterAcceptor(subsection)
639 {
640 gravity_ = 9.81;
641 add_parameter("gravity", gravity_, "Gravitational constant [m/s^2]");
642
643 manning_friction_coefficient_ = 0.;
644 add_parameter("manning friction coefficient",
645 manning_friction_coefficient_,
646 "Roughness coefficient for friction source");
647
648 reference_water_depth_ = 1.;
649 add_parameter("reference water depth",
650 reference_water_depth_,
651 "Problem specific water depth reference");
652
653 dry_state_relaxation_small_ = 1.e2;
654 add_parameter("dry state relaxation small",
655 dry_state_relaxation_small_,
656 "Problem specific dry-state relaxation parameter");
657
658 dry_state_relaxation_large_ = 1.e4;
659 add_parameter("dry state relaxation large",
660 dry_state_relaxation_large_,
661 "Problem specific dry-state relaxation parameter");
662 }
663
664
665 template <int dim, typename Number>
666 template <typename DISPATCH, typename SPARSITY>
667 DEAL_II_ALWAYS_INLINE inline void
669 unsigned int cycle [[maybe_unused]],
670 const DISPATCH &dispatch_check,
671 const SPARSITY &sparsity_simd,
672 StateVector &state_vector,
673 unsigned int left,
674 unsigned int right,
675 const bool skip_constrained_dofs /*= true*/) const
676 {
677 Assert(cycle == 0, dealii::ExcInternalError());
678
679 const auto &U = std::get<0>(state_vector);
680 auto &precomputed = std::get<1>(state_vector);
681
682 /* We are inside a thread parallel context */
683
684 unsigned int stride_size = get_stride_size<Number>;
685
687 for (unsigned int i = left; i < right; i += stride_size) {
688
689 /* Skip constrained degrees of freedom: */
690 const unsigned int row_length = sparsity_simd.row_length(i);
691 if (skip_constrained_dofs && row_length == 1)
692 continue;
693
694 dispatch_check(i);
695
696 const auto U_i = U.template get_tensor<Number>(i);
697
698 const auto eta_m = mathematical_entropy(U_i);
699
700 const auto h_sharp = water_depth_sharp(U_i);
701 const auto h_star = ryujin::pow(h_sharp, ScalarNumber(4. / 3.));
702
703 const precomputed_type prec_i{eta_m, h_star};
704 precomputed.template write_tensor<Number>(prec_i, i);
705 }
706 }
707
708
709 template <int dim, typename Number>
710 DEAL_II_ALWAYS_INLINE inline Number
712 {
713 return U[0];
714 }
715
716
717 template <int dim, typename Number>
718 DEAL_II_ALWAYS_INLINE inline Number
720 const state_type &U) const
721 {
722 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
723
724 const Number h_cutoff_mollified =
725 reference_water_depth() * dry_state_relaxation_large() * Number(eps);
726
727 const Number h = water_depth(U);
728 const Number h_pos = positive_part(water_depth(U));
729 const Number h_max = std::max(h, h_cutoff_mollified);
730 const Number denom = h * h + h_max * h_max;
731 return ScalarNumber(2.) * h_pos / denom;
732 }
733
734
735 template <int dim, typename Number>
736 DEAL_II_ALWAYS_INLINE inline Number
738 const state_type &U) const
739 {
740 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
741
742 const Number h_cutoff_small =
743 reference_water_depth() * dry_state_relaxation_small() * Number(eps);
744
745 const Number h = water_depth(U);
746 const Number h_max = std::max(h, h_cutoff_small);
747 return h_max;
748 }
749
750
751 template <int dim, typename Number>
752 DEAL_II_ALWAYS_INLINE inline Number
754 const state_type &U) const
755 {
756 return ScalarNumber(1.) / water_depth_sharp(U);
757 }
758
759
760 template <int dim, typename Number>
761 DEAL_II_ALWAYS_INLINE inline Number
763 const Number &h) const
764 {
766 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
767
768 const Number h_cutoff_large =
769 reference_water_depth() * dry_state_relaxation_large() * Number(eps);
770
771 return dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
772 std::abs(h), h_cutoff_large, Number(0.), h);
773 }
774
775
776 template <int dim, typename Number>
777 DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim, Number>
779 {
780 dealii::Tensor<1, dim, Number> result;
781
782 for (unsigned int i = 0; i < dim; ++i)
783 result[i] = U[1 + i];
784 return result;
785 }
786
787
788 template <int dim, typename Number>
789 DEAL_II_ALWAYS_INLINE inline Number
791 {
792 const auto h = water_depth(U);
793 const auto vel = momentum(U) * inverse_water_depth_sharp(U);
794
795 /* KE = 1/2 h |v|^2 */
796 return ScalarNumber(0.5) * h * vel.norm_square();
797 }
798
799
800 template <int dim, typename Number>
801 DEAL_II_ALWAYS_INLINE inline Number
803 {
804 const Number h_sqd = U[0] * U[0];
805
806 /* p = 1/2 g h^2 */
807 return ScalarNumber(0.5) * gravity() * h_sqd;
808 }
809
810
811 template <int dim, typename Number>
812 DEAL_II_ALWAYS_INLINE inline Number
814 {
815 /* c^2 = g * h */
816 return std::sqrt(gravity() * U[0]);
817 }
818
819
820 template <int dim, typename Number>
821 DEAL_II_ALWAYS_INLINE inline Number
823 const state_type &U) const
824 {
825 const auto p = pressure(U);
826 const auto k_e = kinetic_energy(U);
827 return p + k_e;
828 }
829
830
831 template <int dim, typename Number>
832 DEAL_II_ALWAYS_INLINE inline auto
834 const state_type &U) const -> state_type
835 {
836 /*
837 * With
838 * eta = 1/2 g h^2 + 1/2 |m|^2 / h
839 *
840 * we get
841 *
842 * eta' = (g h - 1/2 |vel|^2, vel)
843 *
844 * where vel = m / h
845 */
846
847 state_type result;
848
849 const Number &h = U[0];
850 const auto vel = momentum(U) * inverse_water_depth_sharp(U);
851
852 // water depth component
853 result[0] = gravity() * h - ScalarNumber(0.5) * vel.norm_square();
854
855 // momentum components
856 for (unsigned int i = 0; i < dim; ++i) {
857 result[1 + i] = vel[i];
858 }
859
860 return result;
861 }
862
863
864 template <int dim, typename Number>
865 DEAL_II_ALWAYS_INLINE inline bool
867 {
868 const auto h = filter_dry_water_depth(water_depth(U));
869
870 constexpr auto gte = dealii::SIMDComparison::greater_than_or_equal;
871 const auto test = dealii::compare_and_apply_mask<gte>(
872 h, Number(0.), Number(0.), Number(-1.));
873
874#ifdef DEBUG_OUTPUT
875 if (!(test == Number(0.))) {
876 std::cout << std::fixed << std::setprecision(16);
877 std::cout << "Bounds violation: Negative state [h] detected!\n";
878 std::cout << "\t\th: " << h << "\n" << std::endl;
879 __builtin_trap();
880 }
881#endif
882
883 return (test == Number(0.));
884 }
885
886
887 template <int dim, typename Number>
888 template <int component>
889 DEAL_II_ALWAYS_INLINE inline auto
891 const state_type &U,
892 const state_type &U_bar,
893 const dealii::Tensor<1, dim, Number> &normal) const -> state_type
894 {
895 /* Note that U_bar are the dirichlet values that are prescribed */
896 static_assert(component == 1 || component == 2,
897 "component has to be 1 or 2");
898
900
901 const auto m = momentum(U);
902 const auto a = speed_of_sound(U);
903 const auto vn = m * normal * inverse_water_depth_sharp(U);
904
905 const auto m_bar = momentum(U_bar);
906 const auto a_bar = speed_of_sound(U_bar);
907 const auto vn_bar = m_bar * normal * inverse_water_depth_sharp(U_bar);
908
909 /* First Riemann characteristic: v * n - 2 * a */
910
911 const auto R_1 = component == 1 ? vn_bar - ScalarNumber(2.) * a_bar
912 : vn - ScalarNumber(2.) * a;
913
914 /* Second Riemann characteristic: v * n + 2 * a */
915
916 const auto R_2 = component == 2 ? vn_bar + ScalarNumber(2.) * a_bar
917 : vn + ScalarNumber(2.) * a;
918
919 const auto vperp = m * inverse_water_depth_sharp(U) - vn * normal;
920
921 const auto vn_new = ScalarNumber(0.5) * (R_1 + R_2);
922
923 const auto h_new =
924 ryujin::fixed_power<2>((R_2 - R_1) / ScalarNumber(4.)) / gravity();
925
926 state_type U_new;
927 U_new[0] = h_new;
928 for (unsigned int d = 0; d < dim; ++d) {
929 U_new[1 + d] = h_new * (vn_new * normal + vperp)[d];
930 }
931
932 return U_new;
933 }
934
935
936 template <int dim, typename Number>
937 template <typename Lambda>
938 DEAL_II_ALWAYS_INLINE inline auto
940 const dealii::types::boundary_id id,
941 const state_type &U,
942 const dealii::Tensor<1, dim, Number> &normal,
943 const Lambda &get_dirichlet_data) const -> state_type
944 {
945 state_type result = U;
946
947 if (id == Boundary::dirichlet) {
948 result = get_dirichlet_data();
949
950 } else if (id == Boundary::dirichlet_momentum) {
951 /* Only enforce Dirichlet conditions on the momentum: */
952 auto m_dirichlet = momentum(get_dirichlet_data());
953 for (unsigned int k = 0; k < dim; ++k)
954 result[k + 1] = m_dirichlet[k];
955
956 } else if (id == Boundary::slip) {
957 auto m = momentum(U);
958 m -= 1. * (m * normal) * normal;
959 for (unsigned int k = 0; k < dim; ++k)
960 result[k + 1] = m[k];
961
962 } else if (id == Boundary::no_slip) {
963 for (unsigned int k = 0; k < dim; ++k)
964 result[k + 1] = Number(0.);
965
966 } else if (id == Boundary::dynamic) {
967 /*
968 * On dynamic boundary conditions, we distinguish four cases:
969 *
970 * - supersonic inflow: prescribe full state
971 * - subsonic inflow:
972 * decompose into Riemann invariants and leave R_2
973 * characteristic untouched.
974 * - supersonic outflow: do nothing
975 * - subsonic outflow:
976 * decompose into Riemann invariants and prescribe incoming
977 * R_1 characteristic.
978 */
979 const auto m = momentum(U);
980 const auto h_inverse = inverse_water_depth_sharp(U);
981 const auto a = speed_of_sound(U);
982 const auto vn = m * normal * h_inverse;
983
984 /* Supersonic inflow: */
985 if (vn < -a) {
986 result = get_dirichlet_data();
987 }
988
989 /* Subsonic inflow: */
990 if (vn >= -a && vn <= 0.) {
991 const auto U_dirichlet = get_dirichlet_data();
992 result = prescribe_riemann_characteristic<2>(U_dirichlet, U, normal);
993 }
994
995 /* Subsonic outflow: */
996 if (vn > 0. && vn <= a) {
997 const auto U_dirichlet = get_dirichlet_data();
998 result = prescribe_riemann_characteristic<1>(U, U_dirichlet, normal);
999 }
1000
1001 /* Supersonic outflow: do nothing, i.e., keep U as is */
1002
1003 } else {
1004 AssertThrow(false, dealii::ExcNotImplemented());
1005 }
1006
1007 return result;
1008 }
1009
1010
1011 template <int dim, typename Number>
1012 DEAL_II_ALWAYS_INLINE inline auto
1014 {
1015 const auto h_inverse = inverse_water_depth_sharp(U);
1016 const auto m = momentum(U);
1017 const auto p = pressure(U);
1018
1019 flux_type result;
1020
1021 result[0] = (m * h_inverse) * U[0];
1022 for (unsigned int i = 0; i < dim; ++i) {
1023 result[1 + i] = (m * h_inverse) * m[i];
1024 result[1 + i][i] += p;
1025 }
1026 return result;
1027 }
1028
1029
1030 template <int dim, typename Number>
1031 DEAL_II_ALWAYS_INLINE inline auto
1033 {
1034 const auto h_inverse = inverse_water_depth_sharp(U);
1035 const auto m = momentum(U);
1036
1037 flux_type result;
1038
1039 result[0] = (m * h_inverse) * U[0];
1040 for (unsigned int i = 0; i < dim; ++i) {
1041 result[1 + i] = (m * h_inverse) * m[i];
1042 }
1043 return result;
1044 }
1045
1046
1047 template <int dim, typename Number>
1048 DEAL_II_ALWAYS_INLINE inline auto
1050 const Number &Z_left,
1051 const Number &Z_right) const
1052 -> state_type
1053 {
1054 const Number Z_max = std::max(Z_left, Z_right);
1055 const Number h = water_depth(U);
1056 const Number H_star = std::max(Number(0.), h + Z_left - Z_max);
1057
1058 return U * H_star * inverse_water_depth_mollified(U);
1059 }
1060
1061
1062 template <int dim, typename Number>
1063 DEAL_II_ALWAYS_INLINE inline auto
1065 const flux_contribution_type &flux_i,
1066 const flux_contribution_type &flux_j) const -> std::array<state_type, 2>
1067 {
1068 const auto &[U_i, Z_i] = flux_i;
1069 const auto &[U_j, Z_j] = flux_j;
1070
1071 const auto U_star_ij = star_state(U_i, Z_i, Z_j);
1072 const auto U_star_ji = star_state(U_j, Z_j, Z_i);
1073
1074 return {U_star_ij, U_star_ji};
1075 }
1076
1077
1078 template <int dim, typename Number>
1079 DEAL_II_ALWAYS_INLINE inline auto
1081 const PrecomputedVector & /*pv*/,
1082 const InitialPrecomputedVector &piv,
1083 const unsigned int i,
1084 const state_type &U_i) const -> flux_contribution_type
1085 {
1086 const auto Z_i = piv.template get_tensor<Number>(i)[0];
1087 return {U_i, Z_i};
1088 }
1089
1090
1091 template <int dim, typename Number>
1092 DEAL_II_ALWAYS_INLINE inline auto
1094 const PrecomputedVector & /*pv*/,
1095 const InitialPrecomputedVector &piv,
1096 const unsigned int *js,
1097 const state_type &U_j) const -> flux_contribution_type
1098 {
1099 const auto Z_j = piv.template get_tensor<Number>(js)[0];
1100 return {U_j, Z_j};
1101 }
1102
1103
1104 template <int dim, typename Number>
1105 DEAL_II_ALWAYS_INLINE inline auto
1107 const flux_contribution_type &flux_i,
1108 const flux_contribution_type &flux_j,
1109 const dealii::Tensor<1, dim, Number> &c_ij) const -> state_type
1110 {
1111 const auto &[U_i, Z_i] = flux_i;
1112 const auto &[U_star_ij, U_star_ji] = equilibrated_states(flux_i, flux_j);
1113
1114 const auto H_i = water_depth(U_i);
1115 const auto H_star_ij = water_depth(U_star_ij);
1116 const auto H_star_ji = water_depth(U_star_ji);
1117
1118 const auto g_i = g(U_star_ij);
1119 const auto g_j = g(U_star_ji);
1120
1121 auto result = -add(g_i, g_j);
1122
1123 const auto factor =
1124 (ScalarNumber(0.5) * (H_star_ji * H_star_ji - H_star_ij * H_star_ij) +
1125 H_i * H_i) *
1126 gravity();
1127
1128 for (unsigned int i = 0; i < dim; ++i) {
1129 result[1 + i][i] -= factor;
1130 }
1131
1132 return contract(result, c_ij);
1133 }
1134
1135
1136 template <int dim, typename Number>
1137 DEAL_II_ALWAYS_INLINE inline auto
1139 const flux_contribution_type &flux_i,
1140 const flux_contribution_type &flux_j,
1141 const dealii::Tensor<1, dim, Number> &c_ij) const -> state_type
1142 {
1143 const auto &[U_i, Z_i] = flux_i;
1144 const auto &[U_j, Z_j] = flux_j;
1145
1146 const auto H_i = water_depth(U_i);
1147 const auto H_j = water_depth(U_j);
1148
1149 const auto g_i = g(U_i);
1150 const auto g_j = g(U_j);
1151
1152 auto result = -add(g_i, g_j);
1153
1154 const auto factor = gravity() * H_i * (H_j + Z_j - Z_i);
1155 for (unsigned int i = 0; i < dim; ++i) {
1156 result[1 + i][i] -= factor;
1157 }
1158
1159 return contract(result, c_ij);
1160 }
1161
1162
1163 template <int dim, typename Number>
1164 DEAL_II_ALWAYS_INLINE inline auto
1166 const flux_contribution_type &flux_i,
1167 const flux_contribution_type &flux_j,
1168 const dealii::Tensor<1, dim, Number> &c_ij,
1169 const Number &d_ij) const -> state_type
1170 {
1171 const auto &[U_i, Z_i] = flux_i;
1172 const auto &[U_j, Z_j] = flux_j;
1173 const auto U_star_ij = star_state(U_i, Z_i, Z_j);
1174
1175 const auto h_inverse = inverse_water_depth_sharp(U_i);
1176 const auto m = momentum(U_i);
1177 const auto factor = ScalarNumber(2.) * (d_ij + h_inverse * (m * c_ij));
1178
1179 return -factor * (U_star_ij - U_i);
1180 }
1181
1182
1183 template <int dim, typename Number>
1184 DEAL_II_ALWAYS_INLINE inline auto
1186 const state_type &U, const Number &h_star, const ScalarNumber tau) const
1187 -> state_type
1188 {
1189 state_type result;
1190
1191 const auto g = gravity();
1192 const auto n = manning_friction_coefficient();
1193
1194 const auto h_inverse = inverse_water_depth_mollified(U);
1195
1196 const auto m = momentum(U);
1197 const auto v_norm = (m * h_inverse).norm();
1198 const auto factor = ScalarNumber(2.) * g * n * n * v_norm;
1199
1200 const auto denominator = h_star + std::max(h_star, tau * factor);
1201 const auto denominator_inverse = ScalarNumber(1.) / denominator;
1202
1203 for (unsigned int d = 0; d < dim; ++d)
1204 result[d + 1] = -factor * denominator_inverse * m[d];
1205
1206 return result;
1207 }
1208
1209
1210 template <int dim, typename Number>
1211 DEAL_II_ALWAYS_INLINE inline auto
1213 const PrecomputedVector &pv,
1214 const unsigned int i,
1215 const state_type &U_i,
1216 const ScalarNumber tau) const -> state_type
1217 {
1218 const auto &[eta_m, h_star] =
1219 pv.template get_tensor<Number, precomputed_type>(i);
1220
1221 return manning_friction(U_i, h_star, tau);
1222 }
1223
1224
1225 template <int dim, typename Number>
1226 DEAL_II_ALWAYS_INLINE inline auto
1228 const PrecomputedVector &pv,
1229 const unsigned int *js,
1230 const state_type &U_j,
1231 const ScalarNumber tau) const -> state_type
1232 {
1233 const auto &[eta_m, h_star] =
1234 pv.template get_tensor<Number, precomputed_type>(js);
1235
1236 return manning_friction(U_j, h_star, tau);
1237 }
1238
1239
1240 template <int dim, typename Number>
1241 template <typename ST>
1242 DEAL_II_ALWAYS_INLINE inline auto
1244 -> state_type
1245 {
1246 using T = typename ST::value_type;
1247 static_assert(std::is_same_v<Number, T>, "template mismatch");
1248
1249 constexpr auto dim2 = ST::dimension - 1;
1250 static_assert(dim >= dim2,
1251 "the space dimension of the argument state must not be "
1252 "larger than the one of the target state");
1253
1254 state_type result;
1255 result[0] = state[0];
1256 for (unsigned int i = 1; i < dim2 + 1; ++i)
1257 result[i] = state[i];
1258
1259 return result;
1260 }
1261
1262 template <int dim, typename Number>
1263 template <typename ST>
1264 DEAL_II_ALWAYS_INLINE inline auto
1266 const ST &initial_state) const -> state_type
1267 {
1268 const auto primitive_state = expand_state(initial_state);
1269 return from_primitive_state(primitive_state);
1270 }
1271
1272
1273 template <int dim, typename Number>
1274 DEAL_II_ALWAYS_INLINE inline auto
1276 const state_type &primitive_state) const -> state_type
1277 {
1278 const auto &h = primitive_state[0];
1279
1280 auto state = primitive_state;
1281 /* Fix up momentum: */
1282 for (unsigned int i = 1; i < dim + 1; ++i)
1283 state[i] *= h;
1284
1285 return state;
1286 }
1287
1288
1289 template <int dim, typename Number>
1290 DEAL_II_ALWAYS_INLINE inline auto
1292 const state_type &state) const -> state_type
1293 {
1294 const auto h_inverse = inverse_water_depth_sharp(state);
1295
1296 auto primitive_state = state;
1297 /* Fix up velocity: */
1298 for (unsigned int i = 1; i < dim + 1; ++i)
1299 primitive_state[i] *= h_inverse;
1300
1301 return primitive_state;
1302 }
1303
1304
1305 template <int dim, typename Number>
1306 template <typename Lambda>
1307 DEAL_II_ALWAYS_INLINE inline auto
1309 const state_type &state, const Lambda &lambda) const -> state_type
1310 {
1311 auto result = state;
1312 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 ShallowWater
1319} // 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)