11#include <deal.II/base/config.h>
19#include <deal.II/base/parameter_acceptor.h>
20#include <deal.II/base/tensor.h>
26 namespace ShallowWater
28 template <
int dim,
typename Number>
29 class HyperbolicSystemView;
46 static inline const std::string
problem_name =
"Shallow water equations";
59 template <
int dim,
typename Number>
71 double manning_friction_coefficient_;
73 double reference_water_depth_;
74 double dry_state_relaxation_factor_;
75 double dry_state_relaxation_small_;
76 double dry_state_relaxation_large_;
78 template <
int dim,
typename Number>
100 template <
int dim,
typename Number>
109 : hyperbolic_system_(hyperbolic_system)
116 template <
int dim2,
typename Number2>
134 return hyperbolic_system_.gravity_;
140 return hyperbolic_system_.manning_friction_coefficient_;
145 return hyperbolic_system_.reference_water_depth_;
151 return hyperbolic_system_.dry_state_relaxation_factor_;
157 return hyperbolic_system_.dry_state_relaxation_small_;
163 return hyperbolic_system_.dry_state_relaxation_large_;
190 using state_type = dealii::Tensor<1, problem_dimension, Number>;
196 dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim, Number>>;
208 []() -> std::array<std::string, problem_dimension> {
209 if constexpr (dim == 1)
211 else if constexpr (dim == 2)
212 return {
"h",
"m_1",
"m_2"};
213 else if constexpr (dim == 3)
214 return {
"h",
"m_1",
"m_2",
"m_3"};
223 []() -> std::array<std::string, problem_dimension> {
224 if constexpr (dim == 1)
226 else if constexpr (dim == 2)
227 return {
"h",
"v_1",
"v_2"};
228 else if constexpr (dim == 3)
229 return {
"h",
"v_1",
"v_2",
"v_3"};
247 std::array<std::string, n_precomputed_values>{
"eta_m",
"h_star"};
258 std::array<Number, n_initial_precomputed_values>;
264 std::array<std::string, n_initial_precomputed_values>{
"bathymetry"};
270 StateVector<ScalarNumber, problem_dimension, n_precomputed_values>;
301 template <
typename DISPATCH,
typename SPARSITY>
303 const DISPATCH &dispatch_check,
304 const SPARSITY &sparsity_simd,
307 unsigned int right)
const;
416 template <
int component>
420 const dealii::Tensor<1, dim, Number> &normal)
const;
425 template <
typename Lambda>
429 const dealii::Tensor<1, dim, Number> &normal,
430 const Lambda &get_dirichlet_data)
const;
466 const Number &Z_left,
467 const Number &Z_right)
const;
474 std::array<state_type, 2>
501 const unsigned int i,
507 const unsigned int *js,
518 const dealii::Tensor<1, dim, Number> &c_ij)
const;
533 const dealii::Tensor<1, dim, Number> &c_ij)
const;
542 const dealii::Tensor<1, dim, Number> &c_ij,
543 const Number &d_ij)
const;
560 const Number &h_star,
564 const unsigned int i,
569 const unsigned int *js,
590 template <
typename ST>
604 template <
typename ST>
623 template <
typename Lambda>
625 const Lambda &lambda)
const;
638 : ParameterAcceptor(subsection)
641 add_parameter(
"gravity", gravity_,
"Gravitational constant [m/s^2]");
643 manning_friction_coefficient_ = 0.;
644 add_parameter(
"manning friction coefficient",
645 manning_friction_coefficient_,
646 "Roughness coefficient for friction source");
648 reference_water_depth_ = 1.;
649 add_parameter(
"reference water depth",
650 reference_water_depth_,
651 "Problem specific water depth reference");
653 dry_state_relaxation_factor_ = 2.e-1;
654 add_parameter(
"dry state relaxation factor",
655 dry_state_relaxation_factor_,
656 "Problem specific dry-state relaxation parameter");
658 dry_state_relaxation_small_ = 1.e2;
659 add_parameter(
"dry state relaxation small",
660 dry_state_relaxation_small_,
661 "Problem specific dry-state relaxation parameter");
663 dry_state_relaxation_large_ = 1.e4;
664 add_parameter(
"dry state relaxation large",
665 dry_state_relaxation_large_,
666 "Problem specific dry-state relaxation parameter");
670 template <
int dim,
typename Number>
671 template <
typename DISPATCH,
typename SPARSITY>
672 DEAL_II_ALWAYS_INLINE
inline void
674 unsigned int cycle [[maybe_unused]],
675 const DISPATCH &dispatch_check,
676 const SPARSITY &sparsity_simd,
679 unsigned int right)
const
681 Assert(cycle == 0, dealii::ExcInternalError());
683 const auto &U = std::get<0>(state_vector);
684 auto &precomputed = std::get<1>(state_vector);
688 unsigned int stride_size = get_stride_size<Number>;
691 for (
unsigned int i = left; i < right; i += stride_size) {
694 const unsigned int row_length = sparsity_simd.row_length(i);
700 const auto U_i = U.template get_tensor<Number>(i);
702 const auto eta_m = mathematical_entropy(U_i);
704 const auto h_sharp = water_depth_sharp(U_i);
708 precomputed.template write_tensor<Number>(prec_i, i);
713 template <
int dim,
typename Number>
714 DEAL_II_ALWAYS_INLINE
inline Number
721 template <
int dim,
typename Number>
722 DEAL_II_ALWAYS_INLINE
inline Number
726 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
728 const Number h_cutoff_mollified =
729 reference_water_depth() * dry_state_relaxation_large() * Number(eps);
731 const Number h = water_depth(U);
733 const Number h_max = std::max(h, h_cutoff_mollified);
734 const Number denom = h * h + h_max * h_max;
739 template <
int dim,
typename Number>
740 DEAL_II_ALWAYS_INLINE
inline Number
744 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
746 const Number h_cutoff_small =
747 reference_water_depth() * dry_state_relaxation_small() * Number(eps);
749 const Number h = water_depth(U);
750 const Number h_max = std::max(h, h_cutoff_small);
755 template <
int dim,
typename Number>
756 DEAL_II_ALWAYS_INLINE
inline Number
764 template <
int dim,
typename Number>
765 DEAL_II_ALWAYS_INLINE
inline Number
767 const Number &h)
const
770 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
772 const Number h_cutoff_large =
773 reference_water_depth() * dry_state_relaxation_large() * Number(eps);
775 return dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
776 std::abs(h), h_cutoff_large, Number(0.), h);
780 template <
int dim,
typename Number>
781 DEAL_II_ALWAYS_INLINE
inline dealii::Tensor<1, dim, Number>
784 dealii::Tensor<1, dim, Number> result;
786 for (
unsigned int i = 0; i < dim; ++i)
787 result[i] = U[1 + i];
792 template <
int dim,
typename Number>
793 DEAL_II_ALWAYS_INLINE
inline Number
796 const auto h = water_depth(U);
797 const auto vel = momentum(U) * inverse_water_depth_sharp(U);
804 template <
int dim,
typename Number>
805 DEAL_II_ALWAYS_INLINE
inline Number
808 const Number h_sqd = U[0] * U[0];
815 template <
int dim,
typename Number>
816 DEAL_II_ALWAYS_INLINE
inline Number
820 return std::sqrt(gravity() * U[0]);
824 template <
int dim,
typename Number>
825 DEAL_II_ALWAYS_INLINE
inline Number
829 const auto p = pressure(U);
830 const auto k_e = kinetic_energy(U);
835 template <
int dim,
typename Number>
836 DEAL_II_ALWAYS_INLINE
inline auto
853 const Number &h = U[0];
854 const auto vel = momentum(U) * inverse_water_depth_sharp(U);
857 result[0] = gravity() * h -
ScalarNumber(0.5) * vel.norm_square();
860 for (
unsigned int i = 0; i < dim; ++i) {
861 result[1 + i] = vel[i];
868 template <
int dim,
typename Number>
869 DEAL_II_ALWAYS_INLINE
inline bool
872 const auto h = filter_dry_water_depth(water_depth(U));
874 constexpr auto gte = dealii::SIMDComparison::greater_than_or_equal;
875 const auto test = dealii::compare_and_apply_mask<gte>(
876 h, Number(0.), Number(0.), Number(-1.));
879 if (!(test == Number(0.))) {
880 std::cout << std::fixed << std::setprecision(16);
881 std::cout <<
"Bounds violation: Negative state [h] detected!\n";
882 std::cout <<
"\t\th: " << h <<
"\n" << std::endl;
887 return (test == Number(0.));
891 template <
int dim,
typename Number>
892 template <
int component>
893 DEAL_II_ALWAYS_INLINE
inline auto
897 const dealii::Tensor<1, dim, Number> &normal)
const ->
state_type
900 static_assert(component == 1 || component == 2,
901 "component has to be 1 or 2");
905 const auto m = momentum(U);
906 const auto a = speed_of_sound(U);
907 const auto vn = m * normal * inverse_water_depth_sharp(U);
909 const auto m_bar = momentum(U_bar);
910 const auto a_bar = speed_of_sound(U_bar);
911 const auto vn_bar = m_bar * normal * inverse_water_depth_sharp(U_bar);
915 const auto R_1 = component == 1 ? vn_bar -
ScalarNumber(2.) * a_bar
920 const auto R_2 = component == 2 ? vn_bar +
ScalarNumber(2.) * a_bar
923 const auto vperp = m * inverse_water_depth_sharp(U) - vn * normal;
928 ryujin::fixed_power<2>((R_2 - R_1) /
ScalarNumber(4.)) / gravity();
932 for (
unsigned int d = 0; d < dim; ++d) {
933 U_new[1 + d] = h_new * (vn_new * normal + vperp)[d];
940 template <
int dim,
typename Number>
941 template <
typename Lambda>
942 DEAL_II_ALWAYS_INLINE
inline auto
944 const dealii::types::boundary_id
id,
946 const dealii::Tensor<1, dim, Number> &normal,
947 const Lambda &get_dirichlet_data)
const ->
state_type
952 result = get_dirichlet_data();
960 auto m_dir = momentum(get_dirichlet_data());
961 for (
unsigned int k = 0; k < dim; ++k)
962 result[k + 1] = m_dir[k];
965 auto m = momentum(U);
966 m -= 1. * (m * normal) * normal;
967 for (
unsigned int k = 0; k < dim; ++k)
968 result[k + 1] = m[k];
971 for (
unsigned int k = 0; k < dim; ++k)
972 result[k + 1] = Number(0.);
987 const auto m = momentum(U);
988 const auto h_inverse = inverse_water_depth_sharp(U);
989 const auto a = speed_of_sound(U);
990 const auto vn = m * normal * h_inverse;
994 result = get_dirichlet_data();
998 if (vn >= -a && vn <= 0.) {
999 const auto U_dirichlet = get_dirichlet_data();
1000 result = prescribe_riemann_characteristic<2>(U_dirichlet, U, normal);
1004 if (vn > 0. && vn <= a) {
1005 const auto U_dirichlet = get_dirichlet_data();
1006 result = prescribe_riemann_characteristic<1>(U, U_dirichlet, normal);
1016 template <
int dim,
typename Number>
1017 DEAL_II_ALWAYS_INLINE
inline auto
1020 const auto h_inverse = inverse_water_depth_sharp(U);
1021 const auto m = momentum(U);
1022 const auto p = pressure(U);
1026 result[0] = (m * h_inverse) * U[0];
1027 for (
unsigned int i = 0; i < dim; ++i) {
1028 result[1 + i] = (m * h_inverse) * m[i];
1029 result[1 + i][i] += p;
1035 template <
int dim,
typename Number>
1036 DEAL_II_ALWAYS_INLINE
inline auto
1039 const auto h_inverse = inverse_water_depth_sharp(U);
1040 const auto m = momentum(U);
1044 result[0] = (m * h_inverse) * U[0];
1045 for (
unsigned int i = 0; i < dim; ++i) {
1046 result[1 + i] = (m * h_inverse) * m[i];
1052 template <
int dim,
typename Number>
1053 DEAL_II_ALWAYS_INLINE
inline auto
1055 const Number &Z_left,
1056 const Number &Z_right)
const
1059 const Number Z_max = std::max(Z_left, Z_right);
1060 const Number h = water_depth(U);
1061 const Number H_star = std::max(Number(0.), h + Z_left - Z_max);
1063 return U * H_star * inverse_water_depth_mollified(U);
1067 template <
int dim,
typename Number>
1068 DEAL_II_ALWAYS_INLINE
inline auto
1073 const auto &[U_i, Z_i] = flux_i;
1074 const auto &[U_j, Z_j] = flux_j;
1076 const auto U_star_ij = star_state(U_i, Z_i, Z_j);
1077 const auto U_star_ji = star_state(U_j, Z_j, Z_i);
1079 return {U_star_ij, U_star_ji};
1083 template <
int dim,
typename Number>
1084 DEAL_II_ALWAYS_INLINE
inline auto
1088 const unsigned int i,
1091 const auto Z_i = piv.template get_tensor<Number>(i)[0];
1096 template <
int dim,
typename Number>
1097 DEAL_II_ALWAYS_INLINE
inline auto
1101 const unsigned int *js,
1104 const auto Z_j = piv.template get_tensor<Number>(js)[0];
1109 template <
int dim,
typename Number>
1110 DEAL_II_ALWAYS_INLINE
inline auto
1114 const dealii::Tensor<1, dim, Number> &c_ij)
const ->
state_type
1116 const auto &[U_i, Z_i] = flux_i;
1117 const auto &[U_star_ij, U_star_ji] = equilibrated_states(flux_i, flux_j);
1119 const auto H_i = water_depth(U_i);
1120 const auto H_star_ij = water_depth(U_star_ij);
1121 const auto H_star_ji = water_depth(U_star_ji);
1123 const auto g_i = g(U_star_ij);
1124 const auto g_j = g(U_star_ji);
1126 auto result = -
add(g_i, g_j);
1129 (
ScalarNumber(0.5) * (H_star_ji * H_star_ji - H_star_ij * H_star_ij) +
1133 for (
unsigned int i = 0; i < dim; ++i) {
1134 result[1 + i][i] -= factor;
1141 template <
int dim,
typename Number>
1142 DEAL_II_ALWAYS_INLINE
inline auto
1146 const dealii::Tensor<1, dim, Number> &c_ij)
const ->
state_type
1148 const auto &[U_i, Z_i] = flux_i;
1149 const auto &[U_j, Z_j] = flux_j;
1151 const auto H_i = water_depth(U_i);
1152 const auto H_j = water_depth(U_j);
1154 const auto g_i = g(U_i);
1155 const auto g_j = g(U_j);
1157 auto result = -
add(g_i, g_j);
1159 const auto factor = gravity() * H_i * (H_j + Z_j - Z_i);
1160 for (
unsigned int i = 0; i < dim; ++i) {
1161 result[1 + i][i] -= factor;
1168 template <
int dim,
typename Number>
1169 DEAL_II_ALWAYS_INLINE
inline auto
1173 const dealii::Tensor<1, dim, Number> &c_ij,
1176 const auto &[U_i, Z_i] = flux_i;
1177 const auto &[U_j, Z_j] = flux_j;
1178 const auto U_star_ij = star_state(U_i, Z_i, Z_j);
1180 const auto h_inverse = inverse_water_depth_sharp(U_i);
1181 const auto m = momentum(U_i);
1182 const auto factor =
ScalarNumber(2.) * (d_ij + h_inverse * (m * c_ij));
1184 return -factor * (U_star_ij - U_i);
1188 template <
int dim,
typename Number>
1189 DEAL_II_ALWAYS_INLINE
inline auto
1196 const auto g = gravity();
1197 const auto n = manning_friction_coefficient();
1199 const auto h_inverse = inverse_water_depth_mollified(U);
1201 const auto m = momentum(U);
1202 const auto v_norm = (m * h_inverse).norm();
1203 const auto factor =
ScalarNumber(2.) * g * n * n * v_norm;
1205 const auto denominator = h_star + std::max(h_star, tau * factor);
1206 const auto denominator_inverse =
ScalarNumber(1.) / denominator;
1208 for (
unsigned int d = 0; d < dim; ++d)
1209 result[d + 1] = -factor * denominator_inverse * m[d];
1215 template <
int dim,
typename Number>
1216 DEAL_II_ALWAYS_INLINE
inline auto
1219 const unsigned int i,
1223 const auto &[eta_m, h_star] =
1224 pv.template get_tensor<Number, precomputed_type>(i);
1226 return manning_friction(U_i, h_star, tau);
1230 template <
int dim,
typename Number>
1231 DEAL_II_ALWAYS_INLINE
inline auto
1234 const unsigned int *js,
1238 const auto &[eta_m, h_star] =
1239 pv.template get_tensor<Number, precomputed_type>(js);
1241 return manning_friction(U_j, h_star, tau);
1245 template <
int dim,
typename Number>
1246 template <
typename ST>
1247 DEAL_II_ALWAYS_INLINE
inline auto
1251 using T =
typename ST::value_type;
1252 static_assert(std::is_same_v<Number, T>,
"template mismatch");
1254 constexpr auto dim2 = ST::dimension - 1;
1255 static_assert(dim >= dim2,
1256 "the space dimension of the argument state must not be "
1257 "larger than the one of the target state");
1260 result[0] = state[0];
1261 for (
unsigned int i = 1; i < dim2 + 1; ++i)
1262 result[i] = state[i];
1267 template <
int dim,
typename Number>
1268 template <
typename ST>
1269 DEAL_II_ALWAYS_INLINE
inline auto
1273 const auto primitive_state = expand_state(initial_state);
1274 return from_primitive_state(primitive_state);
1278 template <
int dim,
typename Number>
1279 DEAL_II_ALWAYS_INLINE
inline auto
1283 const auto &h = primitive_state[0];
1285 auto state = primitive_state;
1287 for (
unsigned int i = 1; i < dim + 1; ++i)
1294 template <
int dim,
typename Number>
1295 DEAL_II_ALWAYS_INLINE
inline auto
1299 const auto h_inverse = inverse_water_depth_sharp(state);
1301 auto primitive_state = state;
1303 for (
unsigned int i = 1; i < dim + 1; ++i)
1304 primitive_state[i] *= h_inverse;
1306 return primitive_state;
1310 template <
int dim,
typename Number>
1311 template <
typename Lambda>
1312 DEAL_II_ALWAYS_INLINE
inline auto
1316 auto result = state;
1317 auto M = lambda(momentum(state));
1318 for (
unsigned int d = 0; d < dim; ++d)
1319 result[1 + d] = M[d];
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
static const auto precomputed_names
static const auto initial_precomputed_names
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
void precomputation_loop(unsigned int cycle, const DISPATCH &dispatch_check, const SPARSITY &sparsity_simd, StateVector &state_vector, unsigned int left, unsigned int right) 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
static const auto component_names
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
static constexpr bool have_high_order_flux
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
static constexpr bool have_source_terms
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 const auto primitive_component_names
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
DEAL_II_ALWAYS_INLINE ScalarNumber dry_state_relaxation_factor() 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")
static const std::string problem_name
T pow(const T x, const T b)
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)
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)