diff --git a/dumux/discretization/fickslaw.hh b/dumux/discretization/fickslaw.hh
index 9d0bbc0a840c615246e2ebf9186239812aead245..3058e6524105edca9b254904fb30dcdf757d89e5 100644
--- a/dumux/discretization/fickslaw.hh
+++ b/dumux/discretization/fickslaw.hh
@@ -45,5 +45,6 @@ using FicksLaw = FicksLawImplementation<TypeTag, GET_PROP_VALUE(TypeTag, Discret
#include <dumux/discretization/cellcentered/tpfa/fickslaw.hh>
#include <dumux/discretization/cellcentered/mpfa/fickslaw.hh>
#include <dumux/discretization/box/fickslaw.hh>
+#include <dumux/discretization/staggered/freeflow/fickslaw.hh>
#endif
diff --git a/dumux/discretization/staggered/freeflow/fickslaw.hh b/dumux/discretization/staggered/freeflow/fickslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..15c557c349423ef7374872e46334562ab3fe595e
--- /dev/null
+++ b/dumux/discretization/staggered/freeflow/fickslaw.hh
@@ -0,0 +1,195 @@
+// -*- 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 This file contains the data which is required to calculate
+ * diffusive mass fluxes due to molecular diffusion with Fick's law.
+ */
+#ifndef DUMUX_DISCRETIZATION_STAGGERED_FICKS_LAW_HH
+#define DUMUX_DISCRETIZATION_STAGGERED_FICKS_LAW_HH
+
+#include <dune/common/float_cmp.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/common/parameters.hh>
+
+#include <dumux/implicit/properties.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/fluxvariablescaching.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+// forward declaration of properties
+NEW_PROP_TAG(NumPhases);
+NEW_PROP_TAG(ReplaceCompEqIdx);
+NEW_PROP_TAG(UseMoles);
+}
+
+/*!
+ * \ingroup StaggeredFicksLaw
+ * \brief Specialization of Fick's Law for the staggered free flow method.
+ */
+template <class TypeTag>
+class FicksLawImplementation<TypeTag, DiscretizationMethods::Staggered >
+{
+ 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 SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+ using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+ static const int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+
+ using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+
+ 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);
+
+ enum {
+
+ pressureIdx = Indices::pressureIdx,
+ velocityIdx = Indices::velocityIdx,
+
+ massBalanceIdx = Indices::massBalanceIdx,
+ momentumBalanceIdx = Indices::momentumBalanceIdx,
+ conti0EqIdx = Indices::conti0EqIdx
+ };
+
+public:
+ // state the discretization method this implementation belongs to
+ static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::Staggered;
+
+ //! state the type for the corresponding cache and its filler
+ //! We don't cache anything for this law
+ using Cache = FluxVariablesCaching::EmptyDiffusionCache;
+ using CacheFiller = FluxVariablesCaching::EmptyCacheFiller<TypeTag>;
+
+ static CellCenterPrimaryVariables diffusiveFluxForCellCenter(const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvf)
+ {
+ CellCenterPrimaryVariables flux(0.0);
+
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& insideVolVars = elemVolVars[insideScv];
+ const auto& outsideVolVars = scvf.boundary() ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
+
+ const Scalar insideDensity = useMoles ? insideVolVars.molarDensity() : insideVolVars.density();
+
+ for(int compIdx = 1; compIdx < numComponents; ++compIdx)
+ {
+ auto eqIdx = conti0EqIdx + compIdx;
+
+ const Scalar tij = transmissibility_(problem, fvGeometry, elemVolVars, scvf, compIdx);
+ const Scalar insideFraction = useMoles ? insideVolVars.moleFraction(0, compIdx) : insideVolVars.massFraction(0, compIdx);
+
+ if(scvf.boundary())
+ {
+ const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
+ if(bcTypes.isOutflow(eqIdx))
+ return flux;
+ else if(bcTypes.isNeumann(eqIdx))
+ return flux; // TODO: implement neumann
+ else
+ {
+ const Scalar dirichletFraction = problem.dirichletAtPos(scvf.center())[eqIdx];
+ flux[eqIdx] = insideDensity * tij * (insideFraction - dirichletFraction);
+ }
+ }
+ else
+ {
+ const Scalar outsideDensity = useMoles ? outsideVolVars.molarDensity() : outsideVolVars.density();
+ const Scalar avgDensity = 0.5*(insideDensity + outsideDensity);
+ const Scalar outsideFraction = useMoles ? outsideVolVars.moleFraction(0, compIdx) : outsideVolVars.massFraction(0, compIdx);
+ flux[eqIdx] = avgDensity * tij * (insideFraction - outsideFraction);
+ }
+ }
+
+ if (replaceCompEqIdx >= numComponents)
+ {
+ const Scalar cumulativeFlux = std::accumulate(flux.begin(), flux.end(), 0.0);
+ flux[0] = - cumulativeFlux;
+ }
+
+ return flux;
+ }
+
+ static Scalar transmissibility_(const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const int compIdx)
+ {
+ Scalar tij = 0.0;
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
+ const auto& insideVolVars = elemVolVars[insideScv];
+ const auto& outsideVolVars = scvf.boundary() ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
+
+ const Scalar insideDistance = (insideScv.dofPosition() - scvf.ipGlobal()).two_norm();
+ const Scalar insideD = insideVolVars.diffusionCoefficient(0, compIdx);
+ const Scalar ti = calculateOmega_(insideDistance, insideD, 1.0);
+
+ if(scvf.boundary())
+ tij = scvf.area() * ti;
+ else
+ {
+ const Scalar outsideDistance = (outsideScv.dofPosition() - scvf.ipGlobal()).two_norm();
+ const Scalar outsideD = outsideVolVars.diffusionCoefficient(0, compIdx);
+ const Scalar tj = calculateOmega_(outsideDistance, outsideD, 1.0);
+
+ tij = scvf.area()*(ti * tj)/(ti + tj);
+ }
+ return tij;
+ }
+
+ static Scalar calculateOmega_(const Scalar distance,
+ const Scalar D,
+ const Scalar extrusionFactor)
+ {
+ Scalar omega = D / distance;
+ omega *= extrusionFactor;
+
+ return omega;
+ }
+
+};
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/staggered/fluxvariables.hh b/dumux/freeflow/staggered/fluxvariables.hh
index 4679085dac9e968c5ecc1e248c1cf11599a0bea3..f7c2698574a47975fdc5a58f697239f4c0a2e95d 100644
--- a/dumux/freeflow/staggered/fluxvariables.hh
+++ b/dumux/freeflow/staggered/fluxvariables.hh
@@ -367,217 +367,6 @@ private:
}
};
-
-// specialization for miscible, isothermal flow
-template<class TypeTag>
-class FreeFlowFluxVariablesImpl<TypeTag, true, false> : public FreeFlowFluxVariablesImpl<TypeTag, false, 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 GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
- 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>;
-
- static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
- static constexpr auto numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
-
- static constexpr bool useMoles = true;
-
- //! 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, 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
- };
-
-public:
- CellCenterPrimaryVariables computeFluxForCellCenter(const Problem& problem,
- const Element &element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const GlobalFaceVars& globalFaceVars,
- const SubControlVolumeFace &scvf,
- const FluxVariablesCache& fluxVarsCache)
- {
- CellCenterPrimaryVariables flux(0.0);
-
- flux += advectiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, globalFaceVars, scvf);
- flux += diffusiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, scvf);
- return flux;
- }
-
-private:
-
- CellCenterPrimaryVariables advectiveFluxForCellCenter_(const Problem& problem,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const GlobalFaceVars& globalFaceVars,
- const SubControlVolumeFace &scvf)
- {
- CellCenterPrimaryVariables flux(0.0);
-
- const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
- const auto& insideVolVars = elemVolVars[insideScv];
-
- // if we are on an inflow/outflow boundary, use the volVars of the element itself
- const auto& outsideVolVars = scvf.boundary() ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
-
- const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
-
- const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity) ? true : false;
- const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
- const auto& downstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
-
- const Scalar upWindWeight = GET_PROP_VALUE(TypeTag, ImplicitUpwindWeight);
-
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- {
- // get equation index
- auto eqIdx = conti0EqIdx + compIdx;
-
- const Scalar upstreamDensity = useMoles ? upstreamVolVars.molarDensity() : upstreamVolVars.density();
- const Scalar downstreamDensity = useMoles ? downstreamVolVars.molarDensity() : downstreamVolVars.density();
- const Scalar upstreamFraction = useMoles ? upstreamVolVars.moleFraction(0, compIdx) : upstreamVolVars.massFraction(0, compIdx);
- const Scalar downstreamFraction = useMoles ? downstreamVolVars.moleFraction(0, compIdx) : downstreamVolVars.massFraction(0, compIdx);
-
- Scalar advFlux = 0.0;
-
- if(scvf.boundary() && eqIdx > conti0EqIdx)
- {
- const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
- if(bcTypes.isDirichlet(eqIdx))
- advFlux = upstreamDensity * problem.dirichletAtPos(scvf.center())[eqIdx] * velocity;
- if(bcTypes.isOutflow(eqIdx))
- advFlux = upstreamDensity * upstreamFraction * velocity;
- }
- else
- advFlux = (upWindWeight * upstreamDensity * upstreamFraction +
- (1.0 - upWindWeight) * downstreamDensity * downstreamFraction) * velocity;
-
- if (eqIdx != replaceCompEqIdx)
- flux[eqIdx] += advFlux;
-
- // in case one balance is substituted by the total mass balance
- if (replaceCompEqIdx < numComponents)
- flux[replaceCompEqIdx] += advFlux;
- }
-
- flux *= scvf.area() * sign(scvf.outerNormalScalar());
- return flux;
- }
-
-
- CellCenterPrimaryVariables diffusiveFluxForCellCenter_(const Problem& problem,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace &scvf)
- {
- CellCenterPrimaryVariables flux(0.0);
-
- const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
- const auto& insideVolVars = elemVolVars[insideScv];
- const auto& outsideVolVars = scvf.boundary() ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
-
- const Scalar insideDensity = useMoles ? insideVolVars.molarDensity() : insideVolVars.density();
-
- for(int compIdx = 1; compIdx < numComponents; ++compIdx)
- {
- auto eqIdx = conti0EqIdx + compIdx;
-
- const Scalar tij = transmissibility_(problem, fvGeometry, elemVolVars, scvf, compIdx);
- const Scalar insideFraction = useMoles ? insideVolVars.moleFraction(0, compIdx) : insideVolVars.massFraction(0, compIdx);
-
- if(scvf.boundary())
- {
- const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
- if(bcTypes.isOutflow(eqIdx))
- return flux;
- else if(bcTypes.isNeumann(eqIdx))
- return flux; // TODO: implement neumann
- else
- {
- const Scalar dirichletFraction = problem.dirichletAtPos(scvf.center())[eqIdx];
- flux[eqIdx] = insideDensity * tij * (insideFraction - dirichletFraction);
- }
- }
- else
- {
- const Scalar outsideDensity = useMoles ? outsideVolVars.molarDensity() : outsideVolVars.density();
- const Scalar avgDensity = 0.5*(insideDensity + outsideDensity);
- const Scalar outsideFraction = useMoles ? outsideVolVars.moleFraction(0, compIdx) : outsideVolVars.massFraction(0, compIdx);
- flux[eqIdx] = avgDensity * tij * (insideFraction - outsideFraction);
- }
- }
-
- if (replaceCompEqIdx >= numComponents)
- {
- const Scalar cumulativeFlux = std::accumulate(flux.begin(), flux.end(), 0.0);
- flux[0] = - cumulativeFlux;
- }
-
- return flux;
- }
-
- static Scalar transmissibility_(const Problem& problem,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace& scvf,
- const int compIdx)
- {
- Scalar tij = 0.0;
- const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
- const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
- const auto& insideVolVars = elemVolVars[insideScv];
- const auto& outsideVolVars = scvf.boundary() ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
-
- const Scalar insideDistance = (insideScv.dofPosition() - scvf.ipGlobal()).two_norm();
- const Scalar insideD = insideVolVars.diffusionCoefficient(0, compIdx);
- const Scalar ti = calculateOmega_(insideDistance, insideD, 1.0);
-
- if(scvf.boundary())
- tij = scvf.area() * ti;
- else
- {
- const Scalar outsideDistance = (outsideScv.dofPosition() - scvf.ipGlobal()).two_norm();
- const Scalar outsideD = outsideVolVars.diffusionCoefficient(0, compIdx);
- const Scalar tj = calculateOmega_(outsideDistance, outsideD, 1.0);
-
- tij = scvf.area()*(ti * tj)/(ti + tj);
- }
- return tij;
- }
-
- static Scalar calculateOmega_(const Scalar distance,
- const Scalar D,
- const Scalar extrusionFactor)
- {
- Scalar omega = D / distance;
- omega *= extrusionFactor;
-
- return omega;
- }
-};
-
} // end namespace
#endif
diff --git a/dumux/freeflow/staggerednc/fluxvariables.hh b/dumux/freeflow/staggerednc/fluxvariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..44598f6c6081b1e12b008b9430384ae3115d6a7a
--- /dev/null
+++ b/dumux/freeflow/staggerednc/fluxvariables.hh
@@ -0,0 +1,178 @@
+// -*- 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_FREELOW_IMPLICIT_NC_FLUXVARIABLES_HH
+#define DUMUX_FREELOW_IMPLICIT_NC_FLUXVARIABLES_HH
+
+#include <dumux/implicit/properties.hh>
+#include <dumux/discretization/fluxvariablesbase.hh>
+#include "../staggered/fluxvariables.hh"
+
+namespace Dumux
+{
+
+namespace Properties
+{
+// forward declaration
+NEW_PROP_TAG(EnableComponentTransport);
+NEW_PROP_TAG(EnableEnergyBalance);
+NEW_PROP_TAG(EnableInertiaTerms);
+}
+
+// // forward declaration
+// template<class TypeTag, bool enableComponentTransport, bool enableEnergyBalance>
+// class FreeFlowFluxVariablesImpl;
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief The flux variables class
+ * specializations are provided for combinations of physical processes
+ * \note Not all specializations are currently implemented
+ */
+// template<class TypeTag>
+// using FreeFlowFluxVariables = FreeFlowFluxVariablesImpl<TypeTag, GET_PROP_VALUE(TypeTag, EnableComponentTransport),
+// GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>;
+
+// specialization for miscible, isothermal flow
+template<class TypeTag>
+class FreeFlowFluxVariablesImpl<TypeTag, true, false> : public FreeFlowFluxVariablesImpl<TypeTag, false, 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 GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
+ 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, 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
+ };
+
+public:
+ CellCenterPrimaryVariables computeFluxForCellCenter(const Problem& problem,
+ const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const GlobalFaceVars& globalFaceVars,
+ const SubControlVolumeFace &scvf,
+ const FluxVariablesCache& fluxVarsCache)
+ {
+ CellCenterPrimaryVariables flux(0.0);
+
+ flux += advectiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, globalFaceVars, scvf);
+ flux += MolecularDiffusionType::diffusiveFluxForCellCenter(problem, fvGeometry, elemVolVars, scvf);
+ return flux;
+ }
+
+private:
+
+ CellCenterPrimaryVariables advectiveFluxForCellCenter_(const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const GlobalFaceVars& globalFaceVars,
+ const SubControlVolumeFace &scvf)
+ {
+ CellCenterPrimaryVariables flux(0.0);
+
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& insideVolVars = elemVolVars[insideScv];
+
+ // if we are on an inflow/outflow boundary, use the volVars of the element itself
+ const auto& outsideVolVars = scvf.boundary() ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
+
+ const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
+
+ const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity) ? true : false;
+ const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
+ const auto& downstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
+
+ const Scalar upWindWeight = GET_PROP_VALUE(TypeTag, ImplicitUpwindWeight);
+
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ // get equation index
+ auto eqIdx = conti0EqIdx + compIdx;
+
+ const Scalar upstreamDensity = useMoles ? upstreamVolVars.molarDensity() : upstreamVolVars.density();
+ const Scalar downstreamDensity = useMoles ? downstreamVolVars.molarDensity() : downstreamVolVars.density();
+ const Scalar upstreamFraction = useMoles ? upstreamVolVars.moleFraction(0, compIdx) : upstreamVolVars.massFraction(0, compIdx);
+ const Scalar downstreamFraction = useMoles ? downstreamVolVars.moleFraction(0, compIdx) : downstreamVolVars.massFraction(0, compIdx);
+
+ Scalar advFlux = 0.0;
+
+ if(scvf.boundary() && eqIdx > conti0EqIdx)
+ {
+ const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
+ if(bcTypes.isDirichlet(eqIdx))
+ advFlux = upstreamDensity * problem.dirichletAtPos(scvf.center())[eqIdx] * velocity;
+ if(bcTypes.isOutflow(eqIdx))
+ advFlux = upstreamDensity * upstreamFraction * velocity;
+ }
+ else
+ advFlux = (upWindWeight * upstreamDensity * upstreamFraction +
+ (1.0 - upWindWeight) * downstreamDensity * downstreamFraction) * velocity;
+
+ if (eqIdx != replaceCompEqIdx)
+ flux[eqIdx] += advFlux;
+
+ // in case one balance is substituted by the total mass balance
+ if (replaceCompEqIdx < numComponents)
+ flux[replaceCompEqIdx] += advFlux;
+ }
+
+ flux *= scvf.area() * sign(scvf.outerNormalScalar());
+ return flux;
+ }
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/staggerednc/localresidual.hh b/dumux/freeflow/staggerednc/localresidual.hh
index 6b085a835b789f33ee98442291c3671641c52f90..9496f9ca7a2ebbd2c0de967c9461538765c356fa 100644
--- a/dumux/freeflow/staggerednc/localresidual.hh
+++ b/dumux/freeflow/staggerednc/localresidual.hh
@@ -111,6 +111,8 @@ class StaggeredNavierStokesResidualImpl<TypeTag, true, false> : public Staggered
static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+ static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+
static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
//! The index of the component balance equation that gets replaced with the total mass balance
@@ -126,8 +128,7 @@ class StaggeredNavierStokesResidualImpl<TypeTag, true, false> : public Staggered
* the implicit euler time derivative here
*/
CellCenterPrimaryVariables computeStorageForCellCenter(const SubControlVolume& scv,
- const VolumeVariables& volVars,
- bool useMoles = true) const
+ const VolumeVariables& volVars) const
{
CellCenterPrimaryVariables storage(0.0);
diff --git a/dumux/freeflow/staggerednc/model.hh b/dumux/freeflow/staggerednc/model.hh
index e7111ac7a4535d4d1d6a9c409abd3ca446e4f6ed..01198e50f79e7471ba51d0d4a5ca82211c16d6da 100644
--- a/dumux/freeflow/staggerednc/model.hh
+++ b/dumux/freeflow/staggerednc/model.hh
@@ -93,6 +93,7 @@ public:
// register standardized vtk output fields
auto& vtkOutputModule = problem.vtkOutputModule();
+ vtkOutputModule.addSecondaryVariable("rho",[](const VolumeVariables& v){ return v.molarDensity(); });
for (int j = 0; j < numComponents; ++j)
vtkOutputModule.addSecondaryVariable("x" + FluidSystem::componentName(j) + FluidSystem::phaseName(0),
[j](const VolumeVariables& v){ return v.massFraction(0,j); });
diff --git a/dumux/freeflow/staggerednc/properties.hh b/dumux/freeflow/staggerednc/properties.hh
index da3f19c7a7ed8674f4e627d45b4fbf0ce553c872..ef5c28a23f25ec657aeaace2ad3142b3f590eebe 100644
--- a/dumux/freeflow/staggerednc/properties.hh
+++ b/dumux/freeflow/staggerednc/properties.hh
@@ -67,6 +67,7 @@ NEW_PROP_TAG(EnableComponentTransport); //!< Returns whether to consider compone
NEW_PROP_TAG(EnableEnergyTransport); //!< Returns whether to consider energy transport or not
NEW_PROP_TAG(FaceVariables); //!< Returns whether to consider energy transport or not
NEW_PROP_TAG(ReplaceCompEqIdx); //!< Returns whether to consider energy transport or not
+NEW_PROP_TAG(UseMoles); //!< Defines whether molar (true) or mass (false) density is used
// \}
}
diff --git a/dumux/freeflow/staggerednc/propertydefaults.hh b/dumux/freeflow/staggerednc/propertydefaults.hh
index 809511530cb5102bd0e59ddb7e383cf06ba4a485..31d7d5b87a472ad7d18d0f8a91d3c1e4f729ac69 100644
--- a/dumux/freeflow/staggerednc/propertydefaults.hh
+++ b/dumux/freeflow/staggerednc/propertydefaults.hh
@@ -33,6 +33,7 @@
#include "volumevariables.hh"
#include "indices.hh"
#include "localresidual.hh"
+#include "fluxvariables.hh"
#include "../staggered/problem.hh"
// #include "../staggered/model.hh"
#include "../staggered/propertydefaults.hh"
@@ -47,8 +48,6 @@
#include <dumux/material/fluidstates/compositional.hh>
-
-
namespace Dumux
{
@@ -166,6 +165,8 @@ SET_PROP(NavierStokesNC, FluidState)
//
SET_BOOL_PROP(NavierStokesNC, EnableComponentTransport, true);
+SET_BOOL_PROP(NavierStokesNC, UseMoles, false); //!< Defines whether molar (true) or mass (false) density is used
+
//! average is used as default model to compute the effective thermal heat conductivity
// SET_PROP(NavierStokesNI, ThermalConductivityModel)