From c8b68d92b4431791bd5cad832d85ac24d6bf706a Mon Sep 17 00:00:00 2001
From: Melanie Lipp <melanie.lipp@iws.uni-stuttgart.de>
Date: Fri, 26 Jun 2020 16:29:33 +0200
Subject: [PATCH] [staggered] Improve treatment of corner velocities for
advective momentum flux
In advective flux terms, there sometimes exists a corner velocity in a non-dof
position, in which a velocity is needed. This commit aims at using this existing
velocity instead of doing upwinding.
---
.../freeflow/staggeredgeometryhelper.hh | 68 ++++++
.../freeflow/subcontrolvolumeface.hh | 49 +++-
.../staggered/staggeredupwindhelper.hh | 210 +++++++++++++++++-
.../staggered/velocitygradients.hh | 20 ++
4 files changed, 341 insertions(+), 6 deletions(-)
diff --git a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
index 0a731d966f..61948a2ac8 100644
--- a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
+++ b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
@@ -39,6 +39,33 @@ namespace Detail {
/*!
* \ingroup StaggeredDiscretization
* \brief Parallel Data stored per sub face
+ *
+ * ------------
+ * | |
+ * | |
+ * | |
+ * -----------------------
+ * | yyyyyyyy s |
+ * | yyyyyyyy s |
+ * | yyyyyyyy s |
+ * -----------------------
+ * In this corner geometry, hasParallelNeighbor will return true for subcontrolvolumeface s belonging to the
+ * element filled by 'y's, but hasParallelNeighbor will return false for the subcontrolvolumeface that has the
+ * same dofIndex. We name this situation hasHalfParallelNeighbor.
+ *
+ * ------------
+ * | yyyyyyyy s
+ * | yyyyyyyy s
+ * | yyyyyyyy s
+ * -----------------------
+ * | | |
+ * | | |
+ * | | |
+ * -----------------------
+ * In this corner geometry, hasParallelNeighbor will return true for subcontrolvolumeface s belonging to the
+ * element filled by 'y's. However, as there also might be a boundary velocity value known at the corner, which
+ * can be used instead of the standard parallel velocity in some cases, we want to identify this situation. We
+ * name it cornerParallelNeighbor.
*/
template<class GridView, int upwindSchemeOrder>
struct PairData
@@ -49,6 +76,8 @@ struct PairData
using SmallLocalIndexType = typename IndexTraits<GridView>::SmallLocalIndex;
std::bitset<upwindSchemeOrder> hasParallelNeighbor;
+ bool hasHalfParallelNeighbor = false;
+ bool hasCornerParallelNeighbor = false;
std::array<GridIndexType, upwindSchemeOrder> parallelDofs;
std::array<Scalar, upwindSchemeOrder> parallelCellWidths;
bool hasOuterLateral = false;
@@ -378,6 +407,45 @@ private:
// recursively insert parallel neighbor faces into pair data
const auto parallelAxisIdx = directionIndex(intersection);
const auto localLateralIntersectionIndex = intersection.indexInInside();
+
+ /*
+ * ------------
+ * | |
+ * | |
+ * | |
+ * iiiiiiiiiii*bbbbbbbbbbb
+ * | o zzzzzzzz |
+ * | o zzzzzzzz |
+ * | o zzzzzzzz |
+ * -----------------------
+ *
+ * i:intersection,o:intersection_, b: outerIntersection, z: intersection_.outside()
+ */
+ if (intersection_.neighbor())
+ for (const auto& outerIntersection : intersections(gridView_, intersection_.outside()))
+ if (intersection.indexInInside() == outerIntersection.indexInInside())
+ if (!outerIntersection.neighbor())
+ pairData_[numPairParallelIdx].hasHalfParallelNeighbor = true;
+
+ /* ------------
+ * | o
+ * | o
+ * | o
+ * iiiiiiiiiii------------
+ * | zzzzzzzz b |
+ * | zzzzzzzz b |
+ * | zzzzzzzz b |
+ * -----------------------
+ *
+ * i:intersection,o:intersection_, b: outerIntersection, z: intersection.outside()
+ */
+ if (!intersection_.neighbor())
+ for (const auto& outerIntersection : intersections(gridView_, intersection.outside()))
+ if (intersection_.indexInInside() == outerIntersection.indexInInside())
+ if (outerIntersection.neighbor())
+ pairData_[numPairParallelIdx].hasCornerParallelNeighbor = true;
+
+
addParallelNeighborPairData_(intersection.outside(), 0, localLateralIntersectionIndex, parallelLocalIdx, parallelAxisIdx, numPairParallelIdx);
}
diff --git a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
index 29de016a9e..9805aa626a 100644
--- a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
@@ -318,6 +318,53 @@ public:
return pairData(localSubFaceIdx).hasParallelNeighbor[parallelDegreeIdx];
}
+ /*!
+ * \brief Check if the face has a half parallel neighbor
+ *
+ * \param localSubFaceIdx The local index of the subface
+ *
+ * ------------
+ * | |
+ * | |
+ * | |
+ * -----------------------
+ * | yyyyyyyy s |
+ * | yyyyyyyy s |
+ * | yyyyyyyy s |
+ * -----------------------
+ * In this corner geometry, hasParallelNeighbor will return true for subcontrolvolumeface s belonging to the
+ * element filled by 'y's, but hasParallelNeighbor will return false for the subcontrolvolumeface that has the
+ * same dofIndex. We name this situation hasHalfParallelNeighbor.
+ */
+ bool hasHalfParallelNeighbor(const int localSubFaceIdx) const
+ {
+ return pairData(localSubFaceIdx).hasHalfParallelNeighbor;
+ }
+
+ /*!
+ * \brief Check if the face has a corner parallel neighbor
+ *
+ * \param localSubFaceIdx The local index of the subface
+ *
+ * ------------
+ * | yyyyyyyy s
+ * | yyyyyyyy s
+ * | yyyyyyyy s
+ * -----------------------
+ * | | |
+ * | | |
+ * | | |
+ * -----------------------
+ * In this corner geometry, hasParallelNeighbor will return true for subcontrolvolumeface s belonging to the
+ * element filled by 'y's. However, as there also might be a boundary velocity value known at the corner, which
+ * can be used instead of the standard parallel velocity in some cases, we want to identify this situation. We
+ * name it cornerParallelNeighbor.
+ */
+ bool hasCornerParallelNeighbor(const int localSubFaceIdx) const
+ {
+ return pairData(localSubFaceIdx).hasCornerParallelNeighbor;
+ }
+
/*!
* \brief Check if the face has an outer normal neighbor
*
@@ -402,7 +449,7 @@ public:
void setCenter(const GlobalPosition& center)
{ center_ = center; }
- //! set the boudnary flag
+ //! set the boundary flag
void setBoundary(bool boundaryFlag)
{ boundary_ = boundaryFlag; }
diff --git a/dumux/freeflow/navierstokes/staggered/staggeredupwindhelper.hh b/dumux/freeflow/navierstokes/staggered/staggeredupwindhelper.hh
index 27f6961464..9d1a9996bd 100644
--- a/dumux/freeflow/navierstokes/staggered/staggeredupwindhelper.hh
+++ b/dumux/freeflow/navierstokes/staggered/staggeredupwindhelper.hh
@@ -258,6 +258,40 @@ private:
* TODO: In order to get a second order momentum upwind scheme for compressible flow the densities have to be evaluated
* at the same integration points / positions as the velocities. The currently implementation just takes the closest upwind density
* to compute the momentum as a crude approximation.
+ *
+ * ------------
+ * | xxxx o
+ * | xxxx o
+ * | xxxx o
+ * -----------*bbbbbbbbbbb
+ * | yyyy o zzzz |
+ * | yyyy o zzzz |
+ * | yyyy o zzzz |
+ * -----------------------
+ * If scvf_ is touching a corner, at which there is a Dirichlet condition for face b (half-control
+ * volumes x, y and z), the transported velocity at * is given. No upwinding or
+ * higher-order approximation for the velocity at * are required. This also means the transported
+ * velocity at * is the same for the half-control volumes y and z.
+ *
+ * ------------
+ * | |
+ * | |
+ * | |
+ * -----------------------
+ * | xxxx o wwww |
+ * | xxxx o wwww |
+ * | xxxx o wwww |
+ * ------+++++*~~~~~------
+ * | yyyy o zzzz |
+ * | yyyy o zzzz |
+ * | yyyy o zzzz |
+ * -----------------------
+ * If the flux over face + is calculated and a corner occurs a bit further away (here upper right),
+ * no special treatment of the corner geometry is provided. This means, the transported velocity at
+ * star (*) is obtained with an upwind scheme. In particularly, this also means
+ * that while x and y see the same velocity * for the flux over face x (continuity OK), z sees
+ * another velocity * for the flux over face ~ (continuity still OK, as only z and w have to use the
+ * same velocity at *).
*/
std::array<Scalar, 3> getLateralSecondOrderUpwindMomenta_(const Scalar insideDensity,
const Scalar outsideDensity,
@@ -270,11 +304,22 @@ private:
// If the lateral face lies on a boundary, we assume that the parallel velocity on the boundary is actually known,
// thus we always use this value for the computation of the transported momentum.
- if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0))
+ if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0) || scvf_.hasHalfParallelNeighbor(localSubFaceIdx) || scvf_.hasCornerParallelNeighbor(localSubFaceIdx))
{
- const Scalar boundaryVelocity = getParallelVelocityFromBoundary_(element_, scvf_, faceVars_, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
- const Scalar boundaryMomentum = boundaryVelocity*insideDensity;
- return {boundaryMomentum, boundaryMomentum, boundaryMomentum};
+ if ((scvf_.hasHalfParallelNeighbor(localSubFaceIdx) || scvf_.hasCornerParallelNeighbor(localSubFaceIdx)) && dirichletParallelNeighbor_(localSubFaceIdx))
+ {
+ const Scalar boundaryVelocity = getParallelVelocityFromCorner_(localSubFaceIdx);
+ const Scalar boundaryMomentum = boundaryVelocity*insideDensity;
+
+ return {boundaryMomentum, boundaryMomentum, boundaryMomentum};
+ }
+ else if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0))
+ {
+ const Scalar boundaryVelocity = getParallelVelocityFromBoundary_(element_, scvf_, faceVars_, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
+ const Scalar boundaryMomentum = boundaryVelocity*insideDensity;
+
+ return {boundaryMomentum, boundaryMomentum, boundaryMomentum};
+ }
}
if (selfIsUpstream)
@@ -327,7 +372,13 @@ private:
{
// If the lateral face lies on a boundary, we assume that the parallel velocity on the boundary is actually known,
// thus we always use this value for the computation of the transported momentum.
- if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0))
+ if ((scvf_.hasHalfParallelNeighbor(localSubFaceIdx) || scvf_.hasCornerParallelNeighbor(localSubFaceIdx)) && dirichletParallelNeighbor_(localSubFaceIdx))
+ {
+ const Scalar boundaryVelocity = getParallelVelocityFromCorner_(localSubFaceIdx);
+ const Scalar boundaryMomentum = boundaryVelocity*outsideDensity;
+ return {boundaryMomentum, boundaryMomentum};
+ }
+ else if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0))
{
const Scalar boundaryVelocity = getParallelVelocityFromBoundary_(element_, scvf_, faceVars_, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
const Scalar boundaryMomentum = boundaryVelocity*insideDensity;
@@ -450,6 +501,155 @@ private:
return getParallelVelocityFromBoundary_(element, scvf, faceVars, currentScvfBoundaryTypes, lateralOppositeFaceBoundaryTypes, localOppositeSubFaceIdx);
}
+ /*!
+ * \brief Returns the boundary subcontrolvolumeface in a corner geometry.
+ *
+ * ------------
+ * | xxxx o
+ * | xxxx o
+ * | xxxx o
+ * ------------bbbbbbbbbbb
+ * | yyyy o |
+ * | yyyy o |
+ * | yyyy o |
+ * -----------------------
+ *
+ * This function will be entered in such a corner geometry (there is no cell in the upper right, --- and |
+ * stand for the grid cells). The scvf_ will be one of the two ones denoted by o (upper one
+ * hasCornerParallelNeighbor, lower one hasHalfParallelNeighbor). x and y are the two possible corresponding
+ * half-control volumes. In both cases, the returned boundaryScvf is the one marked by b. It needs to be the
+ * same boundaryScvf returned for the sake of flux continuity.
+ */
+ const SubControlVolumeFace& boundaryScvf_(const int localSubFaceIdx) const
+ {
+ if (scvf_.hasHalfParallelNeighbor(localSubFaceIdx))
+ {
+ return fvGeometry_.scvf(scvf_.outsideScvIdx(), scvf_.pairData(localSubFaceIdx).localLateralFaceIdx);
+ }
+ else if (scvf_.hasCornerParallelNeighbor(localSubFaceIdx))
+ {
+ /*
+ * ------------
+ * | xxxxxxxx o
+ * | xxxxxxxx o
+ * | xxxxxxxx o
+ * lllllllllll-bbbbbbbbbbb
+ * | yyyyyyyy p |
+ * | yyyyyyyy p |
+ * | yyyyyyyy p |
+ * -----------------------
+ *
+ * o: scvf_, l: lateralFace, p: parallelFace, b: returned scvf, x: scvf_ inside scv, y: lateralFace
+ * outside scv
+ */
+ const SubControlVolumeFace& lateralFace = fvGeometry_.scvf(scvf_.insideScvIdx(), scvf_.pairData(localSubFaceIdx).localLateralFaceIdx);
+ const SubControlVolumeFace& parallelFace = fvGeometry_.scvf(lateralFace.outsideScvIdx(), scvf_.localFaceIdx());
+
+ const auto& localLateralIdx = scvf_.pairData(localSubFaceIdx).localLateralFaceIdx;
+ const auto& localLateralOppositeIdx = (localLateralIdx % 2) ? (localLateralIdx - 1) : (localLateralIdx + 1);
+
+ return fvGeometry_.scvf(parallelFace.outsideScvIdx(), localLateralOppositeIdx);
+ }
+ else
+ {
+ DUNE_THROW(Dune::InvalidStateException, "The function boundaryScvf_ should only be called when hasHalfParallelNeighbor or hasCornerParallelNeighbor.");
+ }
+ }
+
+ /*!
+ * \brief Gets the boundary element in a corner geometry.
+ *
+ * ------------
+ * | xxxx o
+ * | xxxx o
+ * | xxxx o
+ * -----------------------
+ * | yyyy o bbbbbbbb |
+ * | yyyy o bbbbbbbb |
+ * | yyyy o bbbbbbbb |
+ * -----------------------
+ *
+ * This function will be entered in such a corner geometry (there is no cell in the upper right, --- and |
+ * stand for the grid cells). The scvf_ will be one of the two ones denoted by o (upper one
+ * hasCornerParallelNeighbor, lower one hasHalfParallelNeighbor). x and y are the two possible corresponding
+ * half-control volumes. In both cases, the returned boundaryElement is the one marked by b. It needs to be
+ * the same boundaryScvf returned for the sake of flux continuity.
+ */
+ Element boundaryElement_(const int localSubFaceIdx) const
+ {
+ if (scvf_.hasHalfParallelNeighbor(localSubFaceIdx))
+ {
+ return fvGeometry_.gridGeometry().element(scvf_.outsideScvIdx());
+ }
+ else if (scvf_.hasCornerParallelNeighbor(localSubFaceIdx))
+ {
+ const SubControlVolumeFace& lateralFace = fvGeometry_.scvf(scvf_.insideScvIdx(), scvf_.pairData(localSubFaceIdx).localLateralFaceIdx);
+ const SubControlVolumeFace& parallelFace = fvGeometry_.scvf(lateralFace.outsideScvIdx(), scvf_.localFaceIdx());
+
+ return fvGeometry_.gridGeometry().element(parallelFace.outsideScvIdx());
+ }
+ else
+ {
+ DUNE_THROW(Dune::InvalidStateException, "When entering boundaryElement_ scvf_ should have either hasHalfParallelNeighbor or hasCornerParallelNeighbor true. Not the case here.");
+ }
+ }
+
+ /*!
+ * \brief Sets the bools hasDirichletCornerParallelNeighbor and hasDirichletHalfParallelNeighbor.
+ *
+ * ------------
+ * | xxxx o
+ * | xxxx o
+ * | xxxx o
+ * ------------bbbbbbbbbbb
+ * | yyyy o |
+ * | yyyy o boundary |
+ * | yyyy o element |
+ * -----------------------
+ *
+ * This function will be entered in such a corner geometry (there is no cell in the upper right, --- and |
+ * stand for the grid cells). The scvf_ will be one of the two ones denoted by o (upper one
+ * hasCornerParallelNeighbor, lower one hasHalfParallelNeighbor). x and y are the two possible corresponding
+ * half-control volumes. In both cases, we check if the face bbb, part of the edge of element boundaryElement,
+ * is a Dirichlet boundary.
+ */
+ bool dirichletParallelNeighbor_(const int localSubFaceIdx) const
+ {
+ const auto& problem = elemVolVars_.gridVolVars().problem();
+ const Element& boundaryElement = boundaryElement_(localSubFaceIdx);
+ const SubControlVolumeFace& boundaryScvf = boundaryScvf_(localSubFaceIdx);
+
+ return problem.boundaryTypes(boundaryElement, boundaryScvf).isDirichlet(Indices::velocity(scvf_.directionIndex()));
+ }
+
+ /*!
+ * \brief Gets the parallel velocity from a corner geometry.
+ *
+ * ------------
+ * | xxxx o
+ * | xxxx o
+ * | xxxx o
+ * -----------*-----------
+ * | yyyy o |
+ * | yyyy o |
+ * | yyyy o |
+ * -----------------------
+ *
+ * This function will be entered in such a corner geometry (there is no cell in the upper right, --- and |
+ * stand for the grid cells). The scvf_ will be one of the two ones denoted by o (upper one
+ * hasCornerParallelNeighbor, lower one hasHalfParallelNeighbor). x and y are the two possible corresponding
+ * half-control volumes. In both cases, the returned velocity is situated in the corner (*).
+ */
+ Scalar getParallelVelocityFromCorner_(const int localSubFaceIdx) const
+ {
+ const auto& problem = elemVolVars_.gridVolVars().problem();
+ const Element& boundaryElement = boundaryElement_(localSubFaceIdx);
+ const SubControlVolumeFace& boundaryScvf = boundaryScvf_(localSubFaceIdx);
+ const auto ghostFace = makeStaggeredBoundaryFace(boundaryScvf, scvf_.pairData(localSubFaceIdx).lateralStaggeredFaceCenter);
+
+ return problem.dirichlet(boundaryElement, ghostFace)[Indices::velocity(scvf_.directionIndex())];
+ }
+
const Element& element_;
const FVElementGeometry& fvGeometry_;
const SubControlVolumeFace& scvf_;
diff --git a/dumux/freeflow/navierstokes/staggered/velocitygradients.hh b/dumux/freeflow/navierstokes/staggered/velocitygradients.hh
index b2c9d33a1c..932e1687ac 100644
--- a/dumux/freeflow/navierstokes/staggered/velocitygradients.hh
+++ b/dumux/freeflow/navierstokes/staggered/velocitygradients.hh
@@ -93,6 +93,26 @@ public:
*
* O position at which gradient is evaluated
* \endverbatim
+ *
+ * ------------
+ * | xxxx s
+ * | xxxx a
+ * | xxxx s
+ * -----------O-----------
+ * | yyyy s zzzz |
+ * | yyyy b zzzz |
+ * | yyyy s zzzz |
+ * -----------------------
+ *
+ * In a corner geometry (scvf is sas or sbs), we calculate the velocity gradient at O, by
+ * (velocity(a)-velocity(b))/distance(a,b) for the half-control volumes x and y, but by
+ * (velocity(O)-velocity(b))/distance(O,b) for z. This does not harm flux continuity (x and y use the same
+ * formulation). We do this different formulation for y (approximate gradient by central differncing) and
+ * z (approximate gradient by forward/backward differencing), because it is the natural way of implementing
+ * it and it is not clear which gradient is the better approximation in this case anyway.
+ * Particularly, for the frequent case of no-slip, no-flow boundaries, the velocity would be zero at O and a
+ * and thus, the gradient within the flow domain might be better approximated by velocity(b)/distanc(O,b)
+ * than by velocity(b)/distance(a,b).
*/
template<class Problem, class FaceVariables>
static Scalar velocityGradIJ(const Problem& problem,
--
GitLab