From 938c43160275aae654b29d50929cfaa769fb6406 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Tue, 13 Oct 2020 20:21:01 +0200
Subject: [PATCH] [navierstokes1pnc] Add headers

---
 .../navierstokes/mass/1pnc/advectiveflux.hh   | 104 +++++
 .../navierstokes/mass/1pnc/fluxvariables.hh   | 204 +++++++++
 .../navierstokes/mass/1pnc/indices.hh         |  41 ++
 .../navierstokes/mass/1pnc/iofields.hh        |  87 ++++
 .../navierstokes/mass/1pnc/localresidual.hh   | 131 ++++++
 .../freeflow/navierstokes/mass/1pnc/model.hh  | 411 ++++++++++++++++++
 .../navierstokes/mass/1pnc/volumevariables.hh | 288 ++++++++++++
 .../freeflow/navierstokes/mass/CMakeLists.txt |   1 +
 8 files changed, 1267 insertions(+)
 create mode 100644 dumux/freeflow/navierstokes/mass/1pnc/advectiveflux.hh
 create mode 100644 dumux/freeflow/navierstokes/mass/1pnc/fluxvariables.hh
 create mode 100644 dumux/freeflow/navierstokes/mass/1pnc/indices.hh
 create mode 100644 dumux/freeflow/navierstokes/mass/1pnc/iofields.hh
 create mode 100644 dumux/freeflow/navierstokes/mass/1pnc/localresidual.hh
 create mode 100644 dumux/freeflow/navierstokes/mass/1pnc/model.hh
 create mode 100644 dumux/freeflow/navierstokes/mass/1pnc/volumevariables.hh

