ryujin 2.1.1 revision 075c649df0ea80b3fbb39a7a6e10bf95d4e28390
initial_state_flow_over_bump.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0
3// [LANL Copyright Statement]
4// Copyright (C) 2022 - 2024 by the ryujin authors
5// Copyright (C) 2023 - 2024 by Triad National Security, LLC
6//
7
8#pragma once
9
11
12#include <simd.h>
13
14namespace ryujin
15{
16 namespace ShallowWaterInitialStates
17 {
24 template <typename Description, int dim, typename Number>
25 class FlowOverBump : public InitialState<Description, dim, Number>
26 {
27 public:
29 using View =
30 typename Description::template HyperbolicSystemView<dim, Number>;
31 using state_type = typename View::state_type;
32
33 FlowOverBump(const HyperbolicSystem &hyperbolic_system,
34 const std::string subsection)
35 : InitialState<Description, dim, Number>("flow over bump", subsection)
36 , hyperbolic_system_(hyperbolic_system)
37 {
38 dealii::ParameterAcceptor::parse_parameters_call_back.connect(
40
41 which_case_ = "transcritical";
42 this->add_parameter("flow type",
43 which_case_,
44 "Either 'transcritical' flow with shock "
45 "or 'subsonic' flow.");
46 }
47
49 {
50 AssertThrow(which_case_ == "subsonic" || which_case_ == "transcritical",
51 dealii::ExcMessage("Flow type must be 'transcritical' "
52 "or 'subsonic'. "));
53 }
54
55 state_type compute(const dealii::Point<dim> &point, Number t) final
56 {
57 const auto view = hyperbolic_system_.template view<dim, Number>();
58 const auto g = view.gravity();
59
60 const auto x = point[0];
61
62 /* Define constants for transcritical flow */
63 const Number xM = 10.;
64 const Number xS = 11.7;
65 const Number zM = 0.2;
66 Number h_inflow = 0.28205279813802181;
67 Number q_inflow = 0.18;
68 Number cBer =
69 zM + 1.5 * ryujin::pow(q_inflow * q_inflow / g, Number(1. / 3.));
70
71
72 /* Subsonic flow constants */
73 if (which_case_ == "subsonic") {
74 q_inflow = 4.42;
75 h_inflow = 2.;
76 cBer = std::pow(q_inflow / h_inflow, 2) / (2. * g) + h_inflow;
77 }
78
79 /* General values for Cardano's formula */
80 const Number d = q_inflow * q_inflow / (2. * g);
81 const Number b = compute_bathymetry(point) - cBer;
82 const Number Q = -std::pow(b, 2) / 9.;
83 const Number R = -(27. * d + 2. * std::pow(b, 3)) / 54.;
84 const Number theta = acos(ryujin::pow(-Q, Number(-1.5)) * R);
85
86 /* Define initial and exact solution */
87 const Number h_initial = h_inflow - compute_bathymetry(point);
88
89 /* Explicitly return initial state */
90 if (t < 1e-12) {
91 return state_type{{h_initial, q_inflow}};
92 }
93
94 Number h_exact = 2. * std::sqrt(-Q) * cos(theta / 3.) - b / 3.;
95 if (which_case_ == "transcritical") {
96 if (xM <= x && x < xS) {
97 h_exact = 2. * std::sqrt(-Q) *
98 cos((4. * dealii::numbers::PI + theta) / 3.) -
99 b / 3.;
100 } else if (xS < x) {
101 h_exact = 0.28205279813802181;
102 }
103 }
104
105 return state_type{{h_exact, q_inflow}};
106 }
107
108 auto initial_precomputations(const dealii::Point<dim> &point) ->
110 initial_precomputed_type final
111 {
112 /* Compute bathymetry: */
113 return {compute_bathymetry(point)};
114 }
115
116 private:
117 const HyperbolicSystem &hyperbolic_system_;
118
119 DEAL_II_ALWAYS_INLINE
120 inline Number compute_bathymetry(const dealii::Point<dim> &point) const
121 {
122 const auto x = point[0];
123
124 const Number bump = Number(0.2 / 64.) * std::pow(x - Number(8.), 3) *
125 std::pow(Number(12.) - x, 3);
126
127 Number bath = 0.;
128 if (8. <= x && x <= 12.)
129 bath = bump;
130 return bath;
131 }
132
133 std::string which_case_;
134 };
135
136 } // namespace ShallowWaterInitialStates
137} // namespace ryujin
typename Description::HyperbolicSystem HyperbolicSystem
state_type compute(const dealii::Point< dim > &point, Number t) final
typename Description::template HyperbolicSystemView< dim, Number > View
auto initial_precomputations(const dealii::Point< dim > &point) -> typename InitialState< Description, dim, Number >::initial_precomputed_type final
FlowOverBump(const HyperbolicSystem &hyperbolic_system, const std::string subsection)
T pow(const T x, const T b)
Euler::HyperbolicSystem HyperbolicSystem
Definition: description.h:32