diff --git a/dumux/flux/shallowwaterviscousflux.hh b/dumux/flux/shallowwaterviscousflux.hh
index ce0ce087f8ded89e8c903c5eff2a6ecb36b4e81b..592e4bcd6bd777c2b751815d407d773a50957843 100644
--- a/dumux/flux/shallowwaterviscousflux.hh
+++ b/dumux/flux/shallowwaterviscousflux.hh
@@ -28,6 +28,10 @@
#include <algorithm>
#include <utility>
#include <type_traits>
+#include <array>
+
+#include <dune/common/std/type_traits.hh>
+#include <dune/common/exceptions.hh>
#include <dumux/common/parameters.hh>
#include <dumux/flux/fluxvariablescaching.hh>
@@ -35,6 +39,18 @@
namespace Dumux {
+#ifndef DOXYGEN
+namespace Detail {
+// helper struct detecting if the user-defined spatial params class has a frictionLaw function
+template <typename T, typename ...Ts>
+using FrictionLawDetector = decltype(std::declval<T>().frictionLaw(std::declval<Ts>()...));
+
+template<class T, typename ...Args>
+static constexpr bool implementsFrictionLaw()
+{ return Dune::Std::is_detected<FrictionLawDetector, T, Args...>::value; }
+} // end namespace Detail
+#endif
+
/*!
* \ingroup Flux
* \brief Computes the shallow water viscous momentum flux due to (turbulent) viscosity
@@ -88,7 +104,7 @@ public:
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
- const auto [gradU, gradV] = [&]()
+ const auto gradVelocity = [&]()
{
// The left (inside) and right (outside) states
const auto velocityXLeft = insideVolVars.velocity(0);
@@ -102,10 +118,10 @@ public:
const auto distance = cellCenterToCellCenter.two_norm();
const auto& unitNormal = scvf.unitOuterNormal();
const auto direction = (unitNormal*cellCenterToCellCenter)/distance;
- return std::make_pair(
+ return std::array<Scalar, 2>{
(velocityXRight-velocityXLeft)*direction/distance,
(velocityYRight-velocityYLeft)*direction/distance
- );
+ };
}();
// Use a harmonic average of the depth at the interface.
@@ -114,7 +130,7 @@ public:
const auto averageDepth = 2.0*(waterDepthLeft*waterDepthRight)/(waterDepthLeft + waterDepthRight);
// compute the turbulent viscosity contribution
- const Scalar turbViscosity = [&,gradU=gradU,gradV=gradV]()
+ const Scalar turbViscosity = [&]()
{
// The (constant) background turbulent viscosity
static const auto turbBGViscosity = getParamFromGroup<Scalar>(problem.paramGroup(), "ShallowWater.TurbulentViscosity", 1.0e-6);
@@ -126,74 +142,84 @@ public:
if (!useMixingLengthTurbulenceModel)
return turbBGViscosity;
- // turbulence model based on mixing length
- // 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 = getParamFromGroup<Scalar>(problem.paramGroup(), "ShallowWater.VerticalCoefficientOfMixingLengthModel", 1.0);
- static const auto turbConstH = getParamFromGroup<Scalar>(problem.paramGroup(), "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]
- */
- constexpr Scalar kappa = 0.41;
- // Compute the average shear velocity on the face
- const Scalar ustar = [&]()
+ using SpatialParams = typename Problem::SpatialParams;
+ using E = typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity;
+ using SCV = typename FVElementGeometry::SubControlVolume;
+ if constexpr (!Detail::implementsFrictionLaw<SpatialParams, E, SCV>())
+ DUNE_THROW(Dune::IOError, "Mixing length turbulence model enabled but spatial parameters do not implement the frictionLaw interface!");
+ else
{
- // 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.
- using std::max;
- 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
+ // turbulence model based on mixing length
+ // 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 = getParamFromGroup<Scalar>(problem.paramGroup(), "ShallowWater.VerticalCoefficientOfMixingLengthModel", 1.0);
+ static const auto turbConstH = getParamFromGroup<Scalar>(problem.paramGroup(), "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]
+ */
+ constexpr Scalar kappa = 0.41;
+ // Compute the average shear velocity on the face
+ const Scalar ustar = [&]()
+ {
+ // 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.
+ using std::max;
+ 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
+ using std::sqrt;
+ return sqrt(averageBottomShearStress);
+ }();
+
+ 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).
+ */
using std::sqrt;
- return sqrt(averageBottomShearStress);
- }();
-
- 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).
- */
- using std::sqrt;
- const auto mixingLengthSquared = turbConstH * turbConstH * averageDepth * averageDepth;
- const auto turbViscosityH = mixingLengthSquared * sqrt(2.0*gradU*gradU + 2.0*gradV*gradV);
-
- // Total turbulent viscosity
- return turbBGViscosity + sqrt(turbViscosityV*turbViscosityV + turbViscosityH*turbViscosityH);
+ const auto [gradU, gradV] = gradVelocity;
+ const auto mixingLengthSquared = turbConstH * turbConstH * averageDepth * averageDepth;
+ const auto turbViscosityH = mixingLengthSquared * sqrt(2.0*gradU*gradU + 2.0*gradV*gradV);
+
+ // Total turbulent viscosity
+ return turbBGViscosity + sqrt(turbViscosityV*turbViscosityV + turbViscosityH*turbViscosityH);
+ }
}();
// Compute the viscous momentum fluxes
+ const auto [gradU, gradV] = gradVelocity;
const auto uViscousFlux = turbViscosity * averageDepth * gradU;
const auto vViscousFlux = turbViscosity * averageDepth * gradV;
diff --git a/test/freeflow/shallowwater/poiseuilleflow/spatialparams.hh b/test/freeflow/shallowwater/poiseuilleflow/spatialparams.hh
index eb35e02a1f3dfef7566f3b49250d379054430ad6..5e858f3ae1f88dfc28e53313f2aebc2cbdf8be8d 100644
--- a/test/freeflow/shallowwater/poiseuilleflow/spatialparams.hh
+++ b/test/freeflow/shallowwater/poiseuilleflow/spatialparams.hh
@@ -28,7 +28,6 @@
#include <dumux/common/parameters.hh>
#include <dumux/material/spatialparams/fv.hh>
-#include <dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh>
namespace Dumux {
@@ -65,12 +64,6 @@ public:
Scalar gravity() const
{ return gravity_; }
- const FrictionLaw<VolumeVariables>&
- frictionLaw(const Element& element, const SubControlVolume& scv) const
- {
- DUNE_THROW(Dune::NotImplemented, "No friction law is implemented for this test!");
- }
-
//! Define the bed surface
Scalar bedSurface(const Element& element, const SubControlVolume& scv) const
{