Commit 813cf847 authored by Lars Kaiser's avatar Lars Kaiser Committed by Martin Schneider
Browse files

[newIc] Further incorporate new interface conditions

parent 2265b93a
......@@ -52,7 +52,7 @@ public:
ParentType::resetEq(eqIdx);
boundaryInfo_[eqIdx].isSymmetry = false;
boundaryInfo_[eqIdx].isBeaversJoseph = false;
boundaryInfo_[eqIdx].isSlipCondition = false;
}
/*!
......@@ -77,10 +77,20 @@ public:
* \brief Set a boundary condition for a single equation to
* Beavers-Joseph(-Saffmann) (special case of Dirichlet b.c.).
*/
[[deprecated("Use setSlipCondition(eqIdx) instead! Will be removed after release 3.5.")]]
void setBeaversJoseph(const int eqIdx)
{
setSlipCondition(eqIdx);
}
/*!
* \brief Set a boundary condition for a single equation to
* slip condition: e.g. Beavers-Joseph(-Saffmann) (special case of Dirichlet b.c.).
*/
void setSlipCondition(const int eqIdx)
{
resetEq(eqIdx);
boundaryInfo_[eqIdx].isBeaversJoseph = true;
boundaryInfo_[eqIdx].isSlipCondition = true;
}
/*!
......@@ -89,17 +99,37 @@ public:
*
* \param eqIdx The index of the equation
*/
[[deprecated("Use isSlipCondition(eqIdx) instead! Will be removed after release 3.5.")]]
bool isBeaversJoseph(const int eqIdx) const
{ return boundaryInfo_[eqIdx].isBeaversJoseph; }
/*!
* \brief Returns true if an equation is used to specify a
* slip condition.
*
* \param eqIdx The index of the equation
*/
bool isSlipCondition(const int eqIdx) const
{ return boundaryInfo_[eqIdx].isSlipCondition; }
/*!
* \brief Returns true if some equation is used to specify a
* Beavers-Joseph(-Saffman) boundary condition.
*/
[[deprecated("Use hasSlipCondition(eqIdx) instead! Will be removed after release 3.5.")]]
bool hasBeaversJoseph() const
{
return hasSlipCondition();
}
/*!
* \brief Returns true if some equation is used to specify a
* slip condition.
*/
bool hasSlipCondition() const
{
for (int i = 0; i < numEq; ++i)
if (boundaryInfo_[i].isBeaversJoseph)
if (boundaryInfo_[i].isSlipCondition)
return true;
return false;
}
......@@ -109,7 +139,7 @@ protected:
struct NavierStokesBoundaryInfo
{
bool isSymmetry : 1;
bool isBeaversJoseph : 1;
bool isSlipCondition : 1;
};
std::array<NavierStokesBoundaryInfo, numEq> boundaryInfo_;
......
......@@ -224,58 +224,47 @@ public:
return asImp_().alphaBJ(scvf) / sqrt(asImp_().permeability(element, scvf));
}
/*!
* \brief Returns the scale separation value eps = (macroscopicLength/LengthOfPeriodicityCell).
*/
Scalar epsInterface(const SubControlVolumeFace& scvf) const
{ DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the epsInterface value must be returned in the acutal problem");}
Scalar factorNMomentum(const SubControlVolumeFace& scvf) const
{ DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the factorNMomentum value must be returned in the acutal problem");}
/*!
* \brief Returns the boundary layer constant N^bl computed based on the exact pore geometry.
*/
Scalar factorNTangential(const SubControlVolumeFace& scvf) const
{ DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the factorNTangential value must be returned in the acutal problem");}
Dune::FieldMatrix<Scalar, dim, dim> matrixNTangential(const SubControlVolumeFace& scvf) const
{ DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the matrixNTangential value must be returned in the acutal problem");}
/*!
* \brief Returns the velocity in the porous medium (which is 0 by default according to Saffmann).
*/
VelocityVector porousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const
{
return VelocityVector(0.0);//TODO: -> Dont force implementation?
return VelocityVector(0.0);
}
//! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph condition is used
const Scalar beaversJosephVelocity(const Element& element,
//! helper function to evaluate the slip velocity on the boundary when the slipCondition is used
const Scalar slipVelocity(const Element& element,
const SubControlVolume& scv,
const SubControlVolumeFace& ownScvf,
const SubControlVolumeFace& faceOnPorousBoundary,
const Scalar velocitySelf,
const Scalar tangentialVelocityGradient) const
{
// create a unit normal vector oriented in positive coordinate direction
GlobalPosition orientation = ownScvf.unitOuterNormal();
orientation[ownScvf.directionIndex()] = 1.0;
// du/dy + dv/dx = alpha/sqrt(K) * (u_boundary-uPM)
// beta = alpha/sqrt(K)
const Scalar betaBJ = asImp_().betaBJ(element, faceOnPorousBoundary, orientation);
const Scalar distanceNormalToBoundary = (faceOnPorousBoundary.center() - scv.center()).two_norm();
return (tangentialVelocityGradient*distanceNormalToBoundary
+ asImp_().porousMediumVelocity(element, faceOnPorousBoundary) * orientation * betaBJ * distanceNormalToBoundary
+ velocitySelf) / (betaBJ*distanceNormalToBoundary + 1.0);
}
//TODO: I assume an inverted sign for factorNTangential...
//! helper function to evaluate the slip velocity on the boundary when the new tangential condition is used
const Scalar nTangentialVelocity(const Element& element,
const SubControlVolume& scv,
const SubControlVolumeFace& ownScvf, //stehendes
const SubControlVolumeFace& faceOnPorousBoundary, //liegendes
const Scalar velocitySelf, //vel auf stehendem
const Scalar tangentialVelocityGradient) const //dv/dx=0
const Scalar velocitySelf, //vel perpendicular to tangential vel
const Scalar tangentialVelocityGradient) const //dv/dx (=0)
{
const Scalar factor = 1.0/(asImp_().epsInterface(faceOnPorousBoundary) * asImp_().factorNTangential(faceOnPorousBoundary));
static const bool newIc_ = getParamFromGroup<bool>("Problem", "NewIc", false);
// du/dy + dv/dx = factor * (u_boundary-uPM)
Scalar factor;
if (newIc_){
//return nTangentialVelocity(element, scv,ownScvf, faceOnPorousBoundary, velocitySelf, tangentialVelocityGradient);
factor = -1.0 / asImp_().epsInterface(faceOnPorousBoundary) / asImp_().factorNTangential(faceOnPorousBoundary);
}
else{
//return beaversJosephVelocity(element, scv,ownScvf, faceOnPorousBoundary, velocitySelf, tangentialVelocityGradient);
factor = asImp_().betaBJ(element, faceOnPorousBoundary); //beta = alpha/sqrt(K)
}
const Scalar distanceNormalToBoundary = (faceOnPorousBoundary.center() - scv.center()).two_norm();
// create a unit normal vector oriented in positive coordinate direction
......
......@@ -321,9 +321,9 @@ public:
lateralFaceBoundaryTypes.emplace(problem.boundaryTypes(element, lateralFace));
}
// If the current scvf is on a bounary and if there is a Neumann or Beavers-Joseph-(Saffmann) BC for the stress in tangential direction,
// If the current scvf is on a bounary and if there is a Neumann or SlipCondition BC for the stress in tangential direction,
// assign this value for the lateral momentum flux. No further calculations are required. We assume that all lateral faces
// have the same type of BC (Neumann or Beavers-Joseph-(Saffmann)), but we sample the value at their actual positions.
// have the same type of BC (Neumann or SlipCondition), but we sample the value at their actual positions.
if (currentScvfBoundaryTypes)
{
// Handle Neumann BCs.
......@@ -343,10 +343,9 @@ public:
// It is not clear how to evaluate the BJ condition here.
// For symmetry reasons, our own scvf should then have the same Neumann flux as the lateral face.
// TODO: We should clarify if this is the correct approach.
bool bj = currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex()));
bool nTangential = currentScvfBoundaryTypes->isNTangential(Indices::velocity(lateralScvf.directionIndex()));
bool slipCondition = currentScvfBoundaryTypes->isSlipCondition(Indices::velocity(lateralScvf.directionIndex()));
bool neumann = lateralFaceBoundaryTypes->isNeumann(Indices::velocity(scvf.directionIndex()));
if ( (bj || nTangential) && lateralFaceBoundaryTypes && neumann)
if ( slipCondition && lateralFaceBoundaryTypes && neumann)
{
FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area());
const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
......@@ -389,11 +388,10 @@ public:
// Check the consistency of the boundary conditions, exactly one of the following must be set
if (lateralFaceBoundaryTypes)
{
std::bitset<4> admittableBcTypes;
std::bitset<3> admittableBcTypes;
admittableBcTypes.set(0, lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx));
admittableBcTypes.set(1, lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex())));
admittableBcTypes.set(2, lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())));
admittableBcTypes.set(3, lateralFaceBoundaryTypes->isNTangential(Indices::velocity(scvf.directionIndex())));
admittableBcTypes.set(2, lateralFaceBoundaryTypes->isSlipCondition(Indices::velocity(scvf.directionIndex())));
if (admittableBcTypes.count() != 1)
{
DUNE_THROW(Dune::InvalidStateException, "Invalid boundary conditions for lateral scvf "
......@@ -522,15 +520,11 @@ private:
const auto lateralBoundaryFace = makeStaggeredBoundaryFace(scvf, lateralBoundaryFacePos);
return problem.dirichlet(element, lateralBoundaryFace)[Indices::velocity(lateralFace.directionIndex())];
}
else if (bcTypes.isBeaversJoseph(Indices::velocity(lateralFace.directionIndex())))
else if (bcTypes.isSlipCondition(Indices::velocity(lateralFace.directionIndex())))
{
return VelocityGradients::beaversJosephVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars,
return VelocityGradients::slipVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars,
currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
}
else if (bcTypes.isNTangential(Indices::velocity(lateralFace.directionIndex()))){
return VelocityGradients::nTangentialVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars,
currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
}
else
return faceVars.velocityLateralInside(localSubFaceIdx);
}
......@@ -599,8 +593,7 @@ private:
{
if (!scvf.boundary() ||
currentScvfBoundaryTypes->isDirichlet(Indices::velocity(lateralFace.directionIndex())) ||
currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralFace.directionIndex())) ||
currentScvfBoundaryTypes->isNTangential(Indices::velocity(lateralFace.directionIndex()))
currentScvfBoundaryTypes->isSlipCondition(Indices::velocity(lateralFace.directionIndex()))
)
{
const Scalar velocityGrad_ji = VelocityGradients::velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
......
......@@ -583,8 +583,7 @@ private:
// Find out what boundary type is set on the lateral face
const bool useZeroGradient = lateralFaceBoundaryTypes && (lateralFaceBoundaryTypes->isSymmetry()
|| lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx));
const bool lateralFaceHasBJS = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex()));
const bool lateralFaceHasNTangential = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isNTangential(Indices::velocity(scvf.directionIndex()));
const bool lateralFaceHasSlipCondition = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isSlipCondition(Indices::velocity(scvf.directionIndex()));
const bool lateralFaceHasDirichletVelocity = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex()));
const Scalar velocitySelf = faceVars.velocitySelf();
......@@ -593,12 +592,8 @@ private:
if (useZeroGradient)
return velocitySelf;
if (lateralFaceHasBJS)
return VelocityGradients::beaversJosephVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars,
currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
else if(lateralFaceHasNTangential)
return VelocityGradients::nTangentialVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars,
if (lateralFaceHasSlipCondition)
return VelocityGradients::slipVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars,
currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
else if(lateralFaceHasDirichletVelocity)
{
......
......@@ -142,15 +142,10 @@ public:
const auto lateralBoundaryFace = makeStaggeredBoundaryFace(lateralScvf, lateralBoundaryFacePos);
return problem.dirichlet(element, lateralBoundaryFace)[Indices::velocity(scvf.directionIndex())];
}
else if (lateralFaceBoundaryTypes->isNTangential(Indices::velocity(scvf.directionIndex()))){
return nTangentialVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars,
else if (lateralFaceBoundaryTypes->isSlipCondition(Indices::velocity(scvf.directionIndex()))){
return slipVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars,
currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
}
else if (lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())))
{
return beaversJosephVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars,
currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
}
else
DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary type at " << lateralScvf.center());
}();
......@@ -226,15 +221,10 @@ public:
const auto lateralBoundaryFace = makeStaggeredBoundaryFace(scvf, lateralBoundaryFacePos);
return problem.dirichlet(element, lateralBoundaryFace)[Indices::velocity(lateralScvf.directionIndex())];
}
else if (currentScvfBoundaryTypes->isNTangential(Indices::velocity(lateralScvf.directionIndex()))){
return nTangentialVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars,
else if (currentScvfBoundaryTypes->isSlipCondition(Indices::velocity(lateralScvf.directionIndex()))){
return slipVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars,
currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
}
else if (currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex())))
{
return beaversJosephVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars,
currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
}
else
DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary types at " << lateralScvf.center());
}();
......@@ -248,10 +238,10 @@ public:
}
/*!
* \brief Returns the Beavers-Jospeh slip velocity for a scvf which lies on the boundary itself.
* \brief Returns the slip velocity for a scvf which lies on the boundary itself.
*
* \verbatim
* in.norm. B-J slip
* in.norm. slip
* vel. vel.
* ^ ^
* | |
......@@ -268,7 +258,7 @@ public:
*
*/
template<class Problem, class FaceVariables>
static Scalar beaversJosephVelocityAtCurrentScvf(const Problem& problem,
static Scalar slipVelocityAtCurrentScvf(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
......@@ -283,81 +273,41 @@ public:
const auto tangentialVelocityGradient = [&]()
{
// If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a BJ condition for
// the slip velocity is set there, assume a tangential velocity gradient of zero along the lateral face
// If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a slip condition
// is set there, assume a tangential velocity gradient of zero along the lateral face
// (towards the current scvf).
static const bool unsymmetrizedGradientForBJ = getParamFromGroup<bool>(problem.paramGroup(),
static const bool unsymmetrizedGradientForBeaversJoseph = getParamFromGroup<bool>(problem.paramGroup(),
"FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false);
//TODO: Deprecate and replace unsymmetrizedGradientForBeaversJoseph below by false, when deprecation period expired
static const bool unsymmetrizedGradientForIC = getParamFromGroup<bool>(problem.paramGroup(),
"FreeFlow.EnableUnsymmetrizedVelocityGradientForIC", unsymmetrizedGradientForBeaversJoseph);
if (unsymmetrizedGradientForBJ)
if (unsymmetrizedGradientForIC)
return 0.0;
if (lateralScvf.boundary())
{
if (lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())))
lateralFaceBoundaryTypes->isSlipCondition(Indices::velocity(scvf.directionIndex())))
return 0.0;
}
return velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
}();
return problem.beaversJosephVelocity(element,
fvGeometry.scv(scvf.insideScvIdx()),
lateralScvf,
scvf, /*on boundary*/
innerLateralVelocity,
tangentialVelocityGradient);
}
template<class Problem, class FaceVariables>
static Scalar nTangentialVelocityAtCurrentScvf(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const FaceVariables& faceVars,
const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
const std::size_t localSubFaceIdx)
{
const auto eIdx = scvf.insideScvIdx();
const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
const Scalar innerLateralVelocity = faceVars.velocityLateralInside(localSubFaceIdx);
const auto tangentialVelocityGradient = [&]()
{
// If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a BJ condition for
// the slip velocity is set there, assume a tangential velocity gradient of zero along the lateral face
// (towards the current scvf).
static const bool unsymmetrizedGradientForBJ = getParamFromGroup<bool>(problem.paramGroup(),
"FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false);
if (unsymmetrizedGradientForBJ)
return 0.0;
if (lateralScvf.boundary())
{
if (lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
lateralFaceBoundaryTypes->isNTangential(Indices::velocity(scvf.directionIndex())))
return 0.0;
}
return velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
}();
return problem.nTangentialVelocity(element,
fvGeometry.scv(scvf.insideScvIdx()),
lateralScvf,
scvf, /*on boundary*/
innerLateralVelocity,
tangentialVelocityGradient);
return problem.slipVelocity(element,
fvGeometry.scv(scvf.insideScvIdx()),
lateralScvf,
scvf, /*on boundary*/
innerLateralVelocity,
tangentialVelocityGradient);
}
/*!
* \brief Returns the Beavers-Jospeh slip velocity for a lateral scvf which lies on the boundary.
* \brief Returns the slip velocity for a lateral scvf which lies on the boundary.
*
* \verbatim
* B-J slip * boundary
* slip * boundary
* ************** vel. *****
* scvf ---------##### ~~~~> :::: || and # staggered half-control-volume (own element)
* | || | curr. ::
......@@ -371,7 +321,7 @@ public:
* \endverbatim
*/
template<class Problem, class FaceVariables>
static Scalar beaversJosephVelocityAtLateralScvf(const Problem& problem,
static Scalar slipVelocityAtLateralScvf(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
......@@ -386,74 +336,33 @@ public:
const auto tangentialVelocityGradient = [&]()
{
// If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a BJ condition for
// the slip velocity is set there, assume a tangential velocity gradient of zero along the lateral face
// If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a slip condition
// is set there, assume a tangential velocity gradient of zero along the lateral face
// (towards the current scvf).
static const bool unsymmetrizedGradientForBJ = getParamFromGroup<bool>(problem.paramGroup(),
static const bool unsymmetrizedGradientForBeaversJoseph = getParamFromGroup<bool>(problem.paramGroup(),
"FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false);
if (unsymmetrizedGradientForBJ)
return 0.0;
if (scvf.boundary())
{
if (currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex())))
return 0.0;
}
return velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
}();
return problem.beaversJosephVelocity(element,
fvGeometry.scv(scvf.insideScvIdx()),
scvf,
lateralScvf, /*on boundary*/
innerParallelVelocity,
tangentialVelocityGradient);
}
template<class Problem, class FaceVariables>
static Scalar nTangentialVelocityAtLateralScvf(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const FaceVariables& faceVars,
const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
const std::size_t localSubFaceIdx)
{
const auto eIdx = scvf.insideScvIdx();
const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
const Scalar innerParallelVelocity = faceVars.velocitySelf();
const auto tangentialVelocityGradient = [&]()
{
// If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a BJ condition for
// the slip velocity is set there, assume a tangential velocity gradient of zero along the lateral face
// (towards the current scvf).
static const bool unsymmetrizedGradientForBJ = getParamFromGroup<bool>(problem.paramGroup(),
"FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false);
if (unsymmetrizedGradientForBJ)
// TODO: Deprecate and replace unsymmetrizedGradientForBeaversJoseph below by false, when deprecation period expired
static const bool unsymmetrizedGradientForIC = getParamFromGroup<bool>(problem.paramGroup(),
"FreeFlow.EnableUnsymmetrizedVelocityGradientForIC", unsymmetrizedGradientForBeaversJoseph);
if (unsymmetrizedGradientForIC)
return 0.0;
if (scvf.boundary())
{
if (currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
currentScvfBoundaryTypes->isNTangential(Indices::velocity(lateralScvf.directionIndex())))
currentScvfBoundaryTypes->isSlipCondition(Indices::velocity(lateralScvf.directionIndex())))
return 0.0;
}
return velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
}();
return problem.nTangentialVelocity(element,
fvGeometry.scv(scvf.insideScvIdx()),
scvf,
lateralScvf, /*on boundary*/
innerParallelVelocity,
tangentialVelocityGradient);
return problem.slipVelocity(element,
fvGeometry.scv(scvf.insideScvIdx()),
scvf,
lateralScvf, /*on boundary*/
innerParallelVelocity,
tangentialVelocityGradient);
}
private:
......
......@@ -527,10 +527,10 @@ private:
}
}
// Calculate the BJS-velocity by accounting for all sub faces.
std::vector<int> bjsNumFaces(dim, 0);
std::vector<unsigned int> bjsNeighbor(dim, 0);
DimVector bjsVelocityAverage(0.0);
// Calculate the slip-velocity by accounting for all sub faces.
std::vector<int> slipNumFaces(dim, 0);
std::vector<unsigned int> slipNeighbor(dim, 0);
DimVector slipVelocityAverage(0.0);
DimVector normalNormCoordinate(0.0);
unsigned int velIdx = Indices::velocity(scvfNormDim);
const int numSubFaces = scvf.pairData().size();
......@@ -538,34 +538,34 @@ private:
{
const auto& lateralFace = fvGeometry.scvf(scvf.insideScvIdx(), scvf.pairData()[localSubFaceIdx].localLateralFaceIdx);
// adapt calculations for Beavers-Joseph-Saffman condition
// adapt calculations for slip condition
unsigned int normalNormDim = lateralFace.directionIndex();
if (lateralFace.boundary() && (asImp_().boundaryTypes(element, lateralFace).isBeaversJoseph(Indices::velocity(velIdx))))
if (lateralFace.boundary() && (asImp_().boundaryTypes(element, lateralFace).isSlipCondition(Indices::velocity(velIdx))))
{
unsigned int neighborIdx = neighborIndex(elementIdx, normalNormDim, 0);
if (lateralFace.center()[normalNormDim] < cellCenter(elementIdx)[normalNormDim])
neighborIdx = neighborIndex(elementIdx, normalNormDim, 1);
const SubControlVolume& scv = fvGeometry.scv(scvf.insideScvIdx());
bjsVelocityAverage[normalNormDim] += ParentType::beaversJosephVelocity(element, scv, scvf, lateralFace, ccVelocity(elementIdx, velIdx), 0.0);
if (bjsNumFaces[normalNormDim] > 0 && neighborIdx != bjsNeighbor[normalNormDim])
slipVelocityAverage[normalNormDim] += ParentType::slipVelocity(element, scv, scvf, lateralFace, velocity_[elementIdx][velIdx], 0.0);
if (slipNumFaces[normalNormDim] > 0 && neighborIdx != slipNeighbor[normalNormDim])
DUNE_THROW(Dune::InvalidStateException, "Two different neighborIdx should not occur");
bjsNeighbor[normalNormDim] = neighborIdx;
slipNeighbor[normalNormDim] = neighborIdx;
normalNormCoordinate[normalNormDim] = lateralFace.center()[normalNormDim];
bjsNumFaces[normalNormDim]++;
slipNumFaces[normalNormDim]++;
}
}
for (unsigned dirIdx = 0; dirIdx < dim; ++dirIdx)
{
if (bjsNumFaces[dirIdx] == 0)
if (slipNumFaces[dirIdx] == 0)
continue;
unsigned int neighborIdx = bjsNeighbor[dirIdx];
bjsVelocityAverage[dirIdx] /= bjsNumFaces[dirIdx];
unsigned int neighborIdx = slipNeighbor[dirIdx];
slipVelocityAverage[dirIdx] /= slipNumFaces[dirIdx];
velocityGradients_[elementIdx][velIdx][dirIdx]
= (ccVelocity(neighborIdx, velIdx) - bjsVelocityAverage[dirIdx])
/ (cellCenter(neighborIdx)[dirIdx] - normalNormCoordinate[dirIdx]);
= (velocity_[neighborIdx][velIdx] - slipVelocityAverage[dirIdx])
/ (cellCenter_[neighborIdx][dirIdx] - normalNormCoordinate[dirIdx]);
}
}
......
......@@ -78,6 +78,29 @@ public:
using ParentType::couplingPhaseIdx;
/*!
* \brief Returns the momentum flux across the coupling boundary.
*
* Calls (old/new)MomentumCouplingCondition depending on the value of the parameter "Problem.NewIc"
* Defaults to oldMomentumCouplingCondition
*
*/