From 9115d96e53ceb69417d06d83e203ac81117ed606 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 27 Feb 2021 21:51:47 +0100 Subject: [PATCH] [swe][viscous] Do not require spatial params to implement friction law --- dumux/flux/shallowwaterviscousflux.hh | 160 ++++++++++-------- .../poiseuilleflow/spatialparams.hh | 7 - 2 files changed, 93 insertions(+), 74 deletions(-) diff --git a/dumux/flux/shallowwaterviscousflux.hh b/dumux/flux/shallowwaterviscousflux.hh index ce0ce087f8..592e4bcd6b 100644 --- a/dumux/flux/shallowwaterviscousflux.hh +++ b/dumux/flux/shallowwaterviscousflux.hh @@ -28,6 +28,10 @@ #include #include #include +#include + +#include +#include #include #include @@ -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 +using FrictionLawDetector = decltype(std::declval().frictionLaw(std::declval()...)); + +template +static constexpr bool implementsFrictionLaw() +{ return Dune::Std::is_detected::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{ (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(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(problem.paramGroup(), "ShallowWater.VerticalCoefficientOfMixingLengthModel", 1.0); - static const auto turbConstH = getParamFromGroup(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()) + 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(problem.paramGroup(), "ShallowWater.VerticalCoefficientOfMixingLengthModel", 1.0); + static const auto turbConstH = getParamFromGroup(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 eb35e02a1f..5e858f3ae1 100644 --- a/test/freeflow/shallowwater/poiseuilleflow/spatialparams.hh +++ b/test/freeflow/shallowwater/poiseuilleflow/spatialparams.hh @@ -28,7 +28,6 @@ #include #include -#include namespace Dumux { @@ -65,12 +64,6 @@ public: Scalar gravity() const { return gravity_; } - const FrictionLaw& - 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 { -- GitLab