riemannproblem.hh 7.08 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
21
22
 * \ingroup ShallowWaterFlux
 * \brief This file includes a function to construct the Riemann problem
23
24
25
26
27
28
29
30
31
32
 *
 */
#ifndef DUMUX_FLUX_SHALLOW_WATER_RIEMANN_PROBLEM_HH
#define DUMUX_FLUX_SHALLOW_WATER_RIEMANN_PROBLEM_HH

#include <dumux/flux/shallowwater/fluxlimiterlet.hh>
#include <dumux/flux/shallowwater/exactriemann.hh>

namespace Dumux {
namespace ShallowWater {
33

34
/*!
35
 * \ingroup ShallowWaterFlux
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
 * \brief Construct Riemann Problem and solve it
 *
 *
 * Riemann Problem applies the hydrostatic reconstruction, uses the
 * Riemann invariants to transform the two-dimensional problem to an
 * one-dimensional problem and solves this new problem, and rotates
 * the problem back. Further it applies an flux limiter for the water
 * flux handle drying of elements.
 * The correction of the bed slope surce terme leads to an
 * non-symetric flux term at the interface for the momentum equations
 * since DuMuX computes the fluxes twice from each side this does not
 * matter.
 *
 * So far we have only the exact Riemann solver, and the reconstruction
 * after Audusse but further solvers and reconstructions ca be
 * implemented.
 *
 * \param waterDepthLeft water depth on the left side
 * \param waterDepthRight water depth on the right side
 * \param velocityXLeft veloctiyX on the left side
 * \param velocityXRight velocityX on the right side
 * \param velocityYLeft velocityY on the left side
 * \param velocityYRight velocityY on the right side
 * \param bedSurfaceLeft surface of the bed on the left side
 * \param bedSurfaceRight surface of the bed on the right side
61
 * \param gravity gravity constant
62
 * \param nxy the normal vector
63
64
65
 *
 */
template<class Scalar, class GlobalPosition>
66
67
std::array<Scalar,3> riemannProblem(const Scalar waterDepthLeft,
                                    const Scalar waterDepthRight,
68
69
70
71
                                    Scalar velocityXLeft,
                                    Scalar velocityXRight,
                                    Scalar velocityYLeft,
                                    Scalar velocityYRight,
72
73
74
75
                                    const Scalar bedSurfaceLeft,
                                    const Scalar bedSurfaceRight,
                                    const Scalar gravity,
                                    const GlobalPosition& nxy)
76
77
78
{
    using std::max;

79
80
81
82
83
    // hydrostatic reconstrucion after Audusse
    const Scalar dzl = max(0.0, bedSurfaceRight - bedSurfaceLeft);
    const Scalar waterDepthLeftReconstructed = max(0.0, waterDepthLeft - dzl);
    const Scalar dzr = max(0.0, bedSurfaceLeft - bedSurfaceRight);
    const Scalar waterDepthRightReconstructed = max(0.0, waterDepthRight - dzr);
84

85
    // make rotation of the flux we compute an 1d flux
86
87
88
89
90
91
92
93
    Scalar tempFlux = velocityXLeft;
    velocityXLeft =  nxy[0] * tempFlux + nxy[1] * velocityYLeft;
    velocityYLeft = -nxy[1] * tempFlux + nxy[0] * velocityYLeft;

    tempFlux = velocityXRight;
    velocityXRight =  nxy[0] * tempFlux + nxy[1] * velocityYRight;
    velocityYRight = -nxy[1] * tempFlux + nxy[0] * velocityYRight;

94
95
96
97
98
99
100
    auto riemannResult = ShallowWater::exactRiemann(waterDepthLeftReconstructed,
                                                    waterDepthRightReconstructed,
                                                    velocityXLeft,
                                                    velocityXRight,
                                                    velocityYLeft,
                                                    velocityYRight,
                                                    gravity);
101
102

    //redo rotation
103
104
105
    tempFlux = riemannResult.flux[1];
    riemannResult.flux[1] = nxy[0] * tempFlux - nxy[1] * riemannResult.flux[2];
    riemannResult.flux[2] = nxy[1] * tempFlux + nxy[0] * riemannResult.flux[2];
106

107
108
109
110
    // Add reconstruction flux from Audusse reconstruction
    const Scalar hgzl = 0.5 * (waterDepthLeftReconstructed + waterDepthLeft) * (waterDepthLeftReconstructed - waterDepthLeft);
    const Scalar hdxzl = gravity * nxy[0] * hgzl;
    const Scalar hdyzl = gravity * nxy[1] * hgzl;
111
112
113
114
115
116
117
118
119

    /*Right side is computed from the other side otherwise the
      following "non-symetric" fluxes are needed:

      Scalar hgzr = 0.5 * (waterDepthRightReconstructed + waterDepthRight) * (waterDepthRightReconstructed - waterDepthRight);
      Scalar hdxzr = gravity * nxy[0] * hgzr;
      Scalar hdyzrhdyzr = gravity * nxy[1] * hgzr;
    */

120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
    // compute the mobility of the flux with the fluxlimiter
    static const Scalar upperWaterDepthFluxLimiting = getParam<Scalar>("FluxLimiterLET.UpperWaterDepth", 1e-3);
    static const Scalar lowerWaterDepthFluxLimiting = getParam<Scalar>("FluxLimiterLET.LowerWaterDepth", 1e-5);
    static const bool upwindWaterDepthFluxLimiting = getParam<bool>("FluxLimiterLET.UpwindFluxLimiting", false);

    Scalar limitingDepth = (waterDepthLeftReconstructed + waterDepthRightReconstructed) * 0.5;

    //Using the upwind water depth from the flux direction can improve stability.
    if (upwindWaterDepthFluxLimiting)
    {
        if (riemannResult.flux[0] < 0)
        {
            limitingDepth = waterDepthRightReconstructed;
        }else
        {
            limitingDepth = waterDepthLeftReconstructed;
        }
    }

    const Scalar mobility = ShallowWater::fluxLimiterLET(limitingDepth,
                                                         limitingDepth,
                                                         upperWaterDepthFluxLimiting,
                                                         lowerWaterDepthFluxLimiting);
143
    std::array<Scalar, 3> localFlux;
144
    localFlux[0] = riemannResult.flux[0] * mobility;
145
146
    localFlux[1] = (riemannResult.flux[1] - hdxzl);// * mobility;
    localFlux[2] = (riemannResult.flux[2] - hdyzl);// * mobility;
147
148
149

    return localFlux;
}
150

151
152
153
154
} // end namespace ShallowWater
} // end namespace Dumux

#endif