10#include <compile_time_options.h>
13#include <deal.II/base/config.h>
21#include <deal.II/base/parameter_acceptor.h>
22#include <deal.II/base/tensor.h>
28 namespace ShallowWater
30 template <
int dim,
typename Number>
31 class HyperbolicSystemView;
48 static inline const std::string
problem_name =
"Shallow water equations";
61 template <
int dim,
typename Number>
73 double manning_friction_coefficient_;
75 double reference_water_depth_;
76 double dry_state_relaxation_small_;
77 double dry_state_relaxation_large_;
79 template <
int dim,
typename Number>
101 template <
int dim,
typename Number>
110 : hyperbolic_system_(hyperbolic_system)
117 template <
int dim2,
typename Number2>
135 return hyperbolic_system_.gravity_;
141 return hyperbolic_system_.manning_friction_coefficient_;
146 return hyperbolic_system_.reference_water_depth_;
152 return hyperbolic_system_.dry_state_relaxation_small_;
158 return hyperbolic_system_.dry_state_relaxation_large_;
185 using state_type = dealii::Tensor<1, problem_dimension, Number>;
191 dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim, Number>>;
203 []() -> std::array<std::string, problem_dimension> {
204 if constexpr (dim == 1)
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"};
218 []() -> std::array<std::string, problem_dimension> {
219 if constexpr (dim == 1)
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"};
242 std::array<std::string, n_precomputed_values>{
"eta_m",
"h_star"};
253 std::array<Number, n_initial_precomputed_values>;
259 std::array<std::string, n_initial_precomputed_values>{
"bathymetry"};
265 StateVector<ScalarNumber, problem_dimension, n_precomputed_values>;
302 template <
typename DISPATCH,
typename SPARSITY>
304 const DISPATCH &dispatch_check,
305 const SPARSITY &sparsity_simd,
309 const bool skip_constrained_dofs =
true)
const;
418 template <
int component>
422 const dealii::Tensor<1, dim, Number> &normal)
const;
427 template <
typename Lambda>
431 const dealii::Tensor<1, dim, Number> &normal,
432 const Lambda &get_dirichlet_data)
const;
468 const Number &Z_left,
469 const Number &Z_right)
const;
476 std::array<state_type, 2>
503 const unsigned int i,
509 const unsigned int *js,
520 const dealii::Tensor<1, dim, Number> &c_ij)
const;
535 const dealii::Tensor<1, dim, Number> &c_ij)
const;
544 const dealii::Tensor<1, dim, Number> &c_ij,
545 const Number &d_ij)
const;
562 const Number &h_star,
566 const unsigned int i,
571 const unsigned int *js,
592 template <
typename ST>
606 template <
typename ST>
625 template <
typename Lambda>
627 const Lambda &lambda)
const;
640 : ParameterAcceptor(subsection)
643 add_parameter(
"gravity", gravity_,
"Gravitational constant [m/s^2]");
645 manning_friction_coefficient_ = 0.;
646 add_parameter(
"manning friction coefficient",
647 manning_friction_coefficient_,
648 "Roughness coefficient for friction source");
650 reference_water_depth_ = 1.;
651 add_parameter(
"reference water depth",
652 reference_water_depth_,
653 "Problem specific water depth reference");
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");
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");
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,
677 const bool skip_constrained_dofs )
const
679 Assert(cycle == 0, dealii::ExcInternalError());
681 const auto &U = std::get<0>(state_vector);
682 auto &precomputed = std::get<1>(state_vector);
686 unsigned int stride_size = get_stride_size<Number>;
689 for (
unsigned int i = left; i < right; i += stride_size) {
692 const unsigned int row_length = sparsity_simd.row_length(i);
693 if (skip_constrained_dofs && row_length == 1)
698 const auto U_i = U.template get_tensor<Number>(i);
700 const auto eta_m = mathematical_entropy(U_i);
702 const auto h_sharp = water_depth_sharp(U_i);
706 precomputed.template write_tensor<Number>(prec_i, i);
711 template <
int dim,
typename Number>
712 DEAL_II_ALWAYS_INLINE
inline Number
719 template <
int dim,
typename Number>
720 DEAL_II_ALWAYS_INLINE
inline Number
724 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
726 const Number h_cutoff_mollified =
727 reference_water_depth() * dry_state_relaxation_large() * Number(eps);
729 const Number h = water_depth(U);
731 const Number h_max = std::max(h, h_cutoff_mollified);
732 const Number denom = h * h + h_max * h_max;
737 template <
int dim,
typename Number>
738 DEAL_II_ALWAYS_INLINE
inline Number
742 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
744 const Number h_cutoff_small =
745 reference_water_depth() * dry_state_relaxation_small() * Number(eps);
747 const Number h = water_depth(U);
748 const Number h_max = std::max(h, h_cutoff_small);
753 template <
int dim,
typename Number>
754 DEAL_II_ALWAYS_INLINE
inline Number
762 template <
int dim,
typename Number>
763 DEAL_II_ALWAYS_INLINE
inline Number
765 const Number &h)
const
768 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
770 const Number h_cutoff_large =
771 reference_water_depth() * dry_state_relaxation_large() * Number(eps);
773 return dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
774 std::abs(h), h_cutoff_large, Number(0.), h);
778 template <
int dim,
typename Number>
779 DEAL_II_ALWAYS_INLINE
inline dealii::Tensor<1, dim, Number>
782 dealii::Tensor<1, dim, Number> result;
784 for (
unsigned int i = 0; i < dim; ++i)
785 result[i] = U[1 + i];
790 template <
int dim,
typename Number>
791 DEAL_II_ALWAYS_INLINE
inline Number
794 const auto h = water_depth(U);
795 const auto vel = momentum(U) * inverse_water_depth_sharp(U);
802 template <
int dim,
typename Number>
803 DEAL_II_ALWAYS_INLINE
inline Number
806 const Number h_sqd = U[0] * U[0];
813 template <
int dim,
typename Number>
814 DEAL_II_ALWAYS_INLINE
inline Number
818 return std::sqrt(gravity() * U[0]);
822 template <
int dim,
typename Number>
823 DEAL_II_ALWAYS_INLINE
inline Number
827 const auto p = pressure(U);
828 const auto k_e = kinetic_energy(U);
833 template <
int dim,
typename Number>
834 DEAL_II_ALWAYS_INLINE
inline auto
851 const Number &h = U[0];
852 const auto vel = momentum(U) * inverse_water_depth_sharp(U);
855 result[0] = gravity() * h -
ScalarNumber(0.5) * vel.norm_square();
858 for (
unsigned int i = 0; i < dim; ++i) {
859 result[1 + i] = vel[i];
866 template <
int dim,
typename Number>
867 DEAL_II_ALWAYS_INLINE
inline bool
870 const auto h = filter_dry_water_depth(water_depth(U));
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.));
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;
885 return (test == Number(0.));
889 template <
int dim,
typename Number>
890 template <
int component>
891 DEAL_II_ALWAYS_INLINE
inline auto
895 const dealii::Tensor<1, dim, Number> &normal)
const ->
state_type
898 static_assert(component == 1 || component == 2,
899 "component has to be 1 or 2");
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);
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);
913 const auto R_1 = component == 1 ? vn_bar -
ScalarNumber(2.) * a_bar
918 const auto R_2 = component == 2 ? vn_bar +
ScalarNumber(2.) * a_bar
921 const auto vperp = m * inverse_water_depth_sharp(U) - vn * normal;
926 ryujin::fixed_power<2>((R_2 - R_1) /
ScalarNumber(4.)) / gravity();
930 for (
unsigned int d = 0; d < dim; ++d) {
931 U_new[1 + d] = h_new * (vn_new * normal + vperp)[d];
938 template <
int dim,
typename Number>
939 template <
typename Lambda>
940 DEAL_II_ALWAYS_INLINE
inline auto
942 const dealii::types::boundary_id
id,
944 const dealii::Tensor<1, dim, Number> &normal,
945 const Lambda &get_dirichlet_data)
const ->
state_type
950 result = get_dirichlet_data();
954 auto m_dirichlet = momentum(get_dirichlet_data());
955 for (
unsigned int k = 0; k < dim; ++k)
956 result[k + 1] = m_dirichlet[k];
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];
965 for (
unsigned int k = 0; k < dim; ++k)
966 result[k + 1] = Number(0.);
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;
988 result = get_dirichlet_data();
992 if (vn >= -a && vn <= 0.) {
993 const auto U_dirichlet = get_dirichlet_data();
994 result = prescribe_riemann_characteristic<2>(U_dirichlet, U, normal);
998 if (vn > 0. && vn <= a) {
999 const auto U_dirichlet = get_dirichlet_data();
1000 result = prescribe_riemann_characteristic<1>(U, U_dirichlet, normal);
1006 AssertThrow(
false, dealii::ExcNotImplemented());
1013 template <
int dim,
typename Number>
1014 DEAL_II_ALWAYS_INLINE
inline auto
1017 const auto h_inverse = inverse_water_depth_sharp(U);
1018 const auto m = momentum(U);
1019 const auto p = pressure(U);
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;
1032 template <
int dim,
typename Number>
1033 DEAL_II_ALWAYS_INLINE
inline auto
1036 const auto h_inverse = inverse_water_depth_sharp(U);
1037 const auto m = momentum(U);
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];
1049 template <
int dim,
typename Number>
1050 DEAL_II_ALWAYS_INLINE
inline auto
1052 const Number &Z_left,
1053 const Number &Z_right)
const
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);
1060 return U * H_star * inverse_water_depth_mollified(U);
1064 template <
int dim,
typename Number>
1065 DEAL_II_ALWAYS_INLINE
inline auto
1070 const auto &[U_i, Z_i] = flux_i;
1071 const auto &[U_j, Z_j] = flux_j;
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);
1076 return {U_star_ij, U_star_ji};
1080 template <
int dim,
typename Number>
1081 DEAL_II_ALWAYS_INLINE
inline auto
1085 const unsigned int i,
1088 const auto Z_i = piv.template get_tensor<Number>(i)[0];
1093 template <
int dim,
typename Number>
1094 DEAL_II_ALWAYS_INLINE
inline auto
1098 const unsigned int *js,
1101 const auto Z_j = piv.template get_tensor<Number>(js)[0];
1106 template <
int dim,
typename Number>
1107 DEAL_II_ALWAYS_INLINE
inline auto
1111 const dealii::Tensor<1, dim, Number> &c_ij)
const ->
state_type
1113 const auto &[U_i, Z_i] = flux_i;
1114 const auto &[U_star_ij, U_star_ji] = equilibrated_states(flux_i, flux_j);
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);
1120 const auto g_i = g(U_star_ij);
1121 const auto g_j = g(U_star_ji);
1123 auto result = -
add(g_i, g_j);
1126 (
ScalarNumber(0.5) * (H_star_ji * H_star_ji - H_star_ij * H_star_ij) +
1130 for (
unsigned int i = 0; i < dim; ++i) {
1131 result[1 + i][i] -= factor;
1138 template <
int dim,
typename Number>
1139 DEAL_II_ALWAYS_INLINE
inline auto
1143 const dealii::Tensor<1, dim, Number> &c_ij)
const ->
state_type
1145 const auto &[U_i, Z_i] = flux_i;
1146 const auto &[U_j, Z_j] = flux_j;
1148 const auto H_i = water_depth(U_i);
1149 const auto H_j = water_depth(U_j);
1151 const auto g_i = g(U_i);
1152 const auto g_j = g(U_j);
1154 auto result = -
add(g_i, g_j);
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;
1165 template <
int dim,
typename Number>
1166 DEAL_II_ALWAYS_INLINE
inline auto
1170 const dealii::Tensor<1, dim, Number> &c_ij,
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);
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));
1181 return -factor * (U_star_ij - U_i);
1185 template <
int dim,
typename Number>
1186 DEAL_II_ALWAYS_INLINE
inline auto
1193 const auto g = gravity();
1194 const auto n = manning_friction_coefficient();
1196 const auto h_inverse = inverse_water_depth_mollified(U);
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;
1202 const auto denominator = h_star + std::max(h_star, tau * factor);
1203 const auto denominator_inverse =
ScalarNumber(1.) / denominator;
1205 for (
unsigned int d = 0; d < dim; ++d)
1206 result[d + 1] = -factor * denominator_inverse * m[d];
1212 template <
int dim,
typename Number>
1213 DEAL_II_ALWAYS_INLINE
inline auto
1216 const unsigned int i,
1220 const auto &[eta_m, h_star] =
1221 pv.template get_tensor<Number, precomputed_type>(i);
1223 return manning_friction(U_i, h_star, tau);
1227 template <
int dim,
typename Number>
1228 DEAL_II_ALWAYS_INLINE
inline auto
1231 const unsigned int *js,
1235 const auto &[eta_m, h_star] =
1236 pv.template get_tensor<Number, precomputed_type>(js);
1238 return manning_friction(U_j, h_star, tau);
1242 template <
int dim,
typename Number>
1243 template <
typename ST>
1244 DEAL_II_ALWAYS_INLINE
inline auto
1248 using T =
typename ST::value_type;
1249 static_assert(std::is_same_v<Number, T>,
"template mismatch");
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");
1257 result[0] = state[0];
1258 for (
unsigned int i = 1; i < dim2 + 1; ++i)
1259 result[i] = state[i];
1264 template <
int dim,
typename Number>
1265 template <
typename ST>
1266 DEAL_II_ALWAYS_INLINE
inline auto
1270 const auto primitive_state = expand_state(initial_state);
1271 return from_primitive_state(primitive_state);
1275 template <
int dim,
typename Number>
1276 DEAL_II_ALWAYS_INLINE
inline auto
1280 const auto &h = primitive_state[0];
1282 auto state = primitive_state;
1284 for (
unsigned int i = 1; i < dim + 1; ++i)
1291 template <
int dim,
typename Number>
1292 DEAL_II_ALWAYS_INLINE
inline auto
1296 const auto h_inverse = inverse_water_depth_sharp(state);
1298 auto primitive_state = state;
1300 for (
unsigned int i = 1; i < dim + 1; ++i)
1301 primitive_state[i] *= h_inverse;
1303 return primitive_state;
1307 template <
int dim,
typename Number>
1308 template <
typename Lambda>
1309 DEAL_II_ALWAYS_INLINE
inline auto
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];
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)