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>;
307 template <
typename DISPATCH,
typename SPARSITY>
309 const DISPATCH &dispatch_check,
310 const SPARSITY &sparsity_simd,
313 unsigned int right)
const;
422 template <
int component>
426 const dealii::Tensor<1, dim, Number> &normal)
const;
431 template <
typename Lambda>
435 const dealii::Tensor<1, dim, Number> &normal,
436 const Lambda &get_dirichlet_data)
const;
472 const Number &Z_left,
473 const Number &Z_right)
const;
480 std::array<state_type, 2>
507 const unsigned int i,
513 const unsigned int *js,
524 const dealii::Tensor<1, dim, Number> &c_ij)
const;
539 const dealii::Tensor<1, dim, Number> &c_ij)
const;
548 const dealii::Tensor<1, dim, Number> &c_ij,
549 const Number &d_ij)
const;
566 const Number &h_star,
570 const unsigned int i,
575 const unsigned int *js,
596 template <
typename ST>
610 template <
typename ST>
629 template <
typename Lambda>
631 const Lambda &lambda)
const;
644 : ParameterAcceptor(subsection)
647 add_parameter(
"gravity", gravity_,
"Gravitational constant [m/s^2]");
649 manning_friction_coefficient_ = 0.;
650 add_parameter(
"manning friction coefficient",
651 manning_friction_coefficient_,
652 "Roughness coefficient for friction source");
654 reference_water_depth_ = 1.;
655 add_parameter(
"reference water depth",
656 reference_water_depth_,
657 "Problem specific water depth reference");
659 dry_state_relaxation_factor_ = 2.e-1;
660 add_parameter(
"dry state relaxation factor",
661 dry_state_relaxation_factor_,
662 "Problem specific dry-state relaxation parameter");
664 dry_state_relaxation_small_ = 1.e2;
665 add_parameter(
"dry state relaxation small",
666 dry_state_relaxation_small_,
667 "Problem specific dry-state relaxation parameter");
669 dry_state_relaxation_large_ = 1.e4;
670 add_parameter(
"dry state relaxation large",
671 dry_state_relaxation_large_,
672 "Problem specific dry-state relaxation parameter");
676 template <
int dim,
typename Number>
677 template <
typename DISPATCH,
typename SPARSITY>
678 DEAL_II_ALWAYS_INLINE
inline void
680 unsigned int cycle [[maybe_unused]],
681 const DISPATCH &dispatch_check,
682 const SPARSITY &sparsity_simd,
685 unsigned int right)
const
687 Assert(cycle == 0, dealii::ExcInternalError());
689 const auto &U = std::get<0>(state_vector);
690 auto &precomputed = std::get<1>(state_vector);
694 unsigned int stride_size = get_stride_size<Number>;
697 for (
unsigned int i = left; i < right; i += stride_size) {
700 const unsigned int row_length = sparsity_simd.row_length(i);
706 const auto U_i = U.template get_tensor<Number>(i);
708 const auto eta_m = mathematical_entropy(U_i);
710 const auto h_sharp = water_depth_sharp(U_i);
714 precomputed.template write_tensor<Number>(prec_i, i);
719 template <
int dim,
typename Number>
720 DEAL_II_ALWAYS_INLINE
inline Number
727 template <
int dim,
typename Number>
728 DEAL_II_ALWAYS_INLINE
inline Number
732 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
734 const Number h_cutoff_mollified =
735 reference_water_depth() * dry_state_relaxation_large() * Number(eps);
737 const Number h = water_depth(U);
739 const Number h_max = std::max(h, h_cutoff_mollified);
740 const Number denom = h * h + h_max * h_max;
745 template <
int dim,
typename Number>
746 DEAL_II_ALWAYS_INLINE
inline Number
750 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
752 const Number h_cutoff_small =
753 reference_water_depth() * dry_state_relaxation_small() * Number(eps);
755 const Number h = water_depth(U);
756 const Number h_max = std::max(h, h_cutoff_small);
761 template <
int dim,
typename Number>
762 DEAL_II_ALWAYS_INLINE
inline Number
770 template <
int dim,
typename Number>
771 DEAL_II_ALWAYS_INLINE
inline Number
773 const Number &h)
const
776 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
778 const Number h_cutoff_large =
779 reference_water_depth() * dry_state_relaxation_large() * Number(eps);
781 return dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
782 std::abs(h), h_cutoff_large, Number(0.), h);
786 template <
int dim,
typename Number>
787 DEAL_II_ALWAYS_INLINE
inline dealii::Tensor<1, dim, Number>
790 dealii::Tensor<1, dim, Number> result;
792 for (
unsigned int i = 0; i < dim; ++i)
793 result[i] = U[1 + i];
798 template <
int dim,
typename Number>
799 DEAL_II_ALWAYS_INLINE
inline Number
802 const auto h = water_depth(U);
803 const auto vel = momentum(U) * inverse_water_depth_sharp(U);
810 template <
int dim,
typename Number>
811 DEAL_II_ALWAYS_INLINE
inline Number
814 const Number h_sqd = U[0] * U[0];
821 template <
int dim,
typename Number>
822 DEAL_II_ALWAYS_INLINE
inline Number
826 return std::sqrt(gravity() * U[0]);
830 template <
int dim,
typename Number>
831 DEAL_II_ALWAYS_INLINE
inline Number
835 const auto p = pressure(U);
836 const auto k_e = kinetic_energy(U);
841 template <
int dim,
typename Number>
842 DEAL_II_ALWAYS_INLINE
inline auto
859 const Number &h = U[0];
860 const auto vel = momentum(U) * inverse_water_depth_sharp(U);
863 result[0] = gravity() * h -
ScalarNumber(0.5) * vel.norm_square();
866 for (
unsigned int i = 0; i < dim; ++i) {
867 result[1 + i] = vel[i];
874 template <
int dim,
typename Number>
875 DEAL_II_ALWAYS_INLINE
inline bool
878 const auto h = filter_dry_water_depth(water_depth(U));
880 constexpr auto gte = dealii::SIMDComparison::greater_than_or_equal;
881 const auto test = dealii::compare_and_apply_mask<gte>(
882 h, Number(0.), Number(0.), Number(-1.));
885 if (!(test == Number(0.))) {
886 std::cout << std::fixed << std::setprecision(16);
887 std::cout <<
"Bounds violation: Negative state [h] detected!\n";
888 std::cout <<
"\t\th: " << h <<
"\n" << std::endl;
893 return (test == Number(0.));
897 template <
int dim,
typename Number>
898 template <
int component>
899 DEAL_II_ALWAYS_INLINE
inline auto
903 const dealii::Tensor<1, dim, Number> &normal)
const ->
state_type
906 static_assert(component == 1 || component == 2,
907 "component has to be 1 or 2");
911 const auto m = momentum(U);
912 const auto a = speed_of_sound(U);
913 const auto vn = m * normal * inverse_water_depth_sharp(U);
915 const auto m_bar = momentum(U_bar);
916 const auto a_bar = speed_of_sound(U_bar);
917 const auto vn_bar = m_bar * normal * inverse_water_depth_sharp(U_bar);
921 const auto R_1 = component == 1 ? vn_bar -
ScalarNumber(2.) * a_bar
926 const auto R_2 = component == 2 ? vn_bar +
ScalarNumber(2.) * a_bar
929 const auto vperp = m * inverse_water_depth_sharp(U) - vn * normal;
934 ryujin::fixed_power<2>((R_2 - R_1) /
ScalarNumber(4.)) / gravity();
938 for (
unsigned int d = 0; d < dim; ++d) {
939 U_new[1 + d] = h_new * (vn_new * normal + vperp)[d];
946 template <
int dim,
typename Number>
947 template <
typename Lambda>
948 DEAL_II_ALWAYS_INLINE
inline auto
950 const dealii::types::boundary_id
id,
952 const dealii::Tensor<1, dim, Number> &normal,
953 const Lambda &get_dirichlet_data)
const ->
state_type
958 result = get_dirichlet_data();
966 auto m_dir = momentum(get_dirichlet_data());
967 for (
unsigned int k = 0; k < dim; ++k)
968 result[k + 1] = m_dir[k];
971 auto m = momentum(U);
972 m -= 1. * (m * normal) * normal;
973 for (
unsigned int k = 0; k < dim; ++k)
974 result[k + 1] = m[k];
977 for (
unsigned int k = 0; k < dim; ++k)
978 result[k + 1] = Number(0.);
993 const auto m = momentum(U);
994 const auto h_inverse = inverse_water_depth_sharp(U);
995 const auto a = speed_of_sound(U);
996 const auto vn = m * normal * h_inverse;
1000 result = get_dirichlet_data();
1004 if (vn >= -a && vn <= 0.) {
1005 const auto U_dirichlet = get_dirichlet_data();
1006 result = prescribe_riemann_characteristic<2>(U_dirichlet, U, normal);
1010 if (vn > 0. && vn <= a) {
1011 const auto U_dirichlet = get_dirichlet_data();
1012 result = prescribe_riemann_characteristic<1>(U, U_dirichlet, normal);
1022 template <
int dim,
typename Number>
1023 DEAL_II_ALWAYS_INLINE
inline auto
1026 const auto h_inverse = inverse_water_depth_sharp(U);
1027 const auto m = momentum(U);
1028 const auto p = pressure(U);
1032 result[0] = (m * h_inverse) * U[0];
1033 for (
unsigned int i = 0; i < dim; ++i) {
1034 result[1 + i] = (m * h_inverse) * m[i];
1035 result[1 + i][i] += p;
1041 template <
int dim,
typename Number>
1042 DEAL_II_ALWAYS_INLINE
inline auto
1045 const auto h_inverse = inverse_water_depth_sharp(U);
1046 const auto m = momentum(U);
1050 result[0] = (m * h_inverse) * U[0];
1051 for (
unsigned int i = 0; i < dim; ++i) {
1052 result[1 + i] = (m * h_inverse) * m[i];
1058 template <
int dim,
typename Number>
1059 DEAL_II_ALWAYS_INLINE
inline auto
1061 const Number &Z_left,
1062 const Number &Z_right)
const
1065 const Number Z_max = std::max(Z_left, Z_right);
1066 const Number h = water_depth(U);
1067 const Number H_star = std::max(Number(0.), h + Z_left - Z_max);
1069 return U * H_star * inverse_water_depth_mollified(U);
1073 template <
int dim,
typename Number>
1074 DEAL_II_ALWAYS_INLINE
inline auto
1079 const auto &[U_i, Z_i] = flux_i;
1080 const auto &[U_j, Z_j] = flux_j;
1082 const auto U_star_ij = star_state(U_i, Z_i, Z_j);
1083 const auto U_star_ji = star_state(U_j, Z_j, Z_i);
1085 return {U_star_ij, U_star_ji};
1089 template <
int dim,
typename Number>
1090 DEAL_II_ALWAYS_INLINE
inline auto
1094 const unsigned int i,
1097 const auto Z_i = piv.template get_tensor<Number>(i)[0];
1102 template <
int dim,
typename Number>
1103 DEAL_II_ALWAYS_INLINE
inline auto
1107 const unsigned int *js,
1110 const auto Z_j = piv.template get_tensor<Number>(js)[0];
1115 template <
int dim,
typename Number>
1116 DEAL_II_ALWAYS_INLINE
inline auto
1120 const dealii::Tensor<1, dim, Number> &c_ij)
const ->
state_type
1122 const auto &[U_i, Z_i] = flux_i;
1123 const auto &[U_star_ij, U_star_ji] = equilibrated_states(flux_i, flux_j);
1125 const auto H_i = water_depth(U_i);
1126 const auto H_star_ij = water_depth(U_star_ij);
1127 const auto H_star_ji = water_depth(U_star_ji);
1129 const auto g_i = g(U_star_ij);
1130 const auto g_j = g(U_star_ji);
1132 auto result = -
add(g_i, g_j);
1135 (
ScalarNumber(0.5) * (H_star_ji * H_star_ji - H_star_ij * H_star_ij) +
1139 for (
unsigned int i = 0; i < dim; ++i) {
1140 result[1 + i][i] -= factor;
1147 template <
int dim,
typename Number>
1148 DEAL_II_ALWAYS_INLINE
inline auto
1152 const dealii::Tensor<1, dim, Number> &c_ij)
const ->
state_type
1154 const auto &[U_i, Z_i] = flux_i;
1155 const auto &[U_j, Z_j] = flux_j;
1157 const auto H_i = water_depth(U_i);
1158 const auto H_j = water_depth(U_j);
1160 const auto g_i = g(U_i);
1161 const auto g_j = g(U_j);
1163 auto result = -
add(g_i, g_j);
1165 const auto factor = gravity() * H_i * (H_j + Z_j - Z_i);
1166 for (
unsigned int i = 0; i < dim; ++i) {
1167 result[1 + i][i] -= factor;
1174 template <
int dim,
typename Number>
1175 DEAL_II_ALWAYS_INLINE
inline auto
1179 const dealii::Tensor<1, dim, Number> &c_ij,
1182 const auto &[U_i, Z_i] = flux_i;
1183 const auto &[U_j, Z_j] = flux_j;
1184 const auto U_star_ij = star_state(U_i, Z_i, Z_j);
1186 const auto h_inverse = inverse_water_depth_sharp(U_i);
1187 const auto m = momentum(U_i);
1188 const auto factor =
ScalarNumber(2.) * (d_ij + h_inverse * (m * c_ij));
1190 return -factor * (U_star_ij - U_i);
1194 template <
int dim,
typename Number>
1195 DEAL_II_ALWAYS_INLINE
inline auto
1202 const auto g = gravity();
1203 const auto n = manning_friction_coefficient();
1205 const auto h_inverse = inverse_water_depth_mollified(U);
1207 const auto m = momentum(U);
1208 const auto v_norm = (m * h_inverse).norm();
1209 const auto factor =
ScalarNumber(2.) * g * n * n * v_norm;
1211 const auto denominator = h_star + std::max(h_star, tau * factor);
1212 const auto denominator_inverse =
ScalarNumber(1.) / denominator;
1214 for (
unsigned int d = 0; d < dim; ++d)
1215 result[d + 1] = -factor * denominator_inverse * m[d];
1221 template <
int dim,
typename Number>
1222 DEAL_II_ALWAYS_INLINE
inline auto
1225 const unsigned int i,
1229 const auto &[eta_m, h_star] =
1230 pv.template get_tensor<Number, precomputed_type>(i);
1232 return manning_friction(U_i, h_star, tau);
1236 template <
int dim,
typename Number>
1237 DEAL_II_ALWAYS_INLINE
inline auto
1240 const unsigned int *js,
1244 const auto &[eta_m, h_star] =
1245 pv.template get_tensor<Number, precomputed_type>(js);
1247 return manning_friction(U_j, h_star, tau);
1251 template <
int dim,
typename Number>
1252 template <
typename ST>
1253 DEAL_II_ALWAYS_INLINE
inline auto
1257 using T =
typename ST::value_type;
1258 static_assert(std::is_same_v<Number, T>,
"template mismatch");
1260 constexpr auto dim2 = ST::dimension - 1;
1261 static_assert(dim >= dim2,
1262 "the space dimension of the argument state must not be "
1263 "larger than the one of the target state");
1266 result[0] = state[0];
1267 for (
unsigned int i = 1; i < dim2 + 1; ++i)
1268 result[i] = state[i];
1273 template <
int dim,
typename Number>
1274 template <
typename ST>
1275 DEAL_II_ALWAYS_INLINE
inline auto
1279 const auto primitive_state = expand_state(initial_state);
1280 return from_primitive_state(primitive_state);
1284 template <
int dim,
typename Number>
1285 DEAL_II_ALWAYS_INLINE
inline auto
1289 const auto &h = primitive_state[0];
1291 auto state = primitive_state;
1293 for (
unsigned int i = 1; i < dim + 1; ++i)
1300 template <
int dim,
typename Number>
1301 DEAL_II_ALWAYS_INLINE
inline auto
1305 const auto h_inverse = inverse_water_depth_sharp(state);
1307 auto primitive_state = state;
1309 for (
unsigned int i = 1; i < dim + 1; ++i)
1310 primitive_state[i] *= h_inverse;
1312 return primitive_state;
1316 template <
int dim,
typename Number>
1317 template <
typename Lambda>
1318 DEAL_II_ALWAYS_INLINE
inline auto
1322 auto result = state;
1323 auto M = lambda(momentum(state));
1324 for (
unsigned int d = 0; d < dim; ++d)
1325 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)