diff --git a/dumux/discretization/staggered/globalvolumevariables.hh b/dumux/discretization/staggered/globalvolumevariables.hh
index e9d4ae7446ac1c5433066c03a617a43ff88756e0..d5f3b4868638b34af34e9c44dd9edbb529657c40 100644
--- a/dumux/discretization/staggered/globalvolumevariables.hh
+++ b/dumux/discretization/staggered/globalvolumevariables.hh
@@ -114,7 +114,8 @@ public:
//TODO: this assumes a zero-gradient for e.g. the pressure on the boundary
// could be made more general by allowing a non-zero-gradient, provided in problem file
else
- DUNE_THROW(Dune::InvalidStateException, "Face at: " << scvf.center() << " has neither Dirichlet nor Neumann BC.");
+ if(eqIdx == Indices::pressureIdx)
+ DUNE_THROW(Dune::InvalidStateException, "Face at: " << scvf.center() << " has neither Dirichlet nor Neumann BC.");
}
ElementSolutionVector elemSol{std::move(boundaryPriVars)};
volumeVariables_[scvf.outsideScvIdx()].update(elemSol, problem, element, insideScv);
diff --git a/dumux/freeflow/CMakeLists.txt b/dumux/freeflow/CMakeLists.txt
index aff2016afb1aa6c9440e1850b201fb6bcc4b1c05..fecc1809e5282ceb6568e6a658e0c06d39443b99 100644
--- a/dumux/freeflow/CMakeLists.txt
+++ b/dumux/freeflow/CMakeLists.txt
@@ -1,4 +1,5 @@
add_subdirectory("staggered")
+add_subdirectory("staggerednc")
add_subdirectory("stokes")
add_subdirectory("stokesnc")
add_subdirectory("stokesncni")
diff --git a/dumux/freeflow/staggered/fluxvariables.hh b/dumux/freeflow/staggered/fluxvariables.hh
index edabf142e1d11634d83161763b63144b0d297403..c55b69b50602d0144bc9e48700d9cdcba5fb036a 100644
--- a/dumux/freeflow/staggered/fluxvariables.hh
+++ b/dumux/freeflow/staggered/fluxvariables.hh
@@ -106,17 +106,21 @@ public:
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& insideVolVars = elemVolVars[insideScv];
- const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
-
// if we are on an inflow/outflow boundary, use the volVars of the element itself
const auto& outsideVolVars = scvf.boundary() ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
CellCenterPrimaryVariables flux(0.0);
+ 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);
+
+ flux = (upWindWeight * upstreamVolVars.density() +
+ (1.0 - upWindWeight) * downstreamVolVars.density()) * velocity;
- if(sign(scvf.outerNormalScalar()) == sign(velocity))
- flux[0] = insideVolVars.density() * velocity;
- else
- flux[0] = outsideVolVars.density() * velocity;
return flux * scvf.area() * sign(scvf.outerNormalScalar());
}
@@ -365,8 +369,7 @@ private:
// specialization for miscible, isothermal flow
template<class TypeTag>
-class FreeFlowFluxVariablesImpl<TypeTag, true, false>
-: public FluxVariablesBase<TypeTag>, public FreeFlowFluxVariablesImpl<TypeTag, false, false>
+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);
@@ -386,6 +389,11 @@ class FreeFlowFluxVariablesImpl<TypeTag, true, false>
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 {
@@ -409,10 +417,101 @@ public:
const SubControlVolumeFace &scvf,
const FluxVariablesCache& fluxVarsCache)
{
- CellCenterPrimaryVariables priVars(0.0);
- priVars[0] = ParentType::computeFluxForCellCenter(element, fvGeometry, elemVolVars, globalFaceVars, scvf, fluxVarsCache);
-// priVars[1] = ParentType::computeFluxForCellCenter(element, fvGeometry, elemVolVars, globalFaceVars, scvf, fluxVarsCache);
- return priVars;
+ 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()];
+
+ CellCenterPrimaryVariables flux(0.0);
+ 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 = this->problem().boundaryTypesAtPos(scvf.center());
+ if(bcTypes.isDirichlet(eqIdx))
+ advFlux = upstreamDensity * this->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;
+ }
+
+ //TODO: implement diffusive fluxes
+
+ flux *= scvf.area() * sign(scvf.outerNormalScalar());
+
+ return flux;
+ }
+
+
+ CellCenterPrimaryVariables diffusiveFlux(const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvf)
+ {
+ CellCenterPrimaryVariables flux(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()];
+
+ Scalar distance = (insideScv.dofPosition() - scvf.ipGlobal()).two_norm2();
+
+ Scalar xout = 0.1;
+ Scalar xin = insideVolVars.moleFraction(0,1);
+
+ Scalar Din = insideVolVars.diffusionCoefficient(0, 1);
+ Scalar mean = Din;
+
+
+
+ if(scvf.boundary())
+ return flux;
+
+ // const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
+ // const auto& insideVolVars = elemVolVars[insideScv];
+ // const auto& outsideVolVars = scvf.boundary() ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
+ //
+ // auto distance = (outsideScv.dofPosition() - insideScv.dofPosition()).two_norm();
+ // Scalar xout = outsideVolVars.moleFraction(0,1);
+ // Scalar xin = insideVolVars.moleFraction(0,1);
+ //
+ // Scalar Dout = outsideVolVars.diffusionCoefficient(0, 1);
+ // Scalar Din = insideVolVars.diffusionCoefficient(0, 1);
+ //
+ // Scalar mean = 0.5*(Dout + Din);
+ //
+ // flux[1] = -(xout - xin)/distance *mean;
+
+ return flux;
}
};
diff --git a/dumux/freeflow/staggered/localresidual.hh b/dumux/freeflow/staggered/localresidual.hh
index e6210b546e063233ce84cb7175d56375439d099f..b0a25c8b36166508335a1f86b89d31a719fb3a9d 100644
--- a/dumux/freeflow/staggered/localresidual.hh
+++ b/dumux/freeflow/staggered/localresidual.hh
@@ -39,6 +39,7 @@ namespace Properties
NEW_PROP_TAG(EnableComponentTransport);
NEW_PROP_TAG(EnableEnergyBalance);
NEW_PROP_TAG(EnableInertiaTerms);
+NEW_PROP_TAG(ReplaceCompEqIdx);
}
/*!
@@ -262,8 +263,15 @@ protected:
{
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& insideVolVars = elemVolVars[insideScv];
- this->ccResidual_[pressureIdx] = insideVolVars.pressure() - this->problem().dirichletAtPos(insideScv.dofPosition())[cellCenterIdx][pressureIdx];
+ this->ccResidual_[0] = insideVolVars.pressure() - this->problem().dirichletAtPos(insideScv.dofPosition())[cellCenterIdx][0];
}
+
+// if(scvf.center()[0] < 1e-5)
+// {
+// const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+// const auto& insideVolVars = elemVolVars[insideScv];
+// this->ccResidual_[1] = insideVolVars.moleFraction(0,1) - 0.1;
+// }
}
}
}
@@ -336,91 +344,132 @@ private:
}
};
-// specialization for miscible, isothermal flow
-template<class TypeTag>
-class StaggeredNavierStokesResidualImpl<TypeTag, true, false> : public StaggeredNavierStokesResidualImpl<TypeTag, false, false>
-{
- using ParentType = StaggeredNavierStokesResidualImpl<TypeTag, false, false>;
- friend class StaggeredLocalResidual<TypeTag>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-
- enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
-
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual);
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using Element = typename GridView::template Codim<0>::Entity;
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
- using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
- using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using CellCenterSolutionVector = typename GET_PROP_TYPE(TypeTag, CellCenterSolutionVector);
- using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
- using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
- using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
- using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
- using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
-
-
- using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
- typename DofTypeIndices::CellCenterIdx cellCenterIdx;
- typename DofTypeIndices::FaceIdx faceIdx;
-
- 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
- };
-
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
-
- static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
- static constexpr auto numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
- /*!
- * \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 SubControlVolume& scv,
- const VolumeVariables& volVars)
- {
- CellCenterPrimaryVariables storage;
- // compute storage term of all components within all phases
-
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- {
- auto eqIdx = conti0EqIdx + compIdx;
- auto s = volVars.molarDensity(0)
- * volVars.moleFraction(0, compIdx);
-
-// if (eqIdx != replaceCompEqIdx)
- storage[eqIdx] += s;
-
- // in case one balance is substituted by the total mass balance
-// if (replaceCompEqIdx < numComponents)
-// storage[replaceCompEqIdx] += s;
- }
- }
-};
+// // specialization for miscible, isothermal flow
+// template<class TypeTag>
+// class StaggeredNavierStokesResidualImpl<TypeTag, true, false> : public StaggeredNavierStokesResidualImpl<TypeTag, false, false>
+// {
+// using ParentType = StaggeredNavierStokesResidualImpl<TypeTag, false, false>;
+// friend class StaggeredLocalResidual<TypeTag>;
+// using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+//
+// enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+//
+// using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+// using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual);
+// using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+// using Element = typename GridView::template Codim<0>::Entity;
+// using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+// using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+// using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
+// using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+// using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+// using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+// using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+// using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+// using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+// using CellCenterSolutionVector = typename GET_PROP_TYPE(TypeTag, CellCenterSolutionVector);
+// using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
+// using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+// using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+// using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+// using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+//
+//
+// using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+// typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+// typename DofTypeIndices::FaceIdx faceIdx;
+//
+// 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
+// };
+//
+// using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+// using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
+//
+// static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
+// static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+//
+// static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+//
+// //! The index of the component balance equation that gets replaced with the total mass balance
+// static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
+// /*!
+// * \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 SubControlVolume& scv,
+// const VolumeVariables& volVars,
+// bool useMoles = true) const
+// {
+// PrimaryVariables storage(0.0);
+//
+// // formulation with mole balances
+// if (useMoles)
+// {
+// // compute storage term of all components
+// for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+// {
+// auto eqIdx = conti0EqIdx + compIdx;
+// auto s = volVars.molarDensity(0)
+// * volVars.moleFraction(0, compIdx);
+//
+// if (eqIdx != replaceCompEqIdx)
+// storage[eqIdx] += s;
+//
+// // in case one balance is substituted by the total mass balance
+// if (replaceCompEqIdx < numComponents)
+// storage[replaceCompEqIdx] += s;
+// }
+//
+// //! The energy storage in the fluid phase with index phaseIdx
+// // EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
+//
+// //! The energy storage in the solid matrix
+// // EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
+// }
+// // formulation with mass balances
+// else
+// {
+// for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+// {
+// auto eqIdx = conti0EqIdx + compIdx;
+// auto s = volVars.density(0)
+// * volVars.massFraction(0, compIdx);
+//
+// if (eqIdx != replaceCompEqIdx)
+// storage[eqIdx] += s;
+//
+// // in case one balance is substituted by the total mass balance
+// if (replaceCompEqIdx < numComponents)
+// storage[replaceCompEqIdx] += s;
+// }
+//
+// //! The energy storage in the fluid phase with index phaseIdx
+// // EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
+//
+// //! The energy storage in the solid matrix
+// // EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
+// }
+//
+// return storage;
+// }
+// };
// specialization for immiscible, non-isothermal flow
template<class TypeTag>
diff --git a/dumux/freeflow/staggerednc/CMakeLists.txt b/dumux/freeflow/staggerednc/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8baf812d982cc188659813332319df4a48e74783
--- /dev/null
+++ b/dumux/freeflow/staggerednc/CMakeLists.txt
@@ -0,0 +1,10 @@
+
+#install headers
+install(FILES
+indices.hh
+localresidual.hh
+model.hh
+properties.hh
+propertydefaults.hh
+volumevariables.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/1p/implicit)
diff --git a/dumux/freeflow/staggerednc/indices.hh b/dumux/freeflow/staggerednc/indices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3d2383f480d578e36c76d7c1420b516178f9fc31
--- /dev/null
+++ b/dumux/freeflow/staggerednc/indices.hh
@@ -0,0 +1,47 @@
+// -*- 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 Defines the indices for the one-phase fully implicit model.
+ */
+#ifndef DUMUX_STAGGERED_NAVIERSTOKES_NC_INDICES_HH
+#define DUMUX_STAGGERED_NAVIERSTOKES_NC_INDICES_HH
+
+#include <dumux/freeflow/staggered/indices.hh>
+
+namespace Dumux
+{
+// \{
+/*!
+ * \ingroup NavierStokesModel
+ * \ingroup ImplicitIndices
+ * \brief The common indices for the isothermal stokes model.
+ *
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag, int PVOffset = 0>
+struct NavierStokesNCIndices : public NavierStokesCommonIndices<TypeTag, PVOffset>
+{
+
+};
+
+// \}
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/staggerednc/localresidual.hh b/dumux/freeflow/staggerednc/localresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..2edce91ecc649120e336676e3d5ec8a2b6ae847f
--- /dev/null
+++ b/dumux/freeflow/staggerednc/localresidual.hh
@@ -0,0 +1,187 @@
+// -*- 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 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
+
+#include <dune/istl/matrix.hh>
+
+#include <dumux/common/valgrind.hh>
+#include <dumux/implicit/staggered/localresidual.hh>
+
+#include "properties.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, bool enableEnergyBalance>
+class StaggeredNavierStokesResidualImpl;
+
+// specialization for miscible, isothermal flow
+template<class TypeTag>
+class StaggeredNavierStokesResidualImpl<TypeTag, true, false> : public StaggeredNavierStokesResidualImpl<TypeTag, false, false>
+{
+ using ParentType = StaggeredNavierStokesResidualImpl<TypeTag, false, false>;
+ friend class StaggeredLocalResidual<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using CellCenterSolutionVector = typename GET_PROP_TYPE(TypeTag, CellCenterSolutionVector);
+ using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
+ using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+
+
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+ 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
+ };
+
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
+
+ static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
+ static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+
+ static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+
+ //! The index of the component balance equation that gets replaced with the total mass balance
+ static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
+ /*!
+ * \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 SubControlVolume& scv,
+ const VolumeVariables& volVars,
+ bool useMoles = false) const
+ {
+ CellCenterPrimaryVariables storage(0.0);
+
+ // formulation with mole balances
+ if (useMoles)
+ {
+ // compute storage term of all components
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ auto eqIdx = conti0EqIdx + compIdx;
+ auto s = volVars.molarDensity(0)
+ * volVars.moleFraction(0, compIdx);
+
+ if (eqIdx != replaceCompEqIdx)
+ storage[eqIdx] += s;
+
+ // in case one balance is substituted by the total mass balance
+ if (replaceCompEqIdx < numComponents)
+ storage[replaceCompEqIdx] += s;
+ }
+
+ //! The energy storage in the fluid phase with index phaseIdx
+// EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
+
+ //! The energy storage in the solid matrix
+// EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
+ }
+ // formulation with mass balances
+ else
+ {
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ auto eqIdx = conti0EqIdx + compIdx;
+ auto s = volVars.density(0)
+ * volVars.massFraction(0, compIdx);
+
+ if (eqIdx != replaceCompEqIdx)
+ storage[eqIdx] += s;
+
+ // in case one balance is substituted by the total mass balance
+ if (replaceCompEqIdx < numComponents)
+ storage[replaceCompEqIdx] += s;
+ }
+
+ //! The energy storage in the fluid phase with index phaseIdx
+// EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
+
+ //! The energy storage in the solid matrix
+// EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
+ }
+
+ return storage;
+ }
+};
+}
+
+#endif // DUMUX_CC_LOCAL_RESIDUAL_HH
diff --git a/dumux/freeflow/staggerednc/model.hh b/dumux/freeflow/staggerednc/model.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e7111ac7a4535d4d1d6a9c409abd3ca446e4f6ed
--- /dev/null
+++ b/dumux/freeflow/staggerednc/model.hh
@@ -0,0 +1,107 @@
+// -*- 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 all models which use the one-phase,
+ * fully implicit model.
+ * Adaption of the fully implicit scheme to the one-phase flow model.
+ */
+
+#ifndef DUMUX_NAVIERSTOKES_NC_MODEL_HH
+#define DUMUX_NAVIERSTOKES_NC_MODEL_HH
+
+// #include <dumux/porousmediumflow/implicit/velocityoutput.hh>
+#include "properties.hh"
+#include "../staggered/model.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup NavierStokesModel
+ * \brief A single-phase, isothermal flow model using the fully implicit scheme.
+ *
+ * Single-phase, isothermal flow model, which uses a standard Darcy approach as the
+ * equation for the conservation of momentum:
+ * \f[
+ v = - \frac{\textbf K}{\mu}
+ \left(\textbf{grad}\, p - \varrho {\textbf g} \right)
+ * \f]
+ *
+ * and solves the mass continuity equation:
+ * \f[
+ \phi \frac{\partial \varrho}{\partial t} + \text{div} \left\lbrace
+ - \varrho \frac{\textbf K}{\mu} \left( \textbf{grad}\, p -\varrho {\textbf g} \right) \right\rbrace = q,
+ * \f]
+ * All equations are discretized using a vertex-centered finite volume (box)
+ * or cell-centered finite volume scheme as spatial
+ * and the implicit Euler method as time discretization.
+ * The model supports compressible as well as incompressible fluids.
+ */
+template<class TypeTag >
+class NavierStokesNCModel : public NavierStokesModel<TypeTag>
+{
+ using ParentType = NavierStokesModel<TypeTag>;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GlobalFVGeometry) GlobalFVGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
+
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ enum { dim = GridView::dimension };
+ enum { dimWorld = GridView::dimensionworld };
+
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+
+ static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+public:
+
+ void init(Problem& problem)
+ {
+ ParentType::init(problem);
+
+ // register standardized vtk output fields
+ auto& vtkOutputModule = problem.vtkOutputModule();
+ 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); });
+
+// NonIsothermalModel::maybeAddTemperature(vtkOutputModule);
+ }
+};
+}
+
+#include "propertydefaults.hh"
+
+#endif
diff --git a/dumux/freeflow/staggerednc/properties.hh b/dumux/freeflow/staggerednc/properties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..da3f19c7a7ed8674f4e627d45b4fbf0ce553c872
--- /dev/null
+++ b/dumux/freeflow/staggerednc/properties.hh
@@ -0,0 +1,75 @@
+// -*- 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/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup Properties
+ * \ingroup ImplicitProperties
+ * \ingroup NavierStokesModel
+ * \file
+ *
+ * \brief Defines the properties required for the one-phase fully implicit model.
+ */
+#ifndef DUMUX_NAVIERSTOKES_NC_PROPERTIES_HH
+#define DUMUX_NAVIERSTOKES_NC_PROPERTIES_HH
+
+// #include <dumux/porousmediumflow/nonisothermal/implicit/properties.hh>
+#include <dumux/freeflow/staggered/properties.hh>
+
+namespace Dumux
+{
+// \{
+///////////////////////////////////////////////////////////////////////////
+// properties for the isothermal Navier-Stokes model
+///////////////////////////////////////////////////////////////////////////
+namespace Properties {
+
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tags for the implicit single-phase problems
+NEW_TYPE_TAG(NavierStokesNC, INHERITS_FROM(NavierStokes));
+
+//! The type tags for the corresponding non-isothermal problems
+// NEW_TYPE_TAG(NavierStokesNI, INHERITS_FROM(NavierStokes, NonIsothermal));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system
+NEW_PROP_TAG(Indices); //!< Enumerations for the model
+NEW_PROP_TAG(FluidSystem); //!< The type of the fluid system to use
+NEW_PROP_TAG(Fluid); //!< The fluid used for the default fluid system
+NEW_PROP_TAG(FluidState); //!< The type of the fluid state to use
+NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
+NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< Returns weight of the upwind cell when calculating fluxes
+NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< Weight for the upwind mobility in the velocity calculation
+NEW_PROP_TAG(VtkAddVelocity); //!< Returns whether velocity vectors are written into the vtk output
+NEW_PROP_TAG(EnableInertiaTerms); //!< Returns whether to include inertia terms in the momentum balance eq or not (Stokes / Navier-Stokes)
+NEW_PROP_TAG(BoundaryValues); //!< Type to set values on the boundary
+NEW_PROP_TAG(EnableComponentTransport); //!< Returns whether to consider component transport or not
+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
+// \}
+}
+
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/staggerednc/propertydefaults.hh b/dumux/freeflow/staggerednc/propertydefaults.hh
new file mode 100644
index 0000000000000000000000000000000000000000..809511530cb5102bd0e59ddb7e383cf06ba4a485
--- /dev/null
+++ b/dumux/freeflow/staggerednc/propertydefaults.hh
@@ -0,0 +1,203 @@
+// -*- 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/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup Properties
+ * \ingroup ImplicitProperties
+ * \ingroup OnePModel
+ * \file
+ *
+ * \brief Defines the properties required for the one-phase fully implicit model.
+ */
+#ifndef DUMUX_NAVIER_STOKES_NC_PROPERTY_DEFAULTS_HH
+#define DUMUX_NAVIER_STOKES_NC_PROPERTY_DEFAULTS_HH
+
+#include "properties.hh"
+
+#include "model.hh"
+#include "volumevariables.hh"
+#include "indices.hh"
+#include "localresidual.hh"
+#include "../staggered/problem.hh"
+// #include "../staggered/model.hh"
+#include "../staggered/propertydefaults.hh"
+
+
+#include <dumux/implicit/staggered/localresidual.hh>
+#include <dumux/material/fluidsystems/gasphase.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
+#include <dumux/material/components/nullcomponent.hh>
+#include <dumux/material/fluidsystems/1p.hh>
+
+#include <dumux/material/fluidstates/compositional.hh>
+
+
+
+
+namespace Dumux
+{
+
+namespace Properties
+{
+// forward declaration
+NEW_PROP_TAG(FluxVariables);
+NEW_PROP_TAG(FluxVariablesCache);
+}
+// \{
+
+///////////////////////////////////////////////////////////////////////////
+// default property values for the isothermal single phase model
+///////////////////////////////////////////////////////////////////////////
+namespace Properties {
+// SET_INT_PROP(NavierStokes, NumEqCellCenter, 1); //!< set the number of equations to 1
+// SET_INT_PROP(NavierStokes, NumEqFace, 1); //!< set the number of equations to 1
+// SET_INT_PROP(NavierStokes, NumPhases, 1); //!< The number of phases in the 1p model is 1
+//
+//
+SET_INT_PROP(NavierStokesNC, NumEqCellCenter, 2);
+SET_INT_PROP(NavierStokesNC, ReplaceCompEqIdx, 0);
+
+
+ /*!
+ * \brief Set the property for the number of components.
+ *
+ * We just forward the number from the fluid system
+ *
+ */
+SET_PROP(NavierStokesNC, NumComponents)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ static constexpr int value = FluidSystem::numComponents;
+
+};
+
+
+//! the VolumeVariables property
+SET_TYPE_PROP(NavierStokesNC, VolumeVariables, NavierStokesNCVolumeVariables<TypeTag>);
+SET_TYPE_PROP(NavierStokesNC, Model, NavierStokesNCModel<TypeTag>);
+
+
+/*!
+ * \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.
+ */
+SET_PROP(NavierStokesNC, FluidState)
+{
+ private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ public:
+ typedef CompositionalFluidState<Scalar, FluidSystem> type;
+};
+
+// //! Enable advection
+// SET_BOOL_PROP(NavierStokes, EnableAdvection, true);
+//
+// //! The one-phase model has no molecular diffusion
+// SET_BOOL_PROP(NavierStokes, EnableMolecularDiffusion, false);
+//
+// //! Isothermal model by default
+// SET_BOOL_PROP(NavierStokes, EnableEnergyBalance, false);
+//
+// //! The indices required by the isothermal single-phase model
+// SET_TYPE_PROP(NavierStokes, Indices, NavierStokesCommonIndices<TypeTag>);
+//
+// //! The weight of the upwind control volume when calculating
+// //! fluxes. Use central differences by default.
+// SET_SCALAR_PROP(NavierStokes, ImplicitMassUpwindWeight, 0.5);
+//
+// //! weight for the upwind mobility in the velocity calculation
+// //! fluxes. Use central differences by default.
+// SET_SCALAR_PROP(NavierStokes, ImplicitMobilityUpwindWeight, 0.5);
+//
+// //! The fluid system to use by default
+// SET_TYPE_PROP(NavierStokes, FluidSystem, Dumux::FluidSystems::OneP<typename GET_PROP_TYPE(TypeTag, Scalar), typename GET_PROP_TYPE(TypeTag, Fluid)>);
+//
+// SET_PROP(NavierStokes, Fluid)
+// { private:
+// typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+// public:
+// typedef FluidSystems::LiquidPhase<Scalar, Dumux::NullComponent<Scalar> > type;
+// };
+//
+// /*!
+// * \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.
+// */
+// SET_PROP(NavierStokes, FluidState){
+// private:
+// typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+// typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+// public:
+// typedef Dumux::ImmiscibleFluidState<Scalar, FluidSystem> type;
+// };
+//
+// // disable velocity output by default
+// SET_BOOL_PROP(NavierStokes, VtkAddVelocity, true);
+//
+// // enable gravity by default
+// SET_BOOL_PROP(NavierStokes, ProblemEnableGravity, true);
+//
+// SET_BOOL_PROP(NavierStokes, EnableInertiaTerms, true);
+//
+// SET_BOOL_PROP(NavierStokes, EnableEnergyTransport, false);
+//
+SET_BOOL_PROP(NavierStokesNC, EnableComponentTransport, true);
+
+
+//! average is used as default model to compute the effective thermal heat conductivity
+// SET_PROP(NavierStokesNI, ThermalConductivityModel)
+// { private :
+// typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+// public:
+// typedef ThermalConductivityAverage<Scalar> type;
+// };
+
+//////////////////////////////////////////////////////////////////
+// Property values for isothermal model required for the general non-isothermal model
+//////////////////////////////////////////////////////////////////
+
+// set isothermal Model
+// SET_TYPE_PROP(NavierStokesNI, IsothermalModel, NavierStokesModel<TypeTag>);
+
+//set isothermal VolumeVariables
+// SET_TYPE_PROP(NavierStokesNI, IsothermalVolumeVariables, NavierStokesVolumeVariables<TypeTag>);
+
+//set isothermal LocalResidual
+// SET_TYPE_PROP(NavierStokesNI, IsothermalLocalResidual, ImmiscibleLocalResidual<TypeTag>);
+
+//set isothermal Indices
+// SET_TYPE_PROP(NavierStokesNI, IsothermalIndices, NavierStokesCommonIndices<TypeTag>);
+
+//set isothermal NumEq
+// SET_INT_PROP(NavierStokesNI, IsothermalNumEq, 1);
+
+
+// \}
+} // end namespace Properties
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/freeflow/staggerednc/volumevariables.hh b/dumux/freeflow/staggerednc/volumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9d8515285c2c80c1f0b72eef05ea7ecbe8b923a7
--- /dev/null
+++ b/dumux/freeflow/staggerednc/volumevariables.hh
@@ -0,0 +1,252 @@
+// -*- 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 Quantities required by the one-phase fully implicit model defined on a vertex.
+ */
+#ifndef DUMUX_NAVIER_STOKES_NC_VOLUMEVARIABLES_HH
+#define DUMUX_NAVIER_STOKES_NC_VOLUMEVARIABLES_HH
+
+#include "properties.hh"
+#include <dumux/discretization/volumevariables.hh>
+
+#include <dumux/material/fluidstates/immiscible.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup NavierStokesModel
+ * \ingroup ImplicitVolumeVariables
+ * \brief Contains the quantities which are constant within a
+ * finite volume in the one-phase model.
+ */
+template <class TypeTag>
+class NavierStokesNCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
+{
+ using ParentType = ImplicitVolumeVariables<TypeTag>;
+ using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents)};
+
+public:
+
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+ /*!
+ * \copydoc ImplicitVolumeVariables::update
+ */
+ void update(const ElementSolutionVector &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_);
+ int compIIdx = 0;
+ for (unsigned int compJIdx = 0; compJIdx < numComponents; ++compJIdx)
+ {
+ // binary diffusion coefficents
+ if(compIIdx!= compJIdx)
+ {
+ setDiffusionCoefficient_(0, compJIdx,
+ FluidSystem::binaryDiffusionCoefficient(fluidState_,
+ paramCache,
+ 0,
+ compIIdx,
+ compJIdx));
+ }
+ }
+
+// std::cout << "mole frac 0 is : " << moleFraction(0,0) << std::endl;
+// std::cout << "mole frac 1 is : " << moleFraction(0,1) << std::endl;
+ if(moleFraction(0,0) + moleFraction(0,1) > 1.000001)
+ std::cout << "sum: " << moleFraction(0,0) + moleFraction(0,1) << std::endl;
+ };
+
+ /*!
+ * \copydoc ImplicitModel::completeFluidState
+ */
+ static void completeFluidState(const ElementSolutionVector& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const SubControlVolume& scv,
+ FluidState& fluidState)
+ {
+ Scalar t = ParentType::temperature(elemSol, problem, element, scv);
+ fluidState.setTemperature(t);
+ fluidState.setSaturation(/*phaseIdx=*/0, 1.);
+
+ fluidState.setPressure(/*phaseIdx=*/0, elemSol[0][Indices::pressureIdx]);
+
+ fluidState.setMassFraction(0, 1, elemSol[0][1]);
+ fluidState.setMassFraction(0, 0, 1.0 - elemSol[0][1]);
+
+
+
+ // 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(/*phaseIdx=*/0, 1.0);
+
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState, /*phaseIdx=*/0);
+
+ Scalar value = FluidSystem::density(fluidState, paramCache, /*phaseIdx=*/0);
+ fluidState.setDensity(/*phaseIdx=*/0, value);
+
+ value = FluidSystem::viscosity(fluidState, paramCache, /*phaseIdx=*/0);
+ fluidState.setViscosity(/*phaseIdx=*/0, value);
+
+ // compute and set the enthalpy
+ value = ParentType::enthalpy(fluidState, paramCache, /*phaseIdx=*/0);
+ fluidState.setEnthalpy(/*phaseIdx=*/0, value);
+ }
+
+ /*!
+ * \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 pressure \f$\mathrm{[Pa]}\f$ of a given phase within
+ * the control volume.
+ */
+ Scalar pressure(int phaseIdx = 0) const
+ { return fluidState_.pressure(phaseIdx); }
+
+ /*!
+ * \brief Return 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(phaseIdx); }
+
+ /*!
+ * \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(phaseIdx); }
+
+ /*!
+ * \brief Returns the mass fraction of a component in the phase
+ *
+ * \param phaseIdx the index of the fluid phase
+ * \param compIdx the index of the component
+ */
+ Scalar massFraction(int phaseIdx, int compIdx) const
+ { return fluidState_.massFraction(phaseIdx, compIdx); }
+
+ /*!
+ * \brief Returns the mole fraction of a component in the phase
+ *
+ * \param phaseIdx the index of the fluid phase
+ * \param compIdx the index of the component
+ */
+ Scalar moleFraction(int phaseIdx, int compIdx) const
+ { return fluidState_.moleFraction(phaseIdx, compIdx); }
+
+ /*!
+ * \brief Returns the mass density of a given phase within the
+ * control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar molarDensity(int phaseIdx = 0) const
+ {
+ if (phaseIdx < 1/*numPhases*/)
+ return fluidState_.molarDensity(phaseIdx);
+
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+ /*!
+ * \brief Returns the diffusion coeffiecient
+ */
+ Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
+ {
+ if (compIdx < phaseIdx)
+ return diffCoefficient_[phaseIdx][compIdx];
+ else if (compIdx > phaseIdx)
+ return diffCoefficient_[phaseIdx][compIdx-1];
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Diffusion coeffiecient called for phaseIdx = compIdx");
+ }
+
+ /*!
+ * \brief Return the fluid state of the control volume.
+ */
+ const FluidState &fluidState() const
+ { return fluidState_; }
+
+protected:
+ FluidState fluidState_;
+
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+
+ void setDiffusionCoefficient_(int phaseIdx, int compIdx, Scalar d)
+ {
+ if (compIdx < phaseIdx)
+ diffCoefficient_[phaseIdx][compIdx] = std::move(d);
+ else if (compIdx > phaseIdx)
+ diffCoefficient_[phaseIdx][compIdx-1] = std::move(d);
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Diffusion coeffiecient for phaseIdx = compIdx doesn't exist");
+ }
+
+
+ std::array<std::array<Scalar, numComponents-1>, 1> diffCoefficient_;
+};
+
+}
+
+#endif
diff --git a/dumux/implicit/staggered/localjacobian.hh b/dumux/implicit/staggered/localjacobian.hh
index d2da798d54d4adc7ae07e2eac54e7651a12d8760..739b28f9c5dc7f23e1ab9f54a350f3ec39252515 100644
--- a/dumux/implicit/staggered/localjacobian.hh
+++ b/dumux/implicit/staggered/localjacobian.hh
@@ -295,6 +295,7 @@ private:
// update the global jacobian matrix with the current partial derivatives
this->updateGlobalJacobian_(matrix[cellCenterIdx][cellCenterIdx], ccGlobalI_, globalJ, pvIdx, partialDeriv);
+// printmatrix(std::cout,matrix[cellCenterIdx][cellCenterIdx], "part", "");
// restore the original volVars
curVolVars = origVolVars;
}
diff --git a/test/freeflow/CMakeLists.txt b/test/freeflow/CMakeLists.txt
index f9eb646820fce8066b2aec164a1c436917fdd607..c0a5176e7a88498302dda88012b67fb87dd4d851 100644
--- a/test/freeflow/CMakeLists.txt
+++ b/test/freeflow/CMakeLists.txt
@@ -1,5 +1,6 @@
add_subdirectory("navierstokes")
add_subdirectory("staggered")
+add_subdirectory("staggerednc")
add_subdirectory("stokes")
add_subdirectory("stokes2c")
add_subdirectory("stokes2cni")
diff --git a/test/freeflow/staggerednc/CMakeLists.txt b/test/freeflow/staggerednc/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..9aa0dac746017a67fcab0732d53517f4d98213a7
--- /dev/null
+++ b/test/freeflow/staggerednc/CMakeLists.txt
@@ -0,0 +1,17 @@
+add_input_file_links()
+
+add_dumux_test(test_channel_stokes_nc test_channel_stokes_nc test_channel.cc
+ python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/channel-stokes.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_channel_stokes-00002.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_channel_stokes")
+
+
+#install sources
+install(FILES
+channeltestproblem.hh
+test_channel.cc
+closedsystemtestproblem.hh
+test_closedsystem.cc
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/freeflow/staggered)
diff --git a/test/freeflow/staggerednc/channeltestproblem.hh b/test/freeflow/staggerednc/channeltestproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e4b327ed14b979918cd09f877a39011ad0e09426
--- /dev/null
+++ b/test/freeflow/staggerednc/channeltestproblem.hh
@@ -0,0 +1,300 @@
+// -*- 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 Channel flow test for the staggered grid (Navier-)Stokes model
+ */
+#ifndef DUMUX_CHANNEL_NC_TEST_PROBLEM_HH
+#define DUMUX_CHANNEL_NC_TEST_PROBLEM_HH
+
+#include <dumux/implicit/staggered/properties.hh>
+#include <dumux/freeflow/staggerednc/model.hh>
+#include <dumux/implicit/problem.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
+#include <dumux/material/components/constant.hh>
+
+#include <dumux/material/fluidsystems/h2oair.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class ChannelTestProblem;
+
+namespace Capabilities
+{
+ template<class TypeTag>
+ struct isStationary<ChannelTestProblem<TypeTag>>
+ { static const bool value = false; };
+}
+
+namespace Properties
+{
+NEW_TYPE_TAG(ChannelTestProblem, INHERITS_FROM(StaggeredModel, NavierStokesNC));
+
+// SET_PROP(ChannelTestProblem, Fluid)
+// {
+// private:
+// typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+// public:
+// typedef FluidSystems::LiquidPhase<Scalar, Dumux::Constant<TypeTag, Scalar> > type;
+// };
+// Select the fluid system
+SET_TYPE_PROP(ChannelTestProblem, FluidSystem,
+ FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+
+// Set the grid type
+SET_TYPE_PROP(ChannelTestProblem, Grid, Dune::YaspGrid<2>);
+
+// Set the problem property
+SET_TYPE_PROP(ChannelTestProblem, Problem, Dumux::ChannelTestProblem<TypeTag> );
+
+SET_BOOL_PROP(ChannelTestProblem, EnableGlobalFVGeometryCache, true);
+
+SET_BOOL_PROP(ChannelTestProblem, EnableGlobalFluxVariablesCache, true);
+SET_BOOL_PROP(ChannelTestProblem, EnableGlobalVolumeVariablesCache, true);
+
+// Enable gravity
+SET_BOOL_PROP(ChannelTestProblem, ProblemEnableGravity, true);
+
+#if ENABLE_NAVIERSTOKES
+SET_BOOL_PROP(ChannelTestProblem, EnableInertiaTerms, true);
+#else
+SET_BOOL_PROP(ChannelTestProblem, EnableInertiaTerms, false);
+#endif
+}
+
+/*!
+ * \brief Test problem for the one-phase model:
+ \todo doc me!
+ */
+template <class TypeTag>
+class ChannelTestProblem : public NavierStokesProblem<TypeTag>
+{
+ typedef NavierStokesProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+ // copy some indices for convenience
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ enum {
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+ enum {
+ massBalanceIdx = Indices::massBalanceIdx,
+ momentumBalanceIdx = Indices::momentumBalanceIdx,
+ momentumXBalanceIdx = Indices::momentumXBalanceIdx,
+ momentumYBalanceIdx = Indices::momentumYBalanceIdx,
+ pressureIdx = Indices::pressureIdx,
+ velocityXIdx = Indices::velocityXIdx,
+ velocityYIdx = Indices::velocityYIdx
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
+
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+ using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+
+ using BoundaryValues = typename GET_PROP_TYPE(TypeTag, BoundaryValues);
+ using InitialValues = typename GET_PROP_TYPE(TypeTag, BoundaryValues);
+ using SourceValues = typename GET_PROP_TYPE(TypeTag, BoundaryValues);
+
+public:
+ ChannelTestProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView)
+ {
+ name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag,
+ std::string,
+ Problem,
+ Name);
+
+ inletVelocity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag,
+ Scalar,
+ Problem,
+ InletVelocity);
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ std::string name() const
+ {
+ return name_;
+ }
+
+ bool shouldWriteRestartFile() const
+ {
+ return false;
+ }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature() const
+ { return 273.15 + 10; } // 10C
+
+ /*!
+ * \brief Return the sources within the domain.
+ *
+ * \param globalPos The global position
+ */
+ SourceValues sourceAtPos(const GlobalPosition &globalPos) const
+ {
+ return SourceValues(0.0);
+ }
+ // \}
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary control volume.
+ *
+ * \param globalPos The position of the center of the finite volume
+ */
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+ {
+ BoundaryTypes values;
+
+ // set Dirichlet values for the velocity everywhere
+ values.setDirichlet(momentumBalanceIdx);
+ values.setNeumann(1);
+ // values.setDirichlet(1);
+
+ // set a fixed pressure in one cell
+ if (isOutlet(globalPos))
+ {
+ values.setDirichlet(massBalanceIdx);
+ values.setOutflow(momentumBalanceIdx);
+ values.setOutflow(1);
+ }
+ else
+ values.setNeumann(massBalanceIdx);
+
+ if(isInlet(globalPos))
+ values.setDirichlet(1);
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * control volume.
+ *
+ * \param globalPos The center of the finite volume which ought to be set.
+ */
+ BoundaryValues dirichletAtPos(const GlobalPosition &globalPos) const
+ {
+ BoundaryValues values;
+
+ values[pressureIdx] = 1.1e+5;
+
+ if(isInlet(globalPos))
+ {
+ values[velocityXIdx] = inletVelocity_;
+ values[velocityYIdx] = 0.0;
+ values[pressureIdx+1] =1e-6;
+ }
+ else if(isWall(globalPos))
+ {
+ values[velocityXIdx] = 0.0;
+ values[velocityYIdx] = 0.0;
+ }
+ else if(isOutlet(globalPos))
+ {
+ values[velocityXIdx] = 1.0;
+ values[velocityYIdx]= 0.0;
+ }
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ InitialValues initialAtPos(const GlobalPosition &globalPos) const
+ {
+ InitialValues values;
+ values[pressureIdx] = 1.1e+5;
+ values[velocityXIdx] = 0.0;
+ values[velocityYIdx] = 0.0;
+
+ return values;
+ }
+
+ // \}
+
+private:
+
+ bool isInlet(const GlobalPosition& globalPos) const
+ {
+ return globalPos[0] < 1e-6;
+ }
+
+ bool isOutlet(const GlobalPosition& globalPos) const
+ {
+ return globalPos[0] > 10 - 1e-6;//this->bBoxMax()[0] - 1e-6;
+ }
+
+ bool isWall(const GlobalPosition& globalPos) const
+ {
+ return globalPos[0] > 1e-6 || globalPos[0] < 10 - 1e-6;//this->bBoxMax()[0] - 1e-6;
+ }
+
+ static constexpr Scalar eps_{1e-6}; // TODO: what is wrong with bBoxMax and eps_ ???
+ Scalar inletVelocity_;
+ std::string name_;
+};
+} //end namespace
+
+#endif
diff --git a/test/freeflow/staggerednc/test_channel.cc b/test/freeflow/staggerednc/test_channel.cc
new file mode 100644
index 0000000000000000000000000000000000000000..9e6fef38637fbc07d1e08565410264afab838d7a
--- /dev/null
+++ b/test/freeflow/staggerednc/test_channel.cc
@@ -0,0 +1,64 @@
+// -*- 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 Channel flow test for the staggered grid (Navier-)Stokes model
+ */
+#include <config.h>
+#include "channeltestproblem.hh"
+#include <dumux/common/start.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
+ *
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n"
+ "\t-TimeManager.TEnd End of the simulation [s] \n"
+ "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+ "\t-Grid.File Name of the file containing the grid \n"
+ "\t definition in DGF format\n"
+ "\t-SpatialParams.LensLowerLeftX x-coordinate of the lower left corner of the lens [m] \n"
+ "\t-SpatialParams.LensLowerLeftY y-coordinate of the lower left corner of the lens [m] \n"
+ "\t-SpatialParams.LensUpperRightX x-coordinate of the upper right corner of the lens [m] \n"
+ "\t-SpatialParams.LensUpperRightY y-coordinate of the upper right corner of the lens [m] \n"
+ "\t-SpatialParams.Permeability Permeability of the domain [m^2] \n"
+ "\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(ChannelTestProblem) ProblemTypeTag;
+ return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+}
diff --git a/test/freeflow/staggerednc/test_channel_stokes.input b/test/freeflow/staggerednc/test_channel_stokes.input
new file mode 100644
index 0000000000000000000000000000000000000000..fc4f25ebf94e691d10f8117eb725e146f286ffd0
--- /dev/null
+++ b/test/freeflow/staggerednc/test_channel_stokes.input
@@ -0,0 +1,20 @@
+[TimeManager]
+DtInitial = 1e-3 # [s]
+TEnd = 2 # [s]
+
+[Grid]
+UpperRight = 10 1
+Cells = 10 10
+
+[Problem]
+Name = test_channel_stokes # name passed to the output routines
+InletVelocity = 1
+LiquidDensity = 1
+EnableGravity = false
+
+[ Newton ]
+MaxSteps = 10
+MaxRelativeShift = 1e-5
+
+[Vtk]
+WriteFaceData = false