problem.hh 11 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
// -*- 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
 * \ingroup ShallowWaterTests
 * \brief A test for the Shallow water model (bowl).
 */
#ifndef DUMUX_BOWL_TEST_PROBLEM_HH
#define DUMUX_BOWL_TEST_PROBLEM_HH

#include <dumux/freeflow/shallowwater/problem.hh>
#include <dumux/freeflow/shallowwater/boundaryfluxes.hh>

namespace Dumux {

/*!
 * \ingroup ShallowWaterTests
 * \brief A wetting and drying test with sloshing water in a bowl.
 *
 * The domain is 4 meters long and 4 meters wide. The center of the domain is loacted at
 * x = 0 and y = 0. There is no flow over the boundaries and no friction is considered.
 *
 * This example is demanding for the implicit model if a high mesh resolution is applied
 * (e.g. 150x150 cells) in combination with a large time step size. Using the new limiting
 * (UpwindFluxLimiting = true) will help to improve the convergence for such cases.
 *
 * This test uses a low mesh resoultion and only ensures that UpwindFluxLimiting for the mobility
44
45
 * works. For low mesh resolution the solution is very diffusive so that the oscillation is dampened.
 * This gets better with grid refinement (not tested here).
46
47
48
 *
 * The results are checked against a analytical solution which is based on the "Thacker-Solution"
 * (William Thacker, "Some exact solutions to the nonlinear shallow-water wave equations", Journal
49
50
51
 * of Fluid Mechanics, 107:499–508, 1981, doi: https://doi.org/10.1017/S0022112081001882).
 * This implements the oscillating solution in a circular bowl (Section 4 in the paper).
 * Further examples and details on the solution are given
52
53
54
55
56
57
58
59
60
 * in SWASHES (Shallow Water Analytic Solutions for Hydraulic and Environmental Studies,
 * https://www.idpoisson.fr/swashes/).
 *
 * This problem uses the \ref ShallowWaterModel
 */
template <class TypeTag>
class BowlProblem : public ShallowWaterProblem<TypeTag>
{
    using ParentType = ShallowWaterProblem<TypeTag>;
61

62
63
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
64
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
65
66
67
68
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
    using GridView = GetPropType<TypeTag, Properties::GridView>;
    using Element = typename GridView::template Codim<0>::Entity;
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
69
70

    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
71
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
72
73
74
75
76
77
78
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;

    using NeumannFluxes = GetPropType<TypeTag, Properties::NumEqVector>;
    using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;

79
80
81
82
83
84
85
86
87
88

public:
    BowlProblem(std::shared_ptr<const GridGeometry> gridGeometry)
    : ParentType(gridGeometry)
    {

        name_ = getParam<std::string>("Problem.Name");
        bowlDepthAtCenter_ =  getParam<Scalar>("Problem.BowlDepthAtCenter");
        bowlParaboloidRadius_ =  getParam<Scalar>("Problem.BowlParaboloidRadius");

89
        // Thacker (1981) Eq. (43)
90
        using std::sqrt;
91
92
93
94
95
96
97
98
99
100
101
102
103
        bowlAnalyticParameterOmega_ = sqrt(8.0 * getParam<Scalar>("Problem.Gravity") * bowlDepthAtCenter_) / bowlParaboloidRadius_;
        std::cout << "One full oscillation period of the water table is: "
                  << oscillationPeriodInSeconds() << " seconds." << std::endl;

        // Thacker (1981) Eq. (50)
        const auto D0PlusEta = bowlDepthAtCenter_ + getParam<Scalar>("Problem.BowlInitialWaterElevationAtCenter");
        const auto D0PlusEtaSquared = D0PlusEta*D0PlusEta;
        const auto D0Squared = bowlDepthAtCenter_*bowlDepthAtCenter_;
        bowlAnalyticParameterA_ = (D0PlusEtaSquared - D0Squared)/(D0PlusEtaSquared + D0Squared);

        // check constraint Thacker (1981) text after Eq. (49)
        if (bowlAnalyticParameterA_ >= 1.0)
            DUNE_THROW(Dune::InvalidStateException, "Parameter A has to be smaller than unity!");
104
105
    }

106
107
108
    //! One oscillation period of the water table (analytically this goes on forever)
    Scalar oscillationPeriodInSeconds() const
    { return 2*M_PI/bowlAnalyticParameterOmega_; }
109

110
111
    //! Get the analytical water depth at time t and position pos
    PrimaryVariables analyticalSolution(const Scalar t, const GlobalPosition& pos) const
112
    {
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
        using std::sqrt;
        using std::cos;
        using std::sin;

        // see Thacker (1981) Eq. (51) for formula
        const auto radiusSquared = pos[0]*pos[0] + pos[1]*pos[1];
        const auto LSquared = bowlParaboloidRadius_*bowlParaboloidRadius_;
        const auto A = bowlAnalyticParameterA_;
        const auto omega = bowlAnalyticParameterOmega_;
        const auto D0 = bowlDepthAtCenter_;

        const auto oneMinusASq = 1.0 - A*A;
        const auto oneMinusACosOmegaT = 1.0 - A*cos(omega*t);
        const auto ratioSq = oneMinusASq / (oneMinusACosOmegaT*oneMinusACosOmegaT);
        const auto localRadiusSq = radiusSquared / LSquared;
        // bowl depth function (cf. D in Thacker (1981))
        const auto D = D0*(1.0 - localRadiusSq);
        // height above equilibrium water level (cf. h in Thacker (1981))
        const auto h = D0*(sqrt(ratioSq) - 1.0 - localRadiusSq*(ratioSq - 1.0));
        // see remark about total water depth in Thacker (1981) beginning section 2
        const auto analyticalWaterDepth = h + D;

        const auto halfOmegaASinOmegaT = 0.5*omega*A*sin(omega*t);
        const auto analyticalVelocityX = pos[0]*halfOmegaASinOmegaT/oneMinusACosOmegaT;
        const auto analyticalVelocityY = pos[1]*halfOmegaASinOmegaT/oneMinusACosOmegaT;

        // The radius of the shoreline (where h + D = 0), Eq. (48)
        const auto h0 = D0*(sqrt(ratioSq) - 1.0); // h in the middle of the bowl (r=0)
        const auto shoreLineRadiusSquared = LSquared*(D0/(D0 + h0));

        // outside shoreline the water height and velocity is zero
        if (radiusSquared > shoreLineRadiusSquared)
            return { 0.0, 0.0, 0.0 };
        else
            return { analyticalWaterDepth, analyticalVelocityX, analyticalVelocityY };
148
149
150
151
152
153
154
155
156
157
158
159
    }