diff --git a/dumux/freeflow/navierstokes/mass/1pnc/advectiveflux.hh b/dumux/freeflow/navierstokes/mass/1pnc/advectiveflux.hh
new file mode 100644
index 0000000000..e28440bb9c
--- /dev/null
+++ b/dumux/freeflow/navierstokes/mass/1pnc/advectiveflux.hh
@@ -0,0 +1,104 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+ /*!
+  * \file
+  * \ingroup NavierStokesModel
+  * \brief Helper struct defining the advective fluxes of the single-phase flow multicomponent
+  *        Navier-Stokes mass model
+  */
+#ifndef DUMUX_FREEFLOW_NAVIERSTOKES_MASS_1PNC_ADVECTIVE_FLUX_HH
+#define DUMUX_FREEFLOW_NAVIERSTOKES_MASS_1PNC_ADVECTIVE_FLUX_HH
+
+namespace Dumux {
+
+#ifndef DOXYGEN
+// forward declare
+template<int nComp, bool useM, int repCompEqIdx>
+struct NavierStokesMassOnePNCModelTraits;
+
+template<class IsothermalTraits>
+struct NavierStokesEnergyModelTraits;
+
+template<class ModelTraits>
+struct AdvectiveFlux;
+#endif
+
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief Helper struct defining the advective fluxes of the single-phase flow multicomponent
+ *        Navier-Stokes mass model
+ */
+template<int nComp, bool useM, int repCompEqIdx>
+struct AdvectiveFlux<NavierStokesMassOnePNCModelTraits<nComp, useM, repCompEqIdx>>
+{
+    template<class NumEqVector, class UpwindFunction>
+    static void addAdvectiveFlux(NumEqVector& flux,
+                                 const UpwindFunction& upwind)
+    {
+        using ModelTraits = NavierStokesMassOnePNCModelTraits<nComp, useM, repCompEqIdx>;
+        static constexpr bool useMoles = ModelTraits::useMoles();
+        static constexpr auto numComponents = ModelTraits::numFluidComponents();
+        static constexpr auto replaceCompEqIdx = ModelTraits::replaceCompEqIdx();
+        static constexpr bool useTotalMoleOrMassBalance = replaceCompEqIdx < numComponents;
+
+        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+        {
+            // get equation index
+            const auto eqIdx = ModelTraits::Indices::conti0EqIdx + compIdx;
+
+            if (eqIdx == replaceCompEqIdx)
+                continue;
+
+            auto upwindTerm = [&]()
+            {
+                if constexpr (useMoles)
+                    return [compIdx](const auto& volVars) { return volVars.molarDensity()*volVars.moleFraction(compIdx); };
+                else
+                    return [compIdx](const auto& volVars) { return volVars.density()*volVars.massFraction(compIdx); };
+            }();
+
+            flux[eqIdx] += upwind(upwindTerm);
+        }
+
+        // in case one balance is substituted by the total mole balance
+        if constexpr(useTotalMoleOrMassBalance)
+        {
+            auto upwindTerm = [&]()
+            {
+                if constexpr (useMoles)
+                    return [](const auto& volVars) { return volVars.molarDensity(); };
+                else
+                    return [](const auto& volVars) { return volVars.density(); };
+            }();
+
+            flux[replaceCompEqIdx] += upwind(upwindTerm);
+        }
+    }
+};
+
+// use the same mass flux for the non-isothermal model (heat fluxes are added separately)
+template<int nComp, bool useM, int repCompEqIdx>
+struct AdvectiveFlux<NavierStokesEnergyModelTraits<NavierStokesMassOnePNCModelTraits<nComp, useM, repCompEqIdx>>>
+: public AdvectiveFlux<NavierStokesMassOnePNCModelTraits<nComp, useM, repCompEqIdx>>
+{};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/freeflow/navierstokes/mass/1pnc/fluxvariables.hh b/dumux/freeflow/navierstokes/mass/1pnc/fluxvariables.hh
new file mode 100644
index 0000000000..550a7f87d6
--- /dev/null
+++ b/dumux/freeflow/navierstokes/mass/1pnc/fluxvariables.hh
@@ -0,0 +1,204 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup NavierStokesModel
+ * \copydoc Dumux::NavierStokesMassOnePNCFluxVariables
+ */
+#ifndef DUMUX_NAVIERSTOKES_MASS_1PNC_FLUXVARIABLES_HH
+#define DUMUX_NAVIERSTOKES_MASS_1PNC_FLUXVARIABLES_HH
+
+#include <dumux/common/typetraits/problem.hh>
+#include <dumux/flux/referencesystemformulation.hh>
+#include <dumux/flux/upwindscheme.hh>
+#include <dumux/freeflow/navierstokes/scalarfluxvariables.hh>
+
+#include "advectiveflux.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief The flux variables class for the single-phase flow,
+ *        multi-component Navier-Stokes model.
+ */
+template<class Problem,
+         class ModelTraits,
+         class FluxTs,
+         class ElementVolumeVariables,
+         class ElementFluxVariablesCache,
+         class UpwindScheme = UpwindScheme<typename ProblemTraits<Problem>::GridGeometry>>
+class NavierStokesMassOnePNCFluxVariables
+: public NavierStokesScalarConservationModelFluxVariables<Problem,
+                                                          ModelTraits,
+                                                          FluxTs,
+                                                          ElementVolumeVariables,
+                                                          ElementFluxVariablesCache,
+                                                          UpwindScheme>
+{
+    using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
+    using NumEqVector = typename VolumeVariables::PrimaryVariables;
+    using Scalar = typename VolumeVariables::PrimaryVariables::value_type;
+    using Indices = typename ModelTraits::Indices;
+
+    static constexpr bool enableMolecularDiffusion = ModelTraits::enableMolecularDiffusion();
+    static constexpr auto replaceCompEqIdx = ModelTraits::replaceCompEqIdx();
+    static constexpr bool useTotalMoleOrMassBalance = replaceCompEqIdx < ModelTraits::numFluidComponents();
+
+    using FluidSystem = typename VolumeVariables::FluidSystem;
+
+    using ParentType = NavierStokesScalarConservationModelFluxVariables<Problem,
+                                                                        ModelTraits,
+                                                                        FluxTs,
+                                                                        ElementVolumeVariables,
+                                                                        ElementFluxVariablesCache,
+                                                                        UpwindScheme>;
+
+public:
+
+    static constexpr auto numComponents = ModelTraits::numFluidComponents();
+    static constexpr bool useMoles = ModelTraits::useMoles();
+    using MolecularDiffusionType = typename FluxTs::MolecularDiffusionType;
+
+    /*!
+     * \brief Returns the diffusive fluxes computed by the respective law.
+     */
+    NumEqVector molecularDiffusionFlux(int phaseIdx = 0) const
+    {
+        NumEqVector result(0.0);
+        if constexpr (enableMolecularDiffusion)
+        {
+            const auto diffusiveFluxes = MolecularDiffusionType::flux(this->problem(),
+                                                                      this->element(),
+                                                                      this->fvGeometry(),
+                                                                      this->elemVolVars(),
+                                                                      this->scvFace(),
+                                                                      phaseIdx,
+                                                                      this->elemFluxVarsCache());
+
+            static constexpr auto referenceSystemFormulation = MolecularDiffusionType::referenceSystemFormulation();
+
+            for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+            {
+                // get equation index
+                const auto eqIdx = Indices::conti0EqIdx + compIdx;
+                if (eqIdx == replaceCompEqIdx)
+                    continue;
+
+                //check for the reference system and adapt units of the diffusive flux accordingly.
+                if constexpr (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
+                    result[eqIdx] += useMoles ? diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx)
+                                        : diffusiveFluxes[compIdx];
+                else if constexpr (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged)
+                    result[eqIdx] += useMoles ? diffusiveFluxes[compIdx]
+                                        : diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
+                else
+                    DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
+            }
+
+            // in case one balance is substituted by the total mole balance
+            if constexpr(useTotalMoleOrMassBalance)
+            {
+                for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+                {
+                    //check for the reference system and adapt units of the diffusive flux accordingly.
+                    if constexpr (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
+                        result[replaceCompEqIdx] += useMoles ? diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx) : diffusiveFluxes[compIdx];
+                    else if constexpr (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged)
+                        result[replaceCompEqIdx] += useMoles ? diffusiveFluxes[compIdx]
+                                                : diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
+                    else
+                        DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
+                }
+            }
+        }
+
+        return result;
+    }
+
+    /*!
+     * \brief Returns the flux of enthalpy in J/s carried by diffusing molecules.
+     */
+    Scalar diffusiveEnthalpyFlux(const NumEqVector& diffusiveFlux, int phaseIdx = 0) const
+    {
+        static constexpr auto referenceSystemFormulation = MolecularDiffusionType::referenceSystemFormulation();
+        Scalar flux = 0.0;
+        const auto& scvf = this->scvFace();
+        const auto& elemVolVars = this->elemVolVars();
+
+        const auto componentEnthalpy = [](const auto& volVars, int compIdx)
+        { return FluidSystem::componentEnthalpy(volVars.fluidState(), 0, compIdx); };
+
+        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+        {
+            // do a simple upwinding of enthalpy based on the direction of the diffusive flux
+            using std::signbit;
+            const bool insideIsUpstream = !signbit(diffusiveFlux[compIdx]);
+            const auto& upstreamVolVars = insideIsUpstream ? elemVolVars[scvf.insideScvIdx()] : elemVolVars[scvf.outsideScvIdx()];
+
+            if constexpr (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
+                flux += diffusiveFlux[compIdx] * componentEnthalpy(upstreamVolVars, compIdx);
+            else
+                flux += diffusiveFlux[compIdx] * componentEnthalpy(upstreamVolVars, compIdx)* elemVolVars[scvf.insideScvIdx()].molarMass(compIdx);
+        }
+
+        return flux;
+    }
+
+    /*!
+     * \brief Returns the advective mass flux in kg/s
+     *        or the advective mole flux in mole/s.
+     */
+    NumEqVector advectiveFlux(int phaseIdx = 0) const
+    {
+        NumEqVector result(0.0);
+        // g++ requires to capture 'this' by value
+        const auto upwinding = [this](const auto& term) { return this->getAdvectiveFlux(term); };
+        AdvectiveFlux<ModelTraits>::addAdvectiveFlux(result, upwinding);
+        return result;
+    }
+
+    /*!
+     * \brief Returns all fluxes for the single-phase flow, multi-component
+     *        Navier-Stokes model: the advective mass flux in kg/s
+     *        or the advective mole flux in mole/s and the energy flux
+     *        in J/s (for nonisothermal models).
+     */
+    NumEqVector flux(int phaseIdx = 0) const
+    {
+        const auto diffusiveFlux = molecularDiffusionFlux(phaseIdx);
+        NumEqVector flux = diffusiveFlux;
+        // g++ requires to capture 'this' by value
+        const auto upwinding = [this](const auto& term) { return this->getAdvectiveFlux(term); };
+        AdvectiveFlux<ModelTraits>::addAdvectiveFlux(flux, upwinding);
+
+        if constexpr (ModelTraits::enableEnergyBalance())
+        {
+            ParentType::addHeatFlux(flux);
+            flux[ModelTraits::Indices::energyEqIdx] += diffusiveEnthalpyFlux(diffusiveFlux);
+        }
+
+        return flux;
+    }
+
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/freeflow/navierstokes/mass/1pnc/indices.hh b/dumux/freeflow/navierstokes/mass/1pnc/indices.hh
new file mode 100644
index 0000000000..b059cd6860
--- /dev/null
+++ b/dumux/freeflow/navierstokes/mass/1pnc/indices.hh
@@ -0,0 +1,41 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup NavierStokesModel
+ * \copydoc Dumux::NavierStokesIndices
+ */
+#ifndef DUMUX_NAVIERSTOKES_MASS_1P_INDICES_HH
+#define DUMUX_NAVIERSTOKES_MASS_1P_INDICES_HH
+
+namespace Dumux {
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief The common indices for the isothermal Navier-Stokes mass conservation model.
+ */
+struct NavierStokesMassOnePIndices
+{
+    static constexpr int conti0EqIdx = 0; //!< Index of the first (total for pure-fluid systems) mass balance equation
+    static constexpr int pressureIdx = conti0EqIdx; //!< Index of the pressure
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/freeflow/navierstokes/mass/1pnc/iofields.hh b/dumux/freeflow/navierstokes/mass/1pnc/iofields.hh
new file mode 100644
index 0000000000..6e89bf6912
--- /dev/null
+++ b/dumux/freeflow/navierstokes/mass/1pnc/iofields.hh
@@ -0,0 +1,87 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup FreeflowNCModel
+ * \copydoc Dumux::FreeflowNCIOFields
+ */
+#ifndef DUMUX_NAVIERSTOKES_MASS_1PNC_IO_FIELDS_HH
+#define DUMUX_NAVIERSTOKES_MASS_1PNC_IO_FIELDS_HH
+
+#include <dumux/io/name.hh>
+#include <dumux/freeflow/navierstokes/iofields.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup FreeflowNCModel
+ * \brief Adds I/O fields specific to the FreeflowNC model
+ */
+template<class BaseOutputFields, bool turbulenceModel = false>
+struct NavierStokesMassOnePNCIOFields
+{
+    //! Initialize the FreeflowNC specific output fields.
+    template <class OutputModule>
+    static void initOutputModule(OutputModule& out)
+    {
+        BaseOutputFields::initOutputModule(out);
+
+        // TODO only output molar density if we use mole fractions
+        out.addVolumeVariable([](const auto& v){ return v.molarDensity(); }, IOName::molarDensity());
+
+        using FluidSystem = typename OutputModule::VolumeVariables::FluidSystem;
+        for (int j = 0; j < FluidSystem::numComponents; ++j)
+        {
+            out.addVolumeVariable([j](const auto& v){ return v.massFraction(j); }, IOName::massFraction<FluidSystem>(0, j));
+            out.addVolumeVariable([j](const auto& v){ return v.moleFraction(j); }, IOName::moleFraction<FluidSystem>(0, j));
+
+            if (j != FluidSystem::getMainComponent(0))
+            {
+                out.addVolumeVariable([j](const auto& v){ return v.diffusionCoefficient(0,0, j); }, "D^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(0));
+
+                // the eddy diffusivity is recalculated for an arbitrary component which is not the phase component
+                if (turbulenceModel)
+                    out.addVolumeVariable([j](const auto& v){ return getEffectiveDiffusionCoefficient_(v, 0, j) - v.diffusionCoefficient(0,0, j); }, "D_t^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(0));
+            }
+        }
+    }
+
+    //! return the names of the primary variables
+    template <class ModelTraits, class FluidSystem>
+    static std::string primaryVariableName(int pvIdx = 0, int state = 0)
+    {
+        // priVars: v_0, ..., v_dim-1, p, x_0, ..., x_numComp-1, otherPv ..., T
+        if (pvIdx > ModelTraits::dim() && pvIdx < ModelTraits::dim() + ModelTraits::numFluidComponents())
+            return ModelTraits::useMoles() ? IOName::moleFraction<FluidSystem>(0, pvIdx - ModelTraits::dim())
+                                           : IOName::massFraction<FluidSystem>(0, pvIdx - ModelTraits::dim());
+        else
+            return BaseOutputFields::template primaryVariableName<ModelTraits, FluidSystem>(pvIdx, state);
+
+    }
+
+    template<class VolumeVariables>
+    static double getEffectiveDiffusionCoefficient_(const VolumeVariables& volVars, const int phaseIdx, const int compIdx)
+    {
+        return volVars.effectiveDiffusionCoefficient(phaseIdx, VolumeVariables::FluidSystem::getMainComponent(phaseIdx), compIdx);
+    }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/freeflow/navierstokes/mass/1pnc/localresidual.hh b/dumux/freeflow/navierstokes/mass/1pnc/localresidual.hh
new file mode 100644
index 0000000000..5878df9c16
--- /dev/null
+++ b/dumux/freeflow/navierstokes/mass/1pnc/localresidual.hh
@@ -0,0 +1,131 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup NavierStokesModel
+ * \copydoc Dumux::NavierStokesMassOnePNCLocalResidual
+ */
+#ifndef DUMUX_NAVIERSTOKES_MASS_1PNC_LOCAL_RESIDUAL_HH
+#define DUMUX_NAVIERSTOKES_MASS_1PNC_LOCAL_RESIDUAL_HH
+
+#include <dumux/common/numeqvector.hh>
+#include <dumux/common/properties.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief Element-wise calculation of the Navier-Stokes residual for multicomponent single-phase flow.
+ */
+template<class TypeTag>
+class NavierStokesMassOnePNCLocalResidual : public GetPropType<TypeTag, Properties::BaseLocalResidual>
+{
+    using ParentType = GetPropType<TypeTag, Properties::BaseLocalResidual>;
+    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+
+    using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
+    using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
+    using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
+
+    using GridFluxVariablesCache = typename GridVariables::GridFluxVariablesCache;
+    using ElementFluxVariablesCache = typename GridFluxVariablesCache::LocalView;
+
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+    using FVElementGeometry = typename GridGeometry::LocalView;
+    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+    using GridView = typename GridGeometry::GridView;
+    using Element = typename GridView::template Codim<0>::Entity;
+    using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>;
+    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+    using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
+
+    static_assert(GridGeometry::discMethod == DiscretizationMethod::cctpfa);
+    static constexpr bool useMoles = ModelTraits::useMoles();
+    static constexpr auto numComponents = ModelTraits::numFluidComponents();
+
+    static constexpr int replaceCompEqIdx = ModelTraits::replaceCompEqIdx();
+    static constexpr bool useTotalMoleOrMassBalance = replaceCompEqIdx < numComponents;
+
+public:
+    //! Use the parent type's constructor
+    using ParentType::ParentType;
+
+    /*!
+     * \brief Calculate the storage term of the equation
+     */
+    NumEqVector computeStorage(const Problem& problem,
+                               const SubControlVolume& scv,
+                               const VolumeVariables& volVars) const
+    {
+        NumEqVector 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 = compIdx;
+
+            const Scalar massOrMoleFraction = useMoles ? volVars.moleFraction(compIdx) : volVars.massFraction(compIdx);
+            const Scalar s =  density * massOrMoleFraction;
+
+            if (eqIdx != ModelTraits::replaceCompEqIdx())
+                storage[eqIdx] += s;
+        }
+
+        // in case one balance is substituted by the total mass balance
+        if constexpr (useTotalMoleOrMassBalance)
+            storage[ModelTraits::replaceCompEqIdx()] = density;
+
+        // consider energy storage for non-isothermal models
+        if constexpr (ModelTraits::enableEnergyBalance())
+            storage[ModelTraits::Indices::energyEqIdx] = volVars.density() * volVars.internalEnergy();
+
+        return storage;
+    }
+
+    /*!
+     * \brief Evaluatex the mass or mole flux over a face of a sub control volume.
+     *
+     * \param problem The problem
+     * \param element The element
+     * \param fvGeometry The finite volume geometry context
+     * \param elemVolVars The volume variables for all flux stencil elements
+     * \param scvf The sub control volume face to compute the flux on
+     * \param elemFluxVarsCache The cache related to flux computation
+     */
+    NumEqVector computeFlux(const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const SubControlVolumeFace& scvf,
+                            const ElementFluxVariablesCache& elemFluxVarsCache) const
+    {
+        FluxVariables fluxVars;
+        fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+        return fluxVars.flux(0);
+    }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/freeflow/navierstokes/mass/1pnc/model.hh b/dumux/freeflow/navierstokes/mass/1pnc/model.hh
new file mode 100644
index 0000000000..73fb1d51a9
--- /dev/null
+++ b/dumux/freeflow/navierstokes/mass/1pnc/model.hh
@@ -0,0 +1,411 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup NavierStokesModel
+ *
+ * \brief A single-phase, isothermal Navier-Stokes model
+ *
+ * This model implements a single-phase, isothermal Navier-Stokes model, solving the <B> momentum balance equation </B>
+ * \f[
+ \frac{\partial (\varrho \textbf{v})}{\partial t} + \nabla \cdot (\varrho \textbf{v} \textbf{v}^{\text{T}}) = \nabla \cdot (\mu (\nabla \textbf{v} + \nabla \textbf{v}^{\text{T}}))
+   - \nabla p + \varrho \textbf{g} - \textbf{f}
+ * \f]
+ * By setting the runtime parameter <code>Problem.EnableInertiaTerms</code> to <code>false</code> the Stokes
+ * equation can be solved. In this case the term
+ * \f[
+ *    \nabla \cdot (\varrho \textbf{v} \textbf{v}^{\text{T}})
+ * \f]
+ * is neglected.
+ *
+ * The <B> mass balance equation </B>
+ * \f[
+       \frac{\partial \varrho}{\partial t} + \nabla \cdot (\varrho \textbf{v}) - q = 0
+ * \f]
+ *
+ * closes the system.
+ *
+ *
+ * So far, only the staggered grid spatial discretization (for structured grids) is available.
+ */
+
+#ifndef DUMUX_NAVIERSTOKES_1PNC_MODEL_HH
+#define DUMUX_NAVIERSTOKES_1PNC_MODEL_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/properties/model.hh>
+
+#include <dumux/freeflow/spatialparams.hh>
+#include <dumux/freeflow/turbulencemodel.hh>
+#include "dumux/freeflow/navierstokes/iofields.hh"
+#include <dumux/freeflow/navierstokes/energy/model.hh>
+#include <dumux/freeflow/navierstokes/scalarfluxvariablescachefiller.hh>
+#include <dumux/material/fluidstates/compositional.hh>
+
+#include <dumux/flux/fickslaw.hh>
+#include <dumux/flux/fourierslaw.hh>
+
+#include "localresidual.hh"
+#include "volumevariables.hh"
+#include "fluxvariables.hh"
+#include "indices.hh"
+#include "./../../iofields.hh"
+#include "iofields.hh"
+
+#include <dumux/flux/cctpfa/fourierslaw.hh>
+
+#include <dumux/discretization/method.hh>
+#include <dumux/freeflow/navierstokes/energy/model.hh>
+#include <dumux/flux/fickslaw.hh>
+#include <dumux/freeflow/navierstokes/scalarfluxvariablescachefiller.hh>
+
+
+namespace Dumux {
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief Traits for the Navier-Stokes model
+ *
+ * \tparam dimension The dimension of the problem
+ */
+template<int nComp, bool useM, int repCompEqIdx = nComp>
+struct NavierStokesMassOnePNCModelTraits
+{
+    //! There are as many momentum balance equations as dimensions
+    //! and one mass balance equation.
+    static constexpr int numEq() { return nComp; }
+
+    //! The number of phases is 1
+    static constexpr int numFluidPhases() { return 1; }
+
+    //! The number of components is 1
+    static constexpr int numFluidComponents() { return nComp; }
+
+    //! Use moles or not
+    static constexpr bool useMoles() { return useM; }
+
+    // Enable advection
+    static constexpr bool enableAdvection() { return true; }
+
+    //! The model is isothermal
+    static constexpr bool enableEnergyBalance() { return false; }
+
+    //! The one-phase model has no molecular diffusion
+    static constexpr bool enableMolecularDiffusion() { return true; }
+
+    //! Index of of a component balance eq. to be replaced by a total mass/mole balance
+    static constexpr int replaceCompEqIdx() { return repCompEqIdx; }
+
+    //! The model does not include a turbulence model
+    static constexpr bool usesTurbulenceModel() { return false; }
+
+    //! return the type of turbulence model used
+    static constexpr auto turbulenceModel()
+    { return TurbulenceModel::none; }
+
+    //! the indices
+    using Indices = NavierStokesMassOnePIndices;
+};
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief Traits class for the volume variables of the Navier-Stokes model.
+ *
+ * \tparam PV The type used for primary variables
+ * \tparam FSY The fluid system type
+ * \tparam FST The fluid state type
+ * \tparam MT The model traits
+ */
+template<class PV,
+         class FSY,
+         class FST,
+         class MT>
+struct NavierStokesMassOnePNCVolumeVariablesTraits
+{
+    using PrimaryVariables = PV;
+    using FluidSystem = FSY;
+    using FluidState = FST;
+    using ModelTraits = MT;
+};
+
+// \{
+///////////////////////////////////////////////////////////////////////////
+// properties for the single-phase Navier-Stokes model
+///////////////////////////////////////////////////////////////////////////
+namespace Properties {
+
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+// Create new type tags
+namespace TTag {
+//! The type tag for the single-phase, isothermal Navier-Stokes model
+struct NavierStokesMassOnePNC{ using InheritsFrom = std::tuple<ModelProperties>; };
+struct NavierStokesMassOnePNCNI{ using InheritsFrom = std::tuple<NavierStokesMassOnePNC>; };
+
+}
+
+//! The base model traits. Per default, we use the number of components of the fluid system.
+template<class TypeTag>
+struct BaseModelTraits<TypeTag, TTag::NavierStokesMassOnePNC>
+{
+private:
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+public:
+    using type = NavierStokesMassOnePNCModelTraits<FluidSystem::numComponents, getPropValue<TypeTag, Properties::UseMoles>(), getPropValue<TypeTag, Properties::ReplaceCompEqIdx>()>;
+};
+
+
+//!< states some specifics of the Navier-Stokes model
+template<class TypeTag>
+struct ModelTraits<TypeTag, TTag::NavierStokesMassOnePNC>
+{ using type = GetPropType<TypeTag, Properties::BaseModelTraits>; };
+
+/*!
+ * \brief The fluid state which is used by the volume variables to
+ *        store the thermodynamic state. This should be chosen
+ *        appropriately for the model ((non-)isothermal, equilibrium, ...).
+ *        This can be done in the problem.
+ */
+ template<class TypeTag>
+ struct FluidState<TypeTag, TTag::NavierStokesMassOnePNC>
+ {
+ private:
+     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+     using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+ public:
+     using type = CompositionalFluidState<Scalar, FluidSystem>;
+ };
+
+//! The local residual
+template<class TypeTag>
+struct LocalResidual<TypeTag, TTag::NavierStokesMassOnePNC>
+{ using type = NavierStokesMassOnePNCLocalResidual<TypeTag>; };
+
+//! Set the volume variables property
+template<class TypeTag>
+struct VolumeVariables<TypeTag, TTag::NavierStokesMassOnePNC>
+{
+private:
+    using PV = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using FSY = GetPropType<TypeTag, Properties::FluidSystem>;
+    using FST = GetPropType<TypeTag, Properties::FluidState>;
+    using MT = GetPropType<TypeTag, Properties::ModelTraits>;
+
+    static_assert(FSY::numPhases == MT::numFluidPhases(), "Number of phases mismatch between model and fluid system");
+    static_assert(FST::numPhases == MT::numFluidPhases(), "Number of phases mismatch between model and fluid state");
+    // static_assert(!FSY::isMiscible(), "The Navier-Stokes model only works with immiscible fluid systems.");
+
+    using BaseTraits = NavierStokesMassOnePNCVolumeVariablesTraits<PV, FSY, FST, MT>;
+    using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
+    using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
+
+    template<class BaseTraits, class DT, class EDM>
+    struct NCTraits : public BaseTraits
+    {
+        using DiffusionType = DT;
+        using EffectiveDiffusivityModel = EDM;
+    };
+
+public:
+    using type = NavierStokesMassOnePNCVolumeVariables<NCTraits<BaseTraits, DT, EDM>>;
+};
+
+//! By default, we use fick's law for the diffusive fluxes
+template<class TypeTag>
+struct MolecularDiffusionType<TypeTag, TTag::NavierStokesMassOnePNC> { using type = FicksLaw<TypeTag>; };
+
+//! Use the model after Millington (1961) for the effective diffusivity
+template<class TypeTag>
+struct EffectiveDiffusivityModel<TypeTag, TTag::NavierStokesMassOnePNC>
+{
+    struct type
+    {
+        template<class VolumeVariables>
+        static auto effectiveDiffusionCoefficient(const VolumeVariables& volVars,
+                                                    const int phaseIdx,
+                                                    const int compIdxI,
+                                                    const int compIdxJ)
+        {
+            return volVars.diffusionCoefficient(phaseIdx, compIdxI, compIdxJ);
+        }
+    };
+};
+
+//! The flux variables
+template<class TypeTag>
+struct FluxVariables<TypeTag, TTag::NavierStokesMassOnePNC>
+{
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+
+    struct FluxTypes
+    {
+        using MolecularDiffusionType = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
+    };
+
+    using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+    using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
+
+    using type = NavierStokesMassOnePNCFluxVariables<Problem, ModelTraits, FluxTypes, ElementVolumeVariables, ElementFluxVariablesCache>;
+};
+
+template<class TypeTag>
+struct SolutionDependentMolecularDiffusion<TypeTag, TTag::NavierStokesMassOnePNC> { static constexpr bool value = true; };
+
+
+template<class TypeTag>
+struct FluxVariablesCache<TypeTag, TTag::NavierStokesMassOnePNC>
+{
+    struct type : public GetPropType<TypeTag, Properties::MolecularDiffusionType>::Cache
+    {};
+};
+
+template<class TypeTag>
+struct FluxVariablesCacheFiller<TypeTag, TTag::NavierStokesMassOnePNC>
+{
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+    static constexpr bool diffusionIsSolDependent = getPropValue<TypeTag, Properties::SolutionDependentMolecularDiffusion>();
+    static constexpr bool heatConductionIsSolDependent = false; // no heat conduction
+
+    using type = FreeFlowScalarFluxVariablesCacheFiller<Problem, ModelTraits, diffusionIsSolDependent, heatConductionIsSolDependent>;
+};
+
+// ! The specific I/O fields
+template<class TypeTag>
+struct IOFields<TypeTag, TTag::NavierStokesMassOnePNC> { using type = NavierStokesMassOnePNCIOFields<NavierStokesIOFields>; };
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::NavierStokesMassOnePNC>
+{
+private:
+    struct EmptyCouplingManager {};
+public:
+    using type = EmptyCouplingManager;
+};
+
+// Set the default spatial parameters
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::NavierStokesMassOnePNC>
+{
+    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using type = FreeFlowDefaultSpatialParams<GridGeometry, Scalar>;
+};
+///////////////////////////////////////////////////////////////////////////
+// Properties for the non-isothermal single phase model
+///////////////////////////////////////////////////////////////////////////
+
+//! Add temperature to the output
+template<class TypeTag>
+struct IOFields<TypeTag, TTag::NavierStokesMassOnePNCNI> { using type = NavierStokesEnergyIOFields<NavierStokesMassOnePNCIOFields<NavierStokesIOFields>>; };
+
+//! The model traits of the non-isothermal model
+template<class TypeTag>
+struct ModelTraits<TypeTag, TTag::NavierStokesMassOnePNCNI> { using type = NavierStokesEnergyModelTraits<GetPropType<TypeTag, Properties::BaseModelTraits>>; };
+
+//! Set the volume variables property
+template<class TypeTag>
+struct VolumeVariables<TypeTag, TTag::NavierStokesMassOnePNCNI>
+{
+private:
+    using PV = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using FSY = GetPropType<TypeTag, Properties::FluidSystem>;
+    using FST = GetPropType<TypeTag, Properties::FluidState>;
+    using MT = GetPropType<TypeTag, Properties::ModelTraits>;
+
+    static_assert(FSY::numPhases == MT::numFluidPhases(), "Number of phases mismatch between model and fluid system");
+    static_assert(FST::numPhases == MT::numFluidPhases(), "Number of phases mismatch between model and fluid state");
+    static_assert(!FSY::isMiscible(), "The Navier-Stokes model only works with immiscible fluid systems.");
+
+    using BaseTraits = NavierStokesMassOnePNCVolumeVariablesTraits<PV, FSY, FST, MT>;
+
+    using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
+    using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
+    using ETCM = GetPropType<TypeTag, Properties::ThermalConductivityModel>;
+    using HCT = GetPropType<TypeTag, Properties::HeatConductionType>;
+
+    struct NCNITraits : public BaseTraits
+    {
+        using DiffusionType = DT;
+        using EffectiveDiffusivityModel = EDM;
+        using EffectiveThermalConductivityModel = ETCM;
+        using HeatConductionType = HCT;
+    };
+
+public:
+    using type = NavierStokesMassOnePNCVolumeVariables<NCNITraits>;
+};
+
+//! Use the average for effective conductivities
+template<class TypeTag>
+struct ThermalConductivityModel<TypeTag, TTag::NavierStokesMassOnePNCNI>
+{
+    struct type
+    {
+        template<class VolVars>
+        static auto effectiveThermalConductivity(const VolVars& volVars)
+        {
+            return volVars.fluidThermalConductivity();
+        }
+    };
+};
+
+template<class TypeTag>
+struct HeatConductionType<TypeTag, TTag::NavierStokesMassOnePNCNI>
+{ using type = FouriersLawImplementation<TypeTag, DiscretizationMethod::cctpfa>; };
+
+//! The flux variables
+template<class TypeTag>
+struct FluxVariables<TypeTag, TTag::NavierStokesMassOnePNCNI>
+{
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+
+    struct FluxTypes
+    {
+        using MolecularDiffusionType = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
+        using HeatConductionType = GetPropType<TypeTag, Properties::HeatConductionType>;
+    };
+
+    using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+    using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
+
+    using type = NavierStokesMassOnePNCFluxVariables<Problem, ModelTraits, FluxTypes, ElementVolumeVariables, ElementFluxVariablesCache>;
+};
+
+template<class TypeTag>
+struct FluxVariablesCache<TypeTag, TTag::NavierStokesMassOnePNCNI>
+{
+    struct type
+    : public GetPropType<TypeTag, Properties::MolecularDiffusionType>::Cache
+    , public GetPropType<TypeTag, Properties::HeatConductionType>::Cache
+    {};
+};
+
+template<class TypeTag>
+struct FluxVariablesCacheFiller<TypeTag, TTag::NavierStokesMassOnePNCNI>
+{
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+    static constexpr bool diffusionIsSolDependent = getPropValue<TypeTag, Properties::SolutionDependentMolecularDiffusion>();
+    static constexpr bool heatConductionIsSolDependent = getPropValue<TypeTag, Properties::SolutionDependentHeatConduction>();
+
+    using type = FreeFlowScalarFluxVariablesCacheFiller<Problem, ModelTraits, diffusionIsSolDependent, heatConductionIsSolDependent>;
+};
+
+template<class TypeTag>
+struct SolutionDependentHeatConduction<TypeTag, TTag::NavierStokesMassOnePNCNI> { static constexpr bool value = true; };
+
+} // end namespace Properties
+
+} // end namespace Dumux
+
+#endif // DUMUX_NAVIERSTOKES_MODEL_HH
diff --git a/dumux/freeflow/navierstokes/mass/1pnc/volumevariables.hh b/dumux/freeflow/navierstokes/mass/1pnc/volumevariables.hh
new file mode 100644
index 0000000000..9626e980e2
--- /dev/null
+++ b/dumux/freeflow/navierstokes/mass/1pnc/volumevariables.hh
@@ -0,0 +1,288 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup NavierStokesModel
+ *
+ * \copydoc Dumux::NavierStokesVolumeVariables
+ */
+#ifndef DUMUX_NAVIERSTOKES_MASS_1PNC_VOLUME_VARIABLES_HH
+#define DUMUX_NAVIERSTOKES_MASS_1PNC_VOLUME_VARIABLES_HH
+
+#include <dumux/freeflow/navierstokes/scalarvolumevariables.hh>
+#include <dumux/freeflow/navierstokes/energy/volumevariables.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief Volume variables for the single-phase Navier-Stokes model.
+ */
+template <class Traits>
+class NavierStokesMassOnePNCVolumeVariables
+: public NavierStokesScalarConservationModelVolumeVariables<Traits>
+, public NavierStokesEnergyVolumeVariables<Traits, NavierStokesMassOnePNCVolumeVariables<Traits>>
+{
+    using ParentType = NavierStokesScalarConservationModelVolumeVariables<Traits>;
+    using EnergyVolumeVariables = NavierStokesEnergyVolumeVariables<Traits, NavierStokesMassOnePNCVolumeVariables<Traits>>;
+
+    using Scalar = typename Traits::PrimaryVariables::value_type;
+
+    static constexpr bool useMoles = Traits::ModelTraits::useMoles();
+    using DiffusionCoefficients = typename Traits::DiffusionType::DiffusionCoefficientsContainer;
+
+
+public:
+    //! export the underlying fluid system
+    using FluidSystem = typename Traits::FluidSystem;
+    //! export the fluid state type
+    using FluidState = typename Traits::FluidState;
+    //! export the diffusion type
+    using MolecularDiffusionType = typename Traits::DiffusionType;
+    //! export the indices type
+    using Indices = typename Traits::ModelTraits::Indices;
+
+    /*!
+     * \brief Update all quantities for a given control volume
+     *
+     * \param elemSol A vector containing all primary variables connected to the element
+     * \param problem The object specifying the problem which ought to
+     *                be simulated
+     * \param element An element which contains part of the control volume
+     * \param scv The sub-control volume
+     */
+    template<class ElementSolution, class Problem, class Element, class SubControlVolume>
+    void update(const ElementSolution &elemSol,
+                const Problem &problem,
+                const Element &element,
+                const SubControlVolume& scv)
+    {
+        ParentType::update(elemSol, problem, element, scv);
+
+        completeFluidState(elemSol, problem, element, scv, fluidState_);
+
+        typename FluidSystem::ParameterCache paramCache;
+        paramCache.updateAll(fluidState_);
+
+        auto getDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
+        {
+            return FluidSystem::binaryDiffusionCoefficient(this->fluidState_,
+                                                            paramCache,
+                                                            0,
+                                                            compIIdx,
+                                                            compJIdx);
+        };
+
+        diffCoefficient_.update(getDiffusionCoefficient);
+        EnergyVolumeVariables::updateEffectiveThermalConductivity();
+    }
+
+    /*!
+     * \brief Update the fluid state
+     */
+    template<class ElementSolution, class Problem, class Element, class SubControlVolume>
+    void completeFluidState(const ElementSolution& elemSol,
+                                   const Problem& problem,
+                                   const Element& element,
+                                   const SubControlVolume& scv,
+                                   FluidState& fluidState) const
+    {
+        fluidState.setTemperature(0, EnergyVolumeVariables::getTemperature(elemSol, problem, element, scv));
+        fluidState.setPressure(0, elemSol[0][Indices::pressureIdx]);
+
+        // saturation in a single phase is always 1 and thus redundant
+        // to set. But since we use the fluid state shared by the
+        // immiscible multi-phase models, so we have to set it here...
+        fluidState.setSaturation(0, 1.0);
+
+        Scalar sumFracMinorComp = 0.0;
+
+        for(int compIdx = 1; compIdx < ParentType::numFluidComponents(); ++compIdx)
+        {
+            // temporary add 1.0 to remove spurious differences in mole fractions
+            // which are below the numerical accuracy
+            Scalar moleOrMassFraction = elemSol[0][Indices::conti0EqIdx+compIdx] + 1.0;
+            moleOrMassFraction = moleOrMassFraction - 1.0;
+            if(useMoles)
+                fluidState.setMoleFraction(0, compIdx, moleOrMassFraction);
+            else
+                fluidState.setMassFraction(0, compIdx, moleOrMassFraction);
+            sumFracMinorComp += moleOrMassFraction;
+        }
+        if(useMoles)
+            fluidState.setMoleFraction(0, 0, 1.0 - sumFracMinorComp);
+        else
+            fluidState.setMassFraction(0, 0, 1.0 - sumFracMinorComp);
+
+        typename FluidSystem::ParameterCache paramCache;
+        paramCache.updateAll(fluidState);
+
+        Scalar value = FluidSystem::density(fluidState, paramCache, 0);
+        fluidState.setDensity(0, value);
+
+        value = FluidSystem::molarDensity(fluidState, paramCache, 0);
+        fluidState.setMolarDensity(0, value);
+
+        value = FluidSystem::viscosity(fluidState, paramCache, 0);
+        fluidState.setViscosity(0, value);
+
+        // compute and set the enthalpy
+        const Scalar h = EnergyVolumeVariables::enthalpy(fluidState, paramCache);
+        fluidState.setEnthalpy(0, h);
+    }
+
+    /*!
+     * \brief Return the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
+     *        the control volume.
+     */
+    Scalar pressure(int phaseIdx = 0) const
+    { return fluidState_.pressure(0); }
+
+    /*!
+     * \brief Returns the saturation.
+     */
+    Scalar saturation(int phaseIdx = 0) const
+    { return 1.0; }
+
+    /*!
+     * \brief Return the mass density \f$\mathrm{[kg/m^3]}\f$ of a given phase within the
+     *        control volume.
+     */
+    Scalar density(int phaseIdx = 0) const
+    { return fluidState_.density(0); }
+
+    /*!
+     * \brief Return temperature \f$\mathrm{[K]}\f$ inside the sub-control volume.
+     *
+     * Note that we assume thermodynamic equilibrium, i.e. the
+     * temperatures of the rock matrix and of all fluid phases are
+     * identical.
+     */
+    Scalar temperature() const
+    { return fluidState_.temperature(); }
+
+    /*!
+     * \brief Return the effective dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
+     *        control volume.
+     */
+    Scalar effectiveViscosity() const
+    { return viscosity(); }
+
+    /*!
+     * \brief Return the dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
+     *        control volume.
+     */
+    Scalar viscosity(int phaseIdx = 0) const
+    { return fluidState_.viscosity(0); }
+
+    /*!
+     * \brief Returns the mass fraction of a component in the phase \f$\mathrm{[-]}\f$
+     *
+     * \param phaseIdx the index of the phase
+     * \param compIdx the index of the component
+     */
+    Scalar massFraction(int phaseIdx, int compIdx) const
+    {
+        return fluidState_.massFraction(0, compIdx);
+    }
+
+    /*!
+     * \brief Returns the mass fraction of a component in the phase \f$\mathrm{[-]}\f$
+     *
+     * \param compIdx the index of the component
+     */
+    Scalar massFraction(int compIdx) const
+    {
+        return fluidState_.massFraction(0, compIdx);
+    }
+
+    /*!
+     * \brief Returns the mole fraction of a component in the phase \f$\mathrm{[-]}\f$
+     *
+     * \param phaseIdx the index of the phase
+     * \param compIdx the index of the component
+     */
+    Scalar moleFraction(int phaseIdx, int compIdx) const
+    {
+        return fluidState_.moleFraction(0, compIdx);
+    }
+
+    /*!
+     * \brief Returns the mole fraction of a component in the phase \f$\mathrm{[-]}\f$
+     *
+     * \param compIdx the index of the component
+     */
+    Scalar moleFraction(int compIdx) const
+    {
+        return fluidState_.moleFraction(0, compIdx);
+    }
+
+    /*!
+     * \brief Returns the mass density of a given phase \f$\mathrm{[kg/m^3]}\f$
+     */
+    Scalar molarDensity(int phaseIdx = 0) const
+    {
+        return fluidState_.molarDensity(0);
+    }
+
+    /*!
+     * \brief Returns the molar mass of a given component.
+     *
+     * \param compIdx the index of the component
+     */
+    Scalar molarMass(int compIdx) const
+    {
+        return FluidSystem::molarMass(compIdx);
+    }
+
+    /*!
+     * \brief Returns the average molar mass \f$\mathrm{[kg/mol]}\f$ of the fluid phase.
+     *
+     * \param phaseIdx The phase index
+     */
+    Scalar averageMolarMass(const int phaseIdx = 0) const
+    { return fluidState_.averageMolarMass(phaseIdx); }
+
+    /*!
+     * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
+     */
+    Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
+    { return diffCoefficient_(0, compIIdx, compJIdx); }
+
+    /*!
+     * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
+     */
+    Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
+    { return diffusionCoefficient(0, compIIdx, compJIdx); }
+
+    /*!
+     * \brief Return the fluid state of the control volume.
+     */
+    const FluidState& fluidState() const
+    { return fluidState_; }
+
+protected:
+    FluidState fluidState_;
+    // Binary diffusion coefficient
+    DiffusionCoefficients diffCoefficient_;
+};
+
+} // end namespace Dumux
+
+#endif // DUMUX_NAVIERSTOKES_MOMENTUM_VOLUME_VARIABLES_HH
diff --git a/dumux/freeflow/navierstokes/mass/CMakeLists.txt b/dumux/freeflow/navierstokes/mass/CMakeLists.txt
index 2ec22b780f..13fa6533a8 100644
--- a/dumux/freeflow/navierstokes/mass/CMakeLists.txt
+++ b/dumux/freeflow/navierstokes/mass/CMakeLists.txt
@@ -2,6 +2,7 @@
 # SPDX-License-Identifier: GPL-3.0-or-later
 
 add_subdirectory(1p)
+add_subdirectory(1pnc)
 
 file(GLOB DUMUX_FREEFLOW_NAVIERSTOKES_MASS_HEADERS *.hh *.inc)
 install(FILES ${DUMUX_FREEFLOW_NAVIERSTOKES_MASS_HEADERS}
-- 
GitLab