Skip to content
Snippets Groups Projects
Commit 32d202f2 authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'feature/swe-fix-require-friction' into 'master' 

[swe][viscous] Do not require spatial params to implement friction law

See merge request !2492
parents ae9f7563 9115d96e
No related branches found
No related tags found
1 merge request!2492[swe][viscous] Do not require spatial params to implement friction law
Hide whitespace changes
Inline Side-by-side
dumux/flux/shallowwaterviscousflux.hh
+ 93
− 67
View file @ 32d202f2
......@@ -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;
......
test/freeflow/shallowwater/poiseuilleflow/spatialparams.hh
+ 0
− 7
View file @ 32d202f2
......@@ -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
{
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment