diff --git a/dumux/discretization/staggered/freeflow/connectivitymap.hh b/dumux/discretization/staggered/freeflow/connectivitymap.hh
index 313874e84121d9bf907f3f820f4a566918899a0e..fc0b40e3558d4e9f7b8887dfebc0425dd2877d20 100644
--- a/dumux/discretization/staggered/freeflow/connectivitymap.hh
+++ b/dumux/discretization/staggered/freeflow/connectivitymap.hh
@@ -143,7 +143,7 @@ private:
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
- stencil.push_back(scvf.dofIndex());
+ stencil.push_back(scvf.axisData().selfDof);
}
/*
@@ -163,8 +163,8 @@ private:
auto& normalFace = fvGeometry.scvf(eIdx, data.localNormalFaceIdx);
if(!normalFace.boundary())
{
- const auto outerParallelElementDofIdx = normalFace.outsideScvIdx();
- stencil.push_back(outerParallelElementDofIdx);
+ const auto firstParallelElementDofIdx = normalFace.outsideScvIdx();
+ stencil.push_back(firstParallelElementDofIdx);
}
}
}
@@ -177,21 +177,43 @@ private:
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
- // the first entries are always the face dofIdx itself and the one of the opposing face
if(stencil.empty())
{
- stencil.push_back(scvf.dofIndex());
- stencil.push_back(scvf.dofIndexOpposingFace());
+ for(int i = 0; i < scvf.axisData().inAxisBackwardDofs.size(); i++)
+ {
+ if(scvf.hasBackwardNeighbor(i))
+ {
+ stencil.push_back(scvf.axisData().inAxisBackwardDofs[i]);
+ }
+ }
+
+ stencil.push_back(scvf.axisData().selfDof);
+ stencil.push_back(scvf.axisData().oppositeDof);
+
+ for(int i = 0; i < scvf.axisData().inAxisForwardDofs.size(); i++)
+ {
+ if(scvf.hasForwardNeighbor(i))
+ {
+ stencil.push_back(scvf.axisData().inAxisForwardDofs[i]);
+ }
+ }
}
for(const auto& data : scvf.pairData())
{
+ // add normal dofs
stencil.push_back(data.normalPair.first);
- const auto outerParallelFaceDofIdx = data.outerParallelFaceDofIdx;
- if(outerParallelFaceDofIdx >= 0)
- stencil.push_back(outerParallelFaceDofIdx);
if(!scvf.boundary())
stencil.push_back(data.normalPair.second);
+
+ // add parallel dofs
+ for (int i = 0; i < data.parallelDofs.size(); i++)
+ {
+ if(!(data.parallelDofs[i] < 0))
+ {
+ stencil.push_back(data.parallelDofs[i]);
+ }
+ }
}
}
diff --git a/dumux/discretization/staggered/freeflow/facevariables.hh b/dumux/discretization/staggered/freeflow/facevariables.hh
index 6fd24dc689105446baddd6fb1b78c8d39db3f97c..a4edc1e42683e8f7b05d23f66d6746c89b38bed3 100644
--- a/dumux/discretization/staggered/freeflow/facevariables.hh
+++ b/dumux/discretization/staggered/freeflow/facevariables.hh
@@ -25,6 +25,7 @@
#define DUMUX_DISCRETIZATION_STAGGERED_FREEFLOW_FACEVARIABLES_HH
#include <array>
+#include <vector>
namespace Dumux {
@@ -72,7 +73,30 @@ public:
velocitySelf_ = faceSol[scvf.dofIndex()];
velocityOpposite_ = faceSol[scvf.dofIndexOpposingFace()];
+ // treat the velocity forward of the self face i.e. the face that is
+ // forward wrt the self face by degree i
+ velocityForward_.clear();
+ for (int i = 0; i < scvf.axisData().inAxisForwardDofs.size(); i++)
+ {
+ if(!(scvf.axisData().inAxisForwardDofs[i] < 0))
+ {
+ velocityForward_.push_back(faceSol[scvf.axisData().inAxisForwardDofs[i]]);
+ }
+ }
+
+ // treat the velocity at the first backward face i.e. the face that is
+ // behind the opposite face by degree i
+ velocityBackward_.clear();
+ for (int i = 0; i < scvf.axisData().inAxisBackwardDofs.size(); i++)
+ {
+ if(!(scvf.axisData().inAxisBackwardDofs[i] < 0))
+ {
+ velocityBackward_.push_back(faceSol[scvf.axisData().inAxisBackwardDofs[i]]);
+ }
+ }
+
// handle all sub faces
+ std::vector<Scalar> a;
for(int i = 0; i < scvf.pairData().size(); ++i)
{
const auto& subFaceData = scvf.pairData(i);
@@ -80,12 +104,19 @@ public:
// treat the velocities normal to the face
velocityNormalInside_[i] = faceSol[subFaceData.normalPair.first];
- if(scvf.hasFrontalNeighbor(i))
+ if(scvf.hasOuterNormal(i))
velocityNormalOutside_[i] = faceSol[subFaceData.normalPair.second];
- // treat the velocity parallel to the face
- if(scvf.hasParallelNeighbor(i))
- velocityParallel_[i] = faceSol[subFaceData.outerParallelFaceDofIdx];
+ a.clear();
+ // treat the velocities parallel to the self face
+ for(int j = 0; j < subFaceData.parallelDofs.size(); j++)
+ {
+ if(scvf.hasParallelNeighbor(i,j))
+ {
+ a.push_back(faceSol[subFaceData.parallelDofs[j]]);
+ }
+ }
+ velocityParallel_[i] = a;
}
}
@@ -106,13 +137,34 @@ public:
}
/*!
- * \brief Returns the velocity at the parallel face
+ * \brief Returns the velocity at a backward face
+ *
+ * \param backwardIdx The index describing how many faces backward this dof is from the opposite face
+ */
+ Scalar velocityBackward(const int backwardIdx) const
+ {
+ return velocityBackward_[backwardIdx];
+ }
+
+ /*!
+ * \brief Returns the velocity at a forward face
+ *
+ * \param forwardIdx The index describing how many faces forward this dof is of the self face
+ */
+ Scalar velocityForward(const int forwardIdx) const
+ {
+ return velocityForward_[forwardIdx];
+ }
+
+ /*!
+ * \brief Returns the velocity at a parallel face
*
* \param localSubFaceIdx The local index of the subface
+ * \param parallelDegreeIdx The index describing how many faces parallel this dof is of the parallel face
*/
- Scalar velocityParallel(const int localSubFaceIdx) const
+ Scalar velocityParallel(const int localSubFaceIdx, const int parallelDegreeIdx) const
{
- return velocityParallel_[localSubFaceIdx];
+ return velocityParallel_[localSubFaceIdx][parallelDegreeIdx];
}
/*!
@@ -136,12 +188,14 @@ public:
}
private:
-
Scalar velocitySelf_;
Scalar velocityOpposite_;
- std::array<Scalar, numPairs> velocityParallel_;
- std::array<Scalar, numPairs> velocityNormalInside_;
- std::array<Scalar, numPairs> velocityNormalOutside_;
+ std::vector<Scalar> velocityForward_;
+ std::vector<Scalar> velocityBackward_;
+ std::array<std::vector<Scalar>, numPairs> velocityParallel_;
+ std::array<Scalar, numPairs> velocityNormalInside_;
+ std::array<Scalar, numPairs> velocityNormalOutside_;
+
};
} // end namespace Dumux
diff --git a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
index 7760f173a222308fc2689da2bac8229a316aa0e3..7d87aa882574026435e29e029e2ac2a5a994253d 100644
--- a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
+++ b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
@@ -31,25 +31,41 @@
#include <type_traits>
#include <algorithm>
-namespace Dumux {
+namespace Dumux
+{
/*!
* \ingroup StaggeredDiscretization
- * \brief Data stored per sub face
+ * \brief Parallel Data stored per sub face
*/
template<class Scalar, class GlobalPosition>
struct PairData
{
- int outerParallelFaceDofIdx;
- std::pair<int,int> normalPair;
+ std::vector<int> parallelDofs;
+ std::vector<Scalar> parallelDistances;
+ std::pair<signed int,signed int> normalPair;
int localNormalFaceIdx;
- int globalCommonEntIdx;
- Scalar parallelDistance;
Scalar normalDistance;
- GlobalPosition virtualOuterParallelFaceDofPos;
+ GlobalPosition virtualFirstParallelFaceDofPos;
+};
+
+/*!
+ * \ingroup StaggeredDiscretization
+ * \brief In Axis Data stored per sub face
+ */
+template<class Scalar>
+struct AxisData
+{
+ int selfDof;
+ int oppositeDof;
+ std::vector<int> inAxisForwardDofs;
+ std::vector<int> inAxisBackwardDofs;
+ Scalar selfToOppositeDistance;
+ std::vector<Scalar> inAxisForwardDistances;
+ std::vector<Scalar> inAxisBackwardDistances;
};
- /*!
+/*!
* \ingroup StaggeredDiscretization
* \brief Returns the dirction index of the facet (0 = x, 1 = y, 2 = z)
*/
@@ -65,11 +81,6 @@ inline static unsigned int directionIndex(Vector&& vector)
* \ingroup StaggeredDiscretization
* \brief Helper class constructing the dual grid finite volume geometries
* for the free flow staggered discretization method.
- *
- * A visualization of the variables that are in this document can be found in the following image:
- *
- * \image html staggeredgeometry.png
- *
*/
template<class GridView>
class FreeFlowStaggeredGeometryHelper
@@ -92,7 +103,9 @@ class FreeFlowStaggeredGeometryHelper
static constexpr int codimIntersection = 1;
static constexpr int codimCommonEntity = 2;
static constexpr int numFacetSubEntities = (dim == 2) ? 2 : 4;
+ static constexpr int numfacets = dimWorld * 2;
static constexpr int numPairs = 2 * (dimWorld - 1);
+ static int order_;
public:
@@ -106,95 +119,247 @@ public:
intersection_ = intersection;
innerNormalFacePos_.clear();
fillPairData_();
+ fillAxisData_();
}
/*!
- * \brief Returns the global dofIdx of the intersection itself
+ * \brief Returns the local index of the face (i.e. the intersection)
*/
- int dofIndex() const
+ int localFaceIndex() const
{
- //TODO: use proper intersection mapper!
- const auto inIdx = intersection_.indexInInside();
- return gridView_.indexSet().subIndex(intersection_.inside(), inIdx, codimIntersection);
+ return intersection_.indexInInside();
}
/*!
- * \brief Returns the global dofIdx of the opposing intersection
+ * \brief Returns a copy of the axis data
*/
- int dofIndexOpposingFace() const
+ auto axisData() const
{
- //TODO: use proper intersection mapper!
- const auto inIdx = intersection_.indexInInside();
- return gridView_.indexSet().subIndex(this->intersection_.inside(), localOppositeIdx_(inIdx), codimIntersection);
+ return axisData_;
}
/*!
- * \brief Returns the local index of the face (i.e. the intersection)
+ * \brief Returns a copy of the pair data
*/
- int localFaceIndex() const
+ auto pairData() const
{
- return intersection_.indexInInside();
+ return pairData_;
}
- /*!
- * \brief Returns the distance between dofSelf and dofOpposite
+ /*!
+ * \brief Returns the dirction index of the primary facet (0 = x, 1 = y, 2 = z)
*/
- Scalar selfToOppositeDistance() const
+ int directionIndex() const
{
- const auto inIdx = intersection_.indexInInside();
- const auto self = getFacet_(inIdx, element_);
- const auto opposite = getFacet_(localOppositeIdx_(inIdx), element_);
- return (self.geometry().center() - opposite.geometry().center()).two_norm();
+ return Dumux::directionIndex(std::move(intersection_.centerUnitOuterNormal()));
}
/*!
- * \brief Returns a copy of the pair data
+ * \brief Returns the dirction index of the facet passed as an argument (0 = x, 1 = y, 2 = z)
*/
- auto pairData() const
+ int directionIndex(const Intersection& intersection) const
{
- return pairData_;
+ return Dumux::directionIndex(std::move(intersection.centerUnitOuterNormal()));
}
- /*!
- * \brief Returns the dirction index of the facet (0 = x, 1 = y, 2 = z)
- */
- int directionIndex() const
+ //! \brief Returns the order needed by the scheme
+ static int order()
{
- return Dumux::directionIndex(std::move(intersection_.centerUnitOuterNormal()));
+ return order_;
+ }
+
+ //! \brief Set the order needed by the scheme
+ static void setOrder(const int order)
+ {
+ order_ = order;
}
private:
- /*!
+ /*!
+ * \brief Fills all entries of the in axis data
+ */
+ void fillAxisData_()
+ {
+ unsigned int numForwardBackwardAxisDofs = order_ - 1;
+ const auto inIdx = intersection_.indexInInside();
+ const auto oppIdx = localOppositeIdx_(inIdx);
+
+ // Clear the containers before filling them
+ this->axisData_.inAxisForwardDofs.clear();
+ this->axisData_.inAxisBackwardDofs.clear();
+ this->axisData_.inAxisForwardDistances.clear();
+ this->axisData_.inAxisBackwardDistances.clear();
+
+ // initialize values that could remain unitialized if the intersection lies on a boundary
+ for(int i = 0; i < numForwardBackwardAxisDofs; i++)
+ {
+ this->axisData_.inAxisForwardDofs.push_back(-1);
+ this->axisData_.inAxisBackwardDofs.push_back(-1);
+ this->axisData_.inAxisForwardDistances.push_back(0);
+ this->axisData_.inAxisBackwardDistances.push_back(0);
+ }
+
+ // Set the self Dof
+ this->axisData_.selfDof = gridView_.indexSet().subIndex(this->intersection_.inside(), inIdx, codimIntersection);
+ const auto self = getFacet_(inIdx, element_);
+
+ // Set the opposite Dof
+ this->axisData_.oppositeDof = gridView_.indexSet().subIndex(this->intersection_.inside(), oppIdx, codimIntersection);
+ const auto opposite = getFacet_(oppIdx, element_);
+
+ // Set the Self to Opposite Distance
+ this->axisData_.selfToOppositeDistance = (self.geometry().center() - opposite.geometry().center()).two_norm();
+
+ // Set the Forward Dofs
+ std::stack<Element> inAxisForwardElementStack;
+ auto selfFace = getFacet_(inIdx, element_);
+ if(intersection_.neighbor())
+ {
+ if (order_ > 1)
+ inAxisForwardElementStack.push(intersection_.outside());
+ bool keepStackingForward = (inAxisForwardElementStack.size() < order_-1);
+ while(keepStackingForward)
+ {
+ auto e = inAxisForwardElementStack.top();
+ for(const auto& intersection : intersections(gridView_,e))
+ {
+ if( (intersection.indexInInside() == inIdx ) )
+ {
+ if( intersection.neighbor())
+ {
+ inAxisForwardElementStack.push(intersection.outside());
+ keepStackingForward = (inAxisForwardElementStack.size() < order_-1);
+ }
+ else
+ {
+ keepStackingForward = false;
+ }
+ }
+ }
+ }
+ }
+ std::vector<GlobalPosition> forwardFaceCoordinates(inAxisForwardElementStack.size(), selfFace.geometry().center());
+ while(!inAxisForwardElementStack.empty())
+ {
+ int forwardIdx = inAxisForwardElementStack.size()-1;
+ this->axisData_.inAxisForwardDofs[forwardIdx] = gridView_.indexSet().subIndex(inAxisForwardElementStack.top(), inIdx, codimIntersection);
+ selfFace = getFacet_(inIdx, inAxisForwardElementStack.top());
+ forwardFaceCoordinates[forwardIdx] = selfFace.geometry().center();
+ inAxisForwardElementStack.pop();
+ }
+ for(int i = 0; i< forwardFaceCoordinates.size(); i++)
+ {
+ if(i == 0)
+ {
+ this->axisData_.inAxisForwardDistances[i] = (forwardFaceCoordinates[i] - self.geometry().center()).two_norm();
+ }
+ else
+ {
+ this->axisData_.inAxisForwardDistances[i] = (forwardFaceCoordinates[i]- forwardFaceCoordinates[i-1]).two_norm();
+ }
+ }
+
+ // Set the Backward Dofs
+ std::stack<Element> inAxisBackwardElementStack;
+ auto oppFace = getFacet_(oppIdx, element_);
+ for(const auto& intersection : intersections(gridView_, element_))
+ {
+ if(intersection.indexInInside() == oppIdx)
+ {
+ if(intersection.neighbor())
+ {
+ if (order_ > 1)
+ inAxisBackwardElementStack.push(intersection.outside());
+ bool keepStackingBackward = (inAxisBackwardElementStack.size() < order_-1);
+ while(keepStackingBackward)
+ {
+ auto e = inAxisBackwardElementStack.top();
+ for(const auto& intersectionOut : intersections(gridView_,e))
+ {
+
+ if( (intersectionOut.indexInInside() == oppIdx ) )
+ {
+ if( intersectionOut.neighbor())
+ {
+ inAxisBackwardElementStack.push(intersectionOut.outside());
+ keepStackingBackward = (inAxisBackwardElementStack.size() < order_-1);
+ }
+ else
+ {
+ keepStackingBackward = false;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ std::vector<GlobalPosition> backwardFaceCoordinates(inAxisBackwardElementStack.size(), oppFace.geometry().center());
+ while(!inAxisBackwardElementStack.empty())
+ {
+ int backwardIdx = inAxisBackwardElementStack.size()-1;
+ this->axisData_.inAxisBackwardDofs[backwardIdx] = gridView_.indexSet().subIndex(inAxisBackwardElementStack.top(), oppIdx, codimIntersection);
+ oppFace = getFacet_(oppIdx, inAxisBackwardElementStack.top());
+ backwardFaceCoordinates[backwardIdx] = oppFace.geometry().center();
+ inAxisBackwardElementStack.pop();
+ }
+ for(int i = 0; i< backwardFaceCoordinates.size(); i++)
+ {
+ if(i == 0)
+ {
+ this->axisData_.inAxisBackwardDistances[i] = (backwardFaceCoordinates[i] - opposite.geometry().center()).two_norm();
+ }
+ else
+ {
+ this->axisData_.inAxisBackwardDistances[i] = (backwardFaceCoordinates[i] - backwardFaceCoordinates[i-1]).two_norm();
+ }
+ }
+ }
+
+ /*!
* \brief Fills all entries of the pair data
*/
void fillPairData_()
{
- const int localFacetIdx = intersection_.indexInInside();
-
// initialize values that could remain unitialized if the intersection lies on a boundary
for(auto& data : pairData_)
{
- data.outerParallelFaceDofIdx = -1;
- data.normalPair.second = -1;
- data.normalDistance = -1;
- data.parallelDistance = -1;
- }
+ int numParallelDofs = order_;
+
+ // parallel Dofs
+ data.parallelDofs.clear();
+ data.parallelDofs.resize(numParallelDofs, -1);
- // set the inner parts of the normal pairs
- const auto localIndices = getLocalIndices_(localFacetIdx);
- setInnerIndices_(localIndices);
+ // parallel Distances
+ data.parallelDistances.clear();
+ data.parallelDistances.resize(numParallelDofs + 1, 0.0);
- // get the positions of the faces normal to the intersection within the element
- innerNormalFacePos_.reserve(numPairs);
+ // outer normals
+ data.normalPair.second = -1;
+ }
+ // set basic global positions
+ const auto& elementCenter = this->element_.geometry().center();
+ const auto& FacetCenter = intersection_.geometry().center();
- for(int i = 0; i < numPairs; ++i)
+ // get the inner normal Dof Index
+ int numPairInnerNormalIdx = 0;
+ for(const auto& innerElementIntersection : intersections(gridView_, element_))
{
- const auto& innerNormalFacet = getFacet_(localIndices.normalLocalDofIdx[i], element_);
- innerNormalFacePos_.emplace_back(innerNormalFacet.geometry().center());
+ if(facetIsNormal_(innerElementIntersection.indexInInside(), intersection_.indexInInside()))
+ {
+ const int innerElementIntersectionIdx = innerElementIntersection.indexInInside();
+ setNormalPairFirstInfo_(innerElementIntersectionIdx, element_, numPairInnerNormalIdx);
+ numPairInnerNormalIdx++;
+
+ innerNormalFacePos_.reserve(numPairs);
+ const auto& innerNormalFacet = getFacet_(innerElementIntersection.indexInInside(), element_);
+ innerNormalFacePos_.emplace_back(innerNormalFacet.geometry().center());
+ }
}
- // go into the direct neighbor element
+ // get the outer normal Dof Index
+ int numPairOuterNormalIdx = 0;
if(intersection_.neighbor())
{
// the direct neighbor element and the respective intersection index
@@ -202,60 +367,92 @@ private:
for(const auto& directNeighborElementIntersection : intersections(gridView_, directNeighborElement))
{
-
- const int directNeighborElemIsIdx = directNeighborElementIntersection.indexInInside();
// skip the directly neighboring face itself and its opposing one
- if(facetIsNormal_(directNeighborElemIsIdx, intersection_.indexInOutside()))
+ if(facetIsNormal_(directNeighborElementIntersection.indexInInside(), intersection_.indexInOutside()))
{
- setPairInfo_(directNeighborElemIsIdx, directNeighborElement, false);
-
- // go into the adjacent neighbor element to get parallel dof info
- if(directNeighborElementIntersection.neighbor())
- {
- const auto& diagonalNeighborElement = directNeighborElementIntersection.outside();
- for(const auto& dIs : intersections(gridView_, diagonalNeighborElement))
- {
- if(facetIsNormal_(dIs.indexInInside(), directNeighborElementIntersection.indexInOutside()))
- setPairInfo_(dIs.indexInInside(), diagonalNeighborElement, true);
- }
- }
+ const int directNeighborElemIsIdx = directNeighborElementIntersection.indexInInside();
+ setNormalPairSecondInfo_(directNeighborElemIsIdx, directNeighborElement, numPairOuterNormalIdx);
+ numPairOuterNormalIdx++;
}
}
}
else // intersection is on boundary
{
- // find an intersection normal to the face
- for(const auto& normalIntersection : intersections(gridView_, element_))
+ // fill the normal pair entries
+ for(int pairIdx = 0; pairIdx < numPairs; ++pairIdx)
{
- if(facetIsNormal_(normalIntersection.indexInInside(), intersection_.indexInInside()) && normalIntersection.neighbor())
- {
- const auto& neighborElement = normalIntersection.outside();
+ assert(pairData_[pairIdx].normalPair.second == -1);
+ const auto distance = FacetCenter - elementCenter;
+ this->pairData_[pairIdx].normalDistance = std::move(distance.two_norm());
+ numPairOuterNormalIdx++;
+ }
+ }
+
- for(const auto& neighborElementIs : intersections(gridView_, neighborElement))
+ // get the parallel Dofs
+ const int parallelLocalIdx = intersection_.indexInInside();
+ int numPairParallelIdx = 0;
+ std::stack<Element> parallelElementStack;
+ for(const auto& intersection : intersections(gridView_, element_))
+ {
+ if( facetIsNormal_(intersection.indexInInside(), parallelLocalIdx) )
+ {
+ if( intersection.neighbor() )
+ {
+ auto parallelAxisIdx = directionIndex(intersection);
+ auto localNormalIntersectionIndex = intersection.indexInInside();
+ parallelElementStack.push(element_);
+ bool keepStacking = (parallelElementStack.size() < order_+1);
+ while(keepStacking)
{
- // iterate over facets sub-entities
- if(facetIsNormal_(normalIntersection.indexInInside(), neighborElementIs.indexInInside()))
- setPairInfo_(neighborElementIs.indexInInside(), neighborElement, true);
+ auto e = parallelElementStack.top();
+ for(const auto& normalIntersection : intersections(gridView_, e))
+ {
+ if( normalIntersection.indexInInside() == localNormalIntersectionIndex )
+ {
+ if( normalIntersection.neighbor() )
+ {
+ parallelElementStack.push(normalIntersection.outside());
+ keepStacking = (parallelElementStack.size() < order_+1);
+ }
+ else
+ {
+ keepStacking = false;
+ }
+ }
+ }
+ }
+ while(!parallelElementStack.empty())
+ {
+ if(parallelElementStack.size() > 1)
+ {
+ this->pairData_[numPairParallelIdx].parallelDofs[parallelElementStack.size()-2] = gridView_.indexSet().subIndex(parallelElementStack.top(), parallelLocalIdx, codimIntersection);
+ }
+ this->pairData_[numPairParallelIdx].parallelDistances[parallelElementStack.size()-1] = setParallelPairDistances_(parallelElementStack.top(), parallelAxisIdx);
+ parallelElementStack.pop();
}
}
- }
+ else
+ {
+ // If the intersection has no neighbor we have to deal with the virtual outer parallel dof
+ const auto& elementCenter = this->element_.geometry().center();
+ const auto& boundaryFacetCenter = intersection.geometry().center();
- // fill the normal pair entries
- for(int pairIdx = 0; pairIdx < numPairs; ++pairIdx)
- {
- assert(pairData_[pairIdx].normalPair.second == -1);
- const auto& elementCenter = this->element_.geometry().center();
- const auto& boundaryFacetCenter = intersection_.geometry().center();
- const auto distance = boundaryFacetCenter - elementCenter;
+ const auto distance = boundaryFacetCenter - elementCenter;
+ const auto virtualFirstParallelFaceDofPos = this->intersection_.geometry().center() + distance;
+
+ this->pairData_[numPairParallelIdx].virtualFirstParallelFaceDofPos = std::move(virtualFirstParallelFaceDofPos);
- pairData_[pairIdx].normalDistance = std::move(distance.two_norm());
+ // The distance is saved doubled because with scvf.cellCenteredSelfToFirstParallelDistance
+ // an average between parallelDistances[0] and parallelDistances[1] will be computed
+ this->pairData_[numPairParallelIdx].parallelDistances[0] = std::move(distance.two_norm() * 2);
+ }
+ numPairParallelIdx++;
}
}
-
- treatVirtualOuterParallelFaceDofs_();
}
- /*!
+ /*!
* \brief Returns the local opposing intersection index
*
* \param idx The local index of the intersection itself
@@ -265,7 +462,7 @@ private:
return (idx % 2) ? (idx - 1) : (idx + 1);
}
- /*!
+ /*!
* \brief Returns true if the intersection lies normal to another given intersection
*
* \param selfIdx The local index of the intersection itself
@@ -276,201 +473,85 @@ private:
return !(selfIdx == otherIdx || localOppositeIdx_(selfIdx) == otherIdx);
};
- /*!
- * \brief Returns the global index of the common entity
- *
- * \param localIdx The local index of the common entity
- * \param element The element
- */
- int localToGlobalCommonEntityIdx_(const int localIdx, const Element& element) const
- {
- return this->gridView_.indexSet().subIndex(element, localIdx, codimCommonEntity);
- };
-
auto getFacet_(const int localFacetIdx, const Element& element) const
{
return element.template subEntity <1> (localFacetIdx);
};
- void setPairInfo_(const int isIdx, const Element& element, const bool isParallel)
+ //! Sets the information about the normal faces (within the element)
+ void setNormalPairFirstInfo_(const int isIdx, const Element& element, const int numPairsIdx)
{
- const auto referenceElement = ReferenceElements::general(element_.geometry().type());
+ // store the inner normal dofIdx
+ auto& dofIdx = pairData_[numPairsIdx].normalPair.first;
+ dofIdx = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
- // iterate over facets sub-entities
- for(int i = 0; i < numFacetSubEntities; ++i)
- {
- int localCommonEntIdx = referenceElement.subEntity(isIdx, 1, i, codimCommonEntity);
- int globalCommonEntIdx = localToGlobalCommonEntityIdx_(localCommonEntIdx, element);
-
- // fill the normal pair entries
- for(int pairIdx = 0; pairIdx < numPairs; ++pairIdx)
- {
- if(globalCommonEntIdx == pairData_[pairIdx].globalCommonEntIdx)
- {
- auto& dofIdx = isParallel ? pairData_[pairIdx].outerParallelFaceDofIdx : pairData_[pairIdx].normalPair.second;
- dofIdx = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
- if(isParallel)
- {
- const auto& selfFacet = getFacet_(intersection_.indexInInside(), element_);
- const auto& parallelFacet = getFacet_(isIdx, element);
- pairData_[pairIdx].parallelDistance = (selfFacet.geometry().center() - parallelFacet.geometry().center()).two_norm();
- }
- else
- {
- const auto& outerNormalFacet = getFacet_(isIdx, element);
- const auto outerNormalFacetPos = outerNormalFacet.geometry().center();
- pairData_[pairIdx].normalDistance = (innerNormalFacePos_[pairIdx] - outerNormalFacetPos).two_norm();
- }
- }
- }
- }
+ // store the local normal facet index
+ this->pairData_[numPairsIdx].localNormalFaceIdx = isIdx;
}
- template<class Indices>
- void setInnerIndices_(const Indices& indices)
+ //! Sets the information about the normal faces (in the neighbor element)
+ void setNormalPairSecondInfo_(const int isIdx, const Element& element, const int numPairsIdx)
{
- for(int i = 0; i < numPairs; ++i)
- {
- this->pairData_[i].normalPair.first = this->gridView_.indexSet().subIndex(this->intersection_.inside(), indices.normalLocalDofIdx[i], codimIntersection);
- this->pairData_[i].globalCommonEntIdx = this->gridView_.indexSet().subIndex(this->intersection_.inside(), indices.localCommonEntIdx[i], codimCommonEntity);
- this->pairData_[i].localNormalFaceIdx = indices.normalLocalDofIdx[i];
- }
+ // store the dofIdx
+ auto& dofIdx = pairData_[numPairsIdx].normalPair.second;
+ dofIdx = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
+
+ // store the element distance
+ const auto& outerNormalFacet = getFacet_(isIdx, element);
+ const auto outerNormalFacetPos = outerNormalFacet.geometry().center();
+ const auto& innerNormalFacet = getFacet_(isIdx, element_);
+ const auto innerNormalFacetPos = innerNormalFacet.geometry().center();
+ this->pairData_[numPairsIdx].normalDistance = (innerNormalFacetPos - outerNormalFacetPos).two_norm();
}
- /*!
- * \brief Sets the position of "virtual" dofs on facets which are perpendicular to facets on the domain boundary.
- * This is required to account e.g. for wall friction.
- */
- void treatVirtualOuterParallelFaceDofs_()
+ //! Sets the information about the parallel distances
+ Scalar setParallelPairDistances_(const Element& element, const int parallelAxisIdx) const
{
- // get the local index of the facet we are dealing with in this class
- const int localIntersectionIdx = this->intersection_.indexInInside();
-
- // iterate over all intersections of the element, this facet is part of
- for(const auto& is : intersections(this->gridView_, this->element_))
+ Scalar distance = 0;
+ std::vector<GlobalPosition> faces;
+ faces.reserve(numfacets);
+ for (const auto& intElement : intersections(gridView_, element))
{
- const int otherIsIdx = is.indexInInside();
- // check if any of these intersection, which are normal to our facet, lie on the domain boundary
- if(( !is.neighbor() ) && this->facetIsNormal_(localIntersectionIdx, otherIsIdx) )
- {
- const auto& elementCenter = this->element_.geometry().center();
- const auto& boundaryFacetCenter = is.geometry().center();
-
- const auto distance = boundaryFacetCenter - elementCenter;
- const auto virtualOuterParallelFaceDofPos = this->intersection_.geometry().center() + distance;
-
- auto foundCorrectIdx = [otherIsIdx](const auto& x) { return x.localNormalFaceIdx == otherIsIdx; };
- const int index = std::find_if(this->pairData_.begin(), this->pairData_.end(), foundCorrectIdx) - this->pairData_.begin();
-
- this->pairData_[index].virtualOuterParallelFaceDofPos = std::move(virtualOuterParallelFaceDofPos);
- this->pairData_[index].parallelDistance = std::move(distance.two_norm());
- }
+ faces.push_back(intElement.geometry().center());
}
- }
-
- auto getLocalIndices_(const int localFacetIdx) const
- {
- struct Indices
- {
- std::array<int, numPairs> normalLocalDofIdx;
- std::array<int, numPairs> localCommonEntIdx;
- };
-
- Indices indices;
- if (dim == 1)
- return indices;
-
- switch(localFacetIdx)
+ switch(parallelAxisIdx)
{
case 0:
- indices.normalLocalDofIdx[0] = 3;
- indices.normalLocalDofIdx[1] = 2;
- indices.localCommonEntIdx[0] = 2;
- indices.localCommonEntIdx[1] = 0;
- if(dim == 3)
- {
- indices.normalLocalDofIdx[2] = 4;
- indices.normalLocalDofIdx[3] = 5;
- indices.localCommonEntIdx[2] = 4;
- indices.localCommonEntIdx[3] = 8;
- }
- break;
+ {
+ distance = (faces[1] - faces[0]).two_norm();
+ break;
+ }
case 1:
- indices.normalLocalDofIdx[0] = 2;
- indices.normalLocalDofIdx[1] = 3;
- indices.localCommonEntIdx[0] = 1;
- indices.localCommonEntIdx[1] = 3;
- if(dim == 3)
- {
- indices.normalLocalDofIdx[2] = 4;
- indices.normalLocalDofIdx[3] = 5;
- indices.localCommonEntIdx[2] = 5;
- indices.localCommonEntIdx[3] = 9;
- }
- break;
+ {
+ distance = (faces[3] - faces[2]).two_norm();
+ break;
+ }
case 2:
- indices.normalLocalDofIdx[0] = 0;
- indices.normalLocalDofIdx[1] = 1;
- indices.localCommonEntIdx[0] = 0;
- indices.localCommonEntIdx[1] = 1;
- if(dim == 3)
- {
- indices.normalLocalDofIdx[2] = 4;
- indices.normalLocalDofIdx[3] = 5;
- indices.localCommonEntIdx[2] = 6;
- indices.localCommonEntIdx[3] = 10;
- }
- break;
- case 3:
- indices.normalLocalDofIdx[0] = 1;
- indices.normalLocalDofIdx[1] = 0;
- indices.localCommonEntIdx[0] = 3;
- indices.localCommonEntIdx[1] = 2;
- if(dim == 3)
- {
- indices.normalLocalDofIdx[2] = 4;
- indices.normalLocalDofIdx[3] = 5;
- indices.localCommonEntIdx[2] = 7;
- indices.localCommonEntIdx[3] = 11;
- }
- break;
- case 4:
- assert(dim == 3);
- indices.normalLocalDofIdx[0] = 0;
- indices.normalLocalDofIdx[1] = 1;
- indices.normalLocalDofIdx[2] = 3;
- indices.normalLocalDofIdx[3] = 2;
- indices.localCommonEntIdx[0] = 4;
- indices.localCommonEntIdx[1] = 5;
- indices.localCommonEntIdx[2] = 7;
- indices.localCommonEntIdx[3] = 6;
- break;
- case 5:
+ {
assert(dim == 3);
- indices.normalLocalDofIdx[0] = 0;
- indices.normalLocalDofIdx[1] = 1;
- indices.normalLocalDofIdx[2] = 2;
- indices.normalLocalDofIdx[3] = 3;
- indices.localCommonEntIdx[0] = 8;
- indices.localCommonEntIdx[1] = 9;
- indices.localCommonEntIdx[2] = 10;
- indices.localCommonEntIdx[3] = 11;
- break;
+ distance = (faces[5] - faces[4]).two_norm();
+ break;
+ }
default:
+ {
DUNE_THROW(Dune::InvalidStateException, "Something went terribly wrong");
+ }
}
- return indices;
+ return distance;
}
- // TODO: check whether to use references here or not
Intersection intersection_; //!< The intersection of interest
const Element element_; //!< The respective element
const typename Element::Geometry elementGeometry_; //!< Reference to the element geometry
- const GridView gridView_;
- std::array<PairData<Scalar, GlobalPosition>, numPairs> pairData_; //!< collection of pair information
+ const GridView gridView_; //!< The grid view
+ AxisData<Scalar> axisData_; //!< Data related to forward and backward faces
+ std::array<PairData<Scalar, GlobalPosition>, numPairs> pairData_; //!< Collection of pair information related to normal and parallel faces
std::vector<GlobalPosition> innerNormalFacePos_;
};
+template<class GridView>
+int FreeFlowStaggeredGeometryHelper<GridView>::order_ = 1;
+
} // end namespace Dumux
#endif
diff --git a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
index 924c9b897e6e495dbc25a8a13d4de9a88de0e69c..d9f06de4a8887ef2fe77313ff7269ddcbc680734 100644
--- a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
@@ -73,8 +73,7 @@ public:
const typename Intersection::Geometry& isGeometry,
GridIndexType scvfIndex,
const std::vector<GridIndexType>& scvIndices,
- const GeometryHelper& geometryHelper
- )
+ const GeometryHelper& geometryHelper)
: ParentType(),
geomType_(isGeometry.type()),
area_(isGeometry.volume()),
@@ -83,9 +82,8 @@ public:
scvfIndex_(scvfIndex),
scvIndices_(scvIndices),
boundary_(is.boundary()),
- dofIdx_(geometryHelper.dofIndex()),
- dofIdxOpposingFace_(geometryHelper.dofIndexOpposingFace()),
- selfToOppositeDistance_(geometryHelper.selfToOppositeDistance()),
+
+ axisData_(geometryHelper.axisData()),
pairData_(geometryHelper.pairData()),
localFaceIdx_(geometryHelper.localFaceIndex()),
dirIdx_(geometryHelper.directionIndex()),
@@ -181,18 +179,6 @@ public:
return Geometry(geomType_, corners_);
}
- //! The global index of the dof living on this face
- GridIndexType dofIndex() const
- {
- return dofIdx_;
- }
-
- //! The global index of the dof living on the opposing face
- GridIndexType dofIndexOpposingFace() const
- {
- return dofIdxOpposingFace_;
- }
-
//! The local index of this sub control volume face
LocalIndexType localFaceIdx() const
{
@@ -205,12 +191,6 @@ public:
return dirIdx_;
}
- //! The distance between the position of the dof itself and the one of the oppsing dof
- Scalar selfToOppositeDistance() const
- {
- return selfToOppositeDistance_;
- }
-
//! Returns whether the unitNormal of the face points in positive coordinate direction
bool normalInPosCoordDir() const
{
@@ -235,21 +215,102 @@ public:
return pairData_;
}
+ //! Return an array of all pair data
+ const auto& axisData() const
+ {
+ return axisData_;
+ }
+
+ //! Returns @c true if the face is a ghost face
bool isGhostFace() const
{
return isGhostFace_;
}
- bool hasParallelNeighbor(const int localFaceIdx) const
+ /*!
+ * \brief Check if the face has a parallel neighbor
+ *
+ * \param localSubFaceIdx The local index of the subface
+ * \param parallelDegreeIdx The index describing how many faces away from the self face
+ */
+ bool hasParallelNeighbor(const int localSubFaceIdx, const int parallelDegreeIdx) const
+ {
+ return !(pairData(localSubFaceIdx).parallelDofs[parallelDegreeIdx] < 0);
+ }
+
+ /*!
+ * \brief Check if the face has an outer normal neighbor
+ *
+ * \param localSubFaceIdx The local index of the subface
+ */
+ bool hasOuterNormal(const int localSubFaceIdx) const
+ {
+ return !(pairData_[localSubFaceIdx].normalPair.second < 0);
+ }
+
+ /*!
+ * \brief Check if the face has a backward neighbor
+ *
+ * \param backwardIdx The index describing how many faces backward this dof is from the opposite face
+ */
+ bool hasBackwardNeighbor(const int backwardIdx) const
{
- return pairData_[localFaceIdx].outerParallelFaceDofIdx >= 0;
+ return !(axisData().inAxisBackwardDofs[backwardIdx] < 0);
}
- bool hasFrontalNeighbor(const int localFaceIdx) const
+ /*!
+ * \brief Check if the face has a forward neighbor
+ *
+ * \param forwardIdx The index describing how many faces forward this dof is of the self face
+ */
+ bool hasForwardNeighbor(const int forwardIdx) const
{
- return pairData_[localFaceIdx].normalPair.second >= 0;
+ return !(axisData().inAxisForwardDofs[forwardIdx] < 0);
}
+ //! Returns the dof of the face
+ GridIndexType dofIndex() const
+ {
+ return axisData().selfDof;
+ }
+
+ //! Returns the dof of the opposing face
+ GridIndexType dofIndexOpposingFace() const
+ {
+ return axisData().oppositeDof;
+ }
+
+ //! Returns the dof the first forward face
+ GridIndexType dofIndexForwardFace() const
+ {
+ return axisData().inAxisForwardDofs[0];
+ }
+
+ //! Returns the dof of the first backward face
+ GridIndexType dofIndexBackwardFace() const
+ {
+ return axisData().inAxisBackwardDofs[0];
+ }
+
+ //! Returns the distance between the face and the opposite one
+ Scalar selfToOppositeDistance() const
+ {
+ return axisData().selfToOppositeDistance;
+ }
+
+ /*!
+ * \brief Returns the distance between the parallel dofs
+ *
+ * \param localSubFaceIdx The local index of the subface
+ * \param parallelDegreeIdx The index describing how many faces away from the self
+ */
+ Scalar cellCenteredParallelDistance(const int localSubFaceIdx, const int parallelDegreeIdx) const
+ {
+ return (pairData(localSubFaceIdx).parallelDistances[parallelDegreeIdx] +
+ pairData(localSubFaceIdx).parallelDistances[parallelDegreeIdx+1]) * 0.5;
+ }
+
+
private:
Dune::GeometryType geomType_;
std::vector<GlobalPosition> corners_;
@@ -261,10 +322,11 @@ private:
bool boundary_;
int dofIdx_;
- int dofIdxOpposingFace_;
Scalar selfToOppositeDistance_;
+ AxisData<Scalar> axisData_;
std::array<PairData<Scalar, GlobalPosition>, numPairs> pairData_;
- LocalIndexType localFaceIdx_;
+
+ int localFaceIdx_;
unsigned int dirIdx_;
int outerNormalSign_;
bool isGhostFace_;
@@ -272,4 +334,4 @@ private:
} // end namespace Dumux
-#endif
+#endif // DUMUX_DISCRETIZATION_STAGGERED_FREE_FLOW_SUBCONTROLVOLUMEFACE_HH
diff --git a/dumux/discretization/staggered/fvgridgeometry.hh b/dumux/discretization/staggered/fvgridgeometry.hh
index bdf0c3fd28814c6e195c80555903d28eb61027b4..51b9c01e4474b9ee380eab996ead9ebf9503c3ca 100644
--- a/dumux/discretization/staggered/fvgridgeometry.hh
+++ b/dumux/discretization/staggered/fvgridgeometry.hh
@@ -227,6 +227,8 @@ public:
if (!CheckOverlapSize<DiscretizationMethod::staggered>::isValid(gridView))
DUNE_THROW(Dune::InvalidStateException, "The staggered discretization method needs at least an overlap of 1 for parallel computations. "
<< " Set the parameter \"Grid.Overlap\" in the input file.");
+ if (hasParamInGroup("Discretization", "TvdApproach"))
+ GeometryHelper::setOrder(2);
}
//! The total number of sub control volumes
diff --git a/dumux/freeflow/CMakeLists.txt b/dumux/freeflow/CMakeLists.txt
index 318f762168f1509c66e1f3c7e1b5f8cac466943e..5f16de206e056526ae1a904c3c0261c80257a51b 100644
--- a/dumux/freeflow/CMakeLists.txt
+++ b/dumux/freeflow/CMakeLists.txt
@@ -7,4 +7,5 @@ install(FILES
properties.hh
turbulenceproperties.hh
volumevariables.hh
+higherorderapproximation.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow)
diff --git a/dumux/freeflow/higherorderapproximation.hh b/dumux/freeflow/higherorderapproximation.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0d6da3d71592bb247470bd343902b335da873cca
--- /dev/null
+++ b/dumux/freeflow/higherorderapproximation.hh
@@ -0,0 +1,353 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/** \file
+ * \brief This file contains different higher order methods for approximating the velocity.
+ */
+
+#ifndef DUMUX_HIGHER_ORDER_VELOCITY_APPROXIMATION_HH
+#define DUMUX_HIGHER_ORDER_VELOCITY_APPROXIMATION_HH
+
+#include <cmath>
+#include <functional>
+#include <iostream>
+#include <string>
+
+#include <dumux/common/exceptions.hh>
+
+namespace Dumux {
+
+//! \brief Available Tvd approaches
+enum class TvdApproach
+{
+ none = 0,
+ uniform = 1,
+ li = 2,
+ hou = 3
+};
+
+//! \biraf Available differencing schemes
+enum class DifferencingScheme
+{
+ none = 0,
+ vanleer = 1,
+ vanalbada = 2,
+ minmod = 3,
+ superbee = 4,
+ umist = 5,
+ mclimiter = 6,
+ wahyd = 7
+};
+
+/**
+ * \brief This file contains different higher order methods for approximating the velocity.
+ */
+template<class Scalar>
+class HigherOrderApproximation
+{
+public:
+ HigherOrderApproximation(const std::string& paramGroup = "")
+ {
+ if (hasParamInGroup(paramGroup, "Discretization.TvdApproach"))
+ {
+ // Read the runtime parameters
+ tvdApproach_ = static_cast<TvdApproach>(getParamFromGroup<int>(paramGroup, "Discretization.TvdApproach"));
+ if(tvdApproach_ == TvdApproach::none)
+ {
+ DUNE_THROW(ParameterException, "\nTvd approach 0 is not implemented.\n" <<
+ static_cast<int>(TvdApproach::uniform) << ": Uniform Tvd\n" <<
+ static_cast<int>(TvdApproach::li) << ": Li's approach\n" <<
+ static_cast<int>(TvdApproach::hou) << ": Hou's approach");
+ }
+ differencingScheme_ = static_cast<DifferencingScheme>(getParamFromGroup<int>(paramGroup, "Discretization.DifferencingScheme"));
+
+ // Assign the limiter_ depending on the differencing scheme
+ switch (differencingScheme_)
+ {
+ case DifferencingScheme::vanleer :
+ {
+ limiter_ = this->vanleer;
+ break;
+ }
+ case DifferencingScheme::vanalbada :
+ {
+ if (tvdApproach_ == TvdApproach::hou)
+ DUNE_THROW(ParameterException, "\nDifferencing scheme " << static_cast<int>(DifferencingScheme::vanalbada) <<
+ " (Van Albada) is not implemented for the Tvd approach " << static_cast<int>(TvdApproach::hou) <<" (Hou).");
+ else
+ limiter_ = this->vanalbada;
+ break;
+ }
+ case DifferencingScheme::minmod :
+ {
+ limiter_ = this->minmod;
+ break;
+ }
+ case DifferencingScheme::superbee :
+ {
+ limiter_ = this->superbee;
+ break;
+ }
+ case DifferencingScheme::umist :
+ {
+ limiter_ = this->umist;
+ break;
+ }
+ case DifferencingScheme::mclimiter :
+ {
+ limiter_ = this->mclimiter;
+ break;
+ }
+ case DifferencingScheme::wahyd :
+ {
+ limiter_ = this->wahyd;
+ break;
+ }
+ default:
+ {
+ DUNE_THROW(ParameterException, "\nDifferencing scheme " << static_cast<int>(differencingScheme_) <<
+ " is not implemented.\n");
+ break;
+ }
+ }
+ }
+ else
+ {
+ // If the runtime parameters are not specified we will use upwind
+ tvdApproach_ = TvdApproach::none;
+ differencingScheme_ = DifferencingScheme::none;
+ }
+ }
+
+ /**
+ * \brief Upwind Method
+ */
+ Scalar upwind(const Scalar downstreamVelocity,
+ const Scalar upstreamVelocity,
+ const Scalar density) const
+ {
+ return upstreamVelocity * density;
+ }
+
+ /**
+ * \brief Central Differencing Method
+ */
+ Scalar centralDifference(const Scalar downstreamVelocity,
+ const Scalar upstreamVelocity,
+ const Scalar density) const
+ {
+ return (0.5 * (upstreamVelocity + downstreamVelocity)) * density;
+ }
+
+ /**
+ * \brief Linear Upwind Method
+ */
+ Scalar linearUpwind(const Scalar upstreamVelocity,
+ const Scalar upUpstreamVelocity,
+ const Scalar upstreamToDownstreamDistance,
+ const Scalar upUpstreamToUpstreamDistance,
+ const Scalar density) const
+ {
+ Scalar zeroOrder = upstreamVelocity;
+ Scalar firstOrder = -1.0 * ((upstreamVelocity - upUpstreamVelocity) / upUpstreamToUpstreamDistance) * ( upstreamToDownstreamDistance / -2.0);
+ return (zeroOrder + firstOrder) * density;
+ }
+
+ /**
+ * \brief QUICK upwinding Scheme: Quadratic Upstream Interpolation for Convective Kinematics
+ */
+ Scalar upwindQUICK(const Scalar downstreamVelocity,
+ const Scalar upstreamVelocity,
+ const Scalar upUpstreamVelocity,
+ const Scalar upstreamToDownstreamDistance,
+ const Scalar upUpstreamToUpstreamDistance,
+ const Scalar density) const
+ {
+ Scalar normalDistance = (upUpstreamToUpstreamDistance + upstreamToDownstreamDistance) / 2.0;
+ Scalar zeroOrder = upstreamVelocity;
+ Scalar firstOrder = ((downstreamVelocity - upstreamVelocity) / 2.0);
+ Scalar secondOrder = -(((downstreamVelocity - upstreamVelocity) / upstreamToDownstreamDistance) - ((upstreamVelocity - upUpstreamVelocity) / upUpstreamToUpstreamDistance))
+ * upstreamToDownstreamDistance * upstreamToDownstreamDistance / (8.0 * normalDistance);
+ return (zeroOrder + firstOrder + secondOrder) * density;
+ }
+
+ /**
+ * \brief Tvd Scheme: Total Variation Diminuishing
+ */
+ Scalar tvd(const Scalar downstreamVelocity,
+ const Scalar upstreamVelocity,
+ const Scalar upUpstreamVelocity,
+ const Scalar density) const
+ {
+ const Scalar ratio = (upstreamVelocity - upUpstreamVelocity) / (downstreamVelocity - upstreamVelocity);
+
+ // If the velocity field is uniform (like at the first newton step) we get a NaN
+ if(ratio > 0.0 && std::isfinite(ratio))
+ {
+ const Scalar secondOrderTerm = 0.5 * limiter_(ratio, 2.0) * (downstreamVelocity - upstreamVelocity);
+ return density * (upstreamVelocity + secondOrderTerm);
+ }
+ else
+ return density * upstreamVelocity;
+ }
+
+ /**
+ * \brief Tvd Scheme: Total Variation Diminuishing
+ *
+ * This functions manages the non uniformities of the grid according to [Li, Liao 2007].
+ * It tries to reconstruct the value for the velocity at the upstream-upstream point
+ * if the grid was uniform.
+ */
+ Scalar tvd(const Scalar downstreamVelocity,
+ const Scalar upstreamVelocity,
+ const Scalar upUpstreamVelocity,
+ const Scalar upstreamToDownstreamDistance,
+ const Scalar upUpstreamToUpstreamDistance,
+ const bool selfIsUpstream,
+ const Scalar density) const
+ {
+ // I need the information of selfIsUpstream to get the correct sign because upUpstreamToUpstreamDistance is always positive
+ const Scalar upUpstreamGradient = (upstreamVelocity - upUpstreamVelocity) / upUpstreamToUpstreamDistance * selfIsUpstream;
+
+ // Distance between the upUpstream node and the position where it should be if the grid were uniform.
+ const Scalar correctionDistance = upUpstreamToUpstreamDistance - upstreamToDownstreamDistance;
+ const Scalar reconstrutedUpUpstreamVelocity = upUpstreamVelocity + upUpstreamGradient * correctionDistance;
+ const Scalar ratio = (upstreamVelocity - reconstrutedUpUpstreamVelocity) / (downstreamVelocity - upstreamVelocity);
+
+ // If the velocity field is uniform (like at the first newton step) we get a NaN
+ if(ratio > 0.0 && std::isfinite(ratio))
+ {
+ const Scalar secondOrderTerm = 0.5 * limiter_(ratio, 2.0) * (downstreamVelocity - upstreamVelocity);
+ return density * (upstreamVelocity + secondOrderTerm);
+ }
+ else
+ return density * upstreamVelocity;
+ }
+
+ /**
+ * \brief Tvd Scheme: Total Variation Diminuishing
+ *
+ * This functions manages the non uniformities of the grid according to [Hou, Simons, Hinkelmann 2007].
+ * It should behave better then the Li's version in very stretched grids.
+ */
+ Scalar tvd(const Scalar downstreamVelocity,
+ const Scalar upstreamVelocity,
+ const Scalar upUpstreamVelocity,
+ const Scalar upstreamToDownstreamDistance,
+ const Scalar upUpstreamToUpstreamDistance,
+ const Scalar downstreamStaggeredCellSize,
+ const Scalar density) const
+ {
+ const Scalar ratio = (upstreamVelocity - upUpstreamVelocity) / (downstreamVelocity - upstreamVelocity)
+ * upstreamToDownstreamDistance / upUpstreamToUpstreamDistance;
+
+ // If the velocity field is uniform (like at the first newton step) we get a NaN
+ if(ratio > 0.0 && std::isfinite(ratio))
+ {
+ const Scalar upstreamStaggeredCellSize = 0.5 * (upstreamToDownstreamDistance + upUpstreamToUpstreamDistance);
+ const Scalar R = (upstreamStaggeredCellSize + downstreamStaggeredCellSize) / upstreamStaggeredCellSize;
+ const Scalar secondOrderTerm = limiter_(ratio, R) / R * (downstreamVelocity - upstreamVelocity);
+ return density * (upstreamVelocity + secondOrderTerm);
+ }
+ else
+ return density * upstreamVelocity;
+ }
+
+ /**
+ * \brief Van Leer flux limiter function [Van Leer 1974]
+ *
+ * With R != 2 is the modified Van Leer flux limiter function [Hou, Simons, Hinkelmann 2007]
+ */
+ static Scalar vanleer(const Scalar r, const Scalar R)
+ {
+ return R * r / (R - 1.0 + r);
+ }
+
+ /**
+ * \brief Van Albada flux limiter function [Van Albada et al. 1982]
+ */
+ static Scalar vanalbada(const Scalar r, const Scalar R)
+ {
+ return r * (r + 1.0) / (1.0 + r * r);
+ }
+
+ /**
+ * \brief MinMod flux limiter function [Roe 1985]
+ */
+ static Scalar minmod(const Scalar r, const Scalar R)
+ {
+ return std::min(r, 1.0);
+ }
+
+ /**
+ * \brief SUPERBEE flux limiter function [Roe 1985]
+ *
+ * With R != 2 is the modified SUPERBEE flux limiter function [Hou, Simons, Hinkelmann 2007]
+ */
+ static Scalar superbee(const Scalar r, const Scalar R)
+ {
+ return std::max(std::min(R * r, 1.0), std::min(r, R));
+ }
+
+ /**
+ * \brief UMIST flux limiter function [Lien and Leschziner 1993]
+ */
+ static Scalar umist(const Scalar r, const Scalar R)
+ {
+ return std::min({R * r, (r * (5.0 - R) + R - 1.0) / 4.0, (r * (R - 1.0) + 5.0 - R) / 4.0, R});
+ }
+
+ /*
+ * \brief Monotonized-Central limiter [Van Leer 1977]
+ */
+ static Scalar mclimiter(const Scalar r, const Scalar R)
+ {
+ return std::min({R * r, (r + 1.0) / 2.0, R});
+ }
+
+ /**
+ * \brief WAHYD Scheme [Hou, Simons, Hinkelmann 2007];
+ */
+ static Scalar wahyd(const Scalar r, const Scalar R)
+ {
+ return r > 1 ? std::min((r + R * r * r) / (R + r * r), R)
+ : vanleer(r, R);
+ }
+
+ //! Returns the Tvd approach
+ const TvdApproach& tvdApproach() const
+ {
+ return tvdApproach_;
+ }
+
+ //! Returns the differencing scheme
+ const DifferencingScheme& differencingScheme() const
+ {
+ return differencingScheme_;
+ }
+
+private:
+ TvdApproach tvdApproach_;
+ DifferencingScheme differencingScheme_;
+
+ std::function<Scalar(const Scalar, const Scalar)> limiter_;
+};
+
+} // end namespace Dumux
+
+#endif // DUMUX_HIGHER_ORDER_VELOCITY_APPROXIMATION_HH
diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh
index 58922e3fc8b2cece2feeb7f85a5858cfb7c0832e..d988a03955a7e0471131bbcd694974c0bf31d05a 100644
--- a/dumux/freeflow/navierstokes/problem.hh
+++ b/dumux/freeflow/navierstokes/problem.hh
@@ -28,6 +28,9 @@
#include <dumux/common/properties.hh>
#include <dumux/common/staggeredfvproblem.hh>
#include <dumux/discretization/method.hh>
+
+#include <dumux/freeflow/higherorderapproximation.hh>
+
#include "model.hh"
namespace Dumux {
@@ -92,8 +95,9 @@ public:
* \param paramGroup The parameter group in which to look for runtime parameters first (default is "")
*/
NavierStokesProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, const std::string& paramGroup = "")
- : ParentType(fvGridGeometry, paramGroup)
- , gravity_(0.0)
+ : ParentType(fvGridGeometry, paramGroup),
+ gravity_(0.0),
+ higherOrderApproximation_(paramGroup)
{
if (getParamFromGroup<bool>(paramGroup, "Problem.EnableGravity"))
gravity_[dim-1] = -9.81;
@@ -215,7 +219,13 @@ public:
using std::sqrt;
const Scalar K = asImp_().permeability(element, normalFace);
const Scalar alpha = asImp_().alphaBJ(normalFace);
- return velocitySelf / (alpha / sqrt(K) * scvf.pairData(localSubFaceIdx).parallelDistance + 1.0);
+ return velocitySelf / (alpha / sqrt(K) * scvf.cellCenteredParallelDistance(localSubFaceIdx,0) + 1.0);
+ }
+
+ //! Return the HigherOrderApproximation
+ const HigherOrderApproximation<Scalar>& higherOrderApproximation() const
+ {
+ return higherOrderApproximation_;
}
private:
@@ -229,6 +239,7 @@ private:
{ return *static_cast<const Implementation *>(this); }
GravityVector gravity_;
+ HigherOrderApproximation<Scalar> higherOrderApproximation_;
bool enableInertiaTerms_;
};
diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
index 5bf038741a1b122e14e45804fce996d0d40c228c..46ca4213af4d7a131b1d58cf93728c537641a75a 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
@@ -24,10 +24,15 @@
#ifndef DUMUX_NAVIERSTOKES_STAGGERED_FLUXVARIABLES_HH
#define DUMUX_NAVIERSTOKES_STAGGERED_FLUXVARIABLES_HH
+#include <array>
+
#include <dumux/common/math.hh>
+#include <dumux/common/exceptions.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>
+#include <dumux/freeflow/higherorderapproximation.hh>
+
#include <dumux/flux/fluxvariablesbase.hh>
#include <dumux/discretization/method.hh>
@@ -207,17 +212,61 @@ public:
// Check if the the velocity of the dof at interest lies up- or downstream w.r.t. to the transporting velocity.
const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity);
- // Lamba function to evaluate the transported momentum, regarding an user-specified upwind weight.
- auto computeMomentum = [&insideVolVars, &problem](const Scalar upstreamVelocity, const Scalar downstreamVelocity)
- {
- static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight");
- return (upwindWeight * upstreamVelocity + (1.0 - upwindWeight) * downstreamVelocity) * insideVolVars.density();
- };
+ Scalar momentum = 0.0;
+
+ // Variables that will store the velocities of interest:
+ // velocities[0]: downstream velocity
+ // velocities[1]: upstream velocity
+ // velocities[2]: upstream-upstream velocity
+ std::array<Scalar, 3> velocities{0.0, 0.0, 0.0};
+ // Variables that will store the distances between the dofs of interest:
+ // distances[0]: upstream to downstream distance
+ // distances[1]: upstream-upstream to upstream distance
+ // distances[2]: downstream staggered cell size
+ std::array<Scalar, 3> distances{0.0, 0.0, 0.0};
- // Get the momentum that is advectively transported and account for the flow direction.
- const Scalar momentum = selfIsUpstream ? computeMomentum(velocitySelf, velocityOpposite)
- : computeMomentum(velocityOpposite, velocitySelf);
+ const auto& highOrder = problem.higherOrderApproximation();
+
+ // If a Tvd approach has been specified and I am not too near to the boundary I can use a second order
+ // approximation for the velocity. In this frontal flux I use for the density always the value that I have on the scvf.
+ if (highOrder.tvdApproach() != TvdApproach::none)
+ {
+ if (canFrontalSecondOrder_(scvf, selfIsUpstream, velocities, distances, elemFaceVars[scvf]))
+ {
+ switch (highOrder.tvdApproach())
+ {
+ case TvdApproach::uniform :
+ {
+ momentum = highOrder.tvd(velocities[0], velocities[1], velocities[2], insideVolVars.density());
+ break;
+ }
+ case TvdApproach::li :
+ {
+ momentum = highOrder.tvd(velocities[0], velocities[1], velocities[2], distances[0], distances[1], selfIsUpstream, insideVolVars.density());
+ break;
+ }
+ case TvdApproach::hou :
+ {
+ momentum = highOrder.tvd(velocities[0], velocities[1], velocities[2], distances[0], distances[1], distances[2], insideVolVars.density());
+ break;
+ }
+ default:
+ {
+ DUNE_THROW(ParameterException, "\nTvd approach " << static_cast<int>(highOrder.tvdApproach()) << " is not implemented.\n" <<
+ static_cast<int>(TvdApproach::uniform) << ": Uniform Tvd\n" <<
+ static_cast<int>(TvdApproach::li) << ": Li's approach\n" <<
+ static_cast<int>(TvdApproach::hou) << ": Hou's approach");
+ break;
+ }
+ }
+ }
+ else
+ momentum = highOrder.upwind(velocities[0], velocities[1], insideVolVars.density());
+ }
+ else
+ momentum = selfIsUpstream ? highOrder.upwind(velocityOpposite, velocitySelf, insideVolVars.density())
+ : highOrder.upwind(velocitySelf, velocityOpposite, insideVolVars.density());
// Account for the orientation of the staggered face's normal outer normal vector
// (pointing in opposite direction of the scvf's one).
@@ -290,18 +339,18 @@ public:
{
const auto eIdx = scvf.insideScvIdx();
// Get the face normal to the face the dof lives on. The staggered sub face conincides with half of this normal face.
- const auto& normalFace = fvGeometry.scvf(eIdx, scvf.pairData()[localSubFaceIdx].localNormalFaceIdx);
+ const auto& normalFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localNormalFaceIdx);
bool lateralFaceHasDirichletPressure = false; // check for Dirichlet boundary condition for the pressure
bool lateralFaceHasBJS = false; // check for Beavers-Joseph-Saffman boundary condition
// Check if there is face/element parallel to our face of interest where the dof lives on. If there is no parallel neighbor,
// we are on a boundary where we have to check for boundary conditions.
- if(!scvf.hasParallelNeighbor(localSubFaceIdx))
+ if(!scvf.hasParallelNeighbor(localSubFaceIdx,0))
{
// Construct a temporary scvf which corresponds to the staggered sub face, featuring the location
// the sub faces's center.
- auto localSubFaceCenter = scvf.pairData(localSubFaceIdx).virtualOuterParallelFaceDofPos - normalFace.center();
+ auto localSubFaceCenter = scvf.pairData(localSubFaceIdx).virtualFirstParallelFaceDofPos - normalFace.center();
localSubFaceCenter *= 0.5;
localSubFaceCenter += normalFace.center();
const auto localSubFace = makeGhostFace_(normalFace, localSubFaceCenter);
@@ -384,54 +433,72 @@ private:
const Scalar transportingVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
// Check whether the own or the neighboring element is upstream.
- const bool ownElementIsUpstream = ( normalFace.directionSign() == sign(transportingVelocity) );
-
- // Get the velocities at the current (own) scvf and at the parallel one at the neighboring scvf.
- const Scalar velocitySelf = faceVars.velocitySelf();
+ const bool selfIsUpstream = ( normalFace.directionSign() == sign(transportingVelocity) );
- // Lambda to conveniently get the outer parallel velocity for normal faces that are on the boundary
- // and therefore have no neighbor. Calls the problem to retrieve a fixed value set on the boundary.
- auto getParallelVelocityFromBoundary = [&]()
- {
- // If there is a Dirichlet condition for the pressure we assume zero gradient for the velocity,
- // so the velocity at the boundary equal to that on the scvf.
- if (lateralFaceHasDirichletPressure)
- return velocitySelf;
-
- // Check if we have a Beavers-Joseph-Saffman condition.
- // If yes, the parallel velocity at the boundary is calculated accordingly.
- if (lateralFaceHasBJS)
- return problem.bjsVelocity(element, scvf, normalFace, localSubFaceIdx, velocitySelf);
+ // Get the volume variables of the own and the neighboring element
+ const auto& insideVolVars = elemVolVars[normalFace.insideScvIdx()];
+ const auto& outsideVolVars = elemVolVars[normalFace.outsideScvIdx()];
- const auto ghostFace = makeParallelGhostFace_(scvf, localSubFaceIdx);
+ Scalar momentum = 0.0;
- return problem.dirichlet(element, ghostFace)[Indices::velocity(scvf.directionIndex())];
- };
+ // Variables that will store the velocities of interest:
+ // velocities[0]: downstream velocity
+ // velocities[1]: upsteram velocity
+ // velocities[2]: upstream-upstream velocity
+ std::array<Scalar, 3> velocities{0.0, 0.0, 0.0};
- const Scalar velocityParallel = scvf.hasParallelNeighbor(localSubFaceIdx)
- ? faceVars.velocityParallel(localSubFaceIdx)
- : getParallelVelocityFromBoundary();
+ // Variables that will store the distances between the dofs of interest:
+ // distances[0]: upstream to downstream distance
+ // distances[1]: upstream-upstream to upstream distance
+ // distances[2]: downstream staggered cell size
+ std::array<Scalar, 3> distances{0.0, 0.0, 0.0};
- // Get the volume variables of the own and the neighboring element
- const auto& insideVolVars = elemVolVars[normalFace.insideScvIdx()];
- const auto& outsideVolVars = elemVolVars[normalFace.outsideScvIdx()];
+ const auto& highOrder = problem.higherOrderApproximation();
- // Lamba function to evaluate the transported momentum, regarding an user-specified upwind weight.
- auto computeMomentum = [&problem](const VolumeVariables& upstreamVolVars,
- const VolumeVariables& downstreamVolVars,
- const Scalar upstreamVelocity,
- const Scalar downstreamVelocity)
+ // If a Tvd approach has been specified and I am not too near to the boundary I can use a second order approximation.
+ if (highOrder.tvdApproach() != TvdApproach::none)
{
- static const Scalar upWindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight");
- const Scalar density = upWindWeight * upstreamVolVars.density() + (1.0 - upWindWeight) * downstreamVolVars.density();
- const Scalar transportedVelocity = upWindWeight * upstreamVelocity + (1.0 - upWindWeight) * downstreamVelocity;
- return transportedVelocity * density;
- };
+ if (canLateralSecondOrder_(scvf, selfIsUpstream, localSubFaceIdx, velocities, distances, problem, element, faceVars, lateralFaceHasDirichletPressure, lateralFaceHasBJS))
+ {
+ switch (highOrder.tvdApproach())
+ {
+ case TvdApproach::uniform :
+ {
+ momentum = highOrder.tvd(velocities[0], velocities[1], velocities[2], selfIsUpstream ? insideVolVars.density() : outsideVolVars.density());
+ break;
+ }
+ case TvdApproach::li :
+ {
+ momentum = highOrder.tvd(velocities[0], velocities[1], velocities[2], distances[0], distances[1], selfIsUpstream, selfIsUpstream ? insideVolVars.density() : outsideVolVars.density());
+ break;
+ }
+ case TvdApproach::hou :
+ {
+ momentum = highOrder.tvd(velocities[0], velocities[1], velocities[2], distances[0], distances[1], distances[2], selfIsUpstream ? insideVolVars.density() : outsideVolVars.density());
+ break;
+ }
+ default:
+ {
+ DUNE_THROW(ParameterException, "\nTvd approach " << static_cast<int>(highOrder.tvdApproach()) << " is not implemented.\n" <<
+ static_cast<int>(TvdApproach::uniform) << ": Uniform Tvd\n" <<
+ static_cast<int>(TvdApproach::li) << ": Li's approach\n" <<
+ static_cast<int>(TvdApproach::hou) << ": Hou's approach");
+ break;
+ }
+ }
+ }
+ else
+ momentum = highOrder.upwind(velocities[0], velocities[1], selfIsUpstream ? insideVolVars.density() : outsideVolVars.density());
+ }
+ else
+ {
+ const Scalar velocityFirstParallel = scvf.hasParallelNeighbor(localSubFaceIdx, 0)
+ ? faceVars.velocityParallel(localSubFaceIdx, 0)
+ : getParallelVelocityFromBoundary_(problem, scvf, normalFace, faceVars.velocitySelf(), localSubFaceIdx, element, lateralFaceHasDirichletPressure, lateralFaceHasBJS);
- // Get the momentum that is advectively transported and account for the flow direction.
- const Scalar momentum = ownElementIsUpstream ?
- computeMomentum(insideVolVars, outsideVolVars, velocitySelf, velocityParallel)
- : computeMomentum(outsideVolVars, insideVolVars, velocityParallel, velocitySelf);
+ momentum = selfIsUpstream ? highOrder.upwind(velocityFirstParallel, faceVars.velocitySelf(), insideVolVars.density())
+ : highOrder.upwind(faceVars.velocitySelf(), velocityFirstParallel, outsideVolVars.density());
+ }
// Account for the orientation of the staggered normal face's outer normal vector
// and its area (0.5 of the coinciding scfv).
@@ -519,24 +586,14 @@ private:
// and at the parallel one at the neighboring scvf.
const Scalar innerParallelVelocity = faceVars.velocitySelf();
- // Lambda to conveniently get the outer parallel velocity for normal faces that are on the boundary
- // and therefore have no neighbor. Calls the problem to retrieve a fixed value set on the boundary.
- auto getParallelVelocityFromBoundary = [&]()
- {
- const auto ghostFace = makeParallelGhostFace_(scvf, localSubFaceIdx);
- if (lateralFaceHasBJS)
- return problem.bjsVelocity(element, scvf, normalFace, localSubFaceIdx, innerParallelVelocity);
- return problem.dirichlet(element, ghostFace)[Indices::velocity(scvf.directionIndex())];
- };
-
- const Scalar outerParallelVelocity = scvf.hasParallelNeighbor(localSubFaceIdx)
- ? faceVars.velocityParallel(localSubFaceIdx)
- : getParallelVelocityFromBoundary();
+ const Scalar velocityFirstParallel = scvf.hasParallelNeighbor(localSubFaceIdx,0)
+ ? faceVars.velocityParallel(localSubFaceIdx,0)
+ : getParallelVelocityFromBoundary_(problem, scvf, normalFace, innerParallelVelocity, localSubFaceIdx, element, false, lateralFaceHasBJS);
// The velocity gradient already accounts for the orientation
// of the staggered face's outer normal vector.
- const Scalar parallelGradient = (outerParallelVelocity - innerParallelVelocity)
- / scvf.pairData(localSubFaceIdx).parallelDistance;
+ const Scalar parallelGradient = (velocityFirstParallel - innerParallelVelocity)
+ / scvf.cellCenteredParallelDistance(localSubFaceIdx,0);
normalDiffusiveFlux -= muAvg * parallelGradient;
}
@@ -603,7 +660,7 @@ private:
//! helper function to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
SubControlVolumeFace makeParallelGhostFace_(const SubControlVolumeFace& ownScvf, const int localSubFaceIdx) const
{
- return makeGhostFace_(ownScvf, ownScvf.pairData(localSubFaceIdx).virtualOuterParallelFaceDofPos);
+ return makeGhostFace_(ownScvf, ownScvf.pairData(localSubFaceIdx).virtualFirstParallelFaceDofPos);
};
//! helper function to get the averaged extrusion factor for a face
@@ -613,7 +670,6 @@ private:
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
return harmonicMean(insideVolVars.extrusionFactor(), outsideVolVars.extrusionFactor());
}
-
//! do nothing if no turbulence model is used
template<class ...Args, bool turbulenceModel = ModelTraits::usesTurbulenceModel(), std::enable_if_t<!turbulenceModel, int> = 0>
bool incorporateWallFunction_(Args&&... args) const
@@ -639,9 +695,233 @@ private:
}
else
return false;
- };
+ }
+
+ /*!
+ * \brief Check if a second order approximation for the frontal part of the advective term can be used
+ *
+ * This helper function checks if the scvf of interest is not too near to the
+ * boundary so that a dof upstream with respect to the upstream dof is available.
+ *
+ * \param ownScvf The SubControlVolumeFace we are considering
+ * \param selfIsUpstream @c true if the velocity ownScvf is upstream wrt the transporting velocity
+ * \param velocities Variable that will store the velocities of interest
+ * \param distances Variable that will store the distances of interest
+ * \param faceVars The face variables related to ownScvf
+ */
+ bool canFrontalSecondOrder_(const SubControlVolumeFace& ownScvf,
+ const bool selfIsUpstream,
+ std::array<Scalar, 3>& velocities,
+ std::array<Scalar, 3>& distances,
+ const FaceVariables& faceVars) const
+ {
+ // Depending on selfIsUpstream I can assign the downstream and the upstream velocities,
+ // then I have to check if I have a forward or a backward neighbor to retrieve
+ // an "upstream-upstream velocity" and be able to use a second order scheme.
+ if (selfIsUpstream)
+ {
+ velocities[0] = faceVars.velocityOpposite();
+ velocities[1] = faceVars.velocitySelf();
+
+ if (ownScvf.hasForwardNeighbor(0))
+ {
+ velocities[2] = faceVars.velocityForward(0);
+ distances[0] = ownScvf.selfToOppositeDistance();
+ distances[1] = ownScvf.axisData().inAxisForwardDistances[0];
+ distances[2] = 0.5 * (ownScvf.axisData().selfToOppositeDistance + ownScvf.axisData().inAxisBackwardDistances[0]);
+ return true;
+ }
+ else
+ return false;
+ }
+ else
+ {
+ velocities[0] = faceVars.velocitySelf();
+ velocities[1] = faceVars.velocityOpposite();
+
+ if (ownScvf.hasBackwardNeighbor(0))
+ {
+ velocities[2] = faceVars.velocityBackward(0);
+ distances[0] = ownScvf.selfToOppositeDistance();
+ distances[1] = ownScvf.axisData().inAxisBackwardDistances[0];
+ distances[2] = 0.5 * (ownScvf.axisData().selfToOppositeDistance + ownScvf.axisData().inAxisForwardDistances[0]);
+ return true;
+ }
+ else
+ return false;
+ }
+ }
+
+ /*!
+ * \brief Check if a second order approximation for the lateral part of the advective term can be used
+ *
+ * This helper function checks if the scvf of interest is not too near to the
+ * boundary so that a dof upstream with respect to the upstream dof is available.
+ *
+ * \param ownScvf The SubControlVolumeFace we are considering
+ * \param selfIsUpstream @c true if the velocity ownScvf is upstream wrt the transporting velocity
+ * \param localSubFaceIdx The local subface index
+ * \param velocities Variable that will store the velocities of interest
+ * \param distances Variable that will store the distances of interest
+ */
+ bool canLateralSecondOrder_(const SubControlVolumeFace& ownScvf,
+ const bool selfIsUpstream,
+ const int localSubFaceIdx,
+ std::array<Scalar, 3>& velocities,
+ std::array<Scalar, 3>& distances,
+ const Problem& problem,
+ const Element& element,
+ const FaceVariables& faceVars,
+ const bool lateralFaceHasDirichletPressure,
+ const bool lateralFaceHasBJS) const
+ {
+ const SubControlVolumeFace& normalFace = problem.fvGridGeometry().scvf(ownScvf.insideScvIdx(), ownScvf.pairData(localSubFaceIdx).localNormalFaceIdx);
+
+ // The local index of the faces that is opposite to localSubFaceIdx
+ const int oppositeSubFaceIdx = localSubFaceIdx % 2 ? localSubFaceIdx - 1 : localSubFaceIdx + 1;
+
+ if (selfIsUpstream)
+ {
+ // I can assign the upstream velocity. The downstream velocity can be assigned or retrieved
+ // from the boundary if there is no parallel neighbor.
+ velocities[1] = faceVars.velocitySelf();
+
+ if(ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
+ {
+ velocities[0] = faceVars.velocityParallel(localSubFaceIdx, 0);
+ distances[2] = ownScvf.pairData(localSubFaceIdx).parallelDistances[1];
+ }
+ else
+ {
+ velocities[0] = getParallelVelocityFromBoundary_(problem, ownScvf, normalFace, faceVars.velocitySelf(), localSubFaceIdx, element, lateralFaceHasDirichletPressure, lateralFaceHasBJS);
+ distances[2] = ownScvf.area() / 2.0;
+ }
+
+ // The "upstream-upsteram" velocity is retrieved from the other parallel neighbor
+ // or from the boundary.
+ if (ownScvf.hasParallelNeighbor(oppositeSubFaceIdx, 0))
+ velocities[2] = faceVars.velocityParallel(oppositeSubFaceIdx, 0);
+ else
+ velocities[2] = getParallelVelocityFromOtherBoundary_(problem, ownScvf, oppositeSubFaceIdx, element, velocities[1]);
+
+ distances[0] = ownScvf.cellCenteredParallelDistance(localSubFaceIdx, 0);
+ distances[1] = ownScvf.cellCenteredParallelDistance(oppositeSubFaceIdx, 0);
+
+ return true;
+ }
+ else
+ {
+ // The self velocity is downstream, then if there is no parallel neighbor I can not use
+ // a second order approximation beacuse I have only two velocities.
+ velocities[0] = faceVars.velocitySelf();
+
+ if (!ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
+ {
+ velocities[1] = getParallelVelocityFromBoundary_(problem, ownScvf, normalFace, faceVars.velocitySelf(), localSubFaceIdx, element, lateralFaceHasDirichletPressure, lateralFaceHasBJS);
+ return false;
+ }
+
+ velocities[1] = faceVars.velocityParallel(localSubFaceIdx, 0);
+
+ // If there is another parallel neighbor I can assign the "upstream-upstream"
+ // velocity, otherwise I retrieve it from the boundary.
+ if (ownScvf.hasParallelNeighbor(localSubFaceIdx, 1))
+ velocities[2] = faceVars.velocityParallel(localSubFaceIdx, 1);
+ else
+ {
+ const Element& elementParallel = problem.fvGridGeometry().element(problem.fvGridGeometry().scv(normalFace.outsideScvIdx()));
+ const SubControlVolumeFace& firstParallelScvf = problem.fvGridGeometry().scvf(normalFace.outsideScvIdx(), ownScvf.localFaceIdx());
+ velocities[2] = getParallelVelocityFromOtherBoundary_(problem, firstParallelScvf, localSubFaceIdx, elementParallel, velocities[1]);
+ }
+
+ distances[0] = ownScvf.cellCenteredParallelDistance(localSubFaceIdx, 0);
+ distances[1] = ownScvf.cellCenteredParallelDistance(localSubFaceIdx, 1);
+ distances[2] = ownScvf.area();
+
+ return true;
+ }
+ }
+
+ /*!
+ * \brief Return the outer parallel velocity for normal faces that are on the boundary and therefore have no neighbor.
+ *
+ * Calls the problem to retrieve a fixed value set on the boundary.
+ *
+ * \param problem The problem
+ * \param scvf The SubControlVolumeFace that is normal to the boundary
+ * \param normalFace The face at the boundary
+ * \param velocitySelf the velocity at scvf
+ * \param localSubFaceIdx The local index of the face that is on the boundary
+ * \param element The element that is on the boundary
+ * \param lateralFaceHasDirichletPressure @c true if there is a dirichlet condition for the pressure on the boundary
+ * \param lateralFaceHasBJS @c true if there is a BJS condition fot the velocity on the boundary
+ */
+ Scalar getParallelVelocityFromBoundary_(const Problem& problem,
+ const SubControlVolumeFace& scvf,
+ const SubControlVolumeFace& normalFace,
+ const Scalar velocitySelf,
+ const int localSubFaceIdx,
+ const Element& element,
+ const bool lateralFaceHasDirichletPressure,
+ const bool lateralFaceHasBJS) const
+ {
+ // If there is a Dirichlet condition for the pressure we assume zero gradient for the velocity,
+ // so the velocity at the boundary equal to that on the scvf.
+ if (lateralFaceHasDirichletPressure)
+ return velocitySelf;
+
+ const auto ghostFace = makeParallelGhostFace_(scvf, localSubFaceIdx);
+ if (lateralFaceHasBJS)
+ return problem.bjsVelocity(element, scvf, normalFace, localSubFaceIdx, velocitySelf);
+ return problem.dirichlet(element, ghostFace)[Indices::velocity(scvf.directionIndex())];
+ }
+
+ /*!
+ * \brief Return a velocity value from a boundary for which the boundary conditions have to be checked.
+ *
+ * \param problem The problem
+ * \param scvf The SubControlVolumeFace that is normal to the boundary
+ * \param localIdx The local index of the face that is on the boundary
+ * \param boundaryElement The element that is on the boundary
+ * \param parallelVelocity The velocity over scvf
+ */
+ Scalar getParallelVelocityFromOtherBoundary_(const Problem& problem,
+ const SubControlVolumeFace& scvf,
+ const int localIdx,
+ const Element& boundaryElement,
+ const Scalar parallelVelocity) const
+ {
+ // A ghost subface at the boundary is created, featuring the location of the sub face's center
+ const SubControlVolumeFace& boundaryNormalFace = problem.fvGridGeometry().scvf(scvf.insideScvIdx(), scvf.pairData(localIdx).localNormalFaceIdx);
+ GlobalPosition boundarySubFaceCenter = scvf.pairData(localIdx).virtualFirstParallelFaceDofPos + boundaryNormalFace.center();
+ boundarySubFaceCenter *= 0.5;
+ const SubControlVolumeFace boundarySubFace = makeGhostFace_(boundaryNormalFace, boundarySubFaceCenter);
+
+ // The boundary condition is checked, in case of symmetry or Dirichlet for the pressure
+ // a gradient of zero is assumed in the direction normal to the bounadry, while if there is
+ // Dirichlet of BJS for the velocity the related values are exploited.
+ const auto bcTypes = problem.boundaryTypes(boundaryElement, boundarySubFace);
+
+ if (bcTypes.isDirichlet(Indices::velocity(scvf.directionIndex())))
+ {
+ const SubControlVolumeFace ghostFace = makeParallelGhostFace_(scvf, localIdx);
+ return problem.dirichlet(boundaryElement, ghostFace)[Indices::velocity(scvf.directionIndex())];
+ }
+ else if (bcTypes.isSymmetry() || bcTypes.isDirichlet(Indices::pressureIdx))
+ return parallelVelocity;
+ else if (bcTypes.isBJS(Indices::velocity(scvf.directionIndex())))
+ {
+ const SubControlVolumeFace ghostFace = makeParallelGhostFace_(scvf, localIdx);
+ return problem.bjsVelocity(boundaryElement, scvf, boundaryNormalFace, localIdx, parallelVelocity);
+ }
+ else
+ {
+ // Neumann conditions are not well implemented
+ DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the boundary conditions for the momentum equations at global position " << boundarySubFaceCenter);
+ }
+ }
};
-} // end namespace
+} // end namespace Dumux
-#endif
+#endif // DUMUX_NAVIERSTOKES_STAGGERED_FLUXVARIABLES_HH