ryujin 2.1.1 revision 9dcb748690310d6a540ebb8b066d1a0834fc7604
initial_state_flow_over_bump.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: MIT
3// Copyright (C) 2020 - 2023 by the ryujin authors
4//
5
6#pragma once
7
8#include "hyperbolic_system.h"
9#include <initial_state.h>
10
11namespace ryujin
12{
13 namespace ShallowWater
14 {
21 template <int dim, typename Number, typename state_type>
22 class FlowOverBump : public InitialState<dim, Number, state_type, 1>
23 {
24 public:
25 FlowOverBump(const HyperbolicSystem &hyperbolic_system,
26 const std::string subsec)
27 : InitialState<dim, Number, state_type, 1>("flow over bump", subsec)
28 , hyperbolic_system(hyperbolic_system)
29 {
30 dealii::ParameterAcceptor::parse_parameters_call_back.connect(
32
33 which_case_ = "transcritical";
34 this->add_parameter("flow type",
35 which_case_,
36 "Either 'transcritical' flow with shock "
37 "or 'subsonic' flow.");
38 }
39
41 {
42 AssertThrow(which_case_ == "subsonic" || which_case_ == "transcritical",
43 dealii::ExcMessage("Flow type must be 'transcritical' "
44 "or 'subsonic'. "));
45 }
46
47 state_type compute(const dealii::Point<dim> &point, Number t) final
48 {
49 const auto x = point[0];
50 const Number g = this->hyperbolic_system.gravity();
51
52 /* Define constants for transcritical flow */
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;
58 Number cBer =
59 zM + 1.5 * ryujin::pow(q_inflow * q_inflow / g, Number(1. / 3.));
60
61
62 /* Subsonic flow constants */
63 if (which_case_ == "subsonic") {
64 q_inflow = 4.42;
65 h_inflow = 2.;
66 cBer = std::pow(q_inflow / h_inflow, 2) / (2. * g) + h_inflow;
67 }
68
69 /* General values for Cardano's formula */
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);
75
76 /* Define initial and exact solution */
77 const Number h_initial = h_inflow - compute_bathymetry(point);
78
79 /* Explicitly return initial state */
80 if (t < 1e-12)
81 return hyperbolic_system.template expand_state<dim>(
82 HyperbolicSystem::state_type<1, Number>{{h_initial, q_inflow}});
83
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.) -
89 b / 3.;
90 } else if (xS < x) {
91 h_exact = 0.28205279813802181;
92 }
93 }
94
95 return hyperbolic_system.template expand_state<dim>(
96 HyperbolicSystem::state_type<1, Number>{{h_exact, q_inflow}});
97 }
98
99 auto initial_precomputations(const dealii::Point<dim> &point) ->
101 final
102 {
103 /* Compute bathymetry: */
104 return {compute_bathymetry(point)};
105 }
106
107 private:
108 const HyperbolicSystem &hyperbolic_system;
109
110 DEAL_II_ALWAYS_INLINE
111 inline Number compute_bathymetry(const dealii::Point<dim> &point) const
112 {
113 const auto x = point[0];
114
115 const Number bump = Number(0.2 / 64.) * std::pow(x - Number(8.), 3) *
116 std::pow(Number(12.) - x, 3);
117
118 Number bath = 0.;
119 if (8. <= x && x <= 12.)
120 bath = bump;
121 return bath;
122 }
123
124 std::string which_case_;
125 };
126
127 } // namespace ShallowWater
128} // namespace ryujin
typename HyperbolicSystemView::state_type state_type
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)