From 999c58b63ae9b7ff1af0bff467de4f6aee977518 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Sun, 10 Dec 2017 10:32:42 +0100
Subject: [PATCH] [navierstokesnc] Adapt to NavierStokes structure

---
 .../freeflow/navierstokesnc/fluxvariables.hh  | 134 +--------------
 dumux/freeflow/navierstokesnc/indices.hh      |   2 +-
 .../freeflow/navierstokesnc/localresidual.hh  | 141 +--------------
 dumux/freeflow/navierstokesnc/properties.hh   |  35 ++--
 .../navierstokesnc/staggered/fluxvariables.hh | 129 ++++++++++++++
 .../navierstokesnc/staggered/localresidual.hh | 161 ++++++++++++++++++
 .../navierstokesnc/volumevariables.hh         |   3 +-
 .../navierstokesnc/vtkoutputfields.hh         |   2 +-
 .../staggerednc/channeltestproblem.hh         |  14 +-
 .../staggerednc/densityflowproblem.hh         |  12 +-
 10 files changed, 350 insertions(+), 283 deletions(-)
 create mode 100644 dumux/freeflow/navierstokesnc/staggered/fluxvariables.hh
 create mode 100644 dumux/freeflow/navierstokesnc/staggered/localresidual.hh

diff --git a/dumux/freeflow/navierstokesnc/fluxvariables.hh b/dumux/freeflow/navierstokesnc/fluxvariables.hh
index 8639c26ebf..a68f6e0624 100644
--- a/dumux/freeflow/navierstokesnc/fluxvariables.hh
+++ b/dumux/freeflow/navierstokesnc/fluxvariables.hh
@@ -25,23 +25,16 @@
 
 #include <dumux/common/properties.hh>
 #include <dumux/discretization/fluxvariablesbase.hh>
-#include "../staggered/fluxvariables.hh"
+#include <dumux/freeflow/navierstokes/fluxvariables.hh>
+#include <dumux/freeflow/navierstokesnc/staggered/fluxvariables.hh>
 
 namespace Dumux
 {
 
-namespace Properties
-{
-// forward declaration
-NEW_PROP_TAG(EnableComponentTransport);
-NEW_PROP_TAG(EnableEnergyBalance);
-NEW_PROP_TAG(EnableInertiaTerms);
-NEW_PROP_TAG(ElementFaceVariables);
-}
 
-// // forward declaration
-// template<class TypeTag, bool enableComponentTransport, bool enableEnergyBalance>
-// class FreeFlowFluxVariablesImpl;
+// forward declaration
+template<class TypeTag, DiscretizationMethods Method>
+class NavierStokesNCFluxVariablesImpl;
 
 /*!
  * \ingroup ImplicitModel
@@ -49,124 +42,9 @@ NEW_PROP_TAG(ElementFaceVariables);
  *        specializations are provided for combinations of physical processes
  * \note  Not all specializations are currently implemented
  */
-
-// specialization for miscible, isothermal flow
 template<class TypeTag>
-class FreeFlowFluxVariablesImpl<TypeTag, true> : public FreeFlowFluxVariablesImpl<TypeTag, false>
-{
-    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
-    using Element = typename GridView::template Codim<0>::Entity;
-    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
-    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
-    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
-    using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
-    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
-    using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
-    using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
-    using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
-    using IndexType = typename GridView::IndexSet::IndexType;
-    using Stencil = std::vector<IndexType>;
-
-    using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
-
-    static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
-    static constexpr auto numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
-
-    static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
-
-    //! The index of the component balance equation that gets replaced with the total mass balance
-    static const int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
-
-    using ParentType = FreeFlowFluxVariablesImpl<TypeTag, false>;
-
-    enum {
-         // grid and world dimension
-        dim = GridView::dimension,
-        dimWorld = GridView::dimensionworld,
-
-        pressureIdx = Indices::pressureIdx,
-        velocityIdx = Indices::velocityIdx,
-
-        massBalanceIdx = Indices::massBalanceIdx,
-        momentumBalanceIdx = Indices::momentumBalanceIdx,
-        conti0EqIdx = Indices::conti0EqIdx,
-        phaseIdx = Indices::phaseIdx
-    };
-
-public:
-
-    CellCenterPrimaryVariables computeFluxForCellCenter(const Problem& problem,
-                                                        const Element &element,
-                                                        const FVElementGeometry& fvGeometry,
-                                                        const ElementVolumeVariables& elemVolVars,
-                                                        const ElementFaceVariables& elemFaceVars,
-                                                        const SubControlVolumeFace &scvf,
-                                                        const FluxVariablesCache& fluxVarsCache)
-    {
-        CellCenterPrimaryVariables flux(0.0);
-
-        flux += advectiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, elemFaceVars, scvf);
-        flux += MolecularDiffusionType::diffusiveFluxForCellCenter(problem, fvGeometry, elemVolVars, scvf);
-        return flux;
-    }
-
-private:
-
-    CellCenterPrimaryVariables advectiveFluxForCellCenter_(const Problem& problem,
-                                                          const FVElementGeometry& fvGeometry,
-                                                          const ElementVolumeVariables& elemVolVars,
-                                                          const ElementFaceVariables& elemFaceVars,
-                                                          const SubControlVolumeFace &scvf)
-    {
-        CellCenterPrimaryVariables flux(0.0);
-
-        const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
-        const auto& insideVolVars = elemVolVars[insideScv];
-
-        const Scalar velocity = elemFaceVars[scvf].velocitySelf();
-
-        const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity);
-        const Scalar upWindWeight = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.UpwindWeight");
-
-        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
-        {
-            // get equation index
-            const auto eqIdx = conti0EqIdx + compIdx;
-
-            bool isOutflow = false;
-            if(scvf.boundary())
-            {
-                const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
-                    if(bcTypes.isOutflow(eqIdx))
-                        isOutflow = true;
-            }
-
-            const auto& outsideVolVars = isOutflow ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
-            const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
-            const auto& downstreamVolVars = insideIsUpstream ? outsideVolVars : insideVolVars;
-
-            const Scalar upstreamDensity = useMoles ? upstreamVolVars.molarDensity() : upstreamVolVars.density();
-            const Scalar downstreamDensity = useMoles ? downstreamVolVars.molarDensity() : downstreamVolVars.density();
-            const Scalar upstreamFraction = useMoles ? upstreamVolVars.moleFraction(phaseIdx, compIdx) : upstreamVolVars.massFraction(phaseIdx, compIdx);
-            const Scalar downstreamFraction = useMoles ? downstreamVolVars.moleFraction(phaseIdx, compIdx) : downstreamVolVars.massFraction(phaseIdx, compIdx);
-
-            flux[eqIdx] = (upWindWeight * upstreamDensity * upstreamFraction +
-                          (1.0 - upWindWeight) * downstreamDensity * downstreamFraction)
-                          * velocity;
-        }
-
-        // in case one balance is substituted by the total mass balance
-        if (replaceCompEqIdx < numComponents)
-        {
-            flux[replaceCompEqIdx] = std::accumulate(flux.begin(), flux.end(), 0.0);
-        }
+using NavierStokesNCFluxVariables = NavierStokesNCFluxVariablesImpl<TypeTag, GET_PROP_VALUE(TypeTag, DiscretizationMethod)>;
 
-        flux *= scvf.area() * sign(scvf.outerNormalScalar());
-        return flux;
-    }
-};
 
 } // end namespace
 
diff --git a/dumux/freeflow/navierstokesnc/indices.hh b/dumux/freeflow/navierstokesnc/indices.hh
index d6a83c1eff..f295926bbc 100644
--- a/dumux/freeflow/navierstokesnc/indices.hh
+++ b/dumux/freeflow/navierstokesnc/indices.hh
@@ -23,7 +23,7 @@
 #ifndef DUMUX_STAGGERED_NAVIERSTOKES_NC_INDICES_HH
 #define DUMUX_STAGGERED_NAVIERSTOKES_NC_INDICES_HH
 
-#include <dumux/freeflow/staggered/indices.hh>
+#include <dumux/freeflow/navierstokes/indices.hh>
 #include <dumux/common/properties.hh>
 
 namespace Dumux
diff --git a/dumux/freeflow/navierstokesnc/localresidual.hh b/dumux/freeflow/navierstokesnc/localresidual.hh
index c426ed8aee..d4f6aa5d67 100644
--- a/dumux/freeflow/navierstokesnc/localresidual.hh
+++ b/dumux/freeflow/navierstokesnc/localresidual.hh
@@ -20,152 +20,29 @@
  * \file
  * \brief Calculates the residual of models based on the box scheme element-wise.
  */
-#ifndef DUMUX_STAGGERED_NAVIERSTOKES_NC_LOCAL_RESIDUAL_HH
-#define DUMUX_STAGGERED_NAVIERSTOKES_NC_LOCAL_RESIDUAL_HH
+#ifndef DUMUX_NAVIERSTOKES_NC_LOCAL_RESIDUAL_HH
+#define DUMUX_NAVIERSTOKES_NC_LOCAL_RESIDUAL_HH
 
 #include <dumux/common/properties.hh>
-#include <dumux/common/valgrind.hh>
-#include <dumux/implicit/staggered/localresidual.hh>
-
-#include "properties.hh"
+#include <dumux/discretization/methods.hh>
+#include <dumux/freeflow/navierstokes/localresidual.hh>
+#include <dumux/freeflow/navierstokesnc/staggered/localresidual.hh>
 
 namespace Dumux
 {
 
-namespace Properties
-{
-// forward declaration
-NEW_PROP_TAG(EnableComponentTransport);
-NEW_PROP_TAG(EnableEnergyBalance);
-NEW_PROP_TAG(EnableInertiaTerms);
-NEW_PROP_TAG(ReplaceCompEqIdx);
-}
-
 /*!
- * \ingroup CCModel
- * \ingroup StaggeredLocalResidual
- * \brief Element-wise calculation of the residual for models
- *        based on the fully implicit cell-centered scheme.
  *
  * \todo Please doc me more!
  */
 
-
 // // forward declaration
-template<class TypeTag, bool enableComponentTransport>
-class StaggeredNavierStokesResidualImpl;
+template<class TypeTag, DiscretizationMethods Method>
+class NavierStokesNCResidualImpl;
 
-// specialization for miscible, isothermal flow
 template<class TypeTag>
-class StaggeredNavierStokesResidualImpl<TypeTag, true> : public StaggeredNavierStokesResidualImpl<TypeTag, false>
-{
-    using ParentType = StaggeredNavierStokesResidualImpl<TypeTag, false>;
-    friend class StaggeredLocalResidual<TypeTag>;
-    friend ParentType;
-
-    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
-    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
-    using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
-    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
-    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
-    using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
-
-    using CellCenterResidual = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
-
-    using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
-    typename DofTypeIndices::CellCenterIdx cellCenterIdx;
-
-    enum {
-        conti0EqIdx = Indices::conti0EqIdx,
-        phaseIdx = Indices::phaseIdx,
-
-        // The index of the component balance equation
-        // that gets replaced with the total mass balance
-        replaceCompEqIdx = Indices::replaceCompEqIdx
-    };
-
-    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
-    using EnergyLocalResidual = typename GET_PROP_TYPE(TypeTag, EnergyLocalResidual);
-
-    static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
-    static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
-
-public:
-    using ParentType::ParentType;
-
-     /*!
-     * \brief Evaluate the rate of change of all conservation
-     *        quantites (e.g. phase mass) within a sub-control
-     *        volume of a finite volume element for the immiscible models.
-     * \param scv The sub control volume
-     * \param volVars The current or previous volVars
-     * \note This function should not include the source and sink terms.
-     * \note The volVars can be different to allow computing
-     *       the implicit euler time derivative here
-     */
-    CellCenterPrimaryVariables computeStorageForCellCenter(const Problem& problem,
-                                                           const SubControlVolume& scv,
-                                                           const VolumeVariables& volVars) const
-    {
-        CellCenterPrimaryVariables storage(0.0);
-
-        const Scalar density = useMoles ? volVars.molarDensity() : volVars.density();
-
-        // compute storage term of all components
-        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
-        {
-            const int eqIdx = conti0EqIdx + compIdx;
-
-            const Scalar massOrMoleFraction = useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx);
-            const Scalar s =  density * massOrMoleFraction;
-
-            if (eqIdx != replaceCompEqIdx)
-                storage[eqIdx] += s;
-        }
-            // in case one balance is substituted by the total mass balance
-            if(replaceCompEqIdx < numComponents)
-                storage[replaceCompEqIdx] = density;
-
-        EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars);
-
-        return storage;
-    }
-
-protected:
-
-    /*!
-     * \brief Sets a fixed Dirichlet value for a cell (such as pressure) at the boundary.
-     *        This is a provisional alternative to setting the Dirichlet value on the boundary directly.
-     *
-     * \param insideScv The sub control volume
-     * \param elemVolVars The current or previous element volVars
-     * \param bcTypes The boundary types
-     */
-    void setFixedCell_(CellCenterResidual& residual,
-                       const Problem& problem,
-                       const SubControlVolume& insideScv,
-                       const ElementVolumeVariables& elemVolVars,
-                       const BoundaryTypes& bcTypes) const
-    {
-        ParentType::setFixedCell_(residual, problem, insideScv, elemVolVars, bcTypes);
-
-        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
-        {
-            // get equation index
-            const auto eqIdx = conti0EqIdx + compIdx;
-
-            // set a fixed mole fraction for cells
-            if(eqIdx != conti0EqIdx && bcTypes.isDirichletCell(eqIdx))
-            {
-                const auto& insideVolVars = elemVolVars[insideScv];
-                const Scalar massOrMoleFraction = useMoles ? insideVolVars.moleFraction(phaseIdx, compIdx) : insideVolVars.massFraction(phaseIdx, compIdx);
-                residual[eqIdx] = massOrMoleFraction - problem.dirichletAtPos(insideScv.center())[cellCenterIdx][eqIdx];
-            }
-        }
+using NavierStokesNCResidual = NavierStokesNCResidualImpl<TypeTag, GET_PROP_VALUE(TypeTag, DiscretizationMethod)>;
 
-    }
-};
 }
 
-#endif   // DUMUX_CC_LOCAL_RESIDUAL_HH
+#endif   // DUMUX_NAVIERSTOKES_NC_LOCAL_RESIDUAL_HH
diff --git a/dumux/freeflow/navierstokesnc/properties.hh b/dumux/freeflow/navierstokesnc/properties.hh
index e3cf94efcb..124f041d8c 100644
--- a/dumux/freeflow/navierstokesnc/properties.hh
+++ b/dumux/freeflow/navierstokesnc/properties.hh
@@ -27,8 +27,8 @@
 #ifndef DUMUX_NAVIERSTOKES_NC_PROPERTIES_HH
 #define DUMUX_NAVIERSTOKES_NC_PROPERTIES_HH
 
-#include <dumux/freeflow/staggered/properties.hh>
-#include <dumux/freeflow/staggeredni/properties.hh>
+#include <dumux/freeflow/navierstokes/model.hh>
+#include <dumux/freeflow/nonisothermal/model.hh>
 #include <dumux/discretization/fickslaw.hh>
 
 #include "volumevariables.hh"
@@ -62,12 +62,16 @@ NEW_TYPE_TAG(NavierStokesNCNI, INHERITS_FROM(NavierStokesNC, NavierStokesNonIsot
 ///////////////////////////////////////////////////////////////////////////
 // default property values for the isothermal single phase model
 ///////////////////////////////////////////////////////////////////////////
-SET_PROP(NavierStokesNC, NumEqCellCenter)
+
+//! The number of equations
+SET_PROP(NavierStokesNC, NumEq)
 {
 private:
-    static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    static constexpr auto dim = GridView::dimension;
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
 public:
-    static constexpr int value = numComponents;
+    static constexpr int value = dim + FluidSystem::numComponents;
 };
 
 SET_INT_PROP(NavierStokesNC, ReplaceCompEqIdx, 0);
@@ -81,17 +85,19 @@ SET_INT_PROP(NavierStokesNC, ReplaceCompEqIdx, 0);
 SET_PROP(NavierStokesNC, NumComponents)
 {
 private:
-   typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
-
+   using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
 public:
    static constexpr int value = FluidSystem::numComponents;
-
 };
 
 //! the VolumeVariables property
 SET_TYPE_PROP(NavierStokesNC, VolumeVariables, NavierStokesNCVolumeVariables<TypeTag>);
 SET_TYPE_PROP(NavierStokesNC, Indices, NavierStokesNCIndices<TypeTag>);
 
+SET_TYPE_PROP(NavierStokesNC, FluxVariables, NavierStokesNCFluxVariables<TypeTag>);
+
+SET_TYPE_PROP(NavierStokesNC, LocalResidual, NavierStokesNCResidual<TypeTag>);
+
 /*!
  * \brief The fluid state which is used by the volume variables to
  *        store the thermodynamic state. This should be chosen
@@ -121,11 +127,20 @@ SET_INT_PROP(NavierStokesNC, PhaseIdx, 0); //!< Defines the phaseIdx
 SET_TYPE_PROP(NavierStokesNC, VtkOutputFields, NavierStokesNCVtkOutputFields<TypeTag>); //! the vtk output fields
 
 // non-isothermal properties
-SET_INT_PROP(NavierStokesNCNI, IsothermalNumEqCellCenter, GET_PROP_VALUE(TypeTag, NumComponents)); //!< set the number of equations on the cell centers
-SET_INT_PROP(NavierStokesNCNI, IsothermalNumEqFace, 1); //!< set the number of equations on the faces
 SET_TYPE_PROP(NavierStokesNCNI, IsothermalIndices, NavierStokesNCIndices<TypeTag>); //! the isothermal indices
 SET_TYPE_PROP(NavierStokesNCNI, IsothermalVtkOutputFields, NavierStokesNCVtkOutputFields<TypeTag>); //! the isothermal vtk output fields
 
+//! The number of equations
+SET_PROP(NavierStokesNCNI, IsothermalNumEq)
+{
+private:
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    static constexpr auto dim = GridView::dimension;
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+public:
+    static constexpr int value = dim + FluidSystem::numComponents;
+};
+
 
 // \}
 }
diff --git a/dumux/freeflow/navierstokesnc/staggered/fluxvariables.hh b/dumux/freeflow/navierstokesnc/staggered/fluxvariables.hh
new file mode 100644
index 0000000000..bd3eccdb5f
--- /dev/null
+++ b/dumux/freeflow/navierstokesnc/staggered/fluxvariables.hh
@@ -0,0 +1,129 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for the flux variables
+ */
+#ifndef DUMUX_NAVIERSTOKES_NC_STAGGERED_FLUXVARIABLES_HH
+#define DUMUX_NAVIERSTOKES_NC_STAGGERED_FLUXVARIABLES_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/fluxvariablesbase.hh>
+#include <dumux/discretization/methods.hh>
+
+namespace Dumux
+{
+
+// forward declaration
+template<class TypeTag, DiscretizationMethods Method>
+class NavierStokesNCFluxVariablesImpl;
+
+/*!
+ * \ingroup Discretization
+ * \brief Base class for the flux variables
+ *        Actual flux variables inherit from this class
+ */
+// specialization for immiscible, isothermal flow
+template<class TypeTag>
+class NavierStokesNCFluxVariablesImpl<TypeTag, DiscretizationMethods::Staggered>
+: public NavierStokesFluxVariables<TypeTag>
+{
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+    using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+    using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+    using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+
+    using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
+
+    static constexpr auto numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+
+    static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+
+    //! The index of the component balance equation that gets replaced with the total mass balance
+    static const int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
+    static const int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+
+    using ParentType = NavierStokesFluxVariables<TypeTag>;
+
+    enum {
+
+        pressureIdx = Indices::pressureIdx,
+        velocityIdx = Indices::velocityIdx,
+
+        massBalanceIdx = Indices::massBalanceIdx,
+        momentumBalanceIdx = Indices::momentumBalanceIdx,
+        conti0EqIdx = Indices::conti0EqIdx,
+    };
+
+public:
+
+    CellCenterPrimaryVariables computeFluxForCellCenter(const Problem& problem,
+                                                        const Element &element,
+                                                        const FVElementGeometry& fvGeometry,
+                                                        const ElementVolumeVariables& elemVolVars,
+                                                        const ElementFaceVariables& elemFaceVars,
+                                                        const SubControlVolumeFace &scvf,
+                                                        const FluxVariablesCache& fluxVarsCache)
+    {
+        CellCenterPrimaryVariables flux(0.0);
+
+        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+        {
+            // get equation index
+            const auto eqIdx = conti0EqIdx + compIdx;
+
+            bool isOutflow = false;
+            if(scvf.boundary())
+            {
+                const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
+                    if(bcTypes.isOutflow(eqIdx))
+                        isOutflow = true;
+            }
+
+            auto upwindTerm = [compIdx](const auto& volVars)
+            {
+                const auto density = useMoles ? volVars.molarDensity() : volVars.density();
+                const auto fraction =  useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx);
+                return density * fraction;
+            };
+
+            flux[eqIdx] = ParentType::advectiveFluxForCellCenter(elemVolVars, elemFaceVars, scvf, upwindTerm, isOutflow);
+        }
+
+        // in case one balance is substituted by the total mass balance
+        if (replaceCompEqIdx < numComponents)
+        {
+            flux[replaceCompEqIdx] = std::accumulate(flux.begin(), flux.end(), 0.0);
+        }
+
+        flux += MolecularDiffusionType::diffusiveFluxForCellCenter(problem, fvGeometry, elemVolVars, scvf);
+        return flux;
+    }
+};
+
+} // end namespace
+
+#endif // DUMUX_NAVIERSTOKES_NC_STAGGERED_FLUXVARIABLES_HH
diff --git a/dumux/freeflow/navierstokesnc/staggered/localresidual.hh b/dumux/freeflow/navierstokesnc/staggered/localresidual.hh
new file mode 100644
index 0000000000..d1d585d5f2
--- /dev/null
+++ b/dumux/freeflow/navierstokesnc/staggered/localresidual.hh
@@ -0,0 +1,161 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Element-wise calculation of the residual NavierStokesNC models using the staggered discretization
+ */
+#ifndef DUMUX_STAGGERED_NAVIERSTOKES_NC_LOCAL_RESIDUAL_HH
+#define DUMUX_STAGGERED_NAVIERSTOKES_NC_LOCAL_RESIDUAL_HH
+
+#include <dune/common/hybridutilities.hh>
+#include <dumux/common/properties.hh>
+
+
+namespace Dumux
+{
+
+/*!
+ *
+ * \todo Please doc me more!
+ */
+
+
+// // forward declaration
+template<class TypeTag,  DiscretizationMethods Method>
+class NavierStokesNCResidualImpl;
+
+// specialization for miscible, isothermal flow
+template<class TypeTag>
+class NavierStokesNCResidualImpl<TypeTag, DiscretizationMethods::Staggered> : public NavierStokesResidual<TypeTag>
+{
+    using ParentType = NavierStokesResidual<TypeTag>;
+    friend class StaggeredLocalResidual<TypeTag>;
+    friend ParentType;
+
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+    using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+
+    using CellCenterResidual = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+
+    using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+    typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+
+    enum {
+        conti0EqIdx = Indices::conti0EqIdx,
+        phaseIdx = Indices::phaseIdx,
+
+        // The index of the component balance equation
+        // that gets replaced with the total mass balance
+        replaceCompEqIdx = Indices::replaceCompEqIdx
+    };
+
+    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+
+    static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+    static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+
+public:
+    using ParentType::ParentType;
+
+     /*!
+     * \brief Evaluate the rate of change of all conservation
+     *        quantites (e.g. phase mass) within a sub-control
+     *        volume of a finite volume element for the immiscible models.
+     * \param scv The sub control volume
+     * \param volVars The current or previous volVars
+     * \note This function should not include the source and sink terms.
+     * \note The volVars can be different to allow computing
+     *       the implicit euler time derivative here
+     */
+    CellCenterPrimaryVariables computeStorageForCellCenter(const Problem& problem,
+                                                           const SubControlVolume& scv,
+                                                           const VolumeVariables& volVars) const
+    {
+        CellCenterPrimaryVariables storage(0.0);
+
+        const Scalar density = useMoles ? volVars.molarDensity() : volVars.density();
+
+        // compute storage term of all components
+        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+        {
+            const int eqIdx = conti0EqIdx + compIdx;
+
+            const Scalar massOrMoleFraction = useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx);
+            const Scalar s =  density * massOrMoleFraction;
+
+            if (eqIdx != replaceCompEqIdx)
+                storage[eqIdx] += s;
+        }
+
+        // in case one balance is substituted by the total mass balance
+        if(replaceCompEqIdx < numComponents)
+            storage[replaceCompEqIdx] = density;
+
+        // add energy storage for non-isothermal models
+        Dune::Hybrid::ifElse(std::integral_constant<bool, GET_PROP_VALUE(TypeTag, EnableEnergyBalance) >(),
+        [&](auto IF)
+        {
+            storage[Indices::energyBalanceIdx] = volVars.density() * volVars.internalEnergy();
+        });
+
+        return storage;
+    }
+
+protected:
+
+    /*!
+     * \brief Sets a fixed Dirichlet value for a cell (such as pressure) at the boundary.
+     *        This is a provisional alternative to setting the Dirichlet value on the boundary directly.
+     *
+     * \param insideScv The sub control volume
+     * \param elemVolVars The current or previous element volVars
+     * \param bcTypes The boundary types
+     */
+    void setFixedCell_(CellCenterResidual& residual,
+                       const Problem& problem,
+                       const SubControlVolume& insideScv,
+                       const ElementVolumeVariables& elemVolVars,
+                       const BoundaryTypes& bcTypes) const
+    {
+        ParentType::setFixedCell_(residual, problem, insideScv, elemVolVars, bcTypes);
+
+        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+        {
+            // get equation index
+            const auto eqIdx = conti0EqIdx + compIdx;
+
+            // set a fixed mole fraction for cells
+            if(eqIdx != conti0EqIdx && bcTypes.isDirichletCell(eqIdx))
+            {
+                const auto& insideVolVars = elemVolVars[insideScv];
+                const Scalar massOrMoleFraction = useMoles ? insideVolVars.moleFraction(phaseIdx, compIdx) : insideVolVars.massFraction(phaseIdx, compIdx);
+                residual[eqIdx] = massOrMoleFraction - problem.dirichletAtPos(insideScv.center())[cellCenterIdx][eqIdx];
+            }
+        }
+
+    }
+};
+}
+
+#endif   // DUMUX_CC_LOCAL_RESIDUAL_HH
diff --git a/dumux/freeflow/navierstokesnc/volumevariables.hh b/dumux/freeflow/navierstokesnc/volumevariables.hh
index 463bad3ab2..6444966f68 100644
--- a/dumux/freeflow/navierstokesnc/volumevariables.hh
+++ b/dumux/freeflow/navierstokesnc/volumevariables.hh
@@ -27,7 +27,7 @@
 #include <dumux/common/properties.hh>
 #include "properties.hh"
 
-#include "../staggered/volumevariables.hh"
+#include "../navierstokes/volumevariables.hh"
 
 #include <dumux/material/fluidstates/immiscible.hh>
 
@@ -78,7 +78,6 @@ public:
                 const Element &element,
                 const SubControlVolume& scv)
     {
-        this->priVars_ = this->extractDofPriVars(elemSol, scv);
         this->extrusionFactor_ = problem.extrusionFactor(element, scv, elemSol);
 
         completeFluidState(elemSol, problem, element, scv, this->fluidState_);
diff --git a/dumux/freeflow/navierstokesnc/vtkoutputfields.hh b/dumux/freeflow/navierstokesnc/vtkoutputfields.hh
index 00babf1aa4..4ac91de7f9 100644
--- a/dumux/freeflow/navierstokesnc/vtkoutputfields.hh
+++ b/dumux/freeflow/navierstokesnc/vtkoutputfields.hh
@@ -24,7 +24,7 @@
 #define DUMUX_NAVIER_STOKES_NC_VTK_OUTPUT_FIELDS_HH
 
 #include <dumux/common/properties.hh>
-#include <dumux/freeflow/staggered/vtkoutputfields.hh>
+#include <dumux/freeflow/navierstokes/vtkoutputfields.hh>
 
 namespace Dumux
 {
diff --git a/test/freeflow/staggerednc/channeltestproblem.hh b/test/freeflow/staggerednc/channeltestproblem.hh
index c0d69eeabd..b4f9e1c849 100644
--- a/test/freeflow/staggerednc/channeltestproblem.hh
+++ b/test/freeflow/staggerednc/channeltestproblem.hh
@@ -24,12 +24,16 @@
 #ifndef DUMUX_CHANNEL_NC_TEST_PROBLEM_HH
 #define DUMUX_CHANNEL_NC_TEST_PROBLEM_HH
 
-#include <dumux/freeflow/staggered/problem.hh>
-#include <dumux/discretization/staggered/properties.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
 #include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/components/constant.hh>
 #include <dumux/material/fluidsystems/h2oair.hh>
 
-#include <dumux/freeflow/staggerednc/properties.hh>
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/freeflow/navierstokes/model.hh>
+
+#include <dumux/freeflow/navierstokesnc/properties.hh>
 
 namespace Dumux
 {
@@ -47,9 +51,9 @@ namespace Properties
 {
 
 #if !NONISOTHERMAL
-NEW_TYPE_TAG(ChannelNCTestProblem, INHERITS_FROM(StaggeredModel, NavierStokesNC));
+NEW_TYPE_TAG(ChannelNCTestProblem, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
 #else
-NEW_TYPE_TAG(ChannelNCTestProblem, INHERITS_FROM(StaggeredModel, NavierStokesNCNI));
+NEW_TYPE_TAG(ChannelNCTestProblem, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNCNI));
 #endif
 
 NEW_PROP_TAG(FluidSystem);
diff --git a/test/freeflow/staggerednc/densityflowproblem.hh b/test/freeflow/staggerednc/densityflowproblem.hh
index 31d9da1ab3..431cc6f752 100644
--- a/test/freeflow/staggerednc/densityflowproblem.hh
+++ b/test/freeflow/staggerednc/densityflowproblem.hh
@@ -24,10 +24,14 @@
 #ifndef DUMUX_DENSITY_FLOW_NC_TEST_PROBLEM_HH
 #define DUMUX_DENSITY_FLOW_NC_TEST_PROBLEM_HH
 
-#include <dumux/freeflow/staggered/problem.hh>
-#include <dumux/discretization/staggered/properties.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
 #include <dumux/material/components/simpleh2o.hh>
-#include <dumux/material/fluidsystems/h2oair.hh>
+#include <dumux/material/components/constant.hh>
+#include <dumux/material/fluidsystems/1p.hh>
+
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/freeflow/navierstokes/model.hh>
 
 #include <dumux/freeflow/staggerednc/properties.hh>
 
@@ -45,7 +49,7 @@ namespace Capabilities
 
 namespace Properties
 {
-NEW_TYPE_TAG(DensityDrivenFlowProblem, INHERITS_FROM(StaggeredModel, NavierStokesNC));
+NEW_TYPE_TAG(DensityDrivenFlowProblem, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
 
 NEW_PROP_TAG(FluidSystem);
 
-- 
GitLab