    /*!
     * \name Problem parameters
     */
    // \{

    /*!
     * \brief The problem name
     * This is used as a prefix for files generated by the simulation.
     */
    const std::string& name() const
160
    { return name_; }
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186

    // \}

    /*!
     * \name Boundary conditions
     */
    // \{

    /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary segment.
     */
    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
    {
        BoundaryTypes bcTypes;
        bcTypes.setAllNeumann();
        return bcTypes;
    }

    /*!
     * \brief Specifies the neumann boundary
     *
     *  We need the Riemann invariants to compute the values depending of the boundary type.
     *  Since we use a weak imposition we do not have a dirichlet value. We impose fluxes
     *  based on q, h, etc. computed with the Riemann invariants
     */
187
    template<class ElementFluxVariablesCache>
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
    NeumannFluxes neumann(const Element& element,
                          const FVElementGeometry& fvGeometry,
                          const ElementVolumeVariables& elemVolVars,
                          const ElementFluxVariablesCache& elemFluxVarsCache,
                          const SubControlVolumeFace& scvf) const
    {
        NeumannFluxes values(0.0);

        const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
        const auto& insideVolVars = elemVolVars[insideScv];
        const auto& nxy = scvf.unitOuterNormal();
        const auto gravity = this->spatialParams().gravity(scvf.center());
        std::array<Scalar, 3> boundaryStateVariables;

        //no flow with zero normal velocity and tangential velocity
        const auto vNormalGhost = -(nxy[0] * insideVolVars.velocity(0) +  nxy[1] * insideVolVars.velocity(1));
        const auto vTangentialGhost = -nxy[1] * insideVolVars.velocity(0) + nxy[0] * insideVolVars.velocity(1);

        boundaryStateVariables[0] = insideVolVars.waterDepth();
        boundaryStateVariables[1] =  nxy[0] * vNormalGhost - nxy[1] * vTangentialGhost;
        boundaryStateVariables[2] =  nxy[1] * vNormalGhost + nxy[0] * vTangentialGhost;

210
211
212
213
214
215
        const auto riemannFlux =
            ShallowWater::riemannProblem(insideVolVars.waterDepth(), boundaryStateVariables[0],
                                         insideVolVars.velocity(0), boundaryStateVariables[1],
                                         insideVolVars.velocity(1), boundaryStateVariables[2],
                                         insideVolVars.bedSurface(), insideVolVars.bedSurface(),
                                         gravity, nxy);
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231

        values[Indices::massBalanceIdx] = riemannFlux[0];
        values[Indices::velocityXIdx]   = riemannFlux[1];
        values[Indices::velocityYIdx]   = riemannFlux[2];

        return values;
    }

    // \}

    /*!
     * \name Volume terms
     */
    // \{

    /*!
232
     * \brief Evaluate the initial values at position globalPos
233
     */
234
    PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
235
    {
236
237
        using std::max; // regularize so that we are virtually dry but not completely dry
        return { max(analyticalSolution(0, globalPos)[0], 1e-5), 0.0, 0.0 };
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
    };

    // \}

private:
    Scalar bowlDepthAtCenter_;
    Scalar bowlParaboloidRadius_;
    Scalar bowlAnalyticParameterOmega_;
    Scalar bowlAnalyticParameterA_;
    static constexpr Scalar eps_ = 1.0e-6;
    std::string name_;
};

} //end namespace Dumux

#endif