diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index d3ce6cfa36857a00499795f0246b9e25c741615c..c32c8e1a4dcc42da83d8451e4948a02b00e6c6ab 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -233,7 +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 UpwindSchemeOrder { using type = UndefinedProperty; }; //!< Specifies the order of the upwinding scheme (1 == first order, 2 == second order(tvd methods))
/////////////////////////////////////////////////////////////
// Properties used by the mpnc model
diff --git a/dumux/discretization/staggered.hh b/dumux/discretization/staggered.hh
index 84040d5ec2507a9d497656ebb71a326d30e2f5b9..92696168fc0a94fb3c0da39ce5e6d5f276adf261 100644
--- a/dumux/discretization/staggered.hh
+++ b/dumux/discretization/staggered.hh
@@ -81,12 +81,13 @@ template<class TypeTag>
struct GridFluxVariablesCache<TypeTag, TTag::StaggeredModel>
{
private:
- static constexpr auto enableCache = getPropValue<TypeTag, Properties::EnableGridFluxVariablesCache>();
using Problem = GetPropType<TypeTag, Properties::Problem>;
using FluxVariablesCache = GetPropType<TypeTag, Properties::FluxVariablesCache>;
using FluxVariablesCacheFiller = FluxVariablesCaching::EmptyCacheFiller;
+ static constexpr auto enableCache = getPropValue<TypeTag, Properties::EnableGridFluxVariablesCache>();
+ static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
public:
- using type = StaggeredGridFluxVariablesCache<Problem, FluxVariablesCache, FluxVariablesCacheFiller, enableCache>;
+ using type = StaggeredGridFluxVariablesCache<Problem, FluxVariablesCache, FluxVariablesCacheFiller, enableCache, upwindSchemeOrder>;
};
//! Set the face solution type
diff --git a/dumux/discretization/staggered/freeflow/connectivitymap.hh b/dumux/discretization/staggered/freeflow/connectivitymap.hh
index e43321f1c0c68767972c7212d5b5002febda09c7..1aa82d24fe1e852e2225648e240e1eb4ac79be2b 100644
--- a/dumux/discretization/staggered/freeflow/connectivitymap.hh
+++ b/dumux/discretization/staggered/freeflow/connectivitymap.hh
@@ -52,8 +52,13 @@ class StaggeredFreeFlowConnectivityMap
using FaceToCellCenterMap = std::vector<std::vector<GridIndexType>>;
using FaceToFaceMap = std::vector<std::vector<GridIndexType>>;
+ using SmallLocalIndex = typename IndexTraits<GridView>::SmallLocalIndex;
+
using Stencil = std::vector<GridIndexType>;
+ static constexpr SmallLocalIndex upwindSchemeOrder = FVGridGeometry::upwindSchemeOrder;
+ static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
+
public:
//! Update the map and prepare the stencils
@@ -92,13 +97,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
{
@@ -186,24 +184,9 @@ private:
{
if(stencil.empty())
{
- for(int i = 0; i < stencilOrder_ - 1; i++)
- {
- if(scvf.hasBackwardNeighbor(i))
- {
- stencil.push_back(scvf.axisData().inAxisBackwardDofs[i]);
- }
- }
-
stencil.push_back(scvf.axisData().selfDof);
stencil.push_back(scvf.axisData().oppositeDof);
-
- for(int i = 0; i < stencilOrder_ - 1; i++)
- {
- if(scvf.hasForwardNeighbor(i))
- {
- stencil.push_back(scvf.axisData().inAxisForwardDofs[i]);
- }
- }
+ addHigherOrderInAxisDofs_(scvf, stencil, std::integral_constant<bool, useHigherOrder>{});
}
for(const auto& data : scvf.pairData())
@@ -214,21 +197,33 @@ private:
stencil.push_back(data.normalPair.second);
// add parallel dofs
- for (int i = 0; i < stencilOrder_ /*data.parallelDofs.size()*/; i++)
+ for (SmallLocalIndex i = 0; i < upwindSchemeOrder; i++)
{
- if(!(data.parallelDofs[i] < 0))
- {
+ if (data.hasParallelNeighbor[i])
stencil.push_back(data.parallelDofs[i]);
- }
}
}
}
+ void addHigherOrderInAxisDofs_(const SubControlVolumeFace& scvf, Stencil& stencil, std::false_type) {}
+
+ void addHigherOrderInAxisDofs_(const SubControlVolumeFace& scvf, Stencil& stencil, std::true_type)
+ {
+ for (SmallLocalIndex 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]);
+ }
+ }
+
+
CellCenterToCellCenterMap cellCenterToCellCenterMap_;
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 e227cfdc046c20788df80bdd921032862069b858..4645f66c52013ef06d216e19ec2e33b9b1fdd660 100644
--- a/dumux/discretization/staggered/freeflow/facevariables.hh
+++ b/dumux/discretization/staggered/freeflow/facevariables.hh
@@ -29,17 +29,42 @@
namespace Dumux {
+namespace Detail {
+
+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>
+struct InAxisVelocities<Scalar, 1>
+{
+ Scalar self = 0.0;
+ Scalar opposite = 0.0;
+};
+
+} // 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.
+ * When the upwindSchemeOrder is set to 2, additional velocities located at Dofs
+ * further from the central stencil will be added and used when calculating the
+ * advective term. When the order remains at 1, these velocities will not be provided.
*/
-template<class FacePrimaryVariables, int dim>
+template<class FacePrimaryVariables, int dim, int upwindSchemeOrder>
class StaggeredFaceVariables
{
static constexpr int numPairs = (dim == 2) ? 2 : 4;
- static constexpr int order = 2;
+ static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
+
using Scalar = typename FacePrimaryVariables::block_type;
+ using InAxisVelocities = Detail::InAxisVelocities<Scalar, upwindSchemeOrder>;
public:
@@ -50,7 +75,7 @@ public:
*/
void updateOwnFaceOnly(const FacePrimaryVariables& priVars)
{
- velocitySelf_ = priVars[0];
+ inAxisVelocities_.self = priVars[0];
}
/*!
@@ -71,52 +96,30 @@ public:
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
- velocitySelf_ = faceSol[scvf.dofIndex()];
- velocityOpposite_ = faceSol[scvf.dofIndexOpposingFace()];
+ inAxisVelocities_.self = faceSol[scvf.dofIndex()];
+ inAxisVelocities_.opposite = faceSol[scvf.dofIndexOpposingFace()];
- // treat the velocity forward of the self face i.e. the face that is
- // forward wrt the self face by degree i
- velocityForward_.fill(0.0);
- for (int i = 0; i < scvf.axisData().inAxisForwardDofs.size(); i++)
- {
- if(!(scvf.axisData().inAxisForwardDofs[i] < 0))
- {
- velocityForward_[i]= faceSol[scvf.axisData().inAxisForwardDofs[i]];
- }
- }
-
- // treat the velocity at the first backward face i.e. the face that is
- // behind the opposite face by degree i
- velocityBackward_.fill(0.0);
- for (int i = 0; i < scvf.axisData().inAxisBackwardDofs.size(); i++)
- {
- if(!(scvf.axisData().inAxisBackwardDofs[i] < 0))
- {
- velocityBackward_[i] = faceSol[scvf.axisData().inAxisBackwardDofs[i]];
- }
- }
+ addHigherOrderInAxisVelocities_(faceSol, scvf, std::integral_constant<bool, useHigherOrder>{});
// handle all sub faces
- for(int i = 0; i < velocityParallel_.size(); ++i)
+ for (int i = 0; i < velocityParallel_.size(); ++i)
velocityParallel_[i].fill(0.0);
- for(int i = 0; i < scvf.pairData().size(); ++i)
+ for (int i = 0; i < scvf.pairData().size(); ++i)
{
const auto& subFaceData = scvf.pairData(i);
// treat the velocities normal to the face
velocityNormalInside_[i] = faceSol[subFaceData.normalPair.first];
- if(scvf.hasOuterNormal(i))
+ if (scvf.hasOuterNormal(i))
velocityNormalOutside_[i] = faceSol[subFaceData.normalPair.second];
// treat the velocities parallel to the self face
- for(int j = 0; j < subFaceData.parallelDofs.size(); j++)
+ for (int j = 0; j < upwindSchemeOrder; j++)
{
- if(scvf.hasParallelNeighbor(i,j))
- {
+ if (scvf.hasParallelNeighbor(i,j))
velocityParallel_[i][j] = faceSol[subFaceData.parallelDofs[j]];
- }
}
}
}
@@ -126,7 +129,7 @@ public:
*/
Scalar velocitySelf() const
{
- return velocitySelf_;
+ return inAxisVelocities_.self;
}
/*!
@@ -134,7 +137,7 @@ public:
*/
Scalar velocityOpposite() const
{
- return velocityOpposite_;
+ return inAxisVelocities_.opposite;
}
/*!
@@ -142,9 +145,10 @@ public:
*
* \param backwardIdx The index describing how many faces backward this dof is from the opposite face
*/
+ template<bool enable = useHigherOrder, std::enable_if_t<enable, int> = 0>
Scalar velocityBackward(const int backwardIdx) const
{
- return velocityBackward_[backwardIdx];
+ return inAxisVelocities_.backward[backwardIdx];
}
/*!
@@ -152,9 +156,10 @@ public:
*
* \param forwardIdx The index describing how many faces forward this dof is of the self face
*/
+ template<bool enable = useHigherOrder, std::enable_if_t<enable, int> = 0>
Scalar velocityForward(const int forwardIdx) const
{
- return velocityForward_[forwardIdx];
+ return inAxisVelocities_.forward[forwardIdx];
}
/*!
@@ -163,7 +168,7 @@ public:
* \param localSubFaceIdx The local index of the subface
* \param parallelDegreeIdx The index describing how many faces parallel this dof is of the parallel face
*/
- Scalar velocityParallel(const int localSubFaceIdx, const int parallelDegreeIdx) const
+ Scalar velocityParallel(const int localSubFaceIdx, const int parallelDegreeIdx = 0) const
{
return velocityParallel_[localSubFaceIdx][parallelDegreeIdx];
}
@@ -189,11 +194,33 @@ public:
}
private:
- Scalar velocitySelf_;
- Scalar velocityOpposite_;
- std::array<Scalar, order-1> velocityForward_;
- std::array<Scalar, order-1> velocityBackward_;
- std::array<std::array<Scalar, order>, numPairs> velocityParallel_;
+
+ template<class SolVector, class SubControlVolumeFace>
+ void addHigherOrderInAxisVelocities_(const SolVector& faceSol, const SubControlVolumeFace& scvf, std::false_type) {}
+
+ template<class SolVector, class SubControlVolumeFace>
+ 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
+ // forward wrt the self face by degree i
+ for (int i = 0; i < scvf.axisData().inAxisForwardDofs.size(); i++)
+ {
+ if (scvf.hasForwardNeighbor(i))
+ inAxisVelocities_.forward[i]= faceSol[scvf.axisData().inAxisForwardDofs[i]];
+ }
+
+ // treat the velocity at the first backward face i.e. the face that is
+ // behind the opposite face by degree i
+ for (int i = 0; i < scvf.axisData().inAxisBackwardDofs.size(); i++)
+ {
+ if (scvf.hasBackwardNeighbor(i))
+ inAxisVelocities_.backward[i] = faceSol[scvf.axisData().inAxisBackwardDofs[i]];
+ }
+ }
+
+ InAxisVelocities inAxisVelocities_;
+ 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 4b88982ba9b18d9fd11d11b7e7b0bf032f1c9877..c9e88535066c1de09024b420d6427001c08d74ea 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, class MapperTraits = DefaultMapperTraits<GridView>>
+template<class GridView, int upwOrder, class MapperTraits = DefaultMapperTraits<GridView>>
struct StaggeredFreeFlowDefaultFVGridGeometryTraits
: public MapperTraits
{
using SubControlVolume = CCSubControlVolume<GridView>;
- using SubControlVolumeFace = FreeFlowStaggeredSubControlVolumeFace<GridView>;
+ using SubControlVolumeFace = FreeFlowStaggeredSubControlVolumeFace<GridView, upwOrder>;
using IntersectionMapper = ConformingGridIntersectionMapper<GridView>;
- using GeometryHelper = FreeFlowStaggeredGeometryHelper<GridView>;
+ using GeometryHelper = FreeFlowStaggeredGeometryHelper<GridView, upwOrder>;
+ static constexpr int upwindSchemeOrder = upwOrder;
struct DofTypeIndices
{
diff --git a/dumux/discretization/staggered/freeflow/properties.hh b/dumux/discretization/staggered/freeflow/properties.hh
index 3075fc93c6690e07c2fe0bd4bcbe31c8e367ba74..1588ca21e4783c0eeeee8c980c99fe339ccce426 100644
--- a/dumux/discretization/staggered/freeflow/properties.hh
+++ b/dumux/discretization/staggered/freeflow/properties.hh
@@ -81,9 +81,10 @@ template<class TypeTag>
struct FVGridGeometry<TypeTag, TTag::StaggeredFreeFlowModel>
{
private:
- using GridView = GetPropType<TypeTag, Properties::GridView>;
- using Traits = StaggeredFreeFlowDefaultFVGridGeometryTraits<GridView>;
+ 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, upwindSchemeOrder>;
public:
using type = StaggeredFVGridGeometry<GridView, enableCache, Traits>;
};
@@ -95,8 +96,9 @@ struct FaceVariables<TypeTag, TTag::StaggeredFreeFlowModel>
private:
using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
+ static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
public:
- using type = StaggeredFaceVariables<FacePrimaryVariables, GridView::dimension>;
+ using type = StaggeredFaceVariables<FacePrimaryVariables, GridView::dimension, upwindSchemeOrder>;
};
//! Set the default global volume variables cache vector class
@@ -128,6 +130,12 @@ struct VelocityOutput<TypeTag, TTag::StaggeredFreeFlowModel>
GetPropType<TypeTag, Properties::SolutionVector>>;
};
+/*!
+ * \brief Set the order of the upwinding scheme to 1 by default.
+ */
+template<class TypeTag>
+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 6990c0daed38f81762929ddeedec7f1f85febbbd..071888dfddb0381fe66e3ae4b540fe40dadb856d 100644
--- a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
+++ b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
@@ -25,47 +25,86 @@
#define DUMUX_DISCRETIZATION_STAGGERED_GEOMETRY_HELPER_HH
#include <dune/geometry/multilineargeometry.hh>
-#include <dune/geometry/referenceelements.hh>
#include <dumux/common/math.hh>
+#include <dumux/common/indextraits.hh>
#include <type_traits>
#include <algorithm>
#include <array>
+#include <bitset>
namespace Dumux
{
+namespace Detail {
+
/*!
* \ingroup StaggeredDiscretization
* \brief Parallel Data stored per sub face
*/
-template<class Scalar, class GlobalPosition, int geometryOrder>
+template<class GridView, int upwindSchemeOrder>
struct PairData
{
- std::array<int, geometryOrder> parallelDofs;
- std::array<Scalar, geometryOrder+1> parallelDistances;
- std::pair<signed int,signed int> normalPair;
- int localNormalFaceIdx;
+ using Scalar = typename GridView::ctype;
+ using GlobalPosition = typename GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate;
+ using GridIndexType = typename IndexTraits<GridView>::GridIndex;
+ using SmallLocalIndexType = typename IndexTraits<GridView>::SmallLocalIndex;
+
+ std::bitset<upwindSchemeOrder> hasParallelNeighbor;
+ std::array<GridIndexType, upwindSchemeOrder> parallelDofs;
+ std::array<Scalar, upwindSchemeOrder+1> parallelDistances; // TODO store only two distances, not three.
+ bool hasNormalNeighbor = false;
+ std::pair<GridIndexType, GridIndexType> normalPair;
+ SmallLocalIndexType localNormalFaceIdx;
Scalar normalDistance;
GlobalPosition virtualFirstParallelFaceDofPos;
};
+
+/*!
+ * \ingroup StaggeredDiscretization
+ * \brief In Axis Data stored per sub face
+ */
+template<class GridView, int upwindSchemeOrder>
+struct AxisData;
/*!
* \ingroup StaggeredDiscretization
* \brief In Axis Data stored per sub face
*/
-template<class Scalar, int geometryOrder>
+template<class GridView, int upwindSchemeOrder>
struct AxisData
{
- int selfDof;
- int oppositeDof;
- std::array<int, geometryOrder-1> inAxisForwardDofs;
- std::array<int, geometryOrder-1> inAxisBackwardDofs;
+ using Scalar = typename GridView::ctype;
+ using GridIndexType = typename IndexTraits<GridView>::GridIndex;
+
+ GridIndexType selfDof;
+ GridIndexType oppositeDof;
+ std::bitset<upwindSchemeOrder-1> hasForwardNeighbor;
+ std::bitset<upwindSchemeOrder-1> hasBackwardNeighbor;
+ std::array<GridIndexType, upwindSchemeOrder-1> inAxisForwardDofs;
+ std::array<GridIndexType, 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;
};
+/*!
+ * \ingroup StaggeredDiscretization
+ * \brief In Axis Data stored per sub face
+ */
+template<class GridView>
+struct AxisData<GridView, 1>
+{
+ using Scalar = typename GridView::ctype;
+ using GridIndexType = typename IndexTraits<GridView>::GridIndex;
+
+ GridIndexType selfDof;
+ GridIndexType oppositeDof;
+ Scalar selfToOppositeDistance;
+};
+
+} // namespace Detail
+
/*!
* \ingroup StaggeredDiscretization
* \brief Returns the dirction index of the facet (0 = x, 1 = y, 2 = z)
@@ -74,8 +113,7 @@ template<class Vector>
inline static unsigned int directionIndex(Vector&& vector)
{
const auto eps = 1e-8;
- const int idx = std::find_if(vector.begin(), vector.end(), [eps](const auto& x) { return std::abs(x) > eps; } ) - vector.begin();
- return idx;
+ return std::find_if(vector.begin(), vector.end(), [eps](const auto& x) { return std::abs(x) > eps; } ) - vector.begin();
}
/*!
@@ -83,31 +121,32 @@ 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>
+template<class GridView, int upwindSchemeOrder>
class FreeFlowStaggeredGeometryHelper
{
using Scalar = typename GridView::ctype;
- static constexpr int dim = GridView::dimension;
- static constexpr int dimWorld = GridView::dimensionworld;
-
- using ScvGeometry = Dune::CachedMultiLinearGeometry<Scalar, dim, dimWorld>;
- using ScvfGeometry = Dune::CachedMultiLinearGeometry<Scalar, dim-1, dimWorld>;
-
- using GlobalPosition = typename ScvGeometry::GlobalCoordinate;
+ static constexpr auto dim = GridView::dimension;
+ static constexpr auto dimWorld = GridView::dimensionworld;
using Element = typename GridView::template Codim<0>::Entity;
using Intersection = typename GridView::Intersection;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
- using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
+ using GridIndexType = typename IndexTraits<GridView>::GridIndex;
+ using SmallLocalIndexType = typename IndexTraits<GridView>::SmallLocalIndex;
//TODO include assert that checks for quad geometry
- static constexpr int codimIntersection = 1;
- static constexpr int codimCommonEntity = 2;
- static constexpr int numFacetSubEntities = (dim == 2) ? 2 : 4;
- static constexpr int numfacets = dimWorld * 2;
- static constexpr int numPairs = 2 * (dimWorld - 1);
+ static constexpr auto codimIntersection = 1;
+ static constexpr auto codimCommonEntity = 2;
+ static constexpr auto numFacetSubEntities = (dim == 2) ? 2 : 4;
+ static constexpr auto numfacets = dimWorld * 2;
+ static constexpr auto numPairs = 2 * (dimWorld - 1);
+
+ static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
public:
+ using PairData = Detail::PairData<GridView, upwindSchemeOrder>;
+ using AxisData = Detail::AxisData<GridView, upwindSchemeOrder>;
FreeFlowStaggeredGeometryHelper(const Element& element, const GridView& gridView) : element_(element), elementGeometry_(element.geometry()), gridView_(gridView)
{ }
@@ -117,15 +156,14 @@ public:
void updateLocalFace(const IntersectionMapper& intersectionMapper_, const Intersection& intersection)
{
intersection_ = intersection;
- innerNormalFacePos_.clear();
- fillPairData_();
fillAxisData_();
+ fillPairData_();
}
/*!
* \brief Returns the local index of the face (i.e. the intersection)
*/
- int localFaceIndex() const
+ SmallLocalIndexType localFaceIndex() const
{
return intersection_.indexInInside();
}
@@ -133,7 +171,7 @@ public:
/*!
* \brief Returns a copy of the axis data
*/
- auto axisData() const
+ AxisData axisData() const
{
return axisData_;
}
@@ -141,7 +179,7 @@ public:
/*!
* \brief Returns a copy of the pair data
*/
- auto pairData() const
+ std::array<PairData, numPairs> pairData() const
{
return pairData_;
}
@@ -149,49 +187,65 @@ public:
/*!
* \brief Returns the dirction index of the primary facet (0 = x, 1 = y, 2 = z)
*/
- int directionIndex() const
+ unsigned int directionIndex() const
{
- return Dumux::directionIndex(std::move(intersection_.centerUnitOuterNormal()));
+ return Dumux::directionIndex(intersection_.centerUnitOuterNormal());
}
/*!
* \brief Returns the dirction index of the facet passed as an argument (0 = x, 1 = y, 2 = z)
*/
- int directionIndex(const Intersection& intersection) const
+ unsigned int directionIndex(const Intersection& intersection) const
{
return Dumux::directionIndex(std::move(intersection.centerUnitOuterNormal()));
}
- static constexpr int geometryOrder = 2;
-
private:
+
/*!
* \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_()
{
- unsigned int numForwardBackwardAxisDofs = axisData_.inAxisForwardDofs.size();
+ fillOuterAxisData_(std::integral_constant<bool, useHigherOrder>{});
+
const auto inIdx = intersection_.indexInInside();
const auto oppIdx = localOppositeIdx_(inIdx);
- this->axisData_.inAxisForwardDofs.fill(-1);
- this->axisData_.inAxisBackwardDofs.fill(-1);
- this->axisData_.inAxisForwardDistances.fill(0);
- this->axisData_.inAxisBackwardDistances.fill(0);
-
// Set the self Dof
- this->axisData_.selfDof = gridView_.indexSet().subIndex(this->intersection_.inside(), inIdx, codimIntersection);
- const auto self = getFacet_(inIdx, element_);
+ axisData_.selfDof = gridView_.indexSet().subIndex(intersection_.inside(), inIdx, codimIntersection);
// Set the opposite Dof
- this->axisData_.oppositeDof = gridView_.indexSet().subIndex(this->intersection_.inside(), oppIdx, codimIntersection);
- const auto opposite = getFacet_(oppIdx, element_);
+ axisData_.oppositeDof = gridView_.indexSet().subIndex(intersection_.inside(), oppIdx, codimIntersection);
// Set the Self to Opposite Distance
- this->axisData_.selfToOppositeDistance = (self.geometry().center() - opposite.geometry().center()).two_norm();
+ const auto self = getFacet_(inIdx, element_);
+ const auto opposite = getFacet_(oppIdx, element_);
+ axisData_.selfToOppositeDistance = (self.geometry().center() - opposite.geometry().center()).two_norm();
+
+ }
+
+ /*!
+ * \brief Fills the axis data with the outer Dofs and Distances.
+ * For first order upwind methods no outer information is required.
+ */
+ 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)
+ {
+ // reset the axis data struct
+ axisData_ = {};
+
+ const auto numForwardBackwardAxisDofs = axisData_.inAxisForwardDofs.size();
// Set the Forward Dofs
std::stack<Element> inAxisForwardElementStack;
+ const auto inIdx = intersection_.indexInInside();
auto selfFace = getFacet_(inIdx, element_);
if(intersection_.neighbor())
{
@@ -222,26 +276,25 @@ private:
while(!inAxisForwardElementStack.empty())
{
int forwardIdx = inAxisForwardElementStack.size()-1;
- this->axisData_.inAxisForwardDofs[forwardIdx] = gridView_.indexSet().subIndex(inAxisForwardElementStack.top(), inIdx, codimIntersection);
+ axisData_.hasForwardNeighbor.set(forwardIdx, true);
+ axisData_.inAxisForwardDofs[forwardIdx] = gridView_.indexSet().subIndex(inAxisForwardElementStack.top(), inIdx, codimIntersection);
selfFace = getFacet_(inIdx, inAxisForwardElementStack.top());
forwardFaceCoordinates[forwardIdx] = selfFace.geometry().center();
inAxisForwardElementStack.pop();
}
- for(int i = 0; i< forwardFaceCoordinates.size(); i++)
+ 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();
- }
+ if (i == 0)
+ 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_))
{
@@ -279,78 +332,61 @@ private:
while(!inAxisBackwardElementStack.empty())
{
int backwardIdx = inAxisBackwardElementStack.size()-1;
- this->axisData_.inAxisBackwardDofs[backwardIdx] = gridView_.indexSet().subIndex(inAxisBackwardElementStack.top(), oppIdx, codimIntersection);
+ axisData_.hasBackwardNeighbor.set(backwardIdx, true);
+ axisData_.inAxisBackwardDofs[backwardIdx] = gridView_.indexSet().subIndex(inAxisBackwardElementStack.top(), oppIdx, codimIntersection);
oppFace = getFacet_(oppIdx, inAxisBackwardElementStack.top());
backwardFaceCoordinates[backwardIdx] = oppFace.geometry().center();
inAxisBackwardElementStack.pop();
}
- for(int i = 0; i< backwardFaceCoordinates.size(); i++)
+ const auto opposite = getFacet_(oppIdx, element_);
+ for (int i = 0; i< backwardFaceCoordinates.size(); i++)
{
- if(i == 0)
- {
- this->axisData_.inAxisBackwardDistances[i] = (backwardFaceCoordinates[i] - opposite.geometry().center()).two_norm();
- }
+ if (i == 0)
+ axisData_.inAxisBackwardDistances[i] = (backwardFaceCoordinates[i] - opposite.geometry().center()).two_norm();
else
- {
- this->axisData_.inAxisBackwardDistances[i] = (backwardFaceCoordinates[i] - backwardFaceCoordinates[i-1]).two_norm();
- }
+ axisData_.inAxisBackwardDistances[i] = (backwardFaceCoordinates[i] - backwardFaceCoordinates[i-1]).two_norm();
}
}
/*!
- * \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_()
{
- // 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();
+ // reset the pair data structs
+ pairData_ = {};
// set basic global positions
- const auto& elementCenter = this->element_.geometry().center();
+ const auto& elementCenter = element_.geometry().center();
const auto& FacetCenter = intersection_.geometry().center();
// get the inner normal Dof Index
- int numPairInnerNormalIdx = 0;
- for(const auto& innerElementIntersection : intersections(gridView_, element_))
+ SmallLocalIndexType numPairInnerNormalIdx = 0;
+ for (const auto& innerElementIntersection : intersections(gridView_, element_))
{
- if(facetIsNormal_(innerElementIntersection.indexInInside(), intersection_.indexInInside()))
+ if (facetIsNormal_(innerElementIntersection.indexInInside(), intersection_.indexInInside()))
{
- const int innerElementIntersectionIdx = innerElementIntersection.indexInInside();
+ const auto innerElementIntersectionIdx = innerElementIntersection.indexInInside();
setNormalPairFirstInfo_(innerElementIntersectionIdx, element_, numPairInnerNormalIdx);
numPairInnerNormalIdx++;
-
- innerNormalFacePos_.reserve(numPairs);
- const auto& innerNormalFacet = getFacet_(innerElementIntersection.indexInInside(), element_);
- innerNormalFacePos_.emplace_back(innerNormalFacet.geometry().center());
}
}
// get the outer normal Dof Index
- int numPairOuterNormalIdx = 0;
- if(intersection_.neighbor())
+ SmallLocalIndexType numPairOuterNormalIdx = 0;
+ if (intersection_.neighbor())
{
// the direct neighbor element and the respective intersection index
const auto& directNeighborElement = intersection_.outside();
- for(const auto& directNeighborElementIntersection : intersections(gridView_, directNeighborElement))
+ for (const auto& directNeighborElementIntersection : intersections(gridView_, directNeighborElement))
{
// skip the directly neighboring face itself and its opposing one
- if(facetIsNormal_(directNeighborElementIntersection.indexInInside(), intersection_.indexInOutside()))
+ if (facetIsNormal_(directNeighborElementIntersection.indexInInside(), intersection_.indexInOutside()))
{
- const int directNeighborElemIsIdx = directNeighborElementIntersection.indexInInside();
+ const auto directNeighborElemIsIdx = directNeighborElementIntersection.indexInInside();
setNormalPairSecondInfo_(directNeighborElemIsIdx, directNeighborElement, numPairOuterNormalIdx);
numPairOuterNormalIdx++;
}
@@ -359,18 +395,74 @@ private:
else // intersection is on boundary
{
// fill the normal pair entries
- for(int pairIdx = 0; pairIdx < numPairs; ++pairIdx)
+ for(SmallLocalIndexType pairIdx = 0; pairIdx < numPairs; ++pairIdx)
{
- assert(pairData_[pairIdx].normalPair.second == -1);
+ assert(!pairData_[pairIdx].hasNormalNeighbor);
const auto distance = FacetCenter - elementCenter;
- this->pairData_[pairIdx].normalDistance = std::move(distance.two_norm());
+ pairData_[pairIdx].normalDistance = std::move(distance.two_norm());
numPairOuterNormalIdx++;
}
}
+ 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)
+ {
+ // set basic global positions and stencil size definitions
+ const auto& elementCenter = element_.geometry().center();
+
+ // get the parallel Dofs
+ const auto parallelLocalIdx = intersection_.indexInInside();
+ SmallLocalIndexType 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.
+ const auto& outerElement = intersection.outside();
+ pairData_[numPairParallelIdx].hasParallelNeighbor.set(0, true);
+ 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 = 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)
+ {
+ // set basic global positions and stencil size definitions
+ const auto numParallelFaces = pairData_[0].parallelDistances.size();
+ const auto& elementCenter = element_.geometry().center();
+
// get the parallel Dofs
- const int parallelLocalIdx = intersection_.indexInInside();
- int numPairParallelIdx = 0;
+ const auto parallelLocalIdx = intersection_.indexInInside();
+ SmallLocalIndexType numPairParallelIdx = 0;
std::stack<Element> parallelElementStack;
for(const auto& intersection : intersections(gridView_, element_))
{
@@ -407,9 +499,10 @@ private:
{
if(parallelElementStack.size() > 1)
{
- this->pairData_[numPairParallelIdx].parallelDofs[parallelElementStack.size()-2] = gridView_.indexSet().subIndex(parallelElementStack.top(), parallelLocalIdx, codimIntersection);
+ pairData_[numPairParallelIdx].hasParallelNeighbor.set(parallelElementStack.size()-2, true);
+ pairData_[numPairParallelIdx].parallelDofs[parallelElementStack.size()-2] = gridView_.indexSet().subIndex(parallelElementStack.top(), parallelLocalIdx, codimIntersection);
}
- this->pairData_[numPairParallelIdx].parallelDistances[parallelElementStack.size()-1] = setParallelPairDistances_(parallelElementStack.top(), parallelAxisIdx);
+ pairData_[numPairParallelIdx].parallelDistances[parallelElementStack.size()-1] = setParallelPairDistances_(parallelElementStack.top(), parallelAxisIdx);
parallelElementStack.pop();
}
@@ -417,17 +510,16 @@ 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;
- const auto virtualFirstParallelFaceDofPos = this->intersection_.geometry().center() + distance;
+ const auto virtualFirstParallelFaceDofPos = intersection_.geometry().center() + distance;
- this->pairData_[numPairParallelIdx].virtualFirstParallelFaceDofPos = std::move(virtualFirstParallelFaceDofPos);
+ pairData_[numPairParallelIdx].virtualFirstParallelFaceDofPos = std::move(virtualFirstParallelFaceDofPos);
// The distance is saved doubled because with scvf.cellCenteredSelfToFirstParallelDistance
// an average between parallelDistances[0] and parallelDistances[1] will be computed
- this->pairData_[numPairParallelIdx].parallelDistances[0] = std::move(distance.two_norm() * 2);
+ pairData_[numPairParallelIdx].parallelDistances[0] = std::move(distance.two_norm() * 2);
}
numPairParallelIdx++;
}
@@ -464,26 +556,25 @@ private:
void setNormalPairFirstInfo_(const int isIdx, const Element& element, const int numPairsIdx)
{
// store the inner normal dofIdx
- auto& dofIdx = pairData_[numPairsIdx].normalPair.first;
- dofIdx = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
+ pairData_[numPairsIdx].normalPair.first = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
// store the local normal facet index
- this->pairData_[numPairsIdx].localNormalFaceIdx = isIdx;
+ pairData_[numPairsIdx].localNormalFaceIdx = isIdx;
}
//! Sets the information about the normal faces (in the neighbor element)
void setNormalPairSecondInfo_(const int isIdx, const Element& element, const int numPairsIdx)
{
// store the dofIdx
- auto& dofIdx = pairData_[numPairsIdx].normalPair.second;
- dofIdx = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
+ pairData_[numPairsIdx].hasNormalNeighbor = true;
+ pairData_[numPairsIdx].normalPair.second = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
// store the element distance
const auto& outerNormalFacet = getFacet_(isIdx, element);
const auto outerNormalFacetPos = outerNormalFacet.geometry().center();
const auto& innerNormalFacet = getFacet_(isIdx, element_);
const auto innerNormalFacetPos = innerNormalFacet.geometry().center();
- this->pairData_[numPairsIdx].normalDistance = (innerNormalFacetPos - outerNormalFacetPos).two_norm();
+ pairData_[numPairsIdx].normalDistance = (innerNormalFacetPos - outerNormalFacetPos).two_norm();
}
//! Sets the information about the parallel distances
@@ -526,9 +617,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 axisData_; //!< Data related to forward and backward faces
+ std::array<PairData, 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 0fa6efc60735a797a487304caea533c7781e9032..9f163725c141730cb9a1906fa099adbddbb668d5 100644
--- a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
@@ -36,29 +36,52 @@
namespace Dumux {
+/*!
+ * \ingroup StaggeredDiscretization
+ * \brief Default traits class to be used for the sub-control volume faces
+ * for the free-flow staggered finite volume scheme
+ * \tparam GridView the type of the grid view
+ * \tparam upwindSchemeOrder the order of the upwind scheme
+ */
+template<class GridView, int upwindSchemeOrder>
+struct FreeFlowStaggeredDefaultScvfGeometryTraits
+{
+ using Geometry = typename GridView::template Codim<1>::Geometry;
+ using GridIndexType = typename IndexTraits<GridView>::GridIndex;
+ using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
+ using Scalar = typename GridView::ctype;
+ using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
+ using PairData = typename FreeFlowStaggeredGeometryHelper<GridView, upwindSchemeOrder>::PairData;
+ using AxisData = typename FreeFlowStaggeredGeometryHelper<GridView, upwindSchemeOrder>::AxisData;
+};
+
/*!
* \ingroup StaggeredDiscretization
* \brief Class for a sub control volume face in the staggered method, i.e a part of the boundary
* of a sub control volume we compute fluxes on. This is a specialization for free flow models.
*/
template<class GV,
- class T = StaggeredDefaultScvfGeometryTraits<GV> >
+ int upwindSchemeOrder,
+ class T = FreeFlowStaggeredDefaultScvfGeometryTraits<GV, upwindSchemeOrder>>
class FreeFlowStaggeredSubControlVolumeFace
-: public SubControlVolumeFaceBase<FreeFlowStaggeredSubControlVolumeFace<GV, T>, T>
+: public SubControlVolumeFaceBase<FreeFlowStaggeredSubControlVolumeFace<GV, upwindSchemeOrder, T>, T>
{
- using ThisType = FreeFlowStaggeredSubControlVolumeFace<GV, T>;
+ using ThisType = FreeFlowStaggeredSubControlVolumeFace<GV, upwindSchemeOrder, T>;
using ParentType = SubControlVolumeFaceBase<ThisType, T>;
using Geometry = typename T::Geometry;
using GridIndexType = typename IndexTraits<GV>::GridIndex;
using LocalIndexType = typename IndexTraits<GV>::LocalIndex;
- using GeometryHelper = FreeFlowStaggeredGeometryHelper<GV>;
+
+ using PairData = typename T::PairData;
+ using AxisData = typename T::AxisData;
using Scalar = typename T::Scalar;
static const int dim = Geometry::mydimension;
static const int dimworld = Geometry::coorddimension;
static constexpr int numPairs = 2 * (dimworld - 1);
- static constexpr int geometryOrder = GeometryHelper::geometryOrder;
+
+ static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
public:
using GlobalPosition = typename T::GlobalPosition;
@@ -86,7 +109,7 @@ public:
boundary_(is.boundary()),
axisData_(geometryHelper.axisData()),
- pairData_(geometryHelper.pairData()),
+ pairData_(std::move(geometryHelper.pairData())),
localFaceIdx_(geometryHelper.localFaceIndex()),
dirIdx_(geometryHelper.directionIndex()),
outerNormalSign_(sign(unitOuterNormal_[directionIndex()])),
@@ -206,19 +229,19 @@ public:
}
//! Returns the data for one sub face
- const PairData<Scalar, GlobalPosition, geometryOrder>& pairData(const int idx) const
+ const PairData& pairData(const int idx) const
{
return pairData_[idx];
}
//! Return an array of all pair data
- const auto& pairData() const
+ const std::array<PairData, numPairs>& pairData() const
{
return pairData_;
}
//! Return an array of all pair data
- const auto& axisData() const
+ const AxisData& axisData() const
{
return axisData_;
}
@@ -237,7 +260,7 @@ public:
*/
bool hasParallelNeighbor(const int localSubFaceIdx, const int parallelDegreeIdx) const
{
- return !(pairData(localSubFaceIdx).parallelDofs[parallelDegreeIdx] < 0);
+ return pairData(localSubFaceIdx).hasParallelNeighbor[parallelDegreeIdx];
}
/*!
@@ -247,7 +270,7 @@ public:
*/
bool hasOuterNormal(const int localSubFaceIdx) const
{
- return !(pairData_[localSubFaceIdx].normalPair.second < 0);
+ return pairData(localSubFaceIdx).hasNormalNeighbor;
}
/*!
@@ -255,9 +278,10 @@ public:
*
* \param backwardIdx The index describing how many faces backward this dof is from the opposite face
*/
+ template<bool enable = useHigherOrder, std::enable_if_t<enable, int> = 0>
bool hasBackwardNeighbor(const int backwardIdx) const
{
- return !(axisData().inAxisBackwardDofs[backwardIdx] < 0);
+ return axisData().hasBackwardNeighbor[backwardIdx];
}
/*!
@@ -265,9 +289,10 @@ public:
*
* \param forwardIdx The index describing how many faces forward this dof is of the self face
*/
+ template<bool enable = useHigherOrder, std::enable_if_t<enable, int> = 0>
bool hasForwardNeighbor(const int forwardIdx) const
{
- return !(axisData().inAxisForwardDofs[forwardIdx] < 0);
+ return axisData().hasForwardNeighbor[forwardIdx];
}
//! Returns the dof of the face
@@ -325,8 +350,8 @@ private:
int dofIdx_;
Scalar selfToOppositeDistance_;
- AxisData<Scalar, geometryOrder> axisData_;
- std::array<PairData<Scalar, GlobalPosition, geometryOrder>, numPairs> pairData_;
+ AxisData axisData_;
+ std::array<PairData, numPairs> pairData_;
int localFaceIdx_;
unsigned int dirIdx_;
diff --git a/dumux/discretization/staggered/fvgridgeometry.hh b/dumux/discretization/staggered/fvgridgeometry.hh
index 2b7e9b22f68f18f07a9b1d2ebc413ddc95f69d7b..ccffc03b5d279c3a6fec46d2287146346796b3bb 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 int upwindStencilOrder()
+ { 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 int upwindStencilOrder()
+ { 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..0865bf6e21d2bb9ce611e198c6b34fc5b9df6942 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/staggeredupwindmethods.hh>
namespace Dumux {
@@ -38,14 +38,16 @@ namespace Dumux {
* \tparam P The problem type
* \tparam FVC The flux variables cache type
*/
-template<class P, class FVC, class FVCF>
+template<class P, class FVC, class FVCF, int upwOrder>
struct StaggeredDefaultGridFluxVariablesCacheTraits
{
using Problem = P;
using FluxVariablesCache = FVC;
using FluxVariablesCacheFiller = FVCF;
+
template<class GridFluxVariablesCache, bool cachingEnabled>
using LocalView = StaggeredElementFluxVariablesCache<GridFluxVariablesCache, cachingEnabled>;
+ static constexpr int upwindSchemeOrder = upwOrder;
};
/*!
@@ -56,7 +58,8 @@ template<class Problem,
class FluxVariablesCache,
class FluxVariablesCacheFiller,
bool EnableGridFluxVariablesCache = false,
- class Traits = StaggeredDefaultGridFluxVariablesCacheTraits<Problem, FluxVariablesCache, FluxVariablesCacheFiller> >
+ int upwindSchemeOrder = 1,
+ class Traits = StaggeredDefaultGridFluxVariablesCacheTraits<Problem, FluxVariablesCache, FluxVariablesCacheFiller, upwindSchemeOrder>>
class StaggeredGridFluxVariablesCache;
/*!
@@ -64,17 +67,18 @@ class StaggeredGridFluxVariablesCache;
* \brief Flux variables cache class for staggered models.
Specialization in case of storing the flux cache.
*/
-template<class P, class FVC, class FVCF, class Traits>
-class StaggeredGridFluxVariablesCache<P, FVC, FVCF, true, Traits>
+template<class P, class FVC, class FVCF, int upwindSchemeOrder, class Traits>
+class StaggeredGridFluxVariablesCache<P, FVC, FVCF, true, upwindSchemeOrder, Traits>
{
using Problem = typename Traits::Problem;
- using ThisType = StaggeredGridFluxVariablesCache<P, FVC, FVCF, true, Traits>;
+ using ThisType = StaggeredGridFluxVariablesCache<P, FVC, FVCF, true, upwindSchemeOrder, Traits>;
public:
//! export the flux variable cache type
using FluxVariablesCache = typename Traits::FluxVariablesCache;
using Scalar = typename FluxVariablesCache::Scalar;
+ static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
//! export the flux variable cache filler type
using FluxVariablesCacheFiller = typename Traits::FluxVariablesCacheFiller;
@@ -86,7 +90,7 @@ public:
StaggeredGridFluxVariablesCache(const Problem& problem, const std::string& paramGroup = "")
: problemPtr_(&problem)
- , higherOrderApproximation_(paramGroup)
+ , staggeredUpwindMethods_(paramGroup)
{}
// When global caching is enabled, precompute transmissibilities and stencils for all the scv faces
@@ -123,10 +127,10 @@ public:
}
}
- //! Return the HigherOrderApproximation
- const HigherOrderApproximation<Scalar>& higherOrderApproximation() const
+ //! Return the StaggeredUpwindMethods
+ const StaggeredUpwindMethods<Scalar, upwindSchemeOrder>& staggeredUpwindMethods() const
{
- return higherOrderApproximation_;
+ return staggeredUpwindMethods_;
}
const Problem& problem() const
@@ -141,7 +145,7 @@ public:
private:
const Problem* problemPtr_;
- HigherOrderApproximation<Scalar> higherOrderApproximation_;
+ StaggeredUpwindMethods<Scalar, upwindSchemeOrder> staggeredUpwindMethods_;
std::vector<FluxVariablesCache> fluxVarsCache_;
std::vector<std::size_t> globalScvfIndices_;
@@ -152,11 +156,11 @@ private:
* \brief Flux variables cache class for staggered models.
Specialization in case of not storing the flux cache.
*/
-template<class P, class FVC, class FVCF, class Traits>
-class StaggeredGridFluxVariablesCache<P, FVC, FVCF, false, Traits>
+template<class P, class FVC, class FVCF, int upwindSchemeOrder, class Traits>
+class StaggeredGridFluxVariablesCache<P, FVC, FVCF, false, upwindSchemeOrder, Traits>
{
using Problem = typename Traits::Problem;
- using ThisType = StaggeredGridFluxVariablesCache<P, FVC, FVCF, false, Traits>;
+ using ThisType = StaggeredGridFluxVariablesCache<P, FVC, FVCF, false, upwindSchemeOrder, Traits>;
public:
//! export the flux variable cache type
@@ -174,7 +178,7 @@ public:
StaggeredGridFluxVariablesCache(const Problem& problem, const std::string& paramGroup = "")
: problemPtr_(&problem)
- , higherOrderApproximation_(paramGroup)
+ , staggeredUpwindMethods_(paramGroup)
{}
// When global caching is enabled, precompute transmissibilities and stencils for all the scv faces
@@ -187,15 +191,15 @@ public:
const Problem& problem() const
{ return *problemPtr_; }
- //! Return the HigherOrderApproximation
- const HigherOrderApproximation<Scalar>& higherOrderApproximation() const
+ //! Return the UpwindingMethods
+ const StaggeredUpwindMethods<Scalar, upwindSchemeOrder>& staggeredUpwindMethods() const
{
- return higherOrderApproximation_;
+ return staggeredUpwindMethods_;
}
private:
const Problem* problemPtr_;
- HigherOrderApproximation<Scalar> higherOrderApproximation_;
+ StaggeredUpwindMethods<Scalar, upwindSchemeOrder> staggeredUpwindMethods_;
};
diff --git a/dumux/freeflow/CMakeLists.txt b/dumux/freeflow/CMakeLists.txt
index 5f16de206e056526ae1a904c3c0261c80257a51b..e427be44131d694075407e521615b6a64ff1907d 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
+staggeredupwindmethods.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 a39cd4824f2866c0e120b8c90e98ab13e3781b03..b177c16e251ad13b21fac33e0b95cff24244bbd5 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 "staggeredupwindfluxvariables.hh"
+
namespace Dumux {
// forward declaration
@@ -86,6 +86,9 @@ class NavierStokesFluxVariablesImpl<TypeTag, DiscretizationMethod::staggered>
static constexpr auto cellCenterIdx = FVGridGeometry::cellCenterIdx();
static constexpr auto faceIdx = FVGridGeometry::faceIdx();
+ static constexpr int upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
+ static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
+
public:
using HeatConductionType = GetPropType<TypeTag, Properties::HeatConductionType>;
@@ -200,11 +203,6 @@ public:
// The velocities of the dof at interest and the one of the opposite scvf.
const Scalar velocitySelf = elemFaceVars[scvf].velocitySelf();
const Scalar velocityOpposite = elemFaceVars[scvf].velocityOpposite();
- static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight");
-
- // 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()];
// Advective flux.
if (problem.enableInertiaTerms())
@@ -212,45 +210,14 @@ 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;
- // Check if the the velocity of the dof at interest lies up- or downstream w.r.t. to the transporting velocity.
- const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity);
-
- Scalar momentum = 0.0;
-
- // Variables that will store the velocities of interest:
- // velocities[0]: downstream velocity
- // velocities[1]: upstream velocity
- // velocities[2]: upstream-upstream velocity
- std::array<Scalar, 3> velocities{0.0, 0.0, 0.0};
-
- // Variables that will store the distances between the dofs of interest:
- // distances[0]: upstream to downstream distance
- // distances[1]: upstream-upstream to upstream distance
- // distances[2]: downstream staggered cell size
- std::array<Scalar, 3> distances{0.0, 0.0, 0.0};
-
- 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(), upwindWeight);
- }
- else
- momentum = selfIsUpstream ? highOrder.upwind(velocityOpposite, velocitySelf, insideVolVars.density(), upwindWeight)
- : highOrder.upwind(velocitySelf, velocityOpposite, insideVolVars.density(), upwindWeight);
-
- // 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();
+ frontalFlux += StaggeredUpwindFluxVariables<TypeTag, upwindSchemeOrder>::computeUpwindedFrontalMomentum(scvf, elemFaceVars, elemVolVars, gridFluxVarsCache, transportingVelocity)
+ * transportingVelocity * -1.0 * scvf.directionSign();
}
+ // 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()];
+
// Diffusive flux.
// The velocity gradient already accounts for the orientation
// of the staggered face's outer normal vector.
@@ -310,7 +277,7 @@ public:
const GridFluxVariablesCache& gridFluxVarsCache)
{
FacePrimaryVariables normalFlux(0.0);
- auto& faceVars = elemFaceVars[scvf];
+ const auto& faceVars = elemFaceVars[scvf];
const int numSubFaces = scvf.pairData().size();
// Account for all sub faces.
@@ -367,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;
}
@@ -413,54 +386,10 @@ private:
// of interest is located.
const Scalar transportingVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
- static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight", 1.0);
-
- // Check whether the own or the neighboring element is upstream.
- const bool selfIsUpstream = ( normalFace.directionSign() == sign(transportingVelocity) );
-
- // Get the volume variables of the own and the neighboring element
- const auto& insideVolVars = elemVolVars[normalFace.insideScvIdx()];
- const auto& outsideVolVars = elemVolVars[normalFace.outsideScvIdx()];
-
- Scalar momentum = 0.0;
-
- // 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(), upwindWeight);
- }
- 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(), upwindWeight)
- : highOrder.upwind(faceVars.velocitySelf(), velocityFirstParallel, outsideVolVars.density(), upwindWeight);
- }
-
- // 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);
+ return StaggeredUpwindFluxVariables<TypeTag, upwindSchemeOrder>::computeUpwindedLateralMomentum(problem, fvGeometry, element, scvf, normalFace, elemVolVars, faceVars,
+ gridFluxVarsCache, localSubFaceIdx, lateralFaceHasDirichletPressure,
+ lateralFaceHasBJS)
+ * transportingVelocity * normalFace.directionSign() * normalFace.area() * 0.5 * extrusionFactor_(elemVolVars, normalFace);
}
/*!
@@ -547,7 +476,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.
@@ -613,7 +544,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
@@ -657,151 +589,50 @@ private:
}
/*!
- * \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.
+ * \brief Return a velocity value from a boundary for which the boundary conditions have to be checked.
*
- * \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
+ * \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
*/
- bool canFrontalSecondOrder_(const SubControlVolumeFace& ownScvf,
- const bool selfIsUpstream,
- std::array<Scalar, 3>& velocities,
- std::array<Scalar, 3>& distances,
- const FaceVariables& faceVars) const
+ Scalar getParallelVelocityFromOtherBoundary_(const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf,
+ const int localIdx,
+ const Element& boundaryElement,
+ const Scalar parallelVelocity) 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();
+ // 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);
- 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();
+ // 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 (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;
+ if (bcTypes.isDirichlet(Indices::velocity(scvf.directionIndex())))
+ {
+ const SubControlVolumeFace ghostFace = makeParallelGhostFace_(scvf, localIdx);
+ return problem.dirichlet(boundaryElement, ghostFace)[Indices::velocity(scvf.directionIndex())];
}
- }
-
- /*!
- * \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)
+ else if (bcTypes.isSymmetry() || bcTypes.isDirichlet(Indices::pressureIdx))
+ return parallelVelocity;
+ else if (bcTypes.isBJS(Indices::velocity(scvf.directionIndex())))
{
- // 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;
+ const SubControlVolumeFace ghostFace = makeParallelGhostFace_(scvf, localIdx);
+ return problem.bjsVelocity(boundaryElement, scvf, boundaryNormalFace, localIdx, parallelVelocity);
}
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;
+ // 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);
}
}
-
/*!
* \brief Return the outer parallel velocity for normal faces that are on the boundary and therefore have no neighbor.
*
@@ -836,51 +667,6 @@ private:
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);
- }
- }
};
} // end namespace Dumux
diff --git a/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..39ae0b74f5c856bf8352fcd5d316e51b0499a879
--- /dev/null
+++ b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh
@@ -0,0 +1,664 @@
+// -*- 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/staggeredupwindmethods.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 StaggeredUpwindFluxVariables
+{
+ 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:
+ /*!
+ * \brief Returns the momentum in the frontal directon.
+ *
+ * Checks if the model has higher order methods enabled and if the scvf in
+ * question is far enough from the boundary such that higher order methods can be employed.
+ * Then the corresponding set of momenta are collected and the prescribed
+ * upwinding method is used to calculate the momentum.
+ */
+ static 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 Returns the momentum in the lateral directon.
+ *
+ * Evaluates which face is upstream.
+ * Checks if the model has higher order methods enabled and if the scvf in
+ * question is far enough from the boundary such that higher order methods can be employed.
+ * Then the corresponding set of momenta are collected and the prescribed
+ * upwinding method is used to calculate the momentum.
+ */
+ static 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 Returns false as higher order methods are not enabled.
+ *
+ * If higher order methods are not prescribed, this function will be called and false is returned..
+ *
+ * \param ownScvf the sub control volume face
+ * \param transportingVelocity The average of the self and opposite velocities.
+ * \param false_type False if higher order methods are not enabled.
+ */
+ static bool canFrontalSecondOrder_(const SubControlVolumeFace& ownScvf,
+ const Scalar transportingVelocity,
+ std::false_type)
+ {
+ return false;
+ }
+
+ /*!
+ * \brief Returns whether or not the face in question is far enough from the wall to handle higher order methods.
+ *
+ * Evaluates which face is upstream.
+ * If the face is upstream, and the scvf has a forward neighbor, higher order methods are possible.
+ * If the face is downstream, and the scvf has a backwards neighbor, higher order methods are possible.
+ * Otherwise, higher order methods are not possible.
+ * If higher order methods are not prescribed, this function will not be called.
+ *
+ * \param ownScvf the sub control volume face
+ * \param transportingVelocity The average of the self and opposite velocities.
+ * \param true_type True if higher order methods are enabled.
+ */
+ static bool canFrontalSecondOrder_(const SubControlVolumeFace& ownScvf,
+ const Scalar transportingVelocity,
+ std::true_type)
+ {
+ // Depending on selfIsUpstream I have to check if I have a forward or a backward neighbor to retrieve
+ const bool selfIsUpstream = ownScvf.directionSign() != sign(transportingVelocity);
+
+ if (selfIsUpstream)
+ {
+ if (ownScvf.hasForwardNeighbor(0))
+ return true;
+ else
+ return false;
+ }
+ else
+ {
+ if (ownScvf.hasBackwardNeighbor(0))
+ return true;
+ else
+ return false;
+ }
+ }
+
+ /*!
+ * \brief Returns an array of the two momenta needed for basic upwinding methods.
+ *
+ * \param scvf The sub control volume face
+ * \param elemFaceVars The element face variables
+ * \param density The given density \f$\mathrm{[kg/m^3]}\f$
+ * \param transportingVelocity The average of the self and opposite velocities.
+ * \param false_type False if higher order methods are not enabled.
+ */
+ static std::array<Scalar, 2> getFrontalUpwindingMomenta_(const SubControlVolumeFace& scvf,
+ const ElementFaceVariables& elemFaceVars,
+ const Scalar density,
+ const Scalar transportingVelocity,
+ std::false_type)
+ {
+ 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 Returns an array of the three momenta needed for higher order upwinding methods.
+ *
+ * \param scvf The sub control volume face
+ * \param elemFaceVars The element face variables
+ * \param density The given density \f$\mathrm{[kg/m^3]}\f$
+ * \param transportingVelocity The average of the self and opposite velocities.
+ * \param true_type True if higher order methods are enabled.
+ */
+ static std::array<Scalar, 3> getFrontalUpwindingMomenta_(const SubControlVolumeFace& scvf,
+ const ElementFaceVariables& elemFaceVars,
+ const Scalar density,
+ const Scalar transportingVelocity,
+ std::true_type)
+ {
+ 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(0) * density;
+ }
+ else
+ {
+ momenta[0] = elemFaceVars[scvf].velocitySelf() * density;
+ momenta[1] = elemFaceVars[scvf].velocityOpposite() * density;
+ momenta[2] = elemFaceVars[scvf].velocityBackward(0) * density;
+ }
+
+ return momenta;
+ }
+
+ /*!
+ * \brief Returns the upwinded momentum for basic upwind schemes
+ *
+ * \param scvf The sub control volume face
+ * \param momenta The momenta to be upwinded
+ * \param transportingVelocity The average of the self and opposite velocities.
+ * \param gridFluxVarsCache The grid flux variables cache
+ */
+ static Scalar doFrontalMomentumUpwinding_(const SubControlVolumeFace& scvf,
+ const std::array<Scalar, 2>& momenta,
+ const Scalar transportingVelocity,
+ const GridFluxVariablesCache& gridFluxVarsCache)
+ {
+ const auto& upwindScheme = gridFluxVarsCache.staggeredUpwindMethods();
+ return upwindScheme.upwind(momenta[0], momenta[1]);
+ }
+
+ /*!
+ * \brief Returns the upwinded momentum for higher order upwind schemes
+ *
+ * \param scvf The sub control volume face
+ * \param momenta The momenta to be upwinded
+ * \param transportingVelocity The average of the self and opposite velocities.
+ * \param gridFluxVarsCache The grid flux variables cache
+ */
+ template<bool enable = useHigherOrder, std::enable_if_t<enable, int> = 0>
+ static Scalar doFrontalMomentumUpwinding_(const SubControlVolumeFace& scvf,
+ const std::array<Scalar, 3>& momenta,
+ const Scalar transportingVelocity,
+ const GridFluxVariablesCache& gridFluxVarsCache)
+ {
+ const auto& upwindScheme = gridFluxVarsCache.staggeredUpwindMethods();
+ const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity);
+ std::array<Scalar,3> distances = getFrontalDistances_(scvf, selfIsUpstream);
+ return upwindScheme.tvd(momenta, distances, selfIsUpstream, upwindScheme.tvdApproach());
+ }
+
+ /*!
+ * \brief Returns an array of distances needed for non-uniform higher order upwind schemes
+ *
+ * \param ownScvf The sub control volume face
+ * \param selfIsUpstream bool describing upstream face.
+ */
+ template<bool enable = useHigherOrder, std::enable_if_t<enable, int> = 0>
+ static std::array<Scalar, 3> getFrontalDistances_(const SubControlVolumeFace& ownScvf,
+ const bool selfIsUpstream)
+ {
+ // 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[0];
+ distances[2] = 0.5 * (ownScvf.axisData().selfToOppositeDistance + ownScvf.axisData().inAxisBackwardDistances[0]);
+ }
+ else
+ {
+ distances[0] = ownScvf.selfToOppositeDistance();
+ distances[1] = ownScvf.axisData().inAxisBackwardDistances[0];
+ distances[2] = 0.5 * (ownScvf.axisData().selfToOppositeDistance + ownScvf.axisData().inAxisForwardDistances[0]);
+ }
+ return distances;
+ }
+
+ /*!
+ * \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
+ * \param true_type True when second order is enabled.
+ */
+ static bool canLateralSecondOrder_(const SubControlVolumeFace& ownScvf,
+ const bool selfIsUpstream,
+ const int localSubFaceIdx,
+ std::true_type)
+ {
+ 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.
+ */
+ static bool canLateralSecondOrder_(const SubControlVolumeFace& ownScvf,
+ const bool selfIsUpstream,
+ const int localSubFaceIdx,
+ std::false_type)
+ { return false; }
+
+
+ /*!
+ * \brief Returns an array of the three momenta needed for higher order upwinding methods.
+ *
+ * Only called if higher order methods are enabled and the scvf can use higher order methods.
+ */
+ static 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)
+ {
+ 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();
+ 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 Returns an array of the two momenta needed for basic upwinding methods.
+ *
+ * Called if higher order methods are not enabled of if the scvf can not use higher order methods.
+ */
+ static 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)
+ {
+ // 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 Returns the upwinded momentum for basic upwind schemes
+ *
+ * Fowards to Frontal Momentum Upwinding method
+ */
+ static Scalar doLateralMomentumUpwinding_(const SubControlVolumeFace& scvf,
+ const SubControlVolumeFace& normalScvf,
+ const std::array<Scalar, 2>& momenta,
+ const Scalar transportingVelocity,
+ const int localSubFaceIdx,
+ const GridFluxVariablesCache& gridFluxVarsCache)
+ {
+ return doFrontalMomentumUpwinding_(scvf, momenta, transportingVelocity, gridFluxVarsCache);
+ }
+
+ /*!
+ * \brief Returns the upwinded momentum for higher order upwind schemes
+ *
+ * \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
+ */
+ static Scalar doLateralMomentumUpwinding_(const SubControlVolumeFace& scvf,
+ const SubControlVolumeFace& normalScvf,
+ const std::array<Scalar, 3>& momenta,
+ const Scalar transportingVelocity,
+ const int localSubFaceIdx,
+ const GridFluxVariablesCache& gridFluxVarsCache)
+ {
+ const bool selfIsUpstream = ( normalScvf.directionSign() == sign(transportingVelocity) );
+ const auto& upwindScheme = gridFluxVarsCache.staggeredUpwindMethods();
+ std::array<Scalar,3> distances = getLateralDistances_(scvf, localSubFaceIdx, selfIsUpstream);
+ return upwindScheme.tvd(momenta, distances, selfIsUpstream, upwindScheme.tvdApproach());
+ }
+
+ /*!
+ * \brief Returns an array of distances needed for non-uniform higher order upwind schemes
+ *
+ * \param ownScvf The sub control volume face
+ * \param localSubFaceIdx The local index of the subface
+ * \param selfIsUpstream bool describing upstream face.
+ */
+ template<bool enable = useHigherOrder, std::enable_if_t<enable, int> = 0>
+ static std::array<Scalar, 3> getLateralDistances_(const SubControlVolumeFace& ownScvf,
+ const int localSubFaceIdx,
+ const bool selfIsUpstream)
+ {
+ // 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.pairData(localSubFaceIdx).parallelDistances[1];
+ }
+
+ 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
+ */
+ static 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)
+ {
+ // 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
+ */
+ static Scalar getParallelVelocityFromOtherBoundary_(const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf,
+ const int localIdx,
+ const Element& boundaryElement,
+ const Scalar parallelVelocity)
+ {
+ // 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
+ static SubControlVolumeFace makeGhostFace_(const SubControlVolumeFace& ownScvf, const GlobalPosition& pos)
+ {
+ 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
+ static SubControlVolumeFace makeParallelGhostFace_(const SubControlVolumeFace& ownScvf, const int localSubFaceIdx)
+ {
+ 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/staggeredupwindmethods.hh
similarity index 64%
rename from dumux/freeflow/higherorderapproximation.hh
rename to dumux/freeflow/staggeredupwindmethods.hh
index 601bf115d219eec1d22c79af9b71e02e387feee5..26754bff295067f88a2bf60bf6df7fa8a623f6cb 100644
--- a/dumux/freeflow/higherorderapproximation.hh
+++ b/dumux/freeflow/staggeredupwindmethods.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>
@@ -47,17 +47,19 @@ enum class DifferencingScheme
/**
* \brief This file contains different higher order methods for approximating the velocity.
*/
-template<class Scalar>
-class HigherOrderApproximation
+template<class Scalar, int upwindSchemeOrder>
+class StaggeredUpwindMethods
{
public:
- HigherOrderApproximation(const std::string& paramGroup = "")
+ StaggeredUpwindMethods(const std::string& paramGroup = "")
{
- if (hasParamInGroup(paramGroup, "Discretization.TvdApproach"))
+ upwindWeight_ = getParamFromGroup<Scalar>(paramGroup, "Flux.UpwindWeight");
+
+ if (upwindSchemeOrder > 1)
{
// 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", "Uniform"));
+ differencingScheme_ = differencingSchemeFromString(getParamFromGroup<std::string>(paramGroup, "Flux.DifferencingScheme", "Minmod"));
// Assign the limiter_ depending on the differencing scheme
switch (differencingScheme_)
@@ -107,6 +109,33 @@ public:
break;
}
}
+
+ if (!hasParamInGroup(paramGroup, "Flux.TvdApproach"))
+ {
+ std::cout << "No TvdApproach specified. Defaulting to the Uniform method." << "\n";
+ std::cout << "Other 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 \n";
+ std::cout << "Each approach can be specified as written above in the Flux group under the title TvdApproach in your input file. \n";
+ }
+ if (!hasParamInGroup(paramGroup, "Flux.DifferencingScheme"))
+ {
+
+ std::cout << "No DifferencingScheme specified. Defaulting to the Minmod scheme." << "\n";
+ std::cout << "Other available Differencing Schemes are as follows: \n"
+ << " " << differencingSchemeToString(DifferencingScheme::vanleer) << ": The Vanleer flux limiter\n"
+ << " " << differencingSchemeToString(DifferencingScheme::vanalbada) << ": The Vanalbada flux limiter\n"
+ << " " << differencingSchemeToString(DifferencingScheme::minmod) << ": The Minmod flux limiter\n"
+ << " " << differencingSchemeToString(DifferencingScheme::superbee) << ": The Superbee flux limiter\n"
+ << " " << differencingSchemeToString(DifferencingScheme::umist) << ": The Umist flux limiter\n"
+ << " " << differencingSchemeToString(DifferencingScheme::mclimiter) << ": The Mclimiter flux limiter\n"
+ << " " << differencingSchemeToString(DifferencingScheme::wahyd) << ": The Wahyd flux limiter";
+ std::cout << "Each scheme can be specified as written above in the Flux group under the variable DifferencingScheme in your input file. \n";
+ }
+
+ std::cout << "Using the tvdApproach \"" << tvdApproachToString(tvdApproach_)
+ << "\" and the differencing Scheme \" " << differencingSchemeToString(differencingScheme_) << "\" \n";
}
else
{
@@ -127,9 +156,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");
}
/**
@@ -152,13 +181,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"
@@ -177,13 +206,13 @@ 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
}
}
@@ -191,68 +220,19 @@ public:
/**
* \brief Upwind Method
*/
- Scalar upwind(const Scalar downstreamVelocity,
- const Scalar upstreamVelocity,
- const Scalar density,
- const Scalar upwindWeight) const
- {
- return (upwindWeight * upstreamVelocity + (1.0 - upwindWeight) * downstreamVelocity) * density;
- }
-
- /**
- * \brief Central Differencing Method
- */
- Scalar centralDifference(const Scalar downstreamVelocity,
- const Scalar upstreamVelocity,
- const Scalar density) const
- {
- return (0.5 * (upstreamVelocity + downstreamVelocity)) * density;
- }
-
- /**
- * \brief Linear Upwind Method
- */
- Scalar linearUpwind(const Scalar upstreamVelocity,
- const Scalar upUpstreamVelocity,
- const Scalar upstreamToDownstreamDistance,
- const Scalar upUpstreamToUpstreamDistance,
- const Scalar density) const
+ Scalar upwind(const Scalar downstreamMomentum,
+ const Scalar upstreamMomentum) 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;
+ return (upwindWeight_ * upstreamMomentum + (1.0 - upwindWeight_) * downstreamMomentum);
}
/**
* \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;
@@ -260,17 +240,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:
@@ -287,22 +267,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;
}
/**
@@ -312,31 +294,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;
}
/**
@@ -345,16 +329,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
@@ -362,11 +348,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;
}
/**
@@ -451,10 +437,11 @@ public:
private:
TvdApproach tvdApproach_;
DifferencingScheme differencingScheme_;
+ Scalar upwindWeight_;
std::function<Scalar(const Scalar, const Scalar)> limiter_;
};
} // end namespace Dumux
-#endif // DUMUX_HIGHER_ORDER_VELOCITY_APPROXIMATION_HH
+#endif // DUMUX_UPWINDING_METHODS_HH
diff --git a/test/discretization/staggered/test_staggered_free_flow_geometry.cc b/test/discretization/staggered/test_staggered_free_flow_geometry.cc
index 59c12bf843499cb381d89778a0fbda8a7dc3497a..d2f5387926f24149271b46fbf3d26bc45795860e 100644
--- a/test/discretization/staggered/test_staggered_free_flow_geometry.cc
+++ b/test/discretization/staggered/test_staggered_free_flow_geometry.cc
@@ -53,13 +53,14 @@ public:
} // end namespace Detail
//! the fv grid geometry traits for this test
-template<class GridView>
+template<class GridView, int upwOrder>
struct TestFVGGTraits : public DefaultMapperTraits<GridView>
{
using SubControlVolume = CCSubControlVolume<GridView>;
- using SubControlVolumeFace = FreeFlowStaggeredSubControlVolumeFace<GridView>;
+ using SubControlVolumeFace = FreeFlowStaggeredSubControlVolumeFace<GridView, upwOrder>;
using IntersectionMapper = ConformingGridIntersectionMapper<GridView>;
- using GeometryHelper = FreeFlowStaggeredGeometryHelper<GridView>;
+ using GeometryHelper = FreeFlowStaggeredGeometryHelper<GridView, upwOrder>;
+ static constexpr int upwindSchemeOrder = upwOrder;
struct DofTypeIndices
{
@@ -70,8 +71,8 @@ struct TestFVGGTraits : public DefaultMapperTraits<GridView>
template<class FVGridGeometry>
using ConnectivityMap = StaggeredFreeFlowConnectivityMap<FVGridGeometry>;
- template<class FVGridGeometry, bool enableCache>
- using LocalView = StaggeredFVElementGeometry<FVGridGeometry, enableCache>;
+ template<class FVGridGeometry, bool cachingEnabled>
+ using LocalView = StaggeredFVElementGeometry<FVGridGeometry, cachingEnabled>;
};
} // end namespace Dumux
@@ -92,8 +93,10 @@ int main (int argc, char *argv[]) try
constexpr int dim = Grid::dimension;
+ static constexpr int upwindSchemeOrder = 2;
+
using FVGridGeometry = StaggeredFVGridGeometry<typename Grid::LeafGridView, /*enable caching=*/ true,
- TestFVGGTraits<typename Grid::LeafGridView> >;
+ TestFVGGTraits<typename Grid::LeafGridView, upwindSchemeOrder> >;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
@@ -159,11 +162,11 @@ int main (int argc, char *argv[]) try
std::cout << std::fixed << std::left << std::setprecision(2);
std::cout << " On Axis Dof Index: \n";
- if(fvGridGeometry.order() > 1)
+ if(fvGridGeometry.upwindStencilOrder() > 1)
{std::cout << " | Forward dofIdx : " << std::setw(3) << scvf.axisData().inAxisForwardDofs[0] << "\n";}
std::cout << " | Self dofIdx : " << std::setw(3) << scvf.dofIndex() << "\n";
std::cout << " | Opposite dofIdx : " << std::setw(3) << scvf.dofIndexOpposingFace() << "\n";
- if(fvGridGeometry.order() > 1)
+ if(fvGridGeometry.upwindStencilOrder() > 1)
{std::cout << " | Backward dofIdx : " << std::setw(3) << scvf.axisData().inAxisBackwardDofs[0] << "\n";}
std::cout << " Normal Dof Index: \n";
@@ -176,22 +179,22 @@ int main (int argc, char *argv[]) try
std::cout << " Parallel Dof Index: \n";
for(int i = 0; i < scvf.pairData().size(); i++)
{
- for(int j = 0; j < fvGridGeometry.order(); j++)
+ for(int j = 0; j < fvGridGeometry.upwindStencilOrder(); j++)
{
std::cout << " | Parallel Dof "<< j << " on axis " << i << ": "<< std::setw(3) << scvf.pairData(i).parallelDofs[j] << "\n";
}
}
std::cout << " Distances: \n";
- if(fvGridGeometry.order() > 1)
+ if(fvGridGeometry.upwindStencilOrder() > 1)
{std::cout << " | Opposite To Backwards Face Dist : " << std::setw(3) << scvf.axisData().inAxisBackwardDistances[0] << "\n";}
std::cout << " | self To Opposite Dist : " << std::setw(3) << scvf.selfToOppositeDistance() << "\n";
- if(fvGridGeometry.order() > 1)
- {std::cout << " | Opposite To Backwards Face Dist : " << std::setw(3) << scvf.axisData().inAxisBackwardDistances[0] << "\n";}
+ if(fvGridGeometry.upwindStencilOrder() > 1)
+ {std::cout << " | self To Forwards Face Dist : " << std::setw(3) << scvf.axisData().inAxisForwardDistances[0] << "\n";}
for(int i = 0; i < scvf.pairData().size(); i++)
{
- for(int j = 0; j < fvGridGeometry.order(); j++)
+ for(int j = 0; j < fvGridGeometry.upwindStencilOrder(); j++)
{
std::cout << " | Parallel Distance "<< j << " on axis " << i << ": "<< std::setw(3) << scvf.pairData(i).parallelDistances[j] << "\n";
}
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/problem.hh b/test/freeflow/navierstokes/kovasznay/problem.hh
index bbbb7ab1ec1dab8c282c93977486df2f1cbda0d9..1c5374412fe7afaf30be5a6806a74c5d4296f55a 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::upwindStencilOrder() << "\n";
kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
Scalar reynoldsNumber = 1.0 / kinematicViscosity_;
lambda_ = 0.5 * reynoldsNumber