9#include <initial_state.h>
13 namespace ShallowWater
21 template <
int dim,
typename Number,
typename state_type>
26 const std::string subsec)
28 , hyperbolic_system(hyperbolic_system)
30 dealii::ParameterAcceptor::parse_parameters_call_back.connect(
33 which_case_ =
"transcritical";
34 this->add_parameter(
"flow type",
36 "Either 'transcritical' flow with shock "
37 "or 'subsonic' flow.");
42 AssertThrow(which_case_ ==
"subsonic" || which_case_ ==
"transcritical",
43 dealii::ExcMessage(
"Flow type must be 'transcritical' "
49 const auto x = point[0];
50 const Number g = this->hyperbolic_system.gravity();
53 const Number xM = 10.;
54 const Number xS = 11.7;
55 const Number zM = 0.2;
56 Number h_inflow = 0.28205279813802181;
57 Number q_inflow = 0.18;
59 zM + 1.5 *
ryujin::pow(q_inflow * q_inflow / g, Number(1. / 3.));
63 if (which_case_ ==
"subsonic") {
66 cBer =
std::pow(q_inflow / h_inflow, 2) / (2. * g) + h_inflow;
70 const Number d = q_inflow * q_inflow / (2. * g);
71 const Number b = compute_bathymetry(point) - cBer;
72 const Number Q = -
std::pow(b, 2) / 9.;
73 const Number R = -(27. * d + 2. *
std::pow(b, 3)) / 54.;
74 const Number theta = acos(
ryujin::pow(-Q, Number(-1.5)) * R);
77 const Number h_initial = h_inflow - compute_bathymetry(point);
81 return hyperbolic_system.template expand_state<dim>(
84 Number h_exact = 2. * std::sqrt(-Q) * cos(theta / 3.) - b / 3.;
85 if (which_case_ ==
"transcritical") {
86 if (xM <= x && x < xS) {
87 h_exact = 2. * std::sqrt(-Q) *
88 cos((4. * dealii::numbers::PI + theta) / 3.) -
91 h_exact = 0.28205279813802181;
95 return hyperbolic_system.template expand_state<dim>(
104 return {compute_bathymetry(point)};
110 DEAL_II_ALWAYS_INLINE
111 inline Number compute_bathymetry(
const dealii::Point<dim> &point)
const
113 const auto x = point[0];
115 const Number bump = Number(0.2 / 64.) *
std::pow(x - Number(8.), 3) *
119 if (8. <= x && x <= 12.)
124 std::string which_case_;
typename HyperbolicSystemView::state_type state_type
void parse_parameters_callback()
FlowOverBump(const HyperbolicSystem &hyperbolic_system, const std::string subsec)
auto initial_precomputations(const dealii::Point< dim > &point) -> typename InitialState< dim, Number, state_type, 1 >::precomputed_type final
state_type compute(const dealii::Point< dim > &point, Number t) final
dealii::Tensor< 1, problem_dimension< dim >, Number > state_type
T pow(const T x, const T b)