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_small_;
75 double dry_state_relaxation_large_;
77 template <
int dim,
typename Number>
99 template <
int dim,
typename Number>
108 : hyperbolic_system_(hyperbolic_system)
115 template <
int dim2,
typename Number2>
133 return hyperbolic_system_.gravity_;
139 return hyperbolic_system_.manning_friction_coefficient_;
144 return hyperbolic_system_.reference_water_depth_;
150 return hyperbolic_system_.dry_state_relaxation_small_;
156 return hyperbolic_system_.dry_state_relaxation_large_;
183 using state_type = dealii::Tensor<1, problem_dimension, Number>;
189 dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim, Number>>;
201 []() -> std::array<std::string, problem_dimension> {
202 if constexpr (dim == 1)
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"};
216 []() -> std::array<std::string, problem_dimension> {
217 if constexpr (dim == 1)
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"};
240 std::array<std::string, n_precomputed_values>{
"eta_m",
"h_star"};
251 std::array<Number, n_initial_precomputed_values>;
257 std::array<std::string, n_initial_precomputed_values>{
"bathymetry"};
263 StateVector<ScalarNumber, problem_dimension, n_precomputed_values>;
300 template <
typename DISPATCH,
typename SPARSITY>
302 const DISPATCH &dispatch_check,
303 const SPARSITY &sparsity_simd,
307 const bool skip_constrained_dofs =
true)
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_small_ = 1.e2;
654 add_parameter(
"dry state relaxation small",
655 dry_state_relaxation_small_,
656 "Problem specific dry-state relaxation parameter");
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");
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,
675 const bool skip_constrained_dofs )
const
677 Assert(cycle == 0, dealii::ExcInternalError());
679 const auto &U = std::get<0>(state_vector);
680 auto &precomputed = std::get<1>(state_vector);
684 unsigned int stride_size = get_stride_size<Number>;
687 for (
unsigned int i = left; i < right; i += stride_size) {
690 const unsigned int row_length = sparsity_simd.row_length(i);
691 if (skip_constrained_dofs && row_length == 1)
696 const auto U_i = U.template get_tensor<Number>(i);
698 const auto eta_m = mathematical_entropy(U_i);
700 const auto h_sharp = water_depth_sharp(U_i);
704 precomputed.template write_tensor<Number>(prec_i, i);
709 template <
int dim,
typename Number>
710 DEAL_II_ALWAYS_INLINE
inline Number
717 template <
int dim,
typename Number>
718 DEAL_II_ALWAYS_INLINE
inline Number
722 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
724 const Number h_cutoff_mollified =
725 reference_water_depth() * dry_state_relaxation_large() * Number(eps);
727 const Number h = water_depth(U);
729 const Number h_max = std::max(h, h_cutoff_mollified);
730 const Number denom = h * h + h_max * h_max;
735 template <
int dim,
typename Number>
736 DEAL_II_ALWAYS_INLINE
inline Number
740 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
742 const Number h_cutoff_small =
743 reference_water_depth() * dry_state_relaxation_small() * Number(eps);
745 const Number h = water_depth(U);
746 const Number h_max = std::max(h, h_cutoff_small);
751 template <
int dim,
typename Number>
752 DEAL_II_ALWAYS_INLINE
inline Number
760 template <
int dim,
typename Number>
761 DEAL_II_ALWAYS_INLINE
inline Number
763 const Number &h)
const
766 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
768 const Number h_cutoff_large =
769 reference_water_depth() * dry_state_relaxation_large() * Number(eps);
771 return dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
772 std::abs(h), h_cutoff_large, Number(0.), h);
776 template <
int dim,
typename Number>
777 DEAL_II_ALWAYS_INLINE
inline dealii::Tensor<1, dim, Number>
780 dealii::Tensor<1, dim, Number> result;
782 for (
unsigned int i = 0; i < dim; ++i)
783 result[i] = U[1 + i];
788 template <
int dim,
typename Number>
789 DEAL_II_ALWAYS_INLINE
inline Number
792 const auto h = water_depth(U);
793 const auto vel = momentum(U) * inverse_water_depth_sharp(U);
800 template <
int dim,
typename Number>
801 DEAL_II_ALWAYS_INLINE
inline Number
804 const Number h_sqd = U[0] * U[0];
811 template <
int dim,
typename Number>
812 DEAL_II_ALWAYS_INLINE
inline Number
816 return std::sqrt(gravity() * U[0]);
820 template <
int dim,
typename Number>
821 DEAL_II_ALWAYS_INLINE
inline Number
825 const auto p = pressure(U);
826 const auto k_e = kinetic_energy(U);
831 template <
int dim,
typename Number>
832 DEAL_II_ALWAYS_INLINE
inline auto
849 const Number &h = U[0];
850 const auto vel = momentum(U) * inverse_water_depth_sharp(U);
853 result[0] = gravity() * h -
ScalarNumber(0.5) * vel.norm_square();
856 for (
unsigned int i = 0; i < dim; ++i) {
857 result[1 + i] = vel[i];
864 template <
int dim,
typename Number>
865 DEAL_II_ALWAYS_INLINE
inline bool
868 const auto h = filter_dry_water_depth(water_depth(U));
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.));
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;
883 return (test == Number(0.));
887 template <
int dim,
typename Number>
888 template <
int component>
889 DEAL_II_ALWAYS_INLINE
inline auto
893 const dealii::Tensor<1, dim, Number> &normal)
const ->
state_type
896 static_assert(component == 1 || component == 2,
897 "component has to be 1 or 2");
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);
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);
911 const auto R_1 = component == 1 ? vn_bar -
ScalarNumber(2.) * a_bar
916 const auto R_2 = component == 2 ? vn_bar +
ScalarNumber(2.) * a_bar
919 const auto vperp = m * inverse_water_depth_sharp(U) - vn * normal;
924 ryujin::fixed_power<2>((R_2 - R_1) /
ScalarNumber(4.)) / gravity();
928 for (
unsigned int d = 0; d < dim; ++d) {
929 U_new[1 + d] = h_new * (vn_new * normal + vperp)[d];
936 template <
int dim,
typename Number>
937 template <
typename Lambda>
938 DEAL_II_ALWAYS_INLINE
inline auto
940 const dealii::types::boundary_id
id,
942 const dealii::Tensor<1, dim, Number> &normal,
943 const Lambda &get_dirichlet_data)
const ->
state_type
948 result = get_dirichlet_data();
952 auto m_dirichlet = momentum(get_dirichlet_data());
953 for (
unsigned int k = 0; k < dim; ++k)
954 result[k + 1] = m_dirichlet[k];
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];
963 for (
unsigned int k = 0; k < dim; ++k)
964 result[k + 1] = Number(0.);
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;
986 result = get_dirichlet_data();
990 if (vn >= -a && vn <= 0.) {
991 const auto U_dirichlet = get_dirichlet_data();
992 result = prescribe_riemann_characteristic<2>(U_dirichlet, U, normal);
996 if (vn > 0. && vn <= a) {
997 const auto U_dirichlet = get_dirichlet_data();
998 result = prescribe_riemann_characteristic<1>(U, U_dirichlet, normal);
1004 AssertThrow(
false, dealii::ExcNotImplemented());
1011 template <
int dim,
typename Number>
1012 DEAL_II_ALWAYS_INLINE
inline auto
1015 const auto h_inverse = inverse_water_depth_sharp(U);
1016 const auto m = momentum(U);
1017 const auto p = pressure(U);
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;
1030 template <
int dim,
typename Number>
1031 DEAL_II_ALWAYS_INLINE
inline auto
1034 const auto h_inverse = inverse_water_depth_sharp(U);
1035 const auto m = momentum(U);
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];
1047 template <
int dim,
typename Number>
1048 DEAL_II_ALWAYS_INLINE
inline auto
1050 const Number &Z_left,
1051 const Number &Z_right)
const
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);
1058 return U * H_star * inverse_water_depth_mollified(U);
1062 template <
int dim,
typename Number>
1063 DEAL_II_ALWAYS_INLINE
inline auto
1068 const auto &[U_i, Z_i] = flux_i;
1069 const auto &[U_j, Z_j] = flux_j;
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);
1074 return {U_star_ij, U_star_ji};
1078 template <
int dim,
typename Number>
1079 DEAL_II_ALWAYS_INLINE
inline auto
1083 const unsigned int i,
1086 const auto Z_i = piv.template get_tensor<Number>(i)[0];
1091 template <
int dim,
typename Number>
1092 DEAL_II_ALWAYS_INLINE
inline auto
1096 const unsigned int *js,
1099 const auto Z_j = piv.template get_tensor<Number>(js)[0];
1104 template <
int dim,
typename Number>
1105 DEAL_II_ALWAYS_INLINE
inline auto
1109 const dealii::Tensor<1, dim, Number> &c_ij)
const ->
state_type
1111 const auto &[U_i, Z_i] = flux_i;
1112 const auto &[U_star_ij, U_star_ji] = equilibrated_states(flux_i, flux_j);
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);
1118 const auto g_i = g(U_star_ij);
1119 const auto g_j = g(U_star_ji);
1121 auto result = -
add(g_i, g_j);
1124 (
ScalarNumber(0.5) * (H_star_ji * H_star_ji - H_star_ij * H_star_ij) +
1128 for (
unsigned int i = 0; i < dim; ++i) {
1129 result[1 + i][i] -= factor;
1136 template <
int dim,
typename Number>
1137 DEAL_II_ALWAYS_INLINE
inline auto
1141 const dealii::Tensor<1, dim, Number> &c_ij)
const ->
state_type
1143 const auto &[U_i, Z_i] = flux_i;
1144 const auto &[U_j, Z_j] = flux_j;
1146 const auto H_i = water_depth(U_i);
1147 const auto H_j = water_depth(U_j);
1149 const auto g_i = g(U_i);
1150 const auto g_j = g(U_j);
1152 auto result = -
add(g_i, g_j);
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;
1163 template <
int dim,
typename Number>
1164 DEAL_II_ALWAYS_INLINE
inline auto
1168 const dealii::Tensor<1, dim, Number> &c_ij,
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);
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));
1179 return -factor * (U_star_ij - U_i);
1183 template <
int dim,
typename Number>
1184 DEAL_II_ALWAYS_INLINE
inline auto
1191 const auto g = gravity();
1192 const auto n = manning_friction_coefficient();
1194 const auto h_inverse = inverse_water_depth_mollified(U);
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;
1200 const auto denominator = h_star + std::max(h_star, tau * factor);
1201 const auto denominator_inverse =
ScalarNumber(1.) / denominator;
1203 for (
unsigned int d = 0; d < dim; ++d)
1204 result[d + 1] = -factor * denominator_inverse * m[d];
1210 template <
int dim,
typename Number>
1211 DEAL_II_ALWAYS_INLINE
inline auto
1214 const unsigned int i,
1218 const auto &[eta_m, h_star] =
1219 pv.template get_tensor<Number, precomputed_type>(i);
1221 return manning_friction(U_i, h_star, tau);
1225 template <
int dim,
typename Number>
1226 DEAL_II_ALWAYS_INLINE
inline auto
1229 const unsigned int *js,
1233 const auto &[eta_m, h_star] =
1234 pv.template get_tensor<Number, precomputed_type>(js);
1236 return manning_friction(U_j, h_star, tau);
1240 template <
int dim,
typename Number>
1241 template <
typename ST>
1242 DEAL_II_ALWAYS_INLINE
inline auto
1246 using T =
typename ST::value_type;
1247 static_assert(std::is_same_v<Number, T>,
"template mismatch");
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");
1255 result[0] = state[0];
1256 for (
unsigned int i = 1; i < dim2 + 1; ++i)
1257 result[i] = state[i];
1262 template <
int dim,
typename Number>
1263 template <
typename ST>
1264 DEAL_II_ALWAYS_INLINE
inline auto
1268 const auto primitive_state = expand_state(initial_state);
1269 return from_primitive_state(primitive_state);
1273 template <
int dim,
typename Number>
1274 DEAL_II_ALWAYS_INLINE
inline auto
1278 const auto &h = primitive_state[0];
1280 auto state = primitive_state;
1282 for (
unsigned int i = 1; i < dim + 1; ++i)
1289 template <
int dim,
typename Number>
1290 DEAL_II_ALWAYS_INLINE
inline auto
1294 const auto h_inverse = inverse_water_depth_sharp(state);
1296 auto primitive_state = state;
1298 for (
unsigned int i = 1; i < dim + 1; ++i)
1299 primitive_state[i] *= h_inverse;
1301 return primitive_state;
1305 template <
int dim,
typename Number>
1306 template <
typename Lambda>
1307 DEAL_II_ALWAYS_INLINE
inline auto
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];
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
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
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")
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)