3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
riemannproblem.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 2 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
25#ifndef DUMUX_FLUX_SHALLOW_WATER_RIEMANN_PROBLEM_HH
26#define DUMUX_FLUX_SHALLOW_WATER_RIEMANN_PROBLEM_HH
27
31
32namespace Dumux {
33namespace ShallowWater {
34
71template<class Scalar, class GlobalPosition>
72std::array<Scalar,3> riemannProblem(const Scalar waterDepthLeft,
73 const Scalar waterDepthRight,
74 Scalar velocityXLeft,
75 Scalar velocityXRight,
76 Scalar velocityYLeft,
77 Scalar velocityYRight,
78 const Scalar bedSurfaceLeft,
79 const Scalar bedSurfaceRight,
80 const Scalar gravity,
81 const GlobalPosition& nxy)
82{
83 using std::max;
84
85 // hydrostatic reconstrucion after Audusse
86 const Scalar dzl = max(0.0, bedSurfaceRight - bedSurfaceLeft);
87 const Scalar waterDepthLeftReconstructed = max(0.0, waterDepthLeft - dzl);
88 const Scalar dzr = max(0.0, bedSurfaceLeft - bedSurfaceRight);
89 const Scalar waterDepthRightReconstructed = max(0.0, waterDepthRight - dzr);
90
91 // make rotation of the flux we compute an 1d flux
92 Scalar tempFlux = velocityXLeft;
93 velocityXLeft = nxy[0] * tempFlux + nxy[1] * velocityYLeft;
94 velocityYLeft = -nxy[1] * tempFlux + nxy[0] * velocityYLeft;
95
96 tempFlux = velocityXRight;
97 velocityXRight = nxy[0] * tempFlux + nxy[1] * velocityYRight;
98 velocityYRight = -nxy[1] * tempFlux + nxy[0] * velocityYRight;
99
100 auto riemannResult = ShallowWater::exactRiemann(waterDepthLeftReconstructed,
101 waterDepthRightReconstructed,
102 velocityXLeft,
103 velocityXRight,
104 velocityYLeft,
105 velocityYRight,
106 gravity);
107
108 //redo rotation
109 tempFlux = riemannResult.flux[1];
110 riemannResult.flux[1] = nxy[0] * tempFlux - nxy[1] * riemannResult.flux[2];
111 riemannResult.flux[2] = nxy[1] * tempFlux + nxy[0] * riemannResult.flux[2];
112
113 // Add reconstruction flux from Audusse reconstruction
114 const Scalar hgzl = 0.5 * (waterDepthLeftReconstructed + waterDepthLeft) * (waterDepthLeftReconstructed - waterDepthLeft);
115 const Scalar hdxzl = gravity * nxy[0] * hgzl;
116 const Scalar hdyzl = gravity * nxy[1] * hgzl;
117
118 /*Right side is computed from the other side otherwise the
119 following "non-symetric" fluxes are needed:
120
121 Scalar hgzr = 0.5 * (waterDepthRightReconstructed + waterDepthRight) * (waterDepthRightReconstructed - waterDepthRight);
122 Scalar hdxzr = gravity * nxy[0] * hgzr;
123 Scalar hdyzrhdyzr = gravity * nxy[1] * hgzr;
124 */
125
126 // compute the mobility of the flux with the fluxlimiter
127 static const Scalar upperWaterDepthFluxLimiting = getParam<Scalar>("FluxLimiterLET.UpperWaterDepth", 1e-3);
128 static const Scalar lowerWaterDepthFluxLimiting = getParam<Scalar>("FluxLimiterLET.LowerWaterDepth", 1e-5);
129 static const bool upwindWaterDepthFluxLimiting = getParam<bool>("FluxLimiterLET.UpwindFluxLimiting", false);
130
131 Scalar limitingDepth = (waterDepthLeftReconstructed + waterDepthRightReconstructed) * 0.5;
132
133 //Using the upwind water depth from the flux direction can improve stability.
134 if (upwindWaterDepthFluxLimiting)
135 {
136 if (riemannResult.flux[0] < 0)
137 {
138 limitingDepth = waterDepthRightReconstructed;
139 }else
140 {
141 limitingDepth = waterDepthLeftReconstructed;
142 }
143 }
144
145 const Scalar mobility = ShallowWater::fluxLimiterLET(limitingDepth,
146 limitingDepth,
147 upperWaterDepthFluxLimiting,
148 lowerWaterDepthFluxLimiting);
149 std::array<Scalar, 3> localFlux;
150 localFlux[0] = riemannResult.flux[0] * mobility;
151 localFlux[1] = (riemannResult.flux[1] - hdxzl);
152 localFlux[2] = (riemannResult.flux[2] - hdyzl);
153
154 return localFlux;
155}
156
157} // end namespace ShallowWater
158} // end namespace Dumux
159
160#endif
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
Function to limit the fluxes.
Function to compute the Riemann flux at the interface.
RiemannSolution< Scalar > exactRiemann(const Scalar dl, const Scalar dr, const Scalar ul, const Scalar ur, const Scalar vl, const Scalar vr, const Scalar grav, const Scalar s=0.0)
Exact Riemann solver for the shallow water equations.
Definition: exactriemann.hh:66
static Scalar fluxLimiterLET(const Scalar valueLeft, const Scalar valueRight, const Scalar upperH, const Scalar lowerH)
Flux limiter function to scale fluxes for small water depths.
Definition: fluxlimiterlet.hh:52
std::array< Scalar, 3 > riemannProblem(const Scalar waterDepthLeft, const Scalar waterDepthRight, Scalar velocityXLeft, Scalar velocityXRight, Scalar velocityYLeft, Scalar velocityYRight, const Scalar bedSurfaceLeft, const Scalar bedSurfaceRight, const Scalar gravity, const GlobalPosition &nxy)
Construct a Riemann problem and solve it.
Definition: riemannproblem.hh:72
Definition: adapt.hh:29
std::string mobility(int phaseIdx) noexcept
I/O name of mobility for multiphase systems.
Definition: name.hh:101