[shallowwater][diffusion] add turbulent diffusion to shallow water

dumux/common/properties.hh

@@ -276,6 +276,8 @@ struct SherwoodFormulation { using type = UndefinedProperty; };
 
 template<class TypeTag, class MyTypeTag>
 struct NormalizePressure { using type = UndefinedProperty; }; //!<  Returns whether to normalize the pressure term in the momentum balance or not
+template<class TypeTag, class MyTypeTag>
+struct ViscousFluxType { using type = UndefinedProperty; };              //!< The type for the calculation of the (turbulent) viscous (momentum) fluxes
 
 /////////////////////////////////////////////////////////////
 // Properties used by multidomain simulations

dumux/flux/CMakeLists.txt

@@ -20,6 +20,7 @@ maxwellstefandiffusioncoefficients.hh
 maxwellstefanslaw.hh
 referencesystemformulation.hh
 shallowwaterflux.hh
+shallowwaterviscousflux.hh
 stationaryvelocityfield.hh
 traits.hh
 upwindscheme.hh

dumux/flux/shallowwaterviscousflux.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 3 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 Flux
+ * \copydoc Dumux::ShallowWaterViscousMomentumFlux
+ */
+#ifndef DUMUX_FLUX_SHALLOW_WATER_VISCOUS_FLUX_HH
+#define DUMUX_FLUX_SHALLOW_WATER_VISCOUS_FLUX_HH
+
+#include <dumux/flux/fluxvariablescaching.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Flux
+ * \brief Computes the shallow water viscous momentum flux due to (turbulent) viscosity
+ *        by adding all surrounding shear stresses.
+ *        For now implemented strictly for 2D depth-averaged models (i.e. 3 equations)
+ */
+template<class PrimaryVariables, class NumEqVector, typename std::enable_if_t<NumEqVector::size() == 3, int> = 0>
+class ShallowWaterViscousFlux
+{
+
+public:
+
+    using Cache = FluxVariablesCaching::EmptyDiffusionCache;
+    using CacheFiller = FluxVariablesCaching::EmptyCacheFiller;
+    /*!
+     * \ingroup Flux
+     * \brief Compute the viscous momentum flux contribution from the interface
+     *        shear stress
+     *
+     *        The viscous momentum flux
+     *        \f[
+     *        \int \int_{V} \mathbf{\nabla} \cdot \nu_t h \mathbf{\nabla} \mathbf{u} dV
+     *        \f]
+     *        is re-written using Gauss' divergence theorem to:
+     *        \f[
+     *        \int_{S_f} \nu_t h \mathbf{\nabla} \mathbf{u} \cdot \mathbf{n_f} dS
+     *        \f]
+     *
+     * \todo The implementation now is the simplest one involving
+     *       only direct neighbours. This implementation is not complete/
+     *       correct on non-orthogonal meshes. A more complete implementation
+     *       with a more elaborate stencil that also takes into account
+     *       the non-orthogonal contributions can be considered at a later stage.
+     */
+    template<class Problem, class FVElementGeometry, class ElementVolumeVariables>
+    static NumEqVector flux(const Problem& problem,
+                            const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const typename FVElementGeometry::SubControlVolumeFace& scvf)
+    {
+        using Scalar = typename PrimaryVariables::value_type;
+
+        using std::sqrt;
+        using std::max;
+
+        NumEqVector localFlux(0.0);
+
+        // Get the inside and outside volume variables
+        const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+        const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
+        const auto& unitNormal = scvf.unitOuterNormal();
+
+        // Inside and outside subcontrolvolumes
+        const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+        const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
+
+        // Distance vector between the two cell centers
+        const auto& cellCenterToCellCenter = outsideScv.center() - insideScv.center();
+
+        // Check whether the mixing-length turbulence model is used
+        static const auto useMixingLengthTurbulenceModel = getParam<bool>("ShallowWater.UseMixingLengthTurbulenceModel", false);
+
+        // Compute the average shear velocity on the face
+        const Scalar ustar = [&](){
+            // If the mixing-length turbulence model is used:
+            // compute ustar based on the bottom friction
+            if (useMixingLengthTurbulenceModel)
+            {
+                // Get the bottom shear stress in the two adjacent cells
+                // Note that element is not needed in spatialParams().frictionLaw (should be removed). For now we simply pass the same element twice
+                const auto& bottomShearStressInside = problem.spatialParams().frictionLaw(element, insideScv).shearStress(insideVolVars);
+                const auto& bottomShearStressOutside = problem.spatialParams().frictionLaw(element, outsideScv).shearStress(outsideVolVars);
+                const auto bottomShearStressInsideMag = bottomShearStressInside.two_norm();
+                const auto bottomShearStressOutsideMag = bottomShearStressOutside.two_norm();
+
+                // Use a harmonic average of the viscosity and the depth at the interface.
+                const auto averageBottomShearStress = 2.0*(bottomShearStressInsideMag*bottomShearStressOutsideMag)/max(1.0e-8,bottomShearStressInsideMag + bottomShearStressOutsideMag);
+
+                // Shear velocity possibly needed in mixing-length turbulence model in the computation of the diffusion flux
+                return sqrt(averageBottomShearStress);
+            }
+            else
+                return 0.0;
+        }();
+
+        // The left (inside) and right (outside) states
+        const auto waterDepthLeft  = insideVolVars.waterDepth();
+        const auto velocityXLeft   = insideVolVars.velocity(0);
+        const auto velocityYLeft   = insideVolVars.velocity(1);
+        const auto waterDepthRight = outsideVolVars.waterDepth();
+        const auto velocityXRight  = outsideVolVars.velocity(0);
+        const auto velocityYRight  = outsideVolVars.velocity(1);
+
+        // Distance between the two adjacent cell centres
+        const auto distance = cellCenterToCellCenter.two_norm();
+
+        // Factor that takes the direction of the unit vector into account
+        const auto direction = (unitNormal*cellCenterToCellCenter)/distance;
+
+        // Compute the velocity gradients normal to the face
+        const auto gradU = (velocityXRight-velocityXLeft)*direction/distance;
+        const auto gradV = (velocityYRight-velocityYLeft)*direction/distance;
+
+        // Use a harmonic average of the depth at the interface.
+        const auto averageDepth = 2.0*(waterDepthLeft*waterDepthRight)/(waterDepthLeft + waterDepthRight);
+
+        // Now get the (constant) background turbulent viscosity
+        static const auto turbBGViscosity = getParam<Scalar>("ShallowWater.TurbulentViscosity", 1.0e-6);
+
+        Scalar turbViscosity = 0.0;
+
+        // constant eddy viscosity equal to the prescribed background eddy viscosity
+        if (!useMixingLengthTurbulenceModel)
+            turbViscosity = turbBGViscosity;
+
+        // turbulence model based on mixing length
+        else
+        {
+            // Compute the turbulent viscosity using a combined horizonal/vertical mixing length approach
+            // Turbulence coefficients: vertical (Elder like) and horizontal (Smagorinsky like)
+            static const auto turbConstV = getParam<Scalar>("ShallowWater.VerticalCoefficientOfMixingLengthModel", 1.0);
+            static const auto turbConstH = getParam<Scalar>("ShallowWater.HorizontalCoefficientOfMixingLengthModel", 0.1);
+
+            /** The vertical (Elder-like) contribution to the turbulent viscosity scales with water depth \f[ h \f] and shear velocity \f[ u_{*} \f] :
+            *
+            * \f[
+            * \nu_t^v = c^v \frac{\kappa}{6} u_{*} h
+            * \f]
+            */
+            static const Scalar kappa = 0.41;
+            const auto turbViscosityV = turbConstV * (kappa/6.0) * ustar * averageDepth;
+
+            /** The horizontal (Smagorinsky-like) contribution to the turbulent viscosity scales with the water depth (squared)
+            * and the magnitude of the stress (rate-of-strain) tensor:
+            *
+            * \f[
+            * nu_t^h = (c^h h)^2 \sqrt{ 2\left(\frac{\partial u}{\partial x}\right)^2 + \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right)^2 + 2\left(\frac{\partial v}{\partial y}\right)^2 }
+            * \f]
+            *
+            * However, based on the velocity vectors in the direct neighbours of the volume face, it is not possible to compute all components of the stress tensor.
+            * Only the gradients of u and v in the direction of the vector between the cell centres is available.
+            * To avoid the reconstruction of the full velocity gradient tensor based on a larger stencil,
+            * the horizontal contribution to the eddy viscosity (in the mixing-length model) is computed using only the velocity gradients normal to the face:
+            *
+            * \f[
+            * \frac{\partial u}{\partial n} , \frac{\partial v}{\partial n}
+            * \f]
+            *
+            * In other words, the present approximation of the horizontal contribution to the turbulent viscosity reduces to:
+            *
+            * \f[
+            * nu_t^h = (c^h h)^2 \sqrt{ 2\left(\frac{\partial u}{\partial n}\right)^2 + 2\left(\frac{\partial v}{\partial n}\right)^2 }
+            * \f]
+            *
+            It should be noted that this simplified approach is formally inconsistent and will result in a turbulent viscosity that is dependent on the grid (orientation).
+            */
+            const auto mixingLengthSquared = turbConstH * turbConstH * averageDepth * averageDepth;
+            const auto turbViscosityH      = mixingLengthSquared * sqrt(2.0*gradU*gradU + 2.0*gradV*gradV);
+
+            // Total turbulent viscosity
+            turbViscosity = turbBGViscosity + sqrt(turbViscosityV*turbViscosityV + turbViscosityH*turbViscosityH);
+        }
+
+        // Compute the viscous momentum fluxes
+        const auto uViscousFlux = turbViscosity * averageDepth * gradU;
+        const auto vViscousFlux = turbViscosity * averageDepth * gradV;
+
+        // compute the mobility of the flux with the fluxlimiter
+        static const auto upperWaterDepthFluxLimiting = getParam<double>("FluxLimiterLET.UpperWaterDepth", 1e-3);
+        static const auto lowerWaterDepthFluxLimiting = getParam<double>("FluxLimiterLET.LowerWaterDepth", 1e-5);
+
+        const auto limitingDepth = (waterDepthLeft + waterDepthRight) * 0.5;
+        const auto mobility = ShallowWater::fluxLimiterLET(limitingDepth,
+                                                             limitingDepth,
+                                                             upperWaterDepthFluxLimiting,
+                                                             lowerWaterDepthFluxLimiting);
+
+        localFlux[0] = 0.0;
+        localFlux[1] = -uViscousFlux * mobility * scvf.area();
+        localFlux[2] = -vViscousFlux * mobility * scvf.area();
+
+        return localFlux;
+    }
+};
+
+} // end namespace Dumux
+
+#endif

