From 92caa4f013b77108f087fa1ea5e43cdb2ea25fa2 Mon Sep 17 00:00:00 2001
From: Kilian <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Wed, 27 Feb 2019 17:50:01 +0100
Subject: [PATCH] [freeflow][discretization] Adapt freeflow for compile time
order definition
---
dumux/common/properties.hh | 2 +
.../staggered/freeflow/connectivitymap.hh | 38 +++--
.../staggered/freeflow/facevariables.hh | 113 ++++++++-----
.../freeflow/fvgridgeometrytraits.hh | 6 +-
.../staggered/freeflow/properties.hh | 14 +-
.../freeflow/staggeredgeometryhelper.hh | 91 +++++++---
.../freeflow/subcontrolvolumeface.hh | 11 +-
dumux/freeflow/higherorderapproximation.hh | 12 +-
.../navierstokes/staggered/fluxvariables.hh | 158 ++++++++++++------
9 files changed, 297 insertions(+), 148 deletions(-)
diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index d3ce6cfa36..3dd939e20f 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -234,6 +234,8 @@ 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
+template<class TypeTag, class MyTypeTag>
+struct UpwindGeometryOrder { using type = UndefinedProperty; }; //!< Specifies the order of the upwinding scheme (1 == basic upwinding, 2 == higher order)
/////////////////////////////////////////////////////////////
// Properties used by the mpnc model
diff --git a/dumux/discretization/staggered/freeflow/connectivitymap.hh b/dumux/discretization/staggered/freeflow/connectivitymap.hh
index e43321f1c0..7622726ec6 100644
--- a/dumux/discretization/staggered/freeflow/connectivitymap.hh
+++ b/dumux/discretization/staggered/freeflow/connectivitymap.hh
@@ -54,6 +54,8 @@ class StaggeredFreeFlowConnectivityMap
using Stencil = std::vector<GridIndexType>;
+ static constexpr bool useHigherOrder = false; //TODO
+
public:
//! Update the map and prepare the stencils
@@ -186,24 +188,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(std::integral_constant<bool, useHigherOrder>{}, scvf, stencil);
}
for(const auto& data : scvf.pairData())
@@ -224,6 +211,25 @@ private:
}
}
+ void addHigherOrderInAxisDofs(std::false_type, const SubControlVolumeFace& scvf, Stencil& stencil) {}
+
+ void addHigherOrderInAxisDofs(std::true_type, const SubControlVolumeFace& scvf, Stencil& stencil)
+ {
+ for (int i = 0; i < stencilOrder_ - 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_;
diff --git a/dumux/discretization/staggered/freeflow/facevariables.hh b/dumux/discretization/staggered/freeflow/facevariables.hh
index e227cfdc04..2f1943d9b8 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 geometryOrder>
+struct InAxisVelocities;
+
+template<class Scalar>
+struct InAxisVelocities<Scalar, 1>
+{
+ Scalar self = 0.0;
+ Scalar opposite = 0.0;
+};
+
+template<class Scalar, int geometryOrder>
+struct InAxisVelocities
+{
+ Scalar self = 0.0;
+ Scalar opposite = 0.0;
+ std::array<Scalar, geometryOrder-1> forward{};
+ std::array<Scalar, geometryOrder-1> backward{};
+};
+
+}
+
/*!
* \ingroup StaggeredDiscretization
* \brief The face variables class for free flow staggered grid models.
* Contains all relevant velocities for the assembly of the momentum balance.
*/
-template<class FacePrimaryVariables, int dim>
+template<class FacePrimaryVariables, int dim, int geometryOrder> // TODO doc
class StaggeredFaceVariables
{
static constexpr int numPairs = (dim == 2) ? 2 : 4;
- static constexpr int order = 2;
+ static constexpr bool useHigherOrder = geometryOrder > 1;
+
using Scalar = typename FacePrimaryVariables::block_type;
+ using InAxisVelocities = Detail::InAxisVelocities<Scalar, geometryOrder>;
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(std::integral_constant<bool, useHigherOrder>{}, faceSol, scvf);
// 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 < geometryOrder; 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(std::false_type, const SolVector& faceSol, const SubControlVolumeFace& scvf) {}
+
+ template<class SolVector, class SubControlVolumeFace>
+ void addHigherOrderInAxisVelocities(std::true_type, const SolVector& faceSol, const SubControlVolumeFace& scvf)
+ {
+
+ // 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.axisData().inAxisForwardDofs[i] < 0))
+ 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.axisData().inAxisBackwardDofs[i] < 0))
+ inAxisVelocities_.backward[i] = faceSol[scvf.axisData().inAxisBackwardDofs[i]];
+ }
+ }
+
+ InAxisVelocities inAxisVelocities_;
+ std::array<std::array<Scalar, geometryOrder>, 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 4b88982ba9..c882c55a39 100644
--- a/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh
+++ b/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh
@@ -39,14 +39,14 @@ namespace Dumux {
* \ingroup StaggeredDiscretization
* \brief Default traits for the finite volume grid geometry.
*/
-template<class GridView, class MapperTraits = DefaultMapperTraits<GridView>>
+template<class GridView, int geometryOrder, class MapperTraits = DefaultMapperTraits<GridView>>
struct StaggeredFreeFlowDefaultFVGridGeometryTraits
: public MapperTraits
{
using SubControlVolume = CCSubControlVolume<GridView>;
- using SubControlVolumeFace = FreeFlowStaggeredSubControlVolumeFace<GridView>;
+ using SubControlVolumeFace = FreeFlowStaggeredSubControlVolumeFace<GridView, geometryOrder>;
using IntersectionMapper = ConformingGridIntersectionMapper<GridView>;
- using GeometryHelper = FreeFlowStaggeredGeometryHelper<GridView>;
+ using GeometryHelper = FreeFlowStaggeredGeometryHelper<GridView, geometryOrder>;
struct DofTypeIndices
{
diff --git a/dumux/discretization/staggered/freeflow/properties.hh b/dumux/discretization/staggered/freeflow/properties.hh
index 3075fc93c6..3c568c7d29 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 upwindGeometryOrder = getPropValue<TypeTag, Properties::UpwindGeometryOrder>();
static constexpr bool enableCache = getPropValue<TypeTag, Properties::EnableFVGridGeometryCache>();
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Traits = StaggeredFreeFlowDefaultFVGridGeometryTraits<GridView, upwindGeometryOrder>;
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 upwindGeometryOrder = getPropValue<TypeTag, Properties::UpwindGeometryOrder>();
public:
- using type = StaggeredFaceVariables<FacePrimaryVariables, GridView::dimension>;
+ using type = StaggeredFaceVariables<FacePrimaryVariables, GridView::dimension, upwindGeometryOrder>;
};
//! 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 UpwindGeometryOrder<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 6990c0daed..0872330b0e 100644
--- a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
+++ b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
@@ -43,13 +43,20 @@ template<class Scalar, class GlobalPosition, int geometryOrder>
struct PairData
{
std::array<int, geometryOrder> parallelDofs;
- std::array<Scalar, geometryOrder+1> parallelDistances;
- std::pair<signed int,signed int> normalPair;
+ std::array<Scalar, geometryOrder+1/*geometryOrder*/> parallelDistances; // TODO
+ std::pair<signed int, signed int> normalPair;
int localNormalFaceIdx;
Scalar normalDistance;
GlobalPosition virtualFirstParallelFaceDofPos;
};
+
+/*!
+ * \ingroup StaggeredDiscretization
+ * \brief In Axis Data stored per sub face
+ */
+template<class Scalar, int geometryOrder>
+struct AxisData;
/*!
* \ingroup StaggeredDiscretization
* \brief In Axis Data stored per sub face
@@ -57,6 +64,14 @@ struct PairData
template<class Scalar, int geometryOrder>
struct AxisData
{
+ AxisData()
+ {
+ inAxisForwardDofs.fill(-1);
+ inAxisBackwardDofs.fill(-1);
+ inAxisForwardDistances.fill(0);
+ inAxisBackwardDistances.fill(0);
+ }
+
int selfDof;
int oppositeDof;
std::array<int, geometryOrder-1> inAxisForwardDofs;
@@ -66,6 +81,21 @@ struct AxisData
std::array<Scalar, geometryOrder-1> inAxisBackwardDistances;
};
+/*!
+ * \ingroup StaggeredDiscretization
+ * \brief In Axis Data stored per sub face
+ */
+template<class Scalar>
+struct AxisData<Scalar, 1>
+{
+ int selfDof;
+ int oppositeDof;
+ Scalar selfToOppositeDistance;
+};
+
+// template<class Scalar, int geometryOrder>
+// using AxisData
+
/*!
* \ingroup StaggeredDiscretization
* \brief Returns the dirction index of the facet (0 = x, 1 = y, 2 = z)
@@ -83,7 +113,7 @@ inline static unsigned int directionIndex(Vector&& vector)
* \brief Helper class constructing the dual grid finite volume geometries
* for the free flow staggered discretization method.
*/
-template<class GridView>
+template<class GridView, int geometryOrder>
class FreeFlowStaggeredGeometryHelper
{
using Scalar = typename GridView::ctype;
@@ -107,6 +137,9 @@ class FreeFlowStaggeredGeometryHelper
static constexpr int numfacets = dimWorld * 2;
static constexpr int numPairs = 2 * (dimWorld - 1);
+ static constexpr bool useHigherOrder = geometryOrder > 1;
+
+
public:
FreeFlowStaggeredGeometryHelper(const Element& element, const GridView& gridView) : element_(element), elementGeometry_(element.geometry()), gridView_(gridView)
@@ -118,8 +151,8 @@ public:
{
intersection_ = intersection;
innerNormalFacePos_.clear();
- fillPairData_();
- fillAxisData_();
+ fillPairData_(std::integral_constant<bool, true>{}); //TODO
+ fillAxisData_(std::integral_constant<bool, useHigherOrder>{});
}
/*!
@@ -162,33 +195,44 @@ public:
return Dumux::directionIndex(std::move(intersection.centerUnitOuterNormal()));
}
- static constexpr int geometryOrder = 2;
-
private:
+
+
+ void fillAxisData_(std::false_type)
+ {
+ const auto inIdx = intersection_.indexInInside();
+ const auto oppIdx = localOppositeIdx_(inIdx);
+
+ // Set the self Dof
+ axisData_.selfDof = gridView_.indexSet().subIndex(intersection_.inside(), inIdx, codimIntersection);
+
+ // Set the opposite Dof
+ axisData_.oppositeDof = gridView_.indexSet().subIndex(intersection_.inside(), oppIdx, codimIntersection);
+
+ // Set the Self to Opposite Distance
+ const auto self = getFacet_(inIdx, element_);
+ const auto opposite = getFacet_(oppIdx, element_);
+ axisData_.selfToOppositeDistance = (self.geometry().center() - opposite.geometry().center()).two_norm();
+ }
/*!
* \brief Fills all entries of the in axis data
*/
- void fillAxisData_()
+ void fillAxisData_(std::true_type)
{
- unsigned int numForwardBackwardAxisDofs = axisData_.inAxisForwardDofs.size();
+ const auto numForwardBackwardAxisDofs = axisData_.inAxisForwardDofs.size();
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();
// Set the Forward Dofs
std::stack<Element> inAxisForwardElementStack;
@@ -298,10 +342,17 @@ private:
}
}
+ void fillPairData_(std::false_type)
+ {
+
+ // TODO
+
+ }
+
/*!
* \brief Fills all entries of the pair data
*/
- void fillPairData_()
+ void fillPairData_(std::true_type)
{
// initialize values that could remain unitialized if the intersection lies on a boundary
for(auto& data : pairData_)
diff --git a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
index 0fa6efc607..8dcc9b9327 100644
--- a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
@@ -42,23 +42,24 @@ namespace Dumux {
* of a sub control volume we compute fluxes on. This is a specialization for free flow models.
*/
template<class GV,
+ int geometryOrder,
class T = StaggeredDefaultScvfGeometryTraits<GV> >
class FreeFlowStaggeredSubControlVolumeFace
-: public SubControlVolumeFaceBase<FreeFlowStaggeredSubControlVolumeFace<GV, T>, T>
+: public SubControlVolumeFaceBase<FreeFlowStaggeredSubControlVolumeFace<GV, geometryOrder, T>, T>
{
- using ThisType = FreeFlowStaggeredSubControlVolumeFace<GV, T>;
+ using ThisType = FreeFlowStaggeredSubControlVolumeFace<GV, geometryOrder, 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 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 = geometryOrder > 1;
public:
using GlobalPosition = typename T::GlobalPosition;
@@ -255,6 +256,7 @@ 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);
@@ -265,6 +267,7 @@ 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);
diff --git a/dumux/freeflow/higherorderapproximation.hh b/dumux/freeflow/higherorderapproximation.hh
index 601bf115d2..3e98307415 100644
--- a/dumux/freeflow/higherorderapproximation.hh
+++ b/dumux/freeflow/higherorderapproximation.hh
@@ -53,6 +53,8 @@ class HigherOrderApproximation
public:
HigherOrderApproximation(const std::string& paramGroup = "")
{
+ upwindWeight_ = getParamFromGroup<Scalar>(paramGroup, "Flux.UpwindWeight");
+
if (hasParamInGroup(paramGroup, "Discretization.TvdApproach"))
{
// Read the runtime parameters
@@ -188,15 +190,14 @@ public:
}
}
+
/**
* \brief Upwind Method
*/
- Scalar upwind(const Scalar downstreamVelocity,
- const Scalar upstreamVelocity,
- const Scalar density,
- const Scalar upwindWeight) const
+ Scalar upwind(const Scalar downstreamMomentum,
+ const Scalar upstreamMomentum) const
{
- return (upwindWeight * upstreamVelocity + (1.0 - upwindWeight) * downstreamVelocity) * density;
+ return (upwindWeight_ * upstreamMomentum + (1.0 - upwindWeight_) * downstreamMomentum);
}
/**
@@ -451,6 +452,7 @@ public:
private:
TvdApproach tvdApproach_;
DifferencingScheme differencingScheme_;
+ Scalar upwindWeight_;
std::function<Scalar(const Scalar, const Scalar)> limiter_;
};
diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
index a39cd4824f..37ad3fd3cb 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
@@ -86,6 +86,8 @@ class NavierStokesFluxVariablesImpl<TypeTag, DiscretizationMethod::staggered>
static constexpr auto cellCenterIdx = FVGridGeometry::cellCenterIdx();
static constexpr auto faceIdx = FVGridGeometry::faceIdx();
+ static constexpr bool useHigherOrder = false; // TODO
+
public:
using HeatConductionType = GetPropType<TypeTag, Properties::HeatConductionType>;
@@ -200,7 +202,8 @@ 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");
+
+ // const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity);
// The volume variables within the current element. We only require those (and none of neighboring elements)
// because the fluxes are calculated over the staggered face at the center of the element.
@@ -212,39 +215,37 @@ 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;
+ const auto upwindingMomenta = getFrontalUpwindingMomenta_(scvf, elemFaceVars, insideVolVars.density(), transportingVelocity);
+ const Scalar momentum = doMomentumUpwinding_(scvf, upwindingMomenta, gridFluxVarsCache);
// Variables that will store the velocities of interest:
// velocities[0]: downstream velocity
// velocities[1]: upstream velocity
// velocities[2]: upstream-upstream velocity
- std::array<Scalar, 3> velocities{0.0, 0.0, 0.0};
+ // 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};
+ // std::array<Scalar, 3> distances{0.0, 0.0, 0.0};
- const auto& highOrder = gridFluxVarsCache.higherOrderApproximation();
+ // 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);
+ // if (highOrder.tvdApproach() != TvdApproach::none)
+ // {
+ // if (canFrontalSecondOrder_(scvf, selfIsUpstream, velocities, distances, elemFaceVars[scvf]))
+ // {
+ // momentum = highOrder.tvd(velocities[0], velocities[1], velocities[2], distances[0], distances[1], distances[2], selfIsUpstream, insideVolVars.density(), highOrder.tvdApproach());
+ // }
+ // else
+ // momentum = highOrder.upwind(velocities[0], velocities[1], insideVolVars.density());
+ // }
+ // else
+ // momentum = selfIsUpstream ? highOrder.upwind(velocityOpposite, velocitySelf, insideVolVars.density())
+ // : highOrder.upwind(velocitySelf, velocityOpposite, insideVolVars.density());
// Account for the orientation of the staggered face's normal outer normal vector
// (pointing in opposite direction of the scvf's one).
@@ -310,7 +311,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.
@@ -413,50 +414,45 @@ private:
// of interest is located.
const Scalar transportingVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
- static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight", 1.0);
+ const auto parallelUpwindingMomenta = getParallelUpwindingMomenta_(problem, element, scvf, normalFace, elemVolVars, faceVars,
+ transportingVelocity, localSubFaceIdx,
+ lateralFaceHasDirichletPressure, lateralFaceHasBJS);
- // 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;
+ const Scalar momentum = doMomentumUpwinding_(scvf, parallelUpwindingMomenta, gridFluxVarsCache);
// Variables that will store the velocities of interest:
// velocities[0]: downstream velocity
// velocities[1]: upsteram velocity
// velocities[2]: upstream-upstream velocity
- std::array<Scalar, 3> velocities{0.0, 0.0, 0.0};
+ // 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);
- }
+ // std::array<Scalar, 3> distances{0.0, 0.0, 0.0};
+
+ // const auto& highOrder = gridFluxVarsCache.higherOrderApproximation();
+ //
+ // // If a Tvd approach has been specified and I am not too near to the boundary I can use a second order approximation.
+ // if (highOrder.tvdApproach() != TvdApproach::none)
+ // {
+ // if (canLateralSecondOrder_(scvf, fvGeometry, selfIsUpstream, localSubFaceIdx, velocities, distances, problem, element, faceVars, lateralFaceHasDirichletPressure, lateralFaceHasBJS))
+ // {
+ // momentum = highOrder.tvd(velocities[0], velocities[1], velocities[2], distances[0], distances[1], distances[2], selfIsUpstream, selfIsUpstream ? insideVolVars.density() : outsideVolVars.density(), highOrder.tvdApproach());
+ // }
+ // else
+ // momentum = highOrder.upwind(velocities[0], velocities[1], selfIsUpstream ? insideVolVars.density() : outsideVolVars.density());
+ // }
+ // else
+ // {
+ // const Scalar velocityFirstParallel = scvf.hasParallelNeighbor(localSubFaceIdx, 0)
+ // ? faceVars.velocityParallel(localSubFaceIdx, 0)
+ // : getParallelVelocityFromBoundary_(problem, scvf, normalFace, faceVars.velocitySelf(), localSubFaceIdx, element, lateralFaceHasDirichletPressure, lateralFaceHasBJS);
+ //
+ // momentum = selfIsUpstream ? highOrder.upwind(velocityFirstParallel, faceVars.velocitySelf(), insideVolVars.density())
+ // : highOrder.upwind(faceVars.velocitySelf(), velocityFirstParallel, outsideVolVars.density());
+ // }
// Account for the orientation of the staggered normal face's outer normal vector
// and its area (0.5 of the coinciding scfv).
@@ -881,6 +877,60 @@ private:
DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the boundary conditions for the momentum equations at global position " << boundarySubFaceCenter);
}
}
+
+ Scalar doMomentumUpwinding_(const SubControlVolumeFace& scvf,
+ const std::array<Scalar, 2>& momentum,
+ const GridFluxVariablesCache& gridFluxVarsCache) const
+ {
+ const auto& upwindScheme = gridFluxVarsCache.higherOrderApproximation();
+ return upwindScheme.upwind(momentum[0], momentum[1]);
+ }
+
+ std::array<Scalar, 2> getFrontalUpwindingMomenta_(const SubControlVolumeFace& scvf,
+ const ElementFaceVariables& elemFaceVars,
+ const Scalar density,
+ const Scalar transportingVelocity) const
+ {
+ const Scalar momentumSelf = elemFaceVars[scvf].velocitySelf() * density;
+ const Scalar momentumOpposite = elemFaceVars[scvf].velocityOpposite() * density;
+ const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity);
+
+ return selfIsUpstream ? std::array<Scalar, 2>{momentumOpposite, momentumSelf}
+ : std::array<Scalar, 2>{momentumSelf, momentumOpposite};
+ }
+
+ std::array<Scalar, 2> getParallelUpwindingMomenta_(const Problem& problem,
+ const Element& element,
+ const SubControlVolumeFace& ownFace,
+ const SubControlVolumeFace& lateralFace,
+ const ElementVolumeVariables& elemVolVars,
+ const FaceVariables& faceVars,
+ const Scalar transportingVelocity,
+ const int localSubFaceIdx,
+ bool lateralFaceHasDirichletPressure,
+ bool lateralFaceHasBJS) const
+ {
+
+ // Check whether the own or the neighboring element is upstream.
+ const bool selfIsUpstream = lateralFace.directionSign() == sign(transportingVelocity);
+
+ // Get the volume variables of the own and the neighboring element
+ const auto& insideVolVars = elemVolVars[lateralFace.insideScvIdx()];
+ const auto& outsideVolVars = elemVolVars[lateralFace.outsideScvIdx()];
+
+ const Scalar momentumSelf = faceVars.velocitySelf() * insideVolVars.density();
+
+ const Scalar momentumParallel = ownFace.hasParallelNeighbor(localSubFaceIdx, 0)
+ ? faceVars.velocityParallel(localSubFaceIdx, 0) * outsideVolVars.density()
+ : (getParallelVelocityFromBoundary_(problem, ownFace, lateralFace,
+ faceVars.velocitySelf(), localSubFaceIdx, element,
+ lateralFaceHasDirichletPressure, lateralFaceHasBJS)
+ * insideVolVars.density());
+
+ return selfIsUpstream ? std::array<Scalar, 2>{momentumParallel, momentumSelf}
+ : std::array<Scalar, 2>{momentumSelf, momentumParallel};
+ }
+
};
} // end namespace Dumux
--
GitLab