Merge branch 'feature/upwind-mobility-for-shallow-water-flux' into 'master'

Resolve "[shallowwater] upwind mobility for shallow water flux."

Closes #800

See merge request !1831

dumux/flux/shallowwater/riemannproblem.hh

@@ -82,14 +82,6 @@ std::array<Scalar,3> riemannProblem(const Scalar waterDepthLeft,
     const Scalar dzr = max(0.0, bedSurfaceLeft - bedSurfaceRight);
     const Scalar waterDepthRightReconstructed = max(0.0, waterDepthRight - dzr);
 
-    // 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);
-    const Scalar mobility = ShallowWater::fluxLimiterLET(waterDepthLeftReconstructed,
-                                                         waterDepthRightReconstructed,
-                                                         upperWaterDepthFluxLimiting,
-                                                         lowerWaterDepthFluxLimiting);
-
     // make rotation of the flux we compute an 1d flux
     Scalar tempFlux = velocityXLeft;
     velocityXLeft =  nxy[0] * tempFlux + nxy[1] * velocityYLeft;
@@ -125,10 +117,33 @@ std::array<Scalar,3> riemannProblem(const Scalar waterDepthLeft,
       Scalar hdyzrhdyzr = gravity * nxy[1] * hgzr;
     */
 
+    // 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);
     std::array<Scalar, 3> localFlux;
     localFlux[0] = riemannResult.flux[0] * mobility;
-    localFlux[1] = riemannResult.flux[1] - hdxzl;
-    localFlux[2] = riemannResult.flux[2] - hdyzl;
+    localFlux[1] = (riemannResult.flux[1] - hdxzl);
+    localFlux[2] = (riemannResult.flux[2] - hdyzl);
 
     return localFlux;
 }

test/freeflow/shallowwater/CMakeLists.txt

+add_subdirectory(bowl)
 add_subdirectory(dambreak)
 add_subdirectory(roughchannel)

test/freeflow/shallowwater/bowl/CMakeLists.txt

+dune_symlink_to_source_files(FILES "params.input")
+
+dumux_add_test(NAME test_shallowwater_bowl
+               SOURCES main.cc
+               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+               CMD_ARGS   --script fuzzy
+                          --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_shallowwater_bowl-reference.vtu
+                                  ${CMAKE_CURRENT_BINARY_DIR}/bowl-00013.vtu
+                          --zeroThreshold {"velocityY":1e-14}
+                          --command "${CMAKE_CURRENT_BINARY_DIR}/test_shallowwater_bowl")
+
+dumux_add_test(NAME test_shallowwater_bowl_parallel
+               TARGET test_shallowwater_bowl
+               CMAKE_GUARD MPI_FOUND
+               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+               CMD_ARGS   --script fuzzy
+                          --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_shallowwater_bowl-reference.vtu
+                                  ${CMAKE_CURRENT_BINARY_DIR}/s0002-bowl-parallel-00013.pvtu
+                          --zeroThreshold {"velocityY":1e-14,"process rank":100}
+                          --command "${MPIEXEC} -np 2 ${CMAKE_CURRENT_BINARY_DIR}/test_shallowwater_bowl -Problem.Name bowl-parallel")

test/freeflow/shallowwater/bowl/main.cc

+// -*- 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).
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+
+#include <dumux/io/grid/gridmanager.hh>
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+
+#include "properties.hh"
+
+template<class SolutionVector, int i>
+class SolutionComponent
+{
+    const SolutionVector& sol_;
+public:
+    SolutionComponent(const SolutionVector& sol) : sol_(sol) {}
+    auto operator[] (std::size_t idx) const { return sol_[idx][i]; }
+    auto size() const { return sol_.size(); }
+};
+
+//! Compute the analytical solution at time t
+template<class SolutionVector, class Problem>
+SolutionVector computeAnalyticalSolution(const double t, const Problem& problem)
+{
+    const auto& gg = problem.gridGeometry();
+    SolutionVector exactSolution(gg.numDofs());
+    for (const auto& element : elements(gg.gridView()))
+    {
+        const auto eIdx = gg.elementMapper().index(element);
+        const auto& globalPos = element.geometry().center();
+        exactSolution[eIdx] = problem.analyticalSolution(t, globalPos);
+    }
+    return exactSolution;
+}
+
+//! Compute L2 error wrt the analytical solution at time t
+template<class SolutionVector, class GridGeometry>
+typename SolutionVector::block_type
+computeL2Error(const double t,
+               const SolutionVector& exactSolution,
+               const SolutionVector& curSol,
+               const GridGeometry& gridGeometry)
+{
+    typename SolutionVector::block_type l2Error(0.0);
+    for (const auto& element : elements(gridGeometry.gridView(), Dune::Partitions::interior))
+    {
+        auto fvGeometry = localView(gridGeometry);
+        fvGeometry.bindElement(element);
+
+        using std::pow;
+        for (auto&& scv : scvs(fvGeometry))
+        {
+            auto localDiffSq = exactSolution[scv.dofIndex()] - curSol[scv.dofIndex()];
+            for (int i = 0; i < localDiffSq.size(); ++i)
+                localDiffSq[i] *= localDiffSq[i];
+            l2Error += localDiffSq;
+        }
+    }
+    // sum over processes if we are running with multiple processes in parallel
+    if (gridGeometry.gridView().comm().size() > 0)
+        l2Error = gridGeometry.gridView().comm().sum(l2Error);
+
+    // square root of sum of squared errors is our absolute discrete l2 error
+    for (int i = 0; i < l2Error.size(); ++i)
+        l2Error[i] = std::sqrt(l2Error[i]);
+    return l2Error;
+}
+
+//! Compute and print L2 error wrt the analytical solution at time t
+template<class SolutionVector, class GridGeometry>
+void computeAndPrintL2Error(const double t,
+                            const SolutionVector& exactSolution,
+                            const SolutionVector& curSol,
+                            const GridGeometry& gridGeometry)
+{
+    const auto l2Error = computeL2Error(t, exactSolution, curSol, gridGeometry);
+    auto numElements = gridGeometry.gridView().size(0);
+    // sum over processes if we are running with multiple processes in parallel
+    if (gridGeometry.gridView().comm().size() > 0)
+        numElements = gridGeometry.gridView().comm().sum(numElements);
+
+    if (gridGeometry.gridView().comm().rank() == 0)
+        std::cout << "L2 error in (h, u, v) at t = " <<  t << " seconds for " << numElements << " elements: "
+                  << std::scientific << l2Error << std::endl;
+}
+
+////////////////////////
+// the main function
+////////////////////////
+int main(int argc, char** argv) try
+{
+    using namespace Dumux;
+
+    // define the type tag for this problem
+    using TypeTag = Properties::TTag::Bowl;
+
+    // initialize MPI, finalize is done automatically on exit
+    const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+    // print dumux start message
+    if (mpiHelper.rank() == 0)
+        DumuxMessage::print(/*firstCall=*/true);
+
+    // parse command line arguments and input file
+    Parameters::init(argc, argv);
+
+    // try to create a grid (from the given grid file or the input file)
+    GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
+    gridManager.init();
+
+    ////////////////////////////////////////////////////////////
+    // run instationary non-linear problem on this grid
+    ////////////////////////////////////////////////////////////
+
+    // we compute on the leaf grid view
+    const auto& leafGridView = gridManager.grid().leafGridView();
+
+    // create the finite volume grid geometry
+    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+    auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
+    gridGeometry->update();
+
+    // the problem (initial and boundary conditions)
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    auto problem = std::make_shared<Problem>(gridGeometry);
+
+    // the solution vector
+    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+    SolutionVector x;
+    problem->applyInitialSolution(x);
+    auto xOld = x;
+
+    // the grid variables
+    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+    auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
+    gridVariables->init(x);
+
+    // get some time loop parameters
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    const auto tEnd = 2.0*problem->oscillationPeriodInSeconds();
+    const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+    auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+    // intialize the vtk output module and analytical solution
+    VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables,x, problem->name());
+    auto exactSolution = computeAnalyticalSolution<SolutionVector>(0.0, *problem);
+    const auto exactWaterDepth = SolutionComponent<SolutionVector, 0>(exactSolution);
+    const auto exactVelocityX = SolutionComponent<SolutionVector, 1>(exactSolution);
+    const auto exactVelocityY = SolutionComponent<SolutionVector, 2>(exactSolution);
+    vtkWriter.addField(exactWaterDepth, "exactWaterDepth");
+    vtkWriter.addField(exactVelocityX, "exactVelocityX");
+    vtkWriter.addField(exactVelocityY, "exactVelocityY");
+    using IOFields = GetPropType<TypeTag, Properties::IOFields>;
+    IOFields::initOutputModule(vtkWriter);
+    vtkWriter.write(0.0);
+
+    // instantiate time loop
+    auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, dt, tEnd);
+    timeLoop->setMaxTimeStepSize(maxDt);
+
+    // the assembler with time loop for instationary problem
+    using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+    auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
+
+    // the linear solver
+    using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>>;
+    auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
+
+    // the non-linear solver
+    using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+    NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+    //! Compute L2 error for the initial solution (should be zero because initial solution is exact)
+    computeAndPrintL2Error(0.0, exactSolution, x, *gridGeometry);
+
+    // time loop
+    timeLoop->start(); do
+    {
+        nonLinearSolver.solve(x,*timeLoop);
+
+        // make the new solution the old solution
+        xOld = x;
+        gridVariables->advanceTimeStep();
+
+        // update the analytical solution and print current l2 error
+        exactSolution = computeAnalyticalSolution<SolutionVector>(timeLoop->time(), *problem);
+        computeAndPrintL2Error(timeLoop->time(), exactSolution, x, *gridGeometry);
+
+        // advance to the time loop to the next step
+        timeLoop->advanceTimeStep();
+
+        // write vtk output every 10 time steps
+        if (!(timeLoop->timeStepIndex() % 10))
+            vtkWriter.write(timeLoop->time());
+
+        // report statistics of this time step
+        timeLoop->reportTimeStep();
+
+        // set new dt as suggested by newton controller
+        timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+
+    } while (!timeLoop->finished());
+
+    timeLoop->finalize(leafGridView.comm());
+
+    ////////////////////////////////////////////////////////////
+    // finalize, print dumux message to say goodbye
+    ////////////////////////////////////////////////////////////
+
+    // print dumux end message
+    if (mpiHelper.rank() == 0)
+    {
+        Parameters::print();
+        DumuxMessage::print(/*firstCall=*/false);
+    }
+
+    return 0;
+}
+// Error handling
+catch (const Dumux::ParameterException &e)
+{
+    std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+    return 1;
+}
+catch (const Dune::Exception &e)
+{
+    std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+    return 2;
+}

test/freeflow/shallowwater/bowl/params.input

+[Problem]
+Name = bowl
+BowlDepthAtCenter = 0.03 # [m]
+BowlParaboloidRadius = 1.6 # [m]
+BowlInitialWaterElevationAtCenter = 0.005 # [m]
+Gravity = 9.81 # [m/s^2]
+
+[FluxLimiterLET]
+UpperWaterDepth = 1.0E-5
+LowerWaterDepth = 1.0E-7
+UpwindFluxLimiting = true
+
+[TimeLoop]
+MaxTimeStepSize = 0.1 # [s]
+DtInitial = 0.1 # [s]
+
+[Grid]
+Positions0 = -2.0 2.0
+Positions1 = -2.0 2.0
+Cells0 = 40
+Cells1 = 40
+
+[Newton]
+EnableDynamicOutput = false

test/freeflow/shallowwater/bowl/problem.hh

+// -*- 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
+ * 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).
+ *
+ * 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
+ * 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
+ * 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>;
+
+    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
+    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+    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;
+
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+    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>;
+
+
+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");
+
+        // Thacker (1981) Eq. (43)
+        using std::sqrt;
+        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!");
+    }
+
+    //! One oscillation period of the water table (analytically this goes on forever)
+    Scalar oscillationPeriodInSeconds() const
+    { return 2*M_PI/bowlAnalyticParameterOmega_; }
+
+    //! Get the analytical water depth at time t and position pos
+    PrimaryVariables analyticalSolution(const Scalar t, const GlobalPosition& pos) const
+    {
+        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 };
+    }
+
+    /*!
+     * \name Problem parameters
+     */
+    // \{
+
+    /*!
+     * \brief The problem name
+     * This is used as a prefix for files generated by the simulation.
+     */
+    const std::string& name() const
+    { return name_; }
+
+    // \}
+
+    /*!
+     * \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
+     */
+    template<class ElementFluxVariablesCache>
+    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;
+
+        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);
+
+        values[Indices::massBalanceIdx] = riemannFlux[0];
+        values[Indices::velocityXIdx]   = riemannFlux[1];
+        values[Indices::velocityYIdx]   = riemannFlux[2];
+
+        return values;
+    }
+
+    // \}