From bc83cbdd850d3aa4a1fb61dd44c1ba67760e6206 Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Tue, 5 Mar 2019 15:57:01 +0100
Subject: [PATCH] [staggered][freeflow][fluxvariables] documentation, cleanup,
 and renaming

---
 .../freeflow/fvgridgeometrytraits.hh          |   8 +-
 .../freeflow/subcontrolvolumeface.hh          |  21 +-
 .../staggered/fvgridgeometry.hh               |   8 +-
 .../staggered/gridfluxvariablescache.hh       |  20 +-
 dumux/freeflow/CMakeLists.txt                 |   2 +-
 .../navierstokes/staggered/fluxvariables.hh   |  17 +-
 ...les.hh => staggeredupwindfluxvariables.hh} | 184 +++++++++---------
 ...ngmethods.hh => staggeredupwindmethods.hh} |   4 +-
 .../navierstokes/kovasznay/problem.hh         |   2 +-
 9 files changed, 146 insertions(+), 120 deletions(-)
 rename dumux/freeflow/navierstokes/staggered/{upwindvariables.hh => staggeredupwindfluxvariables.hh} (85%)
 rename dumux/freeflow/{upwindingmethods.hh => staggeredupwindmethods.hh} (99%)

diff --git a/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh b/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh
index d07a0780e0..0032f2cb00 100644
--- a/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh
+++ b/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh
@@ -39,15 +39,15 @@ namespace Dumux {
  * \ingroup StaggeredDiscretization
  * \brief Default traits for the finite volume grid geometry.
  */
-template<class GridView, int upwindSchemeOrder, class MapperTraits = DefaultMapperTraits<GridView>>
+template<class GridView, int USO, class MapperTraits = DefaultMapperTraits<GridView>>
 struct StaggeredFreeFlowDefaultFVGridGeometryTraits
 : public MapperTraits
 {
     using SubControlVolume = CCSubControlVolume<GridView>;
-    using SubControlVolumeFace = FreeFlowStaggeredSubControlVolumeFace<GridView, upwindSchemeOrder>;
+    using SubControlVolumeFace = FreeFlowStaggeredSubControlVolumeFace<GridView, USO>;
     using IntersectionMapper = ConformingGridIntersectionMapper<GridView>;
-    using GeometryHelper = FreeFlowStaggeredGeometryHelper<GridView, upwindSchemeOrder>;
-    static constexpr int UpwindSchemeOrder = upwindSchemeOrder;
+    using GeometryHelper = FreeFlowStaggeredGeometryHelper<GridView, USO>;
+    static constexpr int upwindSchemeOrder = USO;
 
     struct DofTypeIndices
     {
diff --git a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
index a812393899..fe7800a2ca 100644
--- a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
@@ -36,6 +36,25 @@
 
 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
@@ -43,7 +62,7 @@ namespace Dumux {
  */
 template<class GV,
          int upwindSchemeOrder,
-         class T = StaggeredDefaultScvfGeometryTraits<GV> >
+         class T = FreeFlowStaggeredDefaultScvfGeometryTraits<GV, upwindSchemeOrder>>
 class FreeFlowStaggeredSubControlVolumeFace
 : public SubControlVolumeFaceBase<FreeFlowStaggeredSubControlVolumeFace<GV, upwindSchemeOrder, T>, T>
 {
diff --git a/dumux/discretization/staggered/fvgridgeometry.hh b/dumux/discretization/staggered/fvgridgeometry.hh
index c38909a699..ccffc03b5d 100644
--- a/dumux/discretization/staggered/fvgridgeometry.hh
+++ b/dumux/discretization/staggered/fvgridgeometry.hh
@@ -191,7 +191,7 @@ class StaggeredFVGridGeometry<GV, true, Traits>
 public:
     //! export discretization method
     static constexpr DiscretizationMethod discMethod = DiscretizationMethod::staggered;
-    static constexpr int upwindSchemeOrder = Traits::UpwindSchemeOrder;
+    static constexpr int upwindSchemeOrder = Traits::upwindSchemeOrder;
     static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
 
     //! export the type of the fv element geometry (the local view type)
@@ -214,7 +214,7 @@ public:
     { return typename DofTypeIndices::FaceIdx{}; }
 
     //! The order of the stencil built
-    static constexpr std::size_t order()
+    static constexpr int upwindStencilOrder()
     {   return upwindSchemeOrder; }
 
     using CellCenterFVGridGeometryType = CellCenterFVGridGeometry<ThisType>;
@@ -448,7 +448,7 @@ class StaggeredFVGridGeometry<GV, false, Traits>
 public:
     //! export discretization method
     static constexpr DiscretizationMethod discMethod = DiscretizationMethod::staggered;
-    static constexpr int upwindSchemeOrder = Traits::UpwindSchemeOrder;
+    static constexpr int upwindSchemeOrder = Traits::upwindSchemeOrder;
     static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
 
     using GeometryHelper = typename Traits::GeometryHelper;
@@ -473,7 +473,7 @@ public:
     { return typename DofTypeIndices::FaceIdx{}; }
 
     //! The order of the stencil built
-    static constexpr std::size_t order()
+    static constexpr int upwindStencilOrder()
     {   return upwindSchemeOrder; }
 
     using CellCenterFVGridGeometryType = CellCenterFVGridGeometry<ThisType>;
diff --git a/dumux/discretization/staggered/gridfluxvariablescache.hh b/dumux/discretization/staggered/gridfluxvariablescache.hh
index 44d9550018..7a258483d4 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/upwindingmethods.hh>
+#include <dumux/freeflow/staggeredupwindmethods.hh>
 
 namespace Dumux {
 
@@ -86,7 +86,7 @@ public:
 
     StaggeredGridFluxVariablesCache(const Problem& problem, const std::string& paramGroup = "")
     : problemPtr_(&problem)
-    , upwindingMethods_(paramGroup)
+    , staggeredUpwindMethods_(paramGroup)
       {}
 
     // When global caching is enabled, precompute transmissibilities and stencils for all the scv faces
@@ -123,10 +123,10 @@ public:
         }
     }
 
-    //! Return the UpwindingMethods
-    const UpwindingMethods<Scalar>& upwindingMethods() const
+    //! Return the StaggeredUpwindMethods
+    const StaggeredUpwindMethods<Scalar>& staggeredUpwindMethods() const
     {
-        return upwindingMethods_;
+        return staggeredUpwindMethods_;
     }
 
     const Problem& problem() const
@@ -141,7 +141,7 @@ public:
 
 private:
     const Problem* problemPtr_;
-    UpwindingMethods<Scalar> upwindingMethods_;
+    StaggeredUpwindMethods<Scalar> staggeredUpwindMethods_;
 
     std::vector<FluxVariablesCache> fluxVarsCache_;
     std::vector<std::size_t> globalScvfIndices_;
@@ -174,7 +174,7 @@ public:
 
     StaggeredGridFluxVariablesCache(const Problem& problem, const std::string& paramGroup = "")
     : problemPtr_(&problem)
-    , upwindingMethods_(paramGroup)
+    , staggeredUpwindMethods_(paramGroup)
       {}
 
     // When global caching is enabled, precompute transmissibilities and stencils for all the scv faces
@@ -188,14 +188,14 @@ public:
     { return *problemPtr_; }
 
     //! Return the UpwindingMethods
-    const UpwindingMethods<Scalar>& upwindingMethods() const
+    const StaggeredUpwindMethods<Scalar>& staggeredUpwindMethods() const
     {
-        return upwindingMethods_;
+        return staggeredUpwindMethods_;
     }
 
 private:
     const Problem* problemPtr_;
-    UpwindingMethods<Scalar> upwindingMethods_;
+    StaggeredUpwindMethods<Scalar> staggeredUpwindMethods_;
 
 };
 
diff --git a/dumux/freeflow/CMakeLists.txt b/dumux/freeflow/CMakeLists.txt
index d63c7b0937..e427be4413 100644
--- a/dumux/freeflow/CMakeLists.txt
+++ b/dumux/freeflow/CMakeLists.txt
@@ -7,5 +7,5 @@ install(FILES
 properties.hh
 turbulenceproperties.hh
 volumevariables.hh
-upwindingmethods.hh
+staggeredupwindmethods.hh
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow)
diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
index 08d0362184..b177c16e25 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
@@ -34,7 +34,7 @@
 #include <dumux/flux/fluxvariablesbase.hh>
 #include <dumux/discretization/method.hh>
 
-#include "upwindvariables.hh"
+#include "staggeredupwindfluxvariables.hh"
 
 namespace Dumux {
 
@@ -204,20 +204,20 @@ public:
         const Scalar velocitySelf = elemFaceVars[scvf].velocitySelf();
         const Scalar velocityOpposite = elemFaceVars[scvf].velocityOpposite();
 
-        // 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())
         {
             // 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;
-            Dumux::UpwindVariables<TypeTag, upwindSchemeOrder> upwindVariables;
-            frontalFlux += upwindVariables.computeUpwindedFrontalMomentum(scvf, elemFaceVars, elemVolVars, gridFluxVarsCache, transportingVelocity)
+
+            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.
@@ -386,8 +386,7 @@ private:
         // of interest is located.
         const Scalar transportingVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
 
-        Dumux::UpwindVariables<TypeTag, upwindSchemeOrder> upwindVariables;
-        return upwindVariables.computeUpwindedLateralMomentum(problem, fvGeometry, element, scvf, normalFace, elemVolVars, faceVars,
+        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);
diff --git a/dumux/freeflow/navierstokes/staggered/upwindvariables.hh b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh
similarity index 85%
rename from dumux/freeflow/navierstokes/staggered/upwindvariables.hh
rename to dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh
index b5a7bc6d05..62a49aaa3a 100644
--- a/dumux/freeflow/navierstokes/staggered/upwindvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh
@@ -32,7 +32,7 @@
 #include <dumux/common/properties.hh>
 
 #include <dumux/discretization/method.hh>
-#include <dumux/freeflow/upwindingmethods.hh>
+#include <dumux/freeflow/staggeredupwindmethods.hh>
 
 namespace Dumux {
 
@@ -41,7 +41,7 @@ namespace Dumux {
  * \brief The upwinding variables class for the Navier-Stokes model using the staggered grid discretization.
  */
 template<class TypeTag, int upwindSchemeOrder>
-class UpwindVariables
+class StaggeredUpwindFluxVariables
 {
     using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
 
@@ -76,13 +76,13 @@ class UpwindVariables
     static constexpr auto faceIdx = FVGridGeometry::faceIdx();
 
 public:
-    UpwindVariables() {};
-
     /*!
-     * \brief
-     *
-     * TODO: DOC
+     * \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.
      */
     FacePrimaryVariables computeUpwindedFrontalMomentum(const SubControlVolumeFace& scvf,
                                                         const ElementFaceVariables& elemFaceVars,
@@ -110,10 +110,13 @@ public:
     }
 
     /*!
-     * \brief
-     *
-     * TODO: DOC
+     * \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.
      */
     FacePrimaryVariables computeUpwindedLateralMomentum(const Problem& problem,
                                                         const FVElementGeometry& fvGeometry,
@@ -131,7 +134,6 @@ public:
         // 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) );
 
@@ -160,18 +162,42 @@ public:
 
 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.
+     */
+    bool canFrontalSecondOrder_(const SubControlVolumeFace& ownScvf,
+                                const Scalar transportingVelocity,
+                                std::false_type) const
+    {
+        return false;
+    }
+
     /*!
-     * \brief
+     * \brief Returns whether or not the face in question is far enough from the wall to handle higher order methods.
      *
-     * TODO: DOC
+     *        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.
      */
     bool canFrontalSecondOrder_(const SubControlVolumeFace& ownScvf,
                                 const Scalar transportingVelocity,
                                 std::true_type) const
     {
-        const bool selfIsUpstream = ownScvf.directionSign() != sign(transportingVelocity);
         // Depending on selfIsUpstream I have to check if I have a forward or a backward neighbor to retrieve
+        const bool selfIsUpstream = ownScvf.directionSign() != sign(transportingVelocity);
+
         if (selfIsUpstream)
         {
             if (ownScvf.hasForwardNeighbor(0))
@@ -189,23 +215,36 @@ private:
     }
 
     /*!
-     * \brief
-     *
-     * TODO: DOC
+     * \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.
      */
-    bool canFrontalSecondOrder_(const SubControlVolumeFace& ownScvf,
-                                const Scalar transportingVelocity,
-                                std::false_type) const
+    std::array<Scalar, 2> getFrontalUpwindingMomenta_(const SubControlVolumeFace& scvf,
+                                                      const ElementFaceVariables& elemFaceVars,
+                                                      const Scalar density,
+                                                      const Scalar transportingVelocity,
+                                                      std::false_type) const
     {
-        return false;
+        const Scalar momentumSelf = elemFaceVars[scvf].velocitySelf() * density;
+        const Scalar momentumOpposite = elemFaceVars[scvf].velocityOpposite() * density;
+        const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity);
+
+        return selfIsUpstream ? std::array<Scalar, 2>{momentumOpposite, momentumSelf}
+                              : std::array<Scalar, 2>{momentumSelf, momentumOpposite};
     }
 
     /*!
-     * \brief
-     *
-     * TODO: DOC
+     * \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.
      */
     std::array<Scalar, 3> getFrontalUpwindingMomenta_(const SubControlVolumeFace& scvf,
                                                       const ElementFaceVariables& elemFaceVars,
@@ -220,58 +259,42 @@ private:
         {
             momenta[0] = elemFaceVars[scvf].velocityOpposite() * density;
             momenta[1] = elemFaceVars[scvf].velocitySelf() * density;
-            momenta[2] = elemFaceVars[scvf].velocityForward(upwindSchemeOrder-2) * 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(upwindSchemeOrder-2) * density;
+            momenta[2] = elemFaceVars[scvf].velocityBackward(0) * density;
         }
 
         return momenta;
     }
 
     /*!
-     * \brief
-     *
-     * TODO: DOC
-     *
-     */
-    std::array<Scalar, 2> getFrontalUpwindingMomenta_(const SubControlVolumeFace& scvf,
-                                                      const ElementFaceVariables& elemFaceVars,
-                                                      const Scalar density,
-                                                      const Scalar transportingVelocity,
-                                                      std::false_type) const
-    {
-        const Scalar momentumSelf = elemFaceVars[scvf].velocitySelf() * density;
-        const Scalar momentumOpposite = elemFaceVars[scvf].velocityOpposite() * density;
-        const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity);
-
-        return selfIsUpstream ? std::array<Scalar, 2>{momentumOpposite, momentumSelf}
-                              : std::array<Scalar, 2>{momentumSelf, momentumOpposite};
-    }
-
-    /*!
-     * \brief
-     *
-     * TODO: DOC
+     * \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
      */
     Scalar doFrontalMomentumUpwinding_(const SubControlVolumeFace& scvf,
                                 const std::array<Scalar, 2>& momenta,
                                 const Scalar transportingVelocity,
                                 const GridFluxVariablesCache& gridFluxVarsCache) const
     {
-        const auto& upwindScheme = gridFluxVarsCache.upwindingMethods();
+        const auto& upwindScheme = gridFluxVarsCache.staggeredUpwindMethods();
         return upwindScheme.upwind(momenta[0], momenta[1]);
     }
 
     /*!
-     * \brief
-     *
-     * TODO: DOC
+     * \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>
     Scalar doFrontalMomentumUpwinding_(const SubControlVolumeFace& scvf,
@@ -279,17 +302,17 @@ private:
                                        const Scalar transportingVelocity,
                                        const GridFluxVariablesCache& gridFluxVarsCache) const
     {
-        const auto& upwindScheme = gridFluxVarsCache.upwindingMethods();
+        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
-     *
-     * TODO: DOC
+     * \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>
     std::array<Scalar, 3> getFrontalDistances_(const SubControlVolumeFace& ownScvf,
@@ -302,21 +325,19 @@ private:
         if (selfIsUpstream)
         {
             distances[0] = ownScvf.selfToOppositeDistance();
-            distances[1] = ownScvf.axisData().inAxisForwardDistances[upwindSchemeOrder-2];
+            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[upwindSchemeOrder-2];
+            distances[1] = ownScvf.axisData().inAxisBackwardDistances[0];
             distances[2] = 0.5 * (ownScvf.axisData().selfToOppositeDistance + ownScvf.axisData().inAxisForwardDistances[0]);
         }
         return distances;
     }
 
     /*!
-     * TODO: DOC
-     *
      * \brief Check if a second order approximation for the lateral part of the advective term can be used
      *
      * This helper function checks if the scvf of interest is not too near to the
@@ -325,6 +346,7 @@ private:
      * \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.
      */
     bool canLateralSecondOrder_(const SubControlVolumeFace& ownScvf,
                                 const bool selfIsUpstream,
@@ -359,10 +381,9 @@ private:
 
 
     /*!
-     * \brief
-     *
-     * TODO: DOC
+     * \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.
      */
     std::array<Scalar, 3> getLateralUpwindingMomenta_(const Problem& problem,
                                                       const FVElementGeometry& fvGeometry,
@@ -412,17 +433,6 @@ private:
         else
         {
             momenta[0] = faceVars.velocitySelf() * outsideVolVars.density();
-
-            if (!ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
-            {
-                std::cout << "how did we get here \n";
-                momenta[1] = getParallelVelocityFromBoundary_(problem, ownScvf, lateralFace,
-                                                              faceVars.velocitySelf(), localSubFaceIdx,
-                                                              element, lateralFaceHasDirichletPressure,
-                                                              lateralFaceHasBJS) * outsideVolVars.density();
-                return momenta;
-            }
-
             momenta[1] = faceVars.velocityParallel(localSubFaceIdx, 0) * outsideVolVars.density();
 
             // If there is another parallel neighbor I can assign the "upstream-upstream" velocity, otherwise I retrieve it from the boundary.
@@ -442,10 +452,9 @@ private:
     }
 
     /*!
-     * \brief
-     *
-     * TODO: DOC
+     * \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.
      */
     std::array<Scalar, 2> getLateralUpwindingMomenta_(const Problem& problem,
                                                        const FVElementGeometry& fvGeometry,
@@ -481,10 +490,9 @@ private:
     }
 
     /*!
-     * \brief
-     *
-     * TODO: DOC
+     * \brief Returns the upwinded momentum for basic upwind schemes
      *
+     *        Fowards to Frontal Momentum Upwinding method
      */
     Scalar doLateralMomentumUpwinding_(const SubControlVolumeFace& scvf,
                                        const SubControlVolumeFace& normalScvf,
@@ -497,9 +505,7 @@ private:
     }
 
     /*!
-     * \brief
-     *
-     * TODO: DOC
+     * \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
@@ -516,15 +522,17 @@ private:
                                        const GridFluxVariablesCache& gridFluxVarsCache) const
     {
         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
-     *
-     * TODO: DOC
+     * \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>
     std::array<Scalar, 3> getLateralDistances_(const SubControlVolumeFace& ownScvf,
@@ -550,7 +558,7 @@ private:
         {
             distances[0] = ownScvf.cellCenteredParallelDistance(localSubFaceIdx, 0);
             distances[1] = ownScvf.cellCenteredParallelDistance(localSubFaceIdx, 1);
-            distances[2] = ownScvf.area();
+            distances[2] = ownScvf.pairData(localSubFaceIdx).parallelDistances[1];
         }
 
         return distances;
diff --git a/dumux/freeflow/upwindingmethods.hh b/dumux/freeflow/staggeredupwindmethods.hh
similarity index 99%
rename from dumux/freeflow/upwindingmethods.hh
rename to dumux/freeflow/staggeredupwindmethods.hh
index 18f985de6a..98e32b7533 100644
--- a/dumux/freeflow/upwindingmethods.hh
+++ b/dumux/freeflow/staggeredupwindmethods.hh
@@ -48,10 +48,10 @@ enum class DifferencingScheme
   * \brief This file contains different higher order methods for approximating the velocity.
   */
 template<class Scalar>
-class UpwindingMethods
+class StaggeredUpwindMethods
 {
 public:
-    UpwindingMethods(const std::string& paramGroup = "")
+    StaggeredUpwindMethods(const std::string& paramGroup = "")
     {
         upwindWeight_ = getParamFromGroup<Scalar>(paramGroup, "Flux.UpwindWeight");
 
diff --git a/test/freeflow/navierstokes/kovasznay/problem.hh b/test/freeflow/navierstokes/kovasznay/problem.hh
index 8a228736b3..1c5374412f 100644
--- a/test/freeflow/navierstokes/kovasznay/problem.hh
+++ b/test/freeflow/navierstokes/kovasznay/problem.hh
@@ -113,7 +113,7 @@ public:
     : ParentType(fvGridGeometry), eps_(1e-6)
     {
         printL2Error_ = getParam<bool>("Problem.PrintL2Error");
-        std::cout<< "upwindSchemeOrder is: " << FVGridGeometry::order() << "\n";
+        std::cout<< "upwindSchemeOrder is: " << FVGridGeometry::upwindStencilOrder() << "\n";
         kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
         Scalar reynoldsNumber = 1.0 / kinematicViscosity_;
         lambda_ = 0.5 * reynoldsNumber
-- 
GitLab