diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index 3dd939e20fa50aa13515efb2e0a7546a73c87068..22673a693bc40003054f8640d51ea25e3cb03099 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -233,9 +233,9 @@ struct BoundaryValues { using type = UndefinedProperty; };
template<class TypeTag, class MyTypeTag>
struct StaggeredFaceSolution { using type = UndefinedProperty; }; //!< A vector containing the solution for a face (similar to ElementSolution)
template<class TypeTag, class MyTypeTag>
-struct EnableGridFaceVariablesCache { using type = UndefinedProperty; }; //!< Switch on/off caching of face variables
+struct EnableGridFaceVariablesCache { using type = UndefinedProperty; }; //!< Switch on/off caching of face variables
template<class TypeTag, class MyTypeTag>
-struct UpwindGeometryOrder { using type = UndefinedProperty; }; //!< Specifies the order of the upwinding scheme (1 == basic upwinding, 2 == higher order)
+struct UpwindSchemeOrder { using type = UndefinedProperty; }; //!< Specifies the order of the upwinding scheme (1 == basic upwinding, 2 == higher order)
/////////////////////////////////////////////////////////////
// Properties used by the mpnc model
diff --git a/dumux/discretization/staggered/freeflow/connectivitymap.hh b/dumux/discretization/staggered/freeflow/connectivitymap.hh
index 7622726ec6d57549af77c80c9f673438f3020509..3fe85d132fd8ba6577c31d61d2acad16c99285fd 100644
--- a/dumux/discretization/staggered/freeflow/connectivitymap.hh
+++ b/dumux/discretization/staggered/freeflow/connectivitymap.hh
@@ -54,7 +54,8 @@ class StaggeredFreeFlowConnectivityMap
using Stencil = std::vector<GridIndexType>;
- static constexpr bool useHigherOrder = false; //TODO
+ static constexpr int upwindSchemeOrder = FVGridGeometry::upwindSchemeOrder;
+ static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
public:
@@ -94,13 +95,6 @@ public:
}
}
- //! \brief Set the order needed by the scheme
- void setStencilOrder(const int stencilOrder)
- {
- stencilOrder_ = stencilOrder;
- }
-
-
//! Returns the stencil of a cell center dof w.r.t. other cell center dofs
const std::vector<GridIndexType>& operator() (CellCenterIdxType, CellCenterIdxType, const GridIndexType globalI) const
{
@@ -190,7 +184,7 @@ private:
{
stencil.push_back(scvf.axisData().selfDof);
stencil.push_back(scvf.axisData().oppositeDof);
- addHigherOrderInAxisDofs(std::integral_constant<bool, useHigherOrder>{}, scvf, stencil);
+ addHigherOrderInAxisDofs_(scvf, stencil, std::integral_constant<bool, useHigherOrder>{});
}
for(const auto& data : scvf.pairData())
@@ -201,31 +195,25 @@ private:
stencil.push_back(data.normalPair.second);
// add parallel dofs
- for (int i = 0; i < stencilOrder_ /*data.parallelDofs.size()*/; i++)
+ for (int i = 0; i < upwindSchemeOrder; i++)
{
if(!(data.parallelDofs[i] < 0))
- {
stencil.push_back(data.parallelDofs[i]);
- }
}
}
}
- void addHigherOrderInAxisDofs(std::false_type, const SubControlVolumeFace& scvf, Stencil& stencil) {}
+ void addHigherOrderInAxisDofs_(const SubControlVolumeFace& scvf, Stencil& stencil, std::false_type) {}
- void addHigherOrderInAxisDofs(std::true_type, const SubControlVolumeFace& scvf, Stencil& stencil)
+ void addHigherOrderInAxisDofs_(const SubControlVolumeFace& scvf, Stencil& stencil, std::true_type)
{
- for (int i = 0; i < stencilOrder_ - 1; i++)
+ for (int i = 0; i < upwindSchemeOrder - 1; i++)
{
if (scvf.hasBackwardNeighbor(i))
- {
stencil.push_back(scvf.axisData().inAxisBackwardDofs[i]);
- }
if (scvf.hasForwardNeighbor(i))
- {
stencil.push_back(scvf.axisData().inAxisForwardDofs[i]);
- }
}
}
@@ -234,7 +222,6 @@ private:
CellCenterToFaceMap cellCenterToFaceMap_;
FaceToCellCenterMap faceToCellCenterMap_;
FaceToFaceMap faceToFaceMap_;
- int stencilOrder_;
};
} // end namespace Dumux
diff --git a/dumux/discretization/staggered/freeflow/facevariables.hh b/dumux/discretization/staggered/freeflow/facevariables.hh
index 2f1943d9b8681707302d7256eaac675638582f7d..159018d08c95c92185f58858e8d6dce942458346 100644
--- a/dumux/discretization/staggered/freeflow/facevariables.hh
+++ b/dumux/discretization/staggered/freeflow/facevariables.hh
@@ -31,40 +31,37 @@ namespace Dumux {
namespace Detail {
-template<class Scalar, int geometryOrder>
-struct InAxisVelocities;
-
-template<class Scalar>
-struct InAxisVelocities<Scalar, 1>
+template<class Scalar, int upwindSchemeOrder>
+struct InAxisVelocities
{
Scalar self = 0.0;
Scalar opposite = 0.0;
+ std::array<Scalar, upwindSchemeOrder-1> forward{};
+ std::array<Scalar, upwindSchemeOrder-1> backward{};
};
-template<class Scalar, int geometryOrder>
-struct InAxisVelocities
+template<class Scalar>
+struct InAxisVelocities<Scalar, 1>
{
Scalar self = 0.0;
Scalar opposite = 0.0;
- std::array<Scalar, geometryOrder-1> forward{};
- std::array<Scalar, geometryOrder-1> backward{};
};
-}
+} // end namespace Detail
/*!
* \ingroup StaggeredDiscretization
* \brief The face variables class for free flow staggered grid models.
* Contains all relevant velocities for the assembly of the momentum balance.
*/
-template<class FacePrimaryVariables, int dim, int geometryOrder> // TODO doc
+template<class FacePrimaryVariables, int dim, int upwindSchemeOrder> // TODO doc
class StaggeredFaceVariables
{
static constexpr int numPairs = (dim == 2) ? 2 : 4;
- static constexpr bool useHigherOrder = geometryOrder > 1;
+ static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
using Scalar = typename FacePrimaryVariables::block_type;
- using InAxisVelocities = Detail::InAxisVelocities<Scalar, geometryOrder>;
+ using InAxisVelocities = Detail::InAxisVelocities<Scalar, upwindSchemeOrder>;
public:
@@ -99,7 +96,7 @@ public:
inAxisVelocities_.self = faceSol[scvf.dofIndex()];
inAxisVelocities_.opposite = faceSol[scvf.dofIndexOpposingFace()];
- addHigherOrderInAxisVelocities(std::integral_constant<bool, useHigherOrder>{}, faceSol, scvf);
+ addHigherOrderInAxisVelocities_(faceSol, scvf, std::integral_constant<bool, useHigherOrder>{});
// handle all sub faces
for (int i = 0; i < velocityParallel_.size(); ++i)
@@ -116,7 +113,7 @@ public:
velocityNormalOutside_[i] = faceSol[subFaceData.normalPair.second];
// treat the velocities parallel to the self face
- for (int j = 0; j < geometryOrder; j++)
+ for (int j = 0; j < upwindSchemeOrder; j++)
{
if (scvf.hasParallelNeighbor(i,j))
velocityParallel_[i][j] = faceSol[subFaceData.parallelDofs[j]];
@@ -196,10 +193,10 @@ public:
private:
template<class SolVector, class SubControlVolumeFace>
- void addHigherOrderInAxisVelocities(std::false_type, const SolVector& faceSol, const SubControlVolumeFace& scvf) {}
+ void addHigherOrderInAxisVelocities_(const SolVector& faceSol, const SubControlVolumeFace& scvf, std::false_type) {}
template<class SolVector, class SubControlVolumeFace>
- void addHigherOrderInAxisVelocities(std::true_type, const SolVector& faceSol, const SubControlVolumeFace& scvf)
+ void addHigherOrderInAxisVelocities_(const SolVector& faceSol, const SubControlVolumeFace& scvf, std::true_type)
{
// treat the velocity forward of the self face i.e. the face that is
@@ -220,7 +217,7 @@ private:
}
InAxisVelocities inAxisVelocities_;
- std::array<std::array<Scalar, geometryOrder>, numPairs> velocityParallel_;
+ std::array<std::array<Scalar, upwindSchemeOrder>, numPairs> velocityParallel_;
std::array<Scalar, numPairs> velocityNormalInside_;
std::array<Scalar, numPairs> velocityNormalOutside_;
diff --git a/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh b/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh
index c882c55a39329be5f23bc15cebc8817e44351c64..d07a0780e01fdc866ebf2e950ddd2ad4688a0b71 100644
--- a/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh
+++ b/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh
@@ -39,14 +39,15 @@ namespace Dumux {
* \ingroup StaggeredDiscretization
* \brief Default traits for the finite volume grid geometry.
*/
-template<class GridView, int geometryOrder, class MapperTraits = DefaultMapperTraits<GridView>>
+template<class GridView, int upwindSchemeOrder, class MapperTraits = DefaultMapperTraits<GridView>>
struct StaggeredFreeFlowDefaultFVGridGeometryTraits
: public MapperTraits
{
using SubControlVolume = CCSubControlVolume<GridView>;
- using SubControlVolumeFace = FreeFlowStaggeredSubControlVolumeFace<GridView, geometryOrder>;
+ using SubControlVolumeFace = FreeFlowStaggeredSubControlVolumeFace<GridView, upwindSchemeOrder>;
using IntersectionMapper = ConformingGridIntersectionMapper<GridView>;
- using GeometryHelper = FreeFlowStaggeredGeometryHelper<GridView, geometryOrder>;
+ using GeometryHelper = FreeFlowStaggeredGeometryHelper<GridView, upwindSchemeOrder>;
+ static constexpr int UpwindSchemeOrder = upwindSchemeOrder;
struct DofTypeIndices
{
diff --git a/dumux/discretization/staggered/freeflow/properties.hh b/dumux/discretization/staggered/freeflow/properties.hh
index 3c568c7d291cc30181e05c91cc85088ab0633eed..1588ca21e4783c0eeeee8c980c99fe339ccce426 100644
--- a/dumux/discretization/staggered/freeflow/properties.hh
+++ b/dumux/discretization/staggered/freeflow/properties.hh
@@ -81,10 +81,10 @@ template<class TypeTag>
struct FVGridGeometry<TypeTag, TTag::StaggeredFreeFlowModel>
{
private:
- static constexpr auto upwindGeometryOrder = getPropValue<TypeTag, Properties::UpwindGeometryOrder>();
+ static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
static constexpr bool enableCache = getPropValue<TypeTag, Properties::EnableFVGridGeometryCache>();
using GridView = GetPropType<TypeTag, Properties::GridView>;
- using Traits = StaggeredFreeFlowDefaultFVGridGeometryTraits<GridView, upwindGeometryOrder>;
+ using Traits = StaggeredFreeFlowDefaultFVGridGeometryTraits<GridView, upwindSchemeOrder>;
public:
using type = StaggeredFVGridGeometry<GridView, enableCache, Traits>;
};
@@ -96,9 +96,9 @@ struct FaceVariables<TypeTag, TTag::StaggeredFreeFlowModel>
private:
using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
- static constexpr auto upwindGeometryOrder = getPropValue<TypeTag, Properties::UpwindGeometryOrder>();
+ static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
public:
- using type = StaggeredFaceVariables<FacePrimaryVariables, GridView::dimension, upwindGeometryOrder>;
+ using type = StaggeredFaceVariables<FacePrimaryVariables, GridView::dimension, upwindSchemeOrder>;
};
//! Set the default global volume variables cache vector class
@@ -134,7 +134,7 @@ struct VelocityOutput<TypeTag, TTag::StaggeredFreeFlowModel>
* \brief Set the order of the upwinding scheme to 1 by default.
*/
template<class TypeTag>
-struct UpwindGeometryOrder<TypeTag, TTag::StaggeredFreeFlowModel> { static constexpr int value = 1; };
+struct UpwindSchemeOrder<TypeTag, TTag::StaggeredFreeFlowModel> { static constexpr int value = 1; };
} // namespace Properties
} // namespace Dumux
diff --git a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
index 0872330b0ed598713c0941f2a1e42da64242b8e1..09c568e7a904c736034d3f6442fa90d9d5f3a3a0 100644
--- a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
+++ b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
@@ -39,11 +39,11 @@ namespace Dumux
* \ingroup StaggeredDiscretization
* \brief Parallel Data stored per sub face
*/
-template<class Scalar, class GlobalPosition, int geometryOrder>
+template<class Scalar, class GlobalPosition, int upwindSchemeOrder>
struct PairData
{
- std::array<int, geometryOrder> parallelDofs;
- std::array<Scalar, geometryOrder+1/*geometryOrder*/> parallelDistances; // TODO
+ std::array<int, upwindSchemeOrder> parallelDofs;
+ std::array<Scalar, upwindSchemeOrder+1> parallelDistances; // TODO store only two distances, not three.
std::pair<signed int, signed int> normalPair;
int localNormalFaceIdx;
Scalar normalDistance;
@@ -55,13 +55,13 @@ struct PairData
* \ingroup StaggeredDiscretization
* \brief In Axis Data stored per sub face
*/
-template<class Scalar, int geometryOrder>
+template<class Scalar, int upwindSchemeOrder>
struct AxisData;
/*!
* \ingroup StaggeredDiscretization
* \brief In Axis Data stored per sub face
*/
-template<class Scalar, int geometryOrder>
+template<class Scalar, int upwindSchemeOrder>
struct AxisData
{
AxisData()
@@ -74,11 +74,11 @@ struct AxisData
int selfDof;
int oppositeDof;
- std::array<int, geometryOrder-1> inAxisForwardDofs;
- std::array<int, geometryOrder-1> inAxisBackwardDofs;
+ std::array<int, upwindSchemeOrder-1> inAxisForwardDofs;
+ std::array<int, upwindSchemeOrder-1> inAxisBackwardDofs;
Scalar selfToOppositeDistance;
- std::array<Scalar, geometryOrder-1> inAxisForwardDistances;
- std::array<Scalar, geometryOrder-1> inAxisBackwardDistances;
+ std::array<Scalar, upwindSchemeOrder-1> inAxisForwardDistances;
+ std::array<Scalar, upwindSchemeOrder-1> inAxisBackwardDistances;
};
/*!
@@ -93,7 +93,7 @@ struct AxisData<Scalar, 1>
Scalar selfToOppositeDistance;
};
-// template<class Scalar, int geometryOrder>
+// template<class Scalar, int upwindSchemeOrder>
// using AxisData
/*!
@@ -113,7 +113,7 @@ inline static unsigned int directionIndex(Vector&& vector)
* \brief Helper class constructing the dual grid finite volume geometries
* for the free flow staggered discretization method.
*/
-template<class GridView, int geometryOrder>
+template<class GridView, int upwindSchemeOrder>
class FreeFlowStaggeredGeometryHelper
{
using Scalar = typename GridView::ctype;
@@ -137,7 +137,7 @@ class FreeFlowStaggeredGeometryHelper
static constexpr int numfacets = dimWorld * 2;
static constexpr int numPairs = 2 * (dimWorld - 1);
- static constexpr bool useHigherOrder = geometryOrder > 1;
+ static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
public:
@@ -150,9 +150,8 @@ public:
void updateLocalFace(const IntersectionMapper& intersectionMapper_, const Intersection& intersection)
{
intersection_ = intersection;
- innerNormalFacePos_.clear();
- fillPairData_(std::integral_constant<bool, true>{}); //TODO
- fillAxisData_(std::integral_constant<bool, useHigherOrder>{});
+ fillAxisData_();
+ fillPairData_();
}
/*!
@@ -197,8 +196,11 @@ public:
private:
-
- void fillAxisData_(std::false_type)
+ /*!
+ * \brief Fills all entries of the in axis data
+ * Calls a function to extend the axis data for higher order upwind methods if required.
+ */
+ void fillAxisData_()
{
const auto inIdx = intersection_.indexInInside();
const auto oppIdx = localOppositeIdx_(inIdx);
@@ -213,29 +215,27 @@ private:
const auto self = getFacet_(inIdx, element_);
const auto opposite = getFacet_(oppIdx, element_);
axisData_.selfToOppositeDistance = (self.geometry().center() - opposite.geometry().center()).two_norm();
+
+ fillOuterAxisData_(std::integral_constant<bool, useHigherOrder>{});
}
+
/*!
- * \brief Fills all entries of the in axis data
+ * \brief Fills the axis data with the outer Dofs and Distances.
+ * For first order upwind methods no outer information is required.
*/
- void fillAxisData_(std::true_type)
+ void fillOuterAxisData_(std::false_type) {}
+
+ /*!
+ * \brief Fills the axis data with the outer dofs and distances.
+ * This function builds and extended stencil and is therefore only called when higher order upwinding methods are prescribed.
+ */
+ void fillOuterAxisData_(std::true_type)
{
const auto numForwardBackwardAxisDofs = axisData_.inAxisForwardDofs.size();
- const auto inIdx = intersection_.indexInInside();
- const auto oppIdx = localOppositeIdx_(inIdx);
-
- // Set the self Dof
- axisData_.selfDof = gridView_.indexSet().subIndex(intersection_.inside(), inIdx, codimIntersection);
-
- // Set the opposite Dof
- axisData_.oppositeDof = gridView_.indexSet().subIndex(intersection_.inside(), oppIdx, codimIntersection);
-
- // Set the Self to Opposite Distance
- const auto self = getFacet_(inIdx, element_);
- const auto opposite = getFacet_(oppIdx, element_);
- axisData_.selfToOppositeDistance = (self.geometry().center() - opposite.geometry().center()).two_norm();
// Set the Forward Dofs
std::stack<Element> inAxisForwardElementStack;
+ const auto inIdx = intersection_.indexInInside();
auto selfFace = getFacet_(inIdx, element_);
if(intersection_.neighbor())
{
@@ -266,26 +266,28 @@ private:
while(!inAxisForwardElementStack.empty())
{
int forwardIdx = inAxisForwardElementStack.size()-1;
- this->axisData_.inAxisForwardDofs[forwardIdx] = gridView_.indexSet().subIndex(inAxisForwardElementStack.top(), inIdx, codimIntersection);
+ axisData_.inAxisForwardDofs[forwardIdx] = gridView_.indexSet().subIndex(inAxisForwardElementStack.top(), inIdx, codimIntersection);
selfFace = getFacet_(inIdx, inAxisForwardElementStack.top());
forwardFaceCoordinates[forwardIdx] = selfFace.geometry().center();
inAxisForwardElementStack.pop();
}
+ const auto self = getFacet_(inIdx, element_);
for(int i = 0; i< forwardFaceCoordinates.size(); i++)
{
if(i == 0)
{
- this->axisData_.inAxisForwardDistances[i] = (forwardFaceCoordinates[i] - self.geometry().center()).two_norm();
+ axisData_.inAxisForwardDistances[i] = (forwardFaceCoordinates[i] - self.geometry().center()).two_norm();
}
else
{
- this->axisData_.inAxisForwardDistances[i] = (forwardFaceCoordinates[i]- forwardFaceCoordinates[i-1]).two_norm();
+ axisData_.inAxisForwardDistances[i] = (forwardFaceCoordinates[i]- forwardFaceCoordinates[i-1]).two_norm();
}
}
// Set the Backward Dofs
std::stack<Element> inAxisBackwardElementStack;
+ const auto oppIdx = localOppositeIdx_(inIdx);
auto oppFace = getFacet_(oppIdx, element_);
for(const auto& intersection : intersections(gridView_, element_))
{
@@ -329,6 +331,7 @@ private:
inAxisBackwardElementStack.pop();
}
+ const auto opposite = getFacet_(oppIdx, element_);
for(int i = 0; i< backwardFaceCoordinates.size(); i++)
{
if(i == 0)
@@ -342,33 +345,19 @@ private:
}
}
- void fillPairData_(std::false_type)
- {
-
- // TODO
-
- }
-
/*!
- * \brief Fills all entries of the pair data
+ * \brief Fills the pair data with the normal dofs and distances
+ * and calls a further function to collect the parallel dofs and distances
*/
- void fillPairData_(std::true_type)
+ void fillPairData_()
{
// initialize values that could remain unitialized if the intersection lies on a boundary
for(auto& data : pairData_)
{
- // parallel Dofs
- data.parallelDofs.fill(-1);
-
- // parallel Distances
- data.parallelDistances.fill(0.0);
-
// outer normals
data.normalPair.second = -1;
}
- unsigned int numParallelFaces = pairData_[0].parallelDistances.size();
-
// set basic global positions
const auto& elementCenter = this->element_.geometry().center();
const auto& FacetCenter = intersection_.geometry().center();
@@ -382,10 +371,6 @@ private:
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());
}
}
@@ -419,6 +404,76 @@ private:
}
}
+ fillParallelPairData_(std::integral_constant<bool, useHigherOrder>{});
+ }
+
+ /*!
+ * \brief Fills the pair data with the parallel dofs and distances
+ * This function is only called when the simple first order upwind methods are used.
+ */
+ void fillParallelPairData_(std::false_type)
+ {
+ // initialize values that could remain unitialized if the intersection lies on a boundary
+ for(auto& data : pairData_)
+ {
+ // parallel Dofs and Distances
+ data.parallelDofs.fill(-1);
+ data.parallelDistances.fill(0.0);
+ }
+
+ // set basic global positions and stencil size definitions
+ const auto& elementCenter = element_.geometry().center();
+
+ // get the parallel Dofs
+ const int parallelLocalIdx = intersection_.indexInInside();
+ int numPairParallelIdx = 0;
+ for(const auto& intersection : intersections(gridView_, element_))
+ {
+ if( facetIsNormal_(intersection.indexInInside(), parallelLocalIdx) )
+ {
+ // Store the parallel dimension of self cell in the direction of the axis
+ auto parallelAxisIdx = directionIndex(intersection);
+ pairData_[numPairParallelIdx].parallelDistances[0] = setParallelPairDistances_(element_, parallelAxisIdx);
+
+ if( intersection.neighbor() )
+ {
+ // If the normal intersection has a neighboring cell, go in and store the parallel information.
+ auto outerElement = intersection.outside();
+ pairData_[numPairParallelIdx].parallelDofs[0] = gridView_.indexSet().subIndex(outerElement, parallelLocalIdx, codimIntersection);
+ pairData_[numPairParallelIdx].parallelDistances[1] = setParallelPairDistances_(outerElement, parallelAxisIdx);
+ }
+ else // No parallel neighbor available
+ {
+ // If the intersection has no neighbor we have to deal with the virtual outer parallel dof
+ const auto& boundaryFacetCenter = intersection.geometry().center();
+ const auto distance = boundaryFacetCenter - elementCenter;
+ const auto virtualFirstParallelFaceDofPos = this->intersection_.geometry().center() + distance;
+
+ pairData_[numPairParallelIdx].virtualFirstParallelFaceDofPos = std::move(virtualFirstParallelFaceDofPos);
+ }
+ numPairParallelIdx++;
+ }
+ }
+ }
+
+ /*!
+ * \brief Fills the pair data with the parallel dofs and distances.
+ * This function builds and extended stencil and is therefore only called when higher order upwinding methods are prescribed.
+ */
+ void fillParallelPairData_(std::true_type)
+ {
+ // initialize values that could remain unitialized if the intersection lies on a boundary
+ for(auto& data : pairData_)
+ {
+ // parallel Dofs and Distances
+ data.parallelDofs.fill(-1);
+ data.parallelDistances.fill(0.0);
+ }
+
+ // set basic global positions and stencil size definitions
+ unsigned int numParallelFaces = pairData_[0].parallelDistances.size();
+ const auto& elementCenter = this->element_.geometry().center();
+
// get the parallel Dofs
const int parallelLocalIdx = intersection_.indexInInside();
int numPairParallelIdx = 0;
@@ -468,7 +523,6 @@ private:
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();
const auto distance = boundaryFacetCenter - elementCenter;
@@ -577,9 +631,8 @@ private:
const Element element_; //!< The respective element
const typename Element::Geometry elementGeometry_; //!< Reference to the element geometry
const GridView gridView_; //!< The grid view
- AxisData<Scalar, geometryOrder> axisData_; //!< Data related to forward and backward faces
- std::array<PairData<Scalar, GlobalPosition, geometryOrder>, numPairs> pairData_; //!< Collection of pair information related to normal and parallel faces
- std::vector<GlobalPosition> innerNormalFacePos_;
+ AxisData<Scalar, upwindSchemeOrder> axisData_; //!< Data related to forward and backward faces
+ std::array<PairData<Scalar, GlobalPosition, upwindSchemeOrder>, numPairs> pairData_; //!< Collection of pair information related to normal and parallel faces
};
} // end namespace Dumux
diff --git a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
index 8dcc9b93275a7cf3600785e17b11d22fb954e784..a8123938996a7c721e396414e7536e7f758333dc 100644
--- a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
@@ -42,12 +42,12 @@ namespace Dumux {
* of a sub control volume we compute fluxes on. This is a specialization for free flow models.
*/
template<class GV,
- int geometryOrder,
+ int upwindSchemeOrder,
class T = StaggeredDefaultScvfGeometryTraits<GV> >
class FreeFlowStaggeredSubControlVolumeFace
-: public SubControlVolumeFaceBase<FreeFlowStaggeredSubControlVolumeFace<GV, geometryOrder, T>, T>
+: public SubControlVolumeFaceBase<FreeFlowStaggeredSubControlVolumeFace<GV, upwindSchemeOrder, T>, T>
{
- using ThisType = FreeFlowStaggeredSubControlVolumeFace<GV, geometryOrder, T>;
+ using ThisType = FreeFlowStaggeredSubControlVolumeFace<GV, upwindSchemeOrder, T>;
using ParentType = SubControlVolumeFaceBase<ThisType, T>;
using Geometry = typename T::Geometry;
using GridIndexType = typename IndexTraits<GV>::GridIndex;
@@ -59,7 +59,7 @@ class FreeFlowStaggeredSubControlVolumeFace
static constexpr int numPairs = 2 * (dimworld - 1);
- static constexpr bool useHigherOrder = geometryOrder > 1;
+ static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
public:
using GlobalPosition = typename T::GlobalPosition;
@@ -207,7 +207,7 @@ public:
}
//! Returns the data for one sub face
- const PairData<Scalar, GlobalPosition, geometryOrder>& pairData(const int idx) const
+ const PairData<Scalar, GlobalPosition, upwindSchemeOrder>& pairData(const int idx) const
{
return pairData_[idx];
}
@@ -328,8 +328,8 @@ private:
int dofIdx_;
Scalar selfToOppositeDistance_;
- AxisData<Scalar, geometryOrder> axisData_;
- std::array<PairData<Scalar, GlobalPosition, geometryOrder>, numPairs> pairData_;
+ AxisData<Scalar, upwindSchemeOrder> axisData_;
+ std::array<PairData<Scalar, GlobalPosition, upwindSchemeOrder>, numPairs> pairData_;
int localFaceIdx_;
unsigned int dirIdx_;
diff --git a/dumux/discretization/staggered/fvgridgeometry.hh b/dumux/discretization/staggered/fvgridgeometry.hh
index 2b7e9b22f68f18f07a9b1d2ebc413ddc95f69d7b..c38909a6995cbd36e0103c682923f3e1a13344d4 100644
--- a/dumux/discretization/staggered/fvgridgeometry.hh
+++ b/dumux/discretization/staggered/fvgridgeometry.hh
@@ -157,8 +157,6 @@ public:
{ return this->fvGridGeometry_->numFaceDofs(); }
};
-
-
/*!
* \ingroup StaggeredDiscretization
* \brief Base class for the finite volume geometry vector for staggered models
@@ -193,6 +191,8 @@ class StaggeredFVGridGeometry<GV, true, Traits>
public:
//! export discretization method
static constexpr DiscretizationMethod discMethod = DiscretizationMethod::staggered;
+ static constexpr int upwindSchemeOrder = Traits::UpwindSchemeOrder;
+ static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
//! export the type of the fv element geometry (the local view type)
using LocalView = typename Traits::template LocalView<ThisType, true>;
@@ -213,6 +213,10 @@ public:
static constexpr auto faceIdx()
{ return typename DofTypeIndices::FaceIdx{}; }
+ //! The order of the stencil built
+ static constexpr std::size_t order()
+ { return upwindSchemeOrder; }
+
using CellCenterFVGridGeometryType = CellCenterFVGridGeometry<ThisType>;
using FaceFVGridGeometryType = FaceFVGridGeometry<ThisType>;
@@ -227,10 +231,6 @@ 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(paramGroup, "Discretization.TvdApproach"))
- stencilOrder_ = 2;
- else
- stencilOrder_ = 1;
}
//! The total number of sub control volumes
@@ -251,11 +251,6 @@ public:
return numBoundaryScvf_;
}
- //! The order of the stencil built
- std::size_t order() const
- {
- return stencilOrder_;
- }
//! The total number of intersections
std::size_t numIntersections() const
@@ -351,7 +346,6 @@ public:
}
// build the connectivity map for an effecient assembly
- connectivityMap_.setStencilOrder(stencilOrder_);
connectivityMap_.update(*this);
}
@@ -423,7 +417,6 @@ private:
// mappers
ConnectivityMap connectivityMap_;
IntersectionMapper intersectionMapper_;
- int stencilOrder_;
std::vector<SubControlVolume> scvs_;
std::vector<SubControlVolumeFace> scvfs_;
@@ -455,6 +448,8 @@ class StaggeredFVGridGeometry<GV, false, Traits>
public:
//! export discretization method
static constexpr DiscretizationMethod discMethod = DiscretizationMethod::staggered;
+ static constexpr int upwindSchemeOrder = Traits::UpwindSchemeOrder;
+ static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
using GeometryHelper = typename Traits::GeometryHelper;
@@ -477,6 +472,10 @@ public:
static constexpr auto faceIdx()
{ return typename DofTypeIndices::FaceIdx{}; }
+ //! The order of the stencil built
+ static constexpr std::size_t order()
+ { return upwindSchemeOrder; }
+
using CellCenterFVGridGeometryType = CellCenterFVGridGeometry<ThisType>;
using FaceFVGridGeometryType = FaceFVGridGeometry<ThisType>;
@@ -491,10 +490,6 @@ 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(paramGroup, "Discretization.TvdApproach"))
- stencilOrder_ = 2;
- else
- stencilOrder_ = 1;
}
//! update all fvElementGeometries (do this again after grid adaption)
@@ -547,7 +542,6 @@ public:
}
// build the connectivity map for an effecient assembly
- connectivityMap_.setStencilOrder(stencilOrder_);
connectivityMap_.update(*this);
}
@@ -563,12 +557,6 @@ public:
return numScvf_;
}
- //! The order of the stencil built
- std::size_t order() const
- {
- return stencilOrder_;
- }
-
//! The total number of boundary sub control volume faces
std::size_t numBoundaryScvf() const
{
@@ -652,7 +640,6 @@ private:
// mappers
ConnectivityMap connectivityMap_;
IntersectionMapper intersectionMapper_;
- int stencilOrder_;
//! vectors that store the global data
std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_;
diff --git a/dumux/discretization/staggered/gridfluxvariablescache.hh b/dumux/discretization/staggered/gridfluxvariablescache.hh
index 3c87372776d2c95922574fc8730da2ef0249097d..44d9550018bb3870a888bf77e61bf95077ed816d 100644
--- a/dumux/discretization/staggered/gridfluxvariablescache.hh
+++ b/dumux/discretization/staggered/gridfluxvariablescache.hh
@@ -28,7 +28,7 @@
#include <dumux/discretization/localview.hh>
#include <dumux/discretization/staggered/elementfluxvariablescache.hh>
-#include <dumux/freeflow/higherorderapproximation.hh>
+#include <dumux/freeflow/upwindingmethods.hh>
namespace Dumux {
@@ -86,7 +86,7 @@ public:
StaggeredGridFluxVariablesCache(const Problem& problem, const std::string& paramGroup = "")
: problemPtr_(&problem)
- , higherOrderApproximation_(paramGroup)
+ , upwindingMethods_(paramGroup)
{}
// When global caching is enabled, precompute transmissibilities and stencils for all the scv faces
@@ -123,10 +123,10 @@ public:
}
}
- //! Return the HigherOrderApproximation
- const HigherOrderApproximation<Scalar>& higherOrderApproximation() const
+ //! Return the UpwindingMethods
+ const UpwindingMethods<Scalar>& upwindingMethods() const
{
- return higherOrderApproximation_;
+ return upwindingMethods_;
}
const Problem& problem() const
@@ -141,7 +141,7 @@ public:
private:
const Problem* problemPtr_;
- HigherOrderApproximation<Scalar> higherOrderApproximation_;
+ UpwindingMethods<Scalar> upwindingMethods_;
std::vector<FluxVariablesCache> fluxVarsCache_;
std::vector<std::size_t> globalScvfIndices_;
@@ -174,7 +174,7 @@ public:
StaggeredGridFluxVariablesCache(const Problem& problem, const std::string& paramGroup = "")
: problemPtr_(&problem)
- , higherOrderApproximation_(paramGroup)
+ , upwindingMethods_(paramGroup)
{}
// When global caching is enabled, precompute transmissibilities and stencils for all the scv faces
@@ -187,15 +187,15 @@ public:
const Problem& problem() const
{ return *problemPtr_; }
- //! Return the HigherOrderApproximation
- const HigherOrderApproximation<Scalar>& higherOrderApproximation() const
+ //! Return the UpwindingMethods
+ const UpwindingMethods<Scalar>& upwindingMethods() const
{
- return higherOrderApproximation_;
+ return upwindingMethods_;
}
private:
const Problem* problemPtr_;
- HigherOrderApproximation<Scalar> higherOrderApproximation_;
+ UpwindingMethods<Scalar> upwindingMethods_;
};
diff --git a/dumux/freeflow/CMakeLists.txt b/dumux/freeflow/CMakeLists.txt
index 5f16de206e056526ae1a904c3c0261c80257a51b..d63c7b093702e9de86eea7d1dae52ade5f280b7a 100644
--- a/dumux/freeflow/CMakeLists.txt
+++ b/dumux/freeflow/CMakeLists.txt
@@ -7,5 +7,5 @@ install(FILES
properties.hh
turbulenceproperties.hh
volumevariables.hh
-higherorderapproximation.hh
+upwindingmethods.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow)
diff --git a/dumux/freeflow/navierstokes/staggered/CMakeLists.txt b/dumux/freeflow/navierstokes/staggered/CMakeLists.txt
index 5d59b30262acbbf8b4c049b917e591cdc1a74ba7..1f1987ff9f454cac7d604da5aa0205d1822f1077 100644
--- a/dumux/freeflow/navierstokes/staggered/CMakeLists.txt
+++ b/dumux/freeflow/navierstokes/staggered/CMakeLists.txt
@@ -2,4 +2,5 @@ install(FILES
fluxoversurface.hh
fluxvariables.hh
localresidual.hh
+staggeredupwindfluxvariables.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/navierstokes/staggered)
diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
index 37ad3fd3cb934e99cc37bc97477844e4ab205136..08d03621846d337fdc7ffa1e80b9e24a0b2fe4ec 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
@@ -31,11 +31,11 @@
#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>
+#include "upwindvariables.hh"
+
namespace Dumux {
// forward declaration
@@ -86,7 +86,8 @@ class NavierStokesFluxVariablesImpl<TypeTag, DiscretizationMethod::staggered>
static constexpr auto cellCenterIdx = FVGridGeometry::cellCenterIdx();
static constexpr auto faceIdx = FVGridGeometry::faceIdx();
- static constexpr bool useHigherOrder = false; // TODO
+ static constexpr int upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
+ static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
public:
@@ -203,8 +204,6 @@ public:
const Scalar velocitySelf = elemFaceVars[scvf].velocitySelf();
const Scalar velocityOpposite = elemFaceVars[scvf].velocityOpposite();
- // const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity);
-
// The volume variables within the current element. We only require those (and none of neighboring elements)
// because the fluxes are calculated over the staggered face at the center of the element.
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
@@ -214,42 +213,9 @@ public:
{
// Get the average velocity at the center of the element (i.e. the location of the staggered face).
const Scalar transportingVelocity = (velocitySelf + velocityOpposite) * 0.5;
-
- const auto upwindingMomenta = getFrontalUpwindingMomenta_(scvf, elemFaceVars, insideVolVars.density(), transportingVelocity);
- const Scalar momentum = doMomentumUpwinding_(scvf, upwindingMomenta, gridFluxVarsCache);
-
- // 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};
-
- // const auto& highOrder = gridFluxVarsCache.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]))
- // {
- // momentum = highOrder.tvd(velocities[0], velocities[1], velocities[2], distances[0], distances[1], distances[2], selfIsUpstream, insideVolVars.density(), highOrder.tvdApproach());
- // }
- // 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).
- frontalFlux += transportingVelocity * momentum * -1.0 * scvf.directionSign();
+ Dumux::UpwindVariables<TypeTag, upwindSchemeOrder> upwindVariables;
+ frontalFlux += upwindVariables.computeUpwindedFrontalMomentum(scvf, elemFaceVars, elemVolVars, gridFluxVarsCache, transportingVelocity)
+ * transportingVelocity * -1.0 * scvf.directionSign();
}
// Diffusive flux.
@@ -368,9 +334,15 @@ public:
// If there is no symmetry or Neumann boundary condition for the given sub face, proceed to calculate the tangential momentum flux.
if (problem.enableInertiaTerms())
- normalFlux += computeAdvectivePartOfLateralMomentumFlux_(problem, fvGeometry, element, scvf, normalFace, elemVolVars, faceVars, gridFluxVarsCache, localSubFaceIdx, lateralFaceHasDirichletPressure, lateralFaceHasBJS);
-
- normalFlux += computeDiffusivePartOfLateralMomentumFlux_(problem, fvGeometry, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx, lateralFaceHasDirichletPressure, lateralFaceHasBJS);
+ normalFlux += computeAdvectivePartOfLateralMomentumFlux_(problem, fvGeometry, element,
+ scvf, normalFace, elemVolVars, faceVars,
+ gridFluxVarsCache, localSubFaceIdx,
+ lateralFaceHasDirichletPressure, lateralFaceHasBJS);
+
+ normalFlux += computeDiffusivePartOfLateralMomentumFlux_(problem, fvGeometry, element,
+ scvf, normalFace, elemVolVars, faceVars,
+ localSubFaceIdx,
+ lateralFaceHasDirichletPressure, lateralFaceHasBJS);
}
return normalFlux;
}
@@ -414,49 +386,11 @@ private:
// of interest is located.
const Scalar transportingVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
- const auto parallelUpwindingMomenta = getParallelUpwindingMomenta_(problem, element, scvf, normalFace, elemVolVars, faceVars,
- transportingVelocity, localSubFaceIdx,
- lateralFaceHasDirichletPressure, lateralFaceHasBJS);
-
- const Scalar momentum = doMomentumUpwinding_(scvf, parallelUpwindingMomenta, gridFluxVarsCache);
-
- // 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};
-
- // 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};
-
- // const auto& highOrder = gridFluxVarsCache.higherOrderApproximation();
- //
- // // 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)
- // {
- // if (canLateralSecondOrder_(scvf, fvGeometry, selfIsUpstream, localSubFaceIdx, velocities, distances, problem, element, faceVars, lateralFaceHasDirichletPressure, lateralFaceHasBJS))
- // {
- // momentum = highOrder.tvd(velocities[0], velocities[1], velocities[2], distances[0], distances[1], distances[2], selfIsUpstream, selfIsUpstream ? insideVolVars.density() : outsideVolVars.density(), highOrder.tvdApproach());
- // }
- // 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);
- //
- // 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).
- return transportingVelocity * momentum * normalFace.directionSign() * normalFace.area() * 0.5 * extrusionFactor_(elemVolVars, normalFace);
+ Dumux::UpwindVariables<TypeTag, upwindSchemeOrder> upwindVariables;
+ return upwindVariables.computeUpwindedLateralMomentum(problem, fvGeometry, element, scvf, normalFace, elemVolVars, faceVars,
+ gridFluxVarsCache, localSubFaceIdx, lateralFaceHasDirichletPressure,
+ lateralFaceHasBJS)
+ * transportingVelocity * normalFace.directionSign() * normalFace.area() * 0.5 * extrusionFactor_(elemVolVars, normalFace);
}
/*!
@@ -543,7 +477,9 @@ private:
const Scalar velocityFirstParallel = scvf.hasParallelNeighbor(localSubFaceIdx,0)
? faceVars.velocityParallel(localSubFaceIdx,0)
- : getParallelVelocityFromBoundary_(problem, scvf, normalFace, innerParallelVelocity, localSubFaceIdx, element, false, lateralFaceHasBJS);
+ : 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.
@@ -609,7 +545,8 @@ private:
//! helper function to conveniently create a ghost face used to retrieve boundary values from the problem
SubControlVolumeFace makeGhostFace_(const SubControlVolumeFace& ownScvf, const GlobalPosition& pos) const
{
- return SubControlVolumeFace(pos, std::vector<unsigned int>{ownScvf.insideScvIdx(), ownScvf.outsideScvIdx()}, ownScvf.directionIndex(), ownScvf.axisData().selfDof, ownScvf.index());
+ return SubControlVolumeFace(pos, std::vector<unsigned int>{ownScvf.insideScvIdx(), ownScvf.outsideScvIdx()},
+ ownScvf.directionIndex(), ownScvf.axisData().selfDof, ownScvf.index());
};
//! helper function to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
@@ -652,186 +589,6 @@ private:
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 FVElementGeometry& fvGeometry,
- 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 = fvGeometry.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, fvGeometry, 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 = fvGeometry.fvGridGeometry().element(fvGeometry.scv(normalFace.outsideScvIdx()));
- const SubControlVolumeFace& firstParallelScvf = fvGeometry.scvf(normalFace.outsideScvIdx(), ownScvf.localFaceIdx());
- velocities[2] = getParallelVelocityFromOtherBoundary_(problem, fvGeometry, 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.
*
@@ -877,58 +634,38 @@ private:
DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the boundary conditions for the momentum equations at global position " << boundarySubFaceCenter);
}
}
-
- Scalar doMomentumUpwinding_(const SubControlVolumeFace& scvf,
- const std::array<Scalar, 2>& momentum,
- const GridFluxVariablesCache& gridFluxVarsCache) const
- {
- const auto& upwindScheme = gridFluxVarsCache.higherOrderApproximation();
- return upwindScheme.upwind(momentum[0], momentum[1]);
- }
-
- std::array<Scalar, 2> getFrontalUpwindingMomenta_(const SubControlVolumeFace& scvf,
- const ElementFaceVariables& elemFaceVars,
- const Scalar density,
- const Scalar transportingVelocity) const
- {
- const Scalar momentumSelf = elemFaceVars[scvf].velocitySelf() * density;
- const Scalar momentumOpposite = elemFaceVars[scvf].velocityOpposite() * density;
- const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity);
-
- return selfIsUpstream ? std::array<Scalar, 2>{momentumOpposite, momentumSelf}
- : std::array<Scalar, 2>{momentumSelf, momentumOpposite};
- }
-
- std::array<Scalar, 2> getParallelUpwindingMomenta_(const Problem& problem,
- const Element& element,
- const SubControlVolumeFace& ownFace,
- const SubControlVolumeFace& lateralFace,
- const ElementVolumeVariables& elemVolVars,
- const FaceVariables& faceVars,
- const Scalar transportingVelocity,
- const int localSubFaceIdx,
- bool lateralFaceHasDirichletPressure,
- bool lateralFaceHasBJS) const
+ /*!
+ * \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;
- // Check whether the own or the neighboring element is upstream.
- const bool selfIsUpstream = lateralFace.directionSign() == sign(transportingVelocity);
-
- // Get the volume variables of the own and the neighboring element
- const auto& insideVolVars = elemVolVars[lateralFace.insideScvIdx()];
- const auto& outsideVolVars = elemVolVars[lateralFace.outsideScvIdx()];
-
- const Scalar momentumSelf = faceVars.velocitySelf() * insideVolVars.density();
-
- const Scalar momentumParallel = ownFace.hasParallelNeighbor(localSubFaceIdx, 0)
- ? faceVars.velocityParallel(localSubFaceIdx, 0) * outsideVolVars.density()
- : (getParallelVelocityFromBoundary_(problem, ownFace, lateralFace,
- faceVars.velocitySelf(), localSubFaceIdx, element,
- lateralFaceHasDirichletPressure, lateralFaceHasBJS)
- * insideVolVars.density());
-
- return selfIsUpstream ? std::array<Scalar, 2>{momentumParallel, momentumSelf}
- : std::array<Scalar, 2>{momentumSelf, momentumParallel};
+ 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())];
}
};
diff --git a/dumux/freeflow/navierstokes/staggered/upwindvariables.hh b/dumux/freeflow/navierstokes/staggered/upwindvariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b5a7bc6d05dccd55d120c14c2b3e851ef5cd269e
--- /dev/null
+++ b/dumux/freeflow/navierstokes/staggered/upwindvariables.hh
@@ -0,0 +1,656 @@
+// -*- 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 3 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
+ * \ingroup NavierStokesModel
+ * \copydoc Dumux::NavierStokesUpwindVariables
+ */
+#ifndef DUMUX_NAVIERSTOKES_STAGGERED_UPWINDVARIABLES_HH
+#define DUMUX_NAVIERSTOKES_STAGGERED_UPWINDVARIABLES_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/discretization/method.hh>
+#include <dumux/freeflow/upwindingmethods.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief The upwinding variables class for the Navier-Stokes model using the staggered grid discretization.
+ */
+template<class TypeTag, int upwindSchemeOrder>
+class UpwindVariables
+{
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+
+ using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
+ using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
+ using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
+
+ using GridFluxVariablesCache = typename GridVariables::GridFluxVariablesCache;
+ using FluxVariablesCache = typename GridFluxVariablesCache::FluxVariablesCache;
+
+ using GridFaceVariables = typename GridVariables::GridFaceVariables;
+ using ElementFaceVariables = typename GridFaceVariables::LocalView;
+ using FaceVariables = typename GridFaceVariables::FaceVariables;
+
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using GridView = typename FVGridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+ using Indices = typename ModelTraits::Indices;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
+ using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
+
+ static constexpr auto cellCenterIdx = FVGridGeometry::cellCenterIdx();
+ static constexpr auto faceIdx = FVGridGeometry::faceIdx();
+
+public:
+ UpwindVariables() {};
+
+ /*!
+ * \brief
+ *
+ * TODO: DOC
+ *
+ */
+ FacePrimaryVariables computeUpwindedFrontalMomentum(const SubControlVolumeFace& scvf,
+ const ElementFaceVariables& elemFaceVars,
+ const ElementVolumeVariables& elemVolVars,
+ const GridFluxVariablesCache& gridFluxVarsCache,
+ const Scalar transportingVelocity)
+ {
+ Scalar momentum(0.0);
+ const bool canHigherOrder = canFrontalSecondOrder_(scvf, transportingVelocity, std::integral_constant<bool, useHigherOrder>{});
+
+ if (canHigherOrder)
+ {
+ const auto upwindingMomenta = getFrontalUpwindingMomenta_(scvf, elemFaceVars, elemVolVars[scvf.insideScvIdx()].density(),
+ transportingVelocity, std::integral_constant<bool, useHigherOrder>{});
+ momentum = doFrontalMomentumUpwinding_(scvf, upwindingMomenta, transportingVelocity, gridFluxVarsCache);
+ }
+ else
+ {
+ const auto upwindingMomenta = getFrontalUpwindingMomenta_(scvf, elemFaceVars, elemVolVars[scvf.insideScvIdx()].density(),
+ transportingVelocity, std::integral_constant<bool, false>{});
+ momentum = doFrontalMomentumUpwinding_(scvf, upwindingMomenta, transportingVelocity, gridFluxVarsCache);
+ }
+
+ return momentum;
+ }
+
+ /*!
+ * \brief
+ *
+ * TODO: DOC
+ *
+ */
+ FacePrimaryVariables computeUpwindedLateralMomentum(const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const Element& element,
+ const SubControlVolumeFace& scvf,
+ const SubControlVolumeFace& normalFace,
+ const ElementVolumeVariables& elemVolVars,
+ const FaceVariables& faceVars,
+ const GridFluxVariablesCache& gridFluxVarsCache,
+ const int localSubFaceIdx,
+ const bool lateralFaceHasDirichletPressure,
+ const bool lateralFaceHasBJS)
+ {
+ // Get the transporting velocity, located at the scvf perpendicular to the current scvf where the dof
+ // of interest is located.
+ const Scalar transportingVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
+
+
+ // Check whether the own or the neighboring element is upstream.
+ const bool selfIsUpstream = ( normalFace.directionSign() == sign(transportingVelocity) );
+
+ Scalar momentum(0.0);
+ const bool canHigherOrder = canLateralSecondOrder_(scvf, selfIsUpstream, localSubFaceIdx, std::integral_constant<bool, useHigherOrder>{});
+
+ if (canHigherOrder)
+ {
+ const auto parallelUpwindingMomenta = getLateralUpwindingMomenta_(problem, fvGeometry, element, scvf, elemVolVars, faceVars,
+ transportingVelocity, localSubFaceIdx,
+ lateralFaceHasDirichletPressure, lateralFaceHasBJS,
+ std::integral_constant<bool, useHigherOrder>{});
+ momentum = doLateralMomentumUpwinding_(scvf, normalFace, parallelUpwindingMomenta, transportingVelocity, localSubFaceIdx, gridFluxVarsCache);
+ }
+ else
+ {
+ const auto parallelUpwindingMomenta = getLateralUpwindingMomenta_(problem, fvGeometry, element, scvf, elemVolVars, faceVars,
+ transportingVelocity, localSubFaceIdx,
+ lateralFaceHasDirichletPressure, lateralFaceHasBJS,
+ std::integral_constant<bool, false>{});
+ momentum = doLateralMomentumUpwinding_(scvf, normalFace, parallelUpwindingMomenta, transportingVelocity, localSubFaceIdx, gridFluxVarsCache);
+ }
+
+ return momentum;
+ }
+
+private:
+
+ /*!
+ * \brief
+ *
+ * TODO: DOC
+ *
+ */
+ bool canFrontalSecondOrder_(const SubControlVolumeFace& ownScvf,
+ const Scalar transportingVelocity,
+ std::true_type) const
+ {
+ const bool selfIsUpstream = ownScvf.directionSign() != sign(transportingVelocity);
+ // Depending on selfIsUpstream I have to check if I have a forward or a backward neighbor to retrieve
+ if (selfIsUpstream)
+ {
+ if (ownScvf.hasForwardNeighbor(0))
+ return true;
+ else
+ return false;
+ }
+ else
+ {
+ if (ownScvf.hasBackwardNeighbor(0))
+ return true;
+ else
+ return false;
+ }
+ }
+
+ /*!
+ * \brief
+ *
+ * TODO: DOC
+ *
+ */
+ bool canFrontalSecondOrder_(const SubControlVolumeFace& ownScvf,
+ const Scalar transportingVelocity,
+ std::false_type) const
+ {
+ return false;
+ }
+
+ /*!
+ * \brief
+ *
+ * TODO: DOC
+ *
+ */
+ std::array<Scalar, 3> getFrontalUpwindingMomenta_(const SubControlVolumeFace& scvf,
+ const ElementFaceVariables& elemFaceVars,
+ const Scalar density,
+ const Scalar transportingVelocity,
+ std::true_type) const
+ {
+ const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity);
+ std::array<Scalar, 3> momenta{0.0, 0.0, 0.0};
+
+ if (selfIsUpstream)
+ {
+ momenta[0] = elemFaceVars[scvf].velocityOpposite() * density;
+ momenta[1] = elemFaceVars[scvf].velocitySelf() * density;
+ momenta[2] = elemFaceVars[scvf].velocityForward(upwindSchemeOrder-2) * density;
+ }
+ else
+ {
+ momenta[0] = elemFaceVars[scvf].velocitySelf() * density;
+ momenta[1] = elemFaceVars[scvf].velocityOpposite() * density;
+ momenta[2] = elemFaceVars[scvf].velocityBackward(upwindSchemeOrder-2) * density;
+ }
+
+ return momenta;
+ }
+
+ /*!
+ * \brief
+ *
+ * TODO: DOC
+ *
+ */
+ std::array<Scalar, 2> getFrontalUpwindingMomenta_(const SubControlVolumeFace& scvf,
+ const ElementFaceVariables& elemFaceVars,
+ const Scalar density,
+ const Scalar transportingVelocity,
+ std::false_type) const
+ {
+ const Scalar momentumSelf = elemFaceVars[scvf].velocitySelf() * density;
+ const Scalar momentumOpposite = elemFaceVars[scvf].velocityOpposite() * density;
+ const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity);
+
+ return selfIsUpstream ? std::array<Scalar, 2>{momentumOpposite, momentumSelf}
+ : std::array<Scalar, 2>{momentumSelf, momentumOpposite};
+ }
+
+ /*!
+ * \brief
+ *
+ * TODO: DOC
+ *
+ */
+ Scalar doFrontalMomentumUpwinding_(const SubControlVolumeFace& scvf,
+ const std::array<Scalar, 2>& momenta,
+ const Scalar transportingVelocity,
+ const GridFluxVariablesCache& gridFluxVarsCache) const
+ {
+ const auto& upwindScheme = gridFluxVarsCache.upwindingMethods();
+ return upwindScheme.upwind(momenta[0], momenta[1]);
+ }
+
+ /*!
+ * \brief
+ *
+ * TODO: DOC
+ *
+ */
+ template<bool enable = useHigherOrder, std::enable_if_t<enable, int> = 0>
+ Scalar doFrontalMomentumUpwinding_(const SubControlVolumeFace& scvf,
+ const std::array<Scalar, 3>& momenta,
+ const Scalar transportingVelocity,
+ const GridFluxVariablesCache& gridFluxVarsCache) const
+ {
+ const auto& upwindScheme = gridFluxVarsCache.upwindingMethods();
+ const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity);
+ std::array<Scalar,3> distances = getFrontalDistances_(scvf, selfIsUpstream);
+ return upwindScheme.tvd(momenta, distances, selfIsUpstream, upwindScheme.tvdApproach());
+ }
+
+ /*!
+ * \brief
+ *
+ * TODO: DOC
+ *
+ */
+ template<bool enable = useHigherOrder, std::enable_if_t<enable, int> = 0>
+ std::array<Scalar, 3> getFrontalDistances_(const SubControlVolumeFace& ownScvf,
+ const bool selfIsUpstream) const
+ {
+ // Depending on selfIsUpstream the downstream and the (up)upstream distances are saved.
+ // distances {upstream to downstream distance, up-upstream to upstream distance, downstream staggered cell size}
+ std::array<Scalar, 3> distances{0.0, 0.0, 0.0};
+
+ if (selfIsUpstream)
+ {
+ distances[0] = ownScvf.selfToOppositeDistance();
+ distances[1] = ownScvf.axisData().inAxisForwardDistances[upwindSchemeOrder-2];
+ distances[2] = 0.5 * (ownScvf.axisData().selfToOppositeDistance + ownScvf.axisData().inAxisBackwardDistances[0]);
+ }
+ else
+ {
+ distances[0] = ownScvf.selfToOppositeDistance();
+ distances[1] = ownScvf.axisData().inAxisBackwardDistances[upwindSchemeOrder-2];
+ distances[2] = 0.5 * (ownScvf.axisData().selfToOppositeDistance + ownScvf.axisData().inAxisForwardDistances[0]);
+ }
+ return distances;
+ }
+
+ /*!
+ * TODO: DOC
+ *
+ * \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 at ownScvf is upstream wrt the transporting velocity
+ * \param localSubFaceIdx The local subface index
+ */
+ bool canLateralSecondOrder_(const SubControlVolumeFace& ownScvf,
+ const bool selfIsUpstream,
+ const int localSubFaceIdx,
+ std::true_type) const
+ {
+ if (selfIsUpstream)
+ {
+ // The self velocity is upstream. The downstream velocity can be assigned or retrieved
+ // from the boundary, even if there is no parallel neighbor.
+ return true;
+ }
+ else
+ {
+ // The self velocity is downstream. If there is no parallel neighbor I cannot use a second order approximation.
+ if (!ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
+ return false;
+ else
+ return true;
+ }
+ }
+
+ /*!
+ * \brief Check if a second order approximation for the lateral part of the advective term can be used
+ * If higher order methods are not prescribed, this is called and false is returned.
+ */
+ bool canLateralSecondOrder_(const SubControlVolumeFace& ownScvf,
+ const bool selfIsUpstream,
+ const int localSubFaceIdx,
+ std::false_type) const
+ { return false; }
+
+
+ /*!
+ * \brief
+ *
+ * TODO: DOC
+ *
+ */
+ std::array<Scalar, 3> getLateralUpwindingMomenta_(const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const Element& element,
+ const SubControlVolumeFace& ownScvf,
+ const ElementVolumeVariables& elemVolVars,
+ const FaceVariables& faceVars,
+ const Scalar transportingVelocity,
+ const int localSubFaceIdx,
+ bool lateralFaceHasDirichletPressure,
+ bool lateralFaceHasBJS,
+ std::true_type) const
+ {
+ std::array<Scalar, 3> momenta{0.0, 0.0, 0.0};
+ const SubControlVolumeFace& lateralFace = fvGeometry.scvf(ownScvf.insideScvIdx(), ownScvf.pairData(localSubFaceIdx).localNormalFaceIdx);
+
+ // Check whether the own or the neighboring element is upstream.
+ const bool selfIsUpstream = lateralFace.directionSign() == sign(transportingVelocity);
+
+ // Get the volume variables of the own and the neighboring element
+ const auto& insideVolVars = elemVolVars[lateralFace.insideScvIdx()];
+ const auto& outsideVolVars = elemVolVars[lateralFace.outsideScvIdx()];
+
+ // The local index of the faces that is opposite to localSubFaceIdx
+ const int oppositeSubFaceIdx = localSubFaceIdx % 2 ? localSubFaceIdx - 1 : localSubFaceIdx + 1;
+
+ if (selfIsUpstream)
+ {
+ momenta[1] = faceVars.velocitySelf() * insideVolVars.density();
+
+ if(ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
+ momenta[0] = faceVars.velocityParallel(localSubFaceIdx, 0) * insideVolVars.density();
+ else
+ momenta[0] = getParallelVelocityFromBoundary_(problem, ownScvf, lateralFace,
+ faceVars.velocitySelf(), localSubFaceIdx,
+ element, lateralFaceHasDirichletPressure,
+ lateralFaceHasBJS) * insideVolVars.density();
+
+ // The "upstream-upstream" velocity is retrieved from the other parallel neighbor or from the boundary.
+ if (ownScvf.hasParallelNeighbor(oppositeSubFaceIdx, 0))
+ momenta[2] = faceVars.velocityParallel(oppositeSubFaceIdx, 0) * insideVolVars.density();
+ else
+ momenta[2] = getParallelVelocityFromOtherBoundary_(problem, fvGeometry, ownScvf,
+ oppositeSubFaceIdx, element,
+ (momenta[1]/insideVolVars.density())) * insideVolVars.density();
+ }
+ else
+ {
+ momenta[0] = faceVars.velocitySelf() * outsideVolVars.density();
+
+ if (!ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
+ {
+ std::cout << "how did we get here \n";
+ momenta[1] = getParallelVelocityFromBoundary_(problem, ownScvf, lateralFace,
+ faceVars.velocitySelf(), localSubFaceIdx,
+ element, lateralFaceHasDirichletPressure,
+ lateralFaceHasBJS) * outsideVolVars.density();
+ return momenta;
+ }
+
+ momenta[1] = faceVars.velocityParallel(localSubFaceIdx, 0) * outsideVolVars.density();
+
+ // 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))
+ momenta[2] = faceVars.velocityParallel(localSubFaceIdx, 1) * outsideVolVars.density();
+ else
+ {
+ const Element& elementParallel = fvGeometry.fvGridGeometry().element(fvGeometry.scv(lateralFace.outsideScvIdx()));
+ const SubControlVolumeFace& firstParallelScvf = fvGeometry.scvf(lateralFace.outsideScvIdx(), ownScvf.localFaceIdx());
+ momenta[2] = getParallelVelocityFromOtherBoundary_(problem, fvGeometry, firstParallelScvf,
+ localSubFaceIdx, elementParallel,
+ (momenta[1]/outsideVolVars.density())) * outsideVolVars.density();
+ }
+ }
+ return momenta;
+
+ }
+
+ /*!
+ * \brief
+ *
+ * TODO: DOC
+ *
+ */
+ std::array<Scalar, 2> getLateralUpwindingMomenta_(const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const Element& element,
+ const SubControlVolumeFace& ownScvf,
+ const ElementVolumeVariables& elemVolVars,
+ const FaceVariables& faceVars,
+ const Scalar transportingVelocity,
+ const int localSubFaceIdx,
+ bool lateralFaceHasDirichletPressure,
+ bool lateralFaceHasBJS,
+ std::false_type) const
+ {
+ // Check whether the own or the neighboring element is upstream.
+ const SubControlVolumeFace& lateralFace = fvGeometry.scvf(ownScvf.insideScvIdx(), ownScvf.pairData(localSubFaceIdx).localNormalFaceIdx);
+ const bool selfIsUpstream = lateralFace.directionSign() == sign(transportingVelocity);
+
+ // Get the volume variables of the own and the neighboring element
+ const auto& insideVolVars = elemVolVars[lateralFace.insideScvIdx()];
+ const auto& outsideVolVars = elemVolVars[lateralFace.outsideScvIdx()];
+
+ const Scalar momentumSelf = faceVars.velocitySelf() * insideVolVars.density();
+
+ const Scalar momentumParallel = ownScvf.hasParallelNeighbor(localSubFaceIdx, 0)
+ ? faceVars.velocityParallel(localSubFaceIdx, 0) * outsideVolVars.density()
+ : (getParallelVelocityFromBoundary_(problem, ownScvf, lateralFace,
+ faceVars.velocitySelf(), localSubFaceIdx, element,
+ lateralFaceHasDirichletPressure, lateralFaceHasBJS)
+ * insideVolVars.density());
+
+ return selfIsUpstream ? std::array<Scalar, 2>{momentumParallel, momentumSelf}
+ : std::array<Scalar, 2>{momentumSelf, momentumParallel};
+ }
+
+ /*!
+ * \brief
+ *
+ * TODO: DOC
+ *
+ */
+ Scalar doLateralMomentumUpwinding_(const SubControlVolumeFace& scvf,
+ const SubControlVolumeFace& normalScvf,
+ const std::array<Scalar, 2>& momenta,
+ const Scalar transportingVelocity,
+ const int localSubFaceIdx,
+ const GridFluxVariablesCache& gridFluxVarsCache) const
+ {
+ return doFrontalMomentumUpwinding_(scvf, momenta, transportingVelocity, gridFluxVarsCache);
+ }
+
+ /*!
+ * \brief
+ *
+ * TODO: DOC
+ *
+ * \param scvf The sub control volume face
+ * \param normalScvf The normal sub control volume face
+ * \param momenta The momenta to be upwinded
+ * \param transportingVelocity The average of the self and opposite velocities.
+ * \param localSubFaceIdx The local index of the subface
+ * \param gridFluxVarsCache The grid flux variables cache
+ */
+ Scalar doLateralMomentumUpwinding_(const SubControlVolumeFace& scvf,
+ const SubControlVolumeFace& normalScvf,
+ const std::array<Scalar, 3>& momenta,
+ const Scalar transportingVelocity,
+ const int localSubFaceIdx,
+ const GridFluxVariablesCache& gridFluxVarsCache) const
+ {
+ const bool selfIsUpstream = ( normalScvf.directionSign() == sign(transportingVelocity) );
+ std::array<Scalar,3> distances = getLateralDistances_(scvf, localSubFaceIdx, selfIsUpstream);
+ return upwindScheme.tvd(momenta, distances, selfIsUpstream, upwindScheme.tvdApproach());
+ }
+
+ /*!
+ * \brief
+ *
+ * TODO: DOC
+ *
+ */
+ template<bool enable = useHigherOrder, std::enable_if_t<enable, int> = 0>
+ std::array<Scalar, 3> getLateralDistances_(const SubControlVolumeFace& ownScvf,
+ const int localSubFaceIdx,
+ const bool selfIsUpstream) const
+ {
+ // distances {upstream to downstream distance, up-upstream to upstream distance, downstream staggered cell size}
+ std::array<Scalar, 3> distances{0.0, 0.0, 0.0};
+
+ // The local index of the faces that is opposite to localSubFaceIdx
+ const int oppositeSubFaceIdx = localSubFaceIdx % 2 ? localSubFaceIdx - 1 : localSubFaceIdx + 1;
+
+ if (selfIsUpstream)
+ {
+ distances[0] = ownScvf.cellCenteredParallelDistance(localSubFaceIdx, 0);
+ distances[1] = ownScvf.cellCenteredParallelDistance(oppositeSubFaceIdx, 0);
+ if(ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
+ distances[2] = ownScvf.pairData(localSubFaceIdx).parallelDistances[1];
+ else
+ distances[2] = ownScvf.area() / 2.0;
+ }
+ else
+ {
+ distances[0] = ownScvf.cellCenteredParallelDistance(localSubFaceIdx, 0);
+ distances[1] = ownScvf.cellCenteredParallelDistance(localSubFaceIdx, 1);
+ distances[2] = ownScvf.area();
+ }
+
+ return distances;
+ }
+
+ /*!
+ * \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 FVElementGeometry& fvGeometry,
+ 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 = fvGeometry.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);
+ }
+ }
+
+
+ //! helper function to conveniently create a ghost face used to retrieve boundary values from the problem
+ SubControlVolumeFace makeGhostFace_(const SubControlVolumeFace& ownScvf, const GlobalPosition& pos) const
+ {
+ return SubControlVolumeFace(pos, std::vector<unsigned int>{ownScvf.insideScvIdx(), ownScvf.outsideScvIdx()},
+ ownScvf.directionIndex(), ownScvf.axisData().selfDof, ownScvf.index());
+ };
+
+ //! 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).virtualFirstParallelFaceDofPos);
+ };
+};
+
+} // end namespace Dumux
+
+#endif // DUMUX_NAVIERSTOKES_STAGGERED_UPWINDVARIABLES_HH
diff --git a/dumux/freeflow/higherorderapproximation.hh b/dumux/freeflow/upwindingmethods.hh
similarity index 66%
rename from dumux/freeflow/higherorderapproximation.hh
rename to dumux/freeflow/upwindingmethods.hh
index 3e983074159dd35546353aafa47b63029e281bb4..18f985de6ae9d96814bb2dbe7eae4f3f4d4bdae1 100644
--- a/dumux/freeflow/higherorderapproximation.hh
+++ b/dumux/freeflow/upwindingmethods.hh
@@ -20,8 +20,8 @@
* \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
+#ifndef DUMUX_UPWINDING_METHODS_HH
+#define DUMUX_UPWINDING_METHODS_HH
#include <cmath>
#include <functional>
@@ -48,18 +48,18 @@ enum class DifferencingScheme
* \brief This file contains different higher order methods for approximating the velocity.
*/
template<class Scalar>
-class HigherOrderApproximation
+class UpwindingMethods
{
public:
- HigherOrderApproximation(const std::string& paramGroup = "")
+ UpwindingMethods(const std::string& paramGroup = "")
{
upwindWeight_ = getParamFromGroup<Scalar>(paramGroup, "Flux.UpwindWeight");
- if (hasParamInGroup(paramGroup, "Discretization.TvdApproach"))
+ if (hasParamInGroup(paramGroup, "Flux.TvdApproach"))
{
// Read the runtime parameters
- tvdApproach_ = tvdApproachFromString(getParamFromGroup<std::string>(paramGroup, "Discretization.TvdApproach"));
- differencingScheme_ = differencingSchemeFromString(getParamFromGroup<std::string>(paramGroup, "Discretization.DifferencingScheme"));
+ tvdApproach_ = tvdApproachFromString(getParamFromGroup<std::string>(paramGroup, "Flux.TvdApproach"));
+ differencingScheme_ = differencingSchemeFromString(getParamFromGroup<std::string>(paramGroup, "Flux.DifferencingScheme"));
// Assign the limiter_ depending on the differencing scheme
switch (differencingScheme_)
@@ -116,6 +116,8 @@ public:
tvdApproach_ = TvdApproach::none;
differencingScheme_ = DifferencingScheme::none;
}
+ std::cout << "Using the tvdApproach " << tvdApproachToString(tvdApproach_)
+ << " and the differencing Scheme " << differencingSchemeToString(differencingScheme_) << "\n";
}
/**
@@ -129,9 +131,9 @@ public:
if (tvd == "Hou") return TvdApproach::hou;
DUNE_THROW(ParameterException, "\nThis tvd approach : \"" << tvd << "\" is not implemented.\n"
<< "The available TVD approaches for uniform (and nonuniform) grids are as follows: \n"
- << tvdApproachToString(TvdApproach::uniform) << ": Assumes a uniform cell size distribution\n"
- << tvdApproachToString(TvdApproach::li) << ": Li's approach for nonuniform cell sizes\n"
- << tvdApproachToString(TvdApproach::hou) << ": Hou's approach for nonuniform cell sizes");
+ << tvdApproachToString(TvdApproach::uniform) << ": assumes a uniform cell size distribution\n"
+ << tvdApproachToString(TvdApproach::li) << ": li's approach for nonuniform cell sizes\n"
+ << tvdApproachToString(TvdApproach::hou) << ": hou's approach for nonuniform cell sizes");
}
/**
@@ -154,13 +156,13 @@ public:
*/
DifferencingScheme differencingSchemeFromString(const std::string& differencingScheme)
{
- if (differencingScheme == "vanleer") return DifferencingScheme::vanleer;
- if (differencingScheme == "vanalbada") return DifferencingScheme::vanalbada;
- if (differencingScheme == "minmod") return DifferencingScheme::minmod;
- if (differencingScheme == "superbee") return DifferencingScheme::superbee;
- if (differencingScheme == "umist") return DifferencingScheme::umist;
- if (differencingScheme == "mclimiter") return DifferencingScheme::mclimiter;
- if (differencingScheme == "wahyd") return DifferencingScheme::wahyd;
+ if (differencingScheme == "Vanleer") return DifferencingScheme::vanleer;
+ if (differencingScheme == "Vanalbada") return DifferencingScheme::vanalbada;
+ if (differencingScheme == "Minmod") return DifferencingScheme::minmod;
+ if (differencingScheme == "Superbee") return DifferencingScheme::superbee;
+ if (differencingScheme == "Umist") return DifferencingScheme::umist;
+ if (differencingScheme == "Mclimiter") return DifferencingScheme::mclimiter;
+ if (differencingScheme == "Wahyd") return DifferencingScheme::wahyd;
DUNE_THROW(ParameterException, "\nThis differencing scheme: \"" << differencingScheme << "\" is not implemented.\n"
<< "The available differencing schemes are as follows: \n"
<< differencingSchemeToString(DifferencingScheme::vanleer) << ": The Vanleer flux limiter\n"
@@ -179,18 +181,17 @@ public:
{
switch (differencingScheme)
{
- case DifferencingScheme::vanleer: return "vanleer";
- case DifferencingScheme::vanalbada: return "vanalbada";
- case DifferencingScheme::minmod: return "minmod";
- case DifferencingScheme::superbee: return "superbee";
- case DifferencingScheme::umist: return "umist";
- case DifferencingScheme::mclimiter: return "mclimiter";
- case DifferencingScheme::wahyd: return "wahyd";
+ case DifferencingScheme::vanleer: return "Vanleer";
+ case DifferencingScheme::vanalbada: return "Vanalbada";
+ case DifferencingScheme::minmod: return "Minmod";
+ case DifferencingScheme::superbee: return "Superbee";
+ case DifferencingScheme::umist: return "Umist";
+ case DifferencingScheme::mclimiter: return "Mclimiter";
+ case DifferencingScheme::wahyd: return "Wahyd";
default: return "Invalid"; // should never be reached
}
}
-
/**
* \brief Upwind Method
*/
@@ -200,60 +201,13 @@ public:
return (upwindWeight_ * upstreamMomentum + (1.0 - upwindWeight_) * downstreamMomentum);
}
- /**
- * \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 Diminishing
*
*/
- Scalar tvd(const Scalar downstreamVelocity,
- const Scalar upstreamVelocity,
- const Scalar upUpstreamVelocity,
- const Scalar upstreamToDownstreamDistance,
- const Scalar upUpstreamToUpstreamDistance,
- const Scalar downstreamStaggeredCellSize,
+ Scalar tvd(const std::array<Scalar,3>& momenta,
+ const std::array<Scalar,3>& distances,
const bool selfIsUpstream,
- const Scalar density,
const TvdApproach tvdApproach) const
{
Scalar momentum = 0.0;
@@ -261,17 +215,17 @@ public:
{
case TvdApproach::uniform :
{
- momentum += tvdUniform(downstreamVelocity, upstreamVelocity, upUpstreamVelocity, density);
+ momentum += tvdUniform(momenta, distances, selfIsUpstream);
break;
}
case TvdApproach::li :
{
- momentum += tvdLi(downstreamVelocity, upstreamVelocity, upUpstreamVelocity, upstreamToDownstreamDistance, upUpstreamToUpstreamDistance, selfIsUpstream, density);
+ momentum += tvdLi(momenta, distances, selfIsUpstream);
break;
}
case TvdApproach::hou :
{
- momentum += tvdHou(downstreamVelocity, upstreamVelocity, upUpstreamVelocity, upstreamToDownstreamDistance, upUpstreamToUpstreamDistance, downstreamStaggeredCellSize, density);
+ momentum += tvdHou(momenta, distances, selfIsUpstream);
break;
}
default:
@@ -288,22 +242,24 @@ public:
*
* This function assumes the cell size distribution to be uniform.
*/
- Scalar tvdUniform(const Scalar downstreamVelocity,
- const Scalar upstreamVelocity,
- const Scalar upUpstreamVelocity,
- const Scalar density) const
+ Scalar tvdUniform(const std::array<Scalar,3>& momenta,
+ const std::array<Scalar,3>& distances,
+ const bool selfIsUpstream) const
{
using std::isfinite;
- const Scalar ratio = (upstreamVelocity - upUpstreamVelocity) / (downstreamVelocity - upstreamVelocity);
+ const Scalar downstreamMomentum = momenta[0];
+ const Scalar upstreamMomentum = momenta[1];
+ const Scalar upUpstreamMomentum = momenta[2];
+ const Scalar ratio = (upstreamMomentum - upUpstreamMomentum) / (downstreamMomentum - upstreamMomentum);
// If the velocity field is uniform (like at the first newton step) we get a NaN
if(ratio > 0.0 && isfinite(ratio))
{
- const Scalar secondOrderTerm = 0.5 * limiter_(ratio, 2.0) * (downstreamVelocity - upstreamVelocity);
- return density * (upstreamVelocity + secondOrderTerm);
+ const Scalar secondOrderTerm = 0.5 * limiter_(ratio, 2.0) * (downstreamMomentum - upstreamMomentum);
+ return (upstreamMomentum + secondOrderTerm);
}
else
- return density * upstreamVelocity;
+ return upstreamMomentum;
}
/**
@@ -313,31 +269,33 @@ public:
* It tries to reconstruct the value for the velocity at the upstream-upstream point
* if the grid was uniform.
*/
- Scalar tvdLi(const Scalar downstreamVelocity,
- const Scalar upstreamVelocity,
- const Scalar upUpstreamVelocity,
- const Scalar upstreamToDownstreamDistance,
- const Scalar upUpstreamToUpstreamDistance,
- const bool selfIsUpstream,
- const Scalar density) const
+ Scalar tvdLi(const std::array<Scalar,3>& momenta,
+ const std::array<Scalar,3>& distances,
+ const bool selfIsUpstream) const
{
using std::isfinite;
- // I need the information of selfIsUpstream to get the correct sign because upUpstreamToUpstreamDistance is always positive
- const Scalar upUpstreamGradient = (upstreamVelocity - upUpstreamVelocity) / upUpstreamToUpstreamDistance * selfIsUpstream;
+ const Scalar downstreamMomentum = momenta[0];
+ const Scalar upstreamMomentum = momenta[1];
+ const Scalar upUpstreamMomentum = momenta[2];
+ const Scalar upstreamToDownstreamDistance = distances[0];
+ const Scalar upUpstreamToUpstreamDistance = distances[1];
+
+ // selfIsUpstream is required to get the correct sign because upUpstreamToUpstreamDistance is always positive
+ const Scalar upUpstreamGradient = (upstreamMomentum - upUpstreamMomentum) / 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);
+ const Scalar reconstrutedUpUpstreamVelocity = upUpstreamMomentum + upUpstreamGradient * correctionDistance;
+ const Scalar ratio = (upstreamMomentum - reconstrutedUpUpstreamVelocity) / (downstreamMomentum - upstreamMomentum);
// If the velocity field is uniform (like at the first newton step) we get a NaN
if(ratio > 0.0 && isfinite(ratio))
{
- const Scalar secondOrderTerm = 0.5 * limiter_(ratio, 2.0) * (downstreamVelocity - upstreamVelocity);
- return density * (upstreamVelocity + secondOrderTerm);
+ const Scalar secondOrderTerm = 0.5 * limiter_(ratio, 2.0) * (downstreamMomentum - upstreamMomentum);
+ return (upstreamMomentum + secondOrderTerm);
}
else
- return density * upstreamVelocity;
+ return upstreamMomentum;
}
/**
@@ -346,16 +304,18 @@ public:
* This function 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 tvdHou(const Scalar downstreamVelocity,
- const Scalar upstreamVelocity,
- const Scalar upUpstreamVelocity,
- const Scalar upstreamToDownstreamDistance,
- const Scalar upUpstreamToUpstreamDistance,
- const Scalar downstreamStaggeredCellSize,
- const Scalar density) const
+ Scalar tvdHou(const std::array<Scalar,3>& momenta,
+ const std::array<Scalar,3>& distances,
+ const bool selfIsUpstream) const
{
using std::isfinite;
- const Scalar ratio = (upstreamVelocity - upUpstreamVelocity) / (downstreamVelocity - upstreamVelocity)
+ const Scalar downstreamMomentum = momenta[0];
+ const Scalar upstreamMomentum = momenta[1];
+ const Scalar upUpstreamMomentum = momenta[2];
+ const Scalar upstreamToDownstreamDistance = distances[0];
+ const Scalar upUpstreamToUpstreamDistance = distances[1];
+ const Scalar downstreamStaggeredCellSize = distances[2];
+ const Scalar ratio = (upstreamMomentum - upUpstreamMomentum) / (downstreamMomentum - upstreamMomentum)
* upstreamToDownstreamDistance / upUpstreamToUpstreamDistance;
// If the velocity field is uniform (like at the first newton step) we get a NaN
@@ -363,11 +323,11 @@ public:
{
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);
+ const Scalar secondOrderTerm = limiter_(ratio, R) / R * (downstreamMomentum - upstreamMomentum);
+ return (upstreamMomentum + secondOrderTerm);
}
else
- return density * upstreamVelocity;
+ return upstreamMomentum;
}
/**
@@ -459,4 +419,4 @@ private:
} // end namespace Dumux
-#endif // DUMUX_HIGHER_ORDER_VELOCITY_APPROXIMATION_HH
+#endif // DUMUX_UPWINDING_METHODS_HH
diff --git a/test/freeflow/navierstokes/kovasznay/CMakeLists.txt b/test/freeflow/navierstokes/kovasznay/CMakeLists.txt
index f549843b312c1e46c5a53c91d6af8746889a7d06..3b2f353413d0744c1d8f05116eb4fbadeee4c1e0 100644
--- a/test/freeflow/navierstokes/kovasznay/CMakeLists.txt
+++ b/test/freeflow/navierstokes/kovasznay/CMakeLists.txt
@@ -2,6 +2,7 @@ dune_add_test(NAME test_ff_navierstokes_kovasznay
LABELS freeflow navierstokes
SOURCES main.cc
CMAKE_GUARD HAVE_UMFPACK
+ COMPILE_DEFINITIONS UPWINDSCHEMEORDER=1
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_ff_navierstokes_kovasznay-reference.vtu
@@ -10,14 +11,17 @@ dune_add_test(NAME test_ff_navierstokes_kovasznay
-Problem.Name test_ff_navierstokes_kovasznay")
dune_add_test(NAME test_ff_navierstokes_kovasznay_higherorder
- TARGET test_ff_navierstokes_kovasznay
- LABELS freeflow navierstokes
+ SOURCES main.cc
+ LABELS freeflow
+ CMAKE_GUARD HAVE_UMFPACK
+ COMPILE_DEFINITIONS UPWINDSCHEMEORDER=2
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_ff_navierstokes_kovasznay_higherorder-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_kovasznay_higherorder-00001.vtu
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_kovasznay params.input
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_kovasznay_higherorder params.input
-Problem.Name test_ff_navierstokes_kovasznay_higherorder
- -Discretization.TvdApproach Hou
- -Discretization.DifferencingScheme vanleer")
+ -Flux.TvdApproach Hou
+ -Flux.DifferencingScheme Vanleer")
+
dune_symlink_to_source_files(FILES "params.input")
diff --git a/test/freeflow/navierstokes/kovasznay/params.input b/test/freeflow/navierstokes/kovasznay/params.input
index 0f47acc145b5673a089a530996221e5acf63b58a..e3dc1f6b2da8eab3b09ec58d6342113be4c52043 100644
--- a/test/freeflow/navierstokes/kovasznay/params.input
+++ b/test/freeflow/navierstokes/kovasznay/params.input
@@ -19,3 +19,7 @@ MaxRelativeShift = 1e-5
[Vtk]
WriteFaceData = false
+
+[Flux]
+TvdApproach = Hou
+DifferencingScheme = Vanleer
diff --git a/test/freeflow/navierstokes/kovasznay/problem.hh b/test/freeflow/navierstokes/kovasznay/problem.hh
index bbbb7ab1ec1dab8c282c93977486df2f1cbda0d9..8a228736b3221b7e6c934b2e9ea984593a73e970 100644
--- a/test/freeflow/navierstokes/kovasznay/problem.hh
+++ b/test/freeflow/navierstokes/kovasznay/problem.hh
@@ -25,6 +25,10 @@
#ifndef DUMUX_KOVASZNAY_TEST_PROBLEM_HH
#define DUMUX_KOVASZNAY_TEST_PROBLEM_HH
+#ifndef UPWINDSCHEMEORDER
+#define UPWINDSCHEMEORDER 0
+#endif
+
#include <dune/grid/yaspgrid.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
@@ -68,6 +72,9 @@ template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::KovasznayTest> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::KovasznayTest> { static constexpr bool value = true; };
+
+template<class TypeTag>
+struct UpwindSchemeOrder<TypeTag, TTag::KovasznayTest> { static constexpr int value = UPWINDSCHEMEORDER; };
} // end namespace Properties
/*!
@@ -99,12 +106,14 @@ class KovasznayTestProblem : public NavierStokesProblem<TypeTag>
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
+ static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
+
public:
KovasznayTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry), eps_(1e-6)
{
printL2Error_ = getParam<bool>("Problem.PrintL2Error");
-
+ std::cout<< "upwindSchemeOrder is: " << FVGridGeometry::order() << "\n";
kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
Scalar reynoldsNumber = 1.0 / kinematicViscosity_;
lambda_ = 0.5 * reynoldsNumber