dumux/freeflow/shallowwater/boundaryfluxes.hh

@@ -5,7 +5,7 @@
  *                                                                           *
  *   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       *
+ *   the Free Software Foundation, either version 3 of the License, or       *
  *   (at your option) any later version.                                     *
  *                                                                           *
  *   This program is distributed in the hope that it will be useful,         *
@@ -32,14 +32,14 @@
 
 #include <array>
 #include <cmath>
+#include <algorithm>
 
-namespace Dumux {
-namespace ShallowWater {
+namespace Dumux::ShallowWater {
 
 /*!
  * \brief Compute the outer cell state for fixed water depth boundary.
  *
- * \param waterDepthBoundary Discharge per meter at the boundary face [m^2/s]
+ * \param waterDepthBoundary Water depth at the boundary face [m^2/s]
  * \param waterDepthInside Water depth in the inner cell [m]
  * \param velocityXInside Velocity in x-direction in the inner cell [m/s]
  * \param velocityYInside Velocity in y-direction in the inner cell [m/s]
@@ -127,7 +127,98 @@ std::array<Scalar, 3> fixedDischargeBoundary(const Scalar dischargeBoundary,
     return cellStateOutside;
 }
 
-} // end namespace ShallowWater
-} // end namespace Dumux
+/*!
+ * \brief Compute the viscosity/diffusive flux at a rough wall boundary using no-slip formulation.
+ *
+ * \param alphaWall Roughness parameter: alphaWall=0.0 means full slip, alphaWall=1.0 means no slip, 0.0<alphaWall<1.0 means partial slip [-]
+ * \param turbulentViscosity Turbulent viscosity [m^2/s]
+ * \param state Primary variables (water depth, velocities)
+ * \param cellCenterToBoundaryFaceCenter Cell-center to boundary distance
+ * \param unitNormal Normal vector of the boundary face
+ */
+template<class PrimaryVariables, class Scalar, class GlobalPosition>
+std::array<Scalar, 3> noslipWallBoundary(const Scalar alphaWall,
+                                         const Scalar turbulentViscosity,
+                                         const PrimaryVariables& state,
+                                         const GlobalPosition& cellCenterToBoundaryFaceCenter,
+                                         const GlobalPosition& unitNormal)
+{
+    // only impose if abs(alphaWall) > 0
+    using std::abs;
+    if (abs(alphaWall) <= 1.0e-9)
+        return {};
+
+    const auto waterDepth = state[0];
+    // regularization: Set gradients to zero for drying cell
+    // Use LET-limiter instead for differentiability?
+    if (waterDepth <= 0.001)
+        return {};
+
+    const auto xVelocity  = state[1];
+    const auto yVelocity  = state[2];
+    const auto distance = cellCenterToBoundaryFaceCenter.two_norm();
+
+    // Compute the velocity gradients
+    // Outside - inside cell: therefore the minus-sign
+    // Only when cell contains sufficient water.
+    const auto gradU = -alphaWall * xVelocity/distance;
+    const auto gradV = -alphaWall * yVelocity/distance;
+
+    // Factor that takes the direction of the unit vector into account
+    const auto direction = (unitNormal*cellCenterToBoundaryFaceCenter)/distance;
+
+    // Compute the viscosity/diffusive fluxes at the rough wall
+    return {
+        0.0,
+        -turbulentViscosity*waterDepth*gradU*direction,
+        -turbulentViscosity*waterDepth*gradV*direction
+    };
+}
+
+/*!
+ * \brief Compute the viscosity/diffusive flux at a rough wall boundary using Nikuradse formulation.
+ *
+ * \param ksWall Nikuradse roughness height for the wall [m]
+ * \param state the primary variable state (water depth, velocities)
+ * \param cellCenterToBoundaryFaceCenter Cell-center to boundary distance
+ * \param unitNormal Normal vector of the boundary face
+ */
+template<class PrimaryVariables, class Scalar, class GlobalPosition>
+std::array<Scalar, 3> nikuradseWallBoundary(const Scalar ksWall,
+                                            const PrimaryVariables& state,
+                                            const GlobalPosition& cellCenterToBoundaryFaceCenter,
+                                            const GlobalPosition& unitNormal)
+{
+    // only impose if abs(ksWall) > 0
+    using std::abs;
+    if (abs(ksWall) <= 1.0e-9)
+        return {};
+
+    using std::hypot;
+    const Scalar xVelocity = state[1];
+    const Scalar yVelocity = state[2];
+    const Scalar velocityMagnitude = hypot(xVelocity, yVelocity);
+    const Scalar distance = cellCenterToBoundaryFaceCenter.two_norm();
+    const Scalar y0w = ksWall/30.0;
+    constexpr Scalar kappa2 = 0.41*0.41;
+
+    // should distance/y0w be limited to not become too small?
+    using std::log; using std::max;
+    const auto logYPlus = log(distance/y0w+1.0);
+    const auto fac = kappa2*velocityMagnitude / max(1.0e-3,logYPlus*logYPlus);
+
+    // Factor that takes the direction of the unit vector into account
+    const auto direction = (unitNormal*cellCenterToBoundaryFaceCenter)/distance;
+
+    // wall shear stress vector
+    const auto tauWx = direction*fac*xVelocity;
+    const auto tauWy = direction*fac*yVelocity;
+
+    // Compute the viscosity/diffusive fluxes at the rough wall
+    const auto waterDepth = state[0];
+    return {0.0, waterDepth*tauWx, waterDepth*tauWy};
+}
+
+} // end namespace Dumux::ShallowWater
 
 #endif

dumux/freeflow/shallowwater/fluxvariables.hh

@@ -45,7 +45,7 @@ class ShallowWaterFluxVariables
     using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
     using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
     using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>;
-    //using DiffusionType = GetPropType<TypeTag, Properties::DiffusionType>;
+    using ViscousFluxType = GetPropType<TypeTag, Properties::ViscousFluxType>;
 
     using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
     using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
@@ -58,7 +58,6 @@ class ShallowWaterFluxVariables
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
 
     static constexpr bool enableAdvection = ModelTraits::enableAdvection();
-    static constexpr bool enableDiffusion = ModelTraits::enableDiffusion();
 
 public:
 
@@ -79,20 +78,17 @@ public:
     }
 
     /*!
-     * \brief Returns the diffusive flux (e.g. diffusion of tracer)
+     * \brief Returns the viscous momentum flux
      *
      */
-    NumEqVector diffusiveFlux(const Problem& problem,
-                              const Element& element,
-                              const FVElementGeometry& fvGeometry,
-                              const ElementVolumeVariables& elemVolVars,
-                              const SubControlVolumeFace& scvf) const
+    NumEqVector viscousFlux(const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const SubControlVolumeFace& scvf) const
     {
-        // TODO: add diffusive flux (e.g. tracer and viscosity)
-        if (enableDiffusion)
-            return NumEqVector(0.0);
-
-        return NumEqVector(0.0);
+        // Add viscous momentum flux
+        return ViscousFluxType::flux(problem, element, fvGeometry, elemVolVars, scvf);
     }
 };
 

dumux/freeflow/shallowwater/localresidual.hh

@@ -78,7 +78,7 @@ public:
     }
 
        /*!
-     * \brief Evaluate the mass flux over a face of a sub control volume
+     * \brief Evaluate the mass/momentum flux over a face of a sub control volume
      *
      * \param problem The problem
      * \param element The current element.
@@ -97,7 +97,11 @@ public:
         NumEqVector flux(0.0);
         FluxVariables fluxVars;
         flux += fluxVars.advectiveFlux(problem, element, fvGeometry, elemVolVars, scvf);
-        flux += fluxVars.diffusiveFlux(problem, element, fvGeometry, elemVolVars, scvf);
+
+        // Compute viscous momentum flux contribution if required
+        static const bool enableViscousFlux = getParam<bool>("ShallowWater.EnableViscousFlux", false);
+        if (enableViscousFlux)
+            flux += fluxVars.viscousFlux(problem, element, fvGeometry, elemVolVars, scvf);
         return flux;
     }
 };

dumux/freeflow/shallowwater/model.hh

@@ -66,6 +66,7 @@
 #include <dumux/common/properties/model.hh>
 
 #include <dumux/flux/shallowwaterflux.hh>
+#include <dumux/flux/shallowwaterviscousflux.hh>
 #include <dumux/flux/fluxvariablescaching.hh>
 
 #include "localresidual.hh"
@@ -89,9 +90,6 @@ struct ShallowWaterModelTraits
 
     //! Enable advection
     static constexpr bool enableAdvection() { return true; }
-
-    //! Enable diffusion
-    static constexpr bool enableDiffusion() { return false; }
 };
 
 /*!
@@ -160,7 +158,9 @@ template<class TypeTag>
 struct AdvectionType<TypeTag, TTag::ShallowWater>
 { using type = ShallowWaterFlux< GetPropType<TypeTag, Properties::NumEqVector> >; };
 
-//template<class TypeTag> struct DiffusionType<TypeTag, TTag::ShallowWater> {using type = ShallowWaterExactRiemannSolver<TypeTag>;};
+template<class TypeTag>
+struct ViscousFluxType<TypeTag, TTag::ShallowWater>
+{ using type = ShallowWaterViscousFlux< GetPropType<TypeTag, Properties::PrimaryVariables>, GetPropType<TypeTag, Properties::NumEqVector> >; };
 
 } // end properties
 } // end namespace Dumux

dumux/freeflow/shallowwater/volumevariables.hh

@@ -57,6 +57,10 @@ public:
     Scalar extrusionFactor() const
     { return 1.0; }
 
+    //! Return the vector of primary variables
+    const PrimaryVariables& priVars() const
+    { return priVars_; }
+
      /*!
      * \brief Return water detph h inside the sub-control volume.
      *

test/freeflow/shallowwater/CMakeLists.txt

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

test/freeflow/shallowwater/poiseuilleflow/CMakeLists.txt

+dune_symlink_to_source_files(FILES "params.input")
+
+dumux_add_test(NAME test_shallowwater_poiseuilleflow
+               SOURCES main.cc
+               LABELS shallowwater
+               COMPILE_DEFINITIONS GRIDTYPE=Dune::YaspGrid<2,Dune::EquidistantOffsetCoordinates<double,2>>
+               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+               CMD_ARGS       --script fuzzy
+                              --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_shallowwater_poiseuilleflow-reference.vtu
+                                      ${CMAKE_CURRENT_BINARY_DIR}/poiseuilleflow-00007.vtu
+                              --zeroThreshold {"velocityY":1e-14}
+                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_shallowwater_poiseuilleflow params.input")
+
+dumux_add_test(NAME test_shallowwater_poiseuilleflow_parallel
+               TARGET test_shallowwater_poiseuilleflow
+               LABELS shallowwater
+               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_poiseuilleflow-reference.vtu
+                                  ${CMAKE_CURRENT_BINARY_DIR}/s0002-poiseuilleflow-parallel-00007.pvtu
+                          --zeroThreshold {"velocityY":1e-14,"process rank":100}
+                          --command "${MPIEXEC} -np 2 ${CMAKE_CURRENT_BINARY_DIR}/test_shallowwater_poiseuilleflow params.input -Problem.Name poiseuilleflow-parallel")
+
+dumux_add_test(NAME test_shallowwater_poiseuilleflow_unstructured
+               SOURCES main.cc
+               LABELS shallowwater
+               CMAKE_GUARD dune-uggrid_FOUND
+               COMPILE_DEFINITIONS GRIDTYPE=Dune::UGGrid<2>
+               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+               CMD_ARGS       --script fuzzy
+                              --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_shallowwater_poiseuilleflow_unstructured-reference.vtu
+                                      ${CMAKE_CURRENT_BINARY_DIR}/poiseuilleflow-unstructured-00007.vtu
+                              --zeroThreshold {"velocityY":1e-14}
+                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_shallowwater_poiseuilleflow_unstructured params.input -Problem.Name poiseuilleflow-unstructured -Grid.File grids/irregular_grid_10m.dgf")
+
+dumux_add_test(NAME test_shallowwater_poiseuilleflow_unstructured_parallel
+               TARGET test_shallowwater_poiseuilleflow_unstructured
+               LABELS shallowwater
+               CMAKE_GUARD "( MPI_FOUND AND dune-uggrid_FOUND )"
+               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+               CMD_ARGS   --script fuzzy
+                          --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_shallowwater_poiseuilleflow_unstructured-reference.vtu
+                                  ${CMAKE_CURRENT_BINARY_DIR}/s0002-poiseuilleflow-unstructured-parallel-00007.pvtu
+                          --zeroThreshold {"velocityY":1e-14,"process rank":100}
+                          --command "${MPIEXEC} -np 2 ${CMAKE_CURRENT_BINARY_DIR}/test_shallowwater_poiseuilleflow_unstructured params.input -Problem.Name poiseuilleflow-unstructured-parallel -Grid.File grids/irregular_grid_10m.dgf")
+
+dune_symlink_to_source_files(FILES "grids")

test/freeflow/shallowwater/poiseuilleflow/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 3 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/>.   *
+ *****************************************************************************/
+
+#include <config.h>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+
+#include <dumux/linear/linearsolvertraits.hh>
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+#include <dumux/assembly/fvassembler.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+
+#include <dumux/io/grid/gridmanager_yasp.hh>
+#include <dumux/io/grid/gridmanager_ug.hh>
+
+#include "properties.hh"
+
+int main(int argc, char** argv)
+{
+    using namespace Dumux;
+
+    const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+    // parse command line arguments and input file
+    Parameters::init(argc, argv);
+
+    // the property tag
+    using TypeTag = Properties::TTag::PoiseuilleFlow;
+
+    // create grid
+    GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
+    gridManager.init();
+    const auto& leafGridView = gridManager.grid().leafGridView();
+
+    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+    auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
+    gridGeometry->update();
+
+    // problem, in which we define the boundary and initial conditions.
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    auto problem = std::make_shared<Problem>(gridGeometry);
+
+    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+    SolutionVector x;