diff --git a/dumux/freeflow/navierstokes/staggered/fluxoversurface.hh b/dumux/freeflow/navierstokes/staggered/fluxoversurface.hh
index d12d0ceb4efe4e5f5348906350d0abc9a5a65d1e..5191502f46c45a7058b63e393a589a8123230c86 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxoversurface.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxoversurface.hh
@@ -227,7 +227,11 @@ public:
const auto& elemFluxVarsCache)
{
LocalResidual localResidual(&problem_());
- return localResidual.computeFluxForCellCenter(problem_(), element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
+
+ if (scvf.boundary())
+ return localResidual.computeBoundaryFluxForCellCenter(problem_(), element, fvGeometry, scvf, elemVolVars, elemFaceVars, /*elemBcTypes*/{}, elemFluxVarsCache);
+ else
+ return localResidual.computeFluxForCellCenter(problem_(), element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
};
calculateFluxes(fluxType);
@@ -238,8 +242,19 @@ public:
*/
void calculateVolumeFluxes()
{
- const auto isCompositional = std::integral_constant<bool, (ModelTraits::numFluidComponents() > 1) >();
- calculateVolumeFluxesImpl_(isCompositional);
+ auto fluxType = [this](const auto& element,
+ const auto& fvGeometry,
+ const auto& elemVolVars,
+ const auto& elemFaceVars,
+ const auto& scvf,
+ const auto& elemFluxVarsCache)
+ {
+ CellCenterPrimaryVariables result(0.0);
+ result[0] = elemFaceVars[scvf].velocitySelf() * scvf.area() * scvf.directionSign();
+ return result;
+ };
+
+ calculateFluxes(fluxType);
}
/*!
@@ -339,88 +354,6 @@ private:
const auto& problem_() const { return gridVariables_.curGridVolVars().problem(); }
- /*!
- * \brief Calculate the volume fluxes over all surfaces for compositional models.
- * This method simply averages the densities between two adjacent cells.
- */
- void calculateVolumeFluxesImpl_(std::true_type)
- {
- auto fluxType = [this](const auto& element,
- const auto& fvGeometry,
- const auto& elemVolVars,
- const auto& elemFaceVars,
- const auto& scvf,
- const auto& elemFluxVarsCache)
- {
- LocalResidual localResidual(&problem_());
- const auto massOrMoleFlux = localResidual.computeFluxForCellCenter(problem_(), element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
-
- const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
- const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
-
- constexpr bool useMoles = false;//getPropValue<TypeTag, Properties::UseMoles>();
- const auto density = [useMoles](const auto& volVars)
- {
- return useMoles ? volVars.molarDensity() : volVars.density() ;
- };
-
- const auto avgDensity = 0.5*density(insideVolVars) + 0.5*density(outsideVolVars);
-
- constexpr auto replaceCompEqIdx = ModelTraits::ReplaceCompEqIdx();
- constexpr auto numComponents = ModelTraits::numFluidComponents();
-
- const Scalar cumulativeFlux = [replaceCompEqIdx, numComponents, &massOrMoleFlux]()
- {
- Scalar result = 0.0;
- if(replaceCompEqIdx < numComponents)
- result = massOrMoleFlux[replaceCompEqIdx];
- else
- {
- for(int i = 0; i < numComponents; ++i)
- result += massOrMoleFlux[i];
- }
-
- return result;
- }();
-
- CellCenterPrimaryVariables tmp(0.0);
- tmp[0] = cumulativeFlux / avgDensity;
- return tmp;
- };
-
- calculateFluxes(fluxType);
- }
-
- /*!
- * \brief Calculate the volume fluxes over all surfaces for non-compositional models.
- * This method simply averages the densities between two adjacent cells.
- */
- void calculateVolumeFluxesImpl_(std::false_type)
- {
- auto fluxType = [this](const auto& element,
- const auto& fvGeometry,
- const auto& elemVolVars,
- const auto& elemFaceVars,
- const auto& scvf,
- const auto& elemFluxVarsCache)
- {
- LocalResidual localResidual(&problem_());
- const Scalar totalMassFlux = localResidual.computeFluxForCellCenter(problem_(), element, fvGeometry, elemVolVars,
- elemFaceVars, scvf, elemFluxVarsCache)[0];
-
- const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
- const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
-
- const auto avgDensity = 0.5*insideVolVars.density() + 0.5*outsideVolVars.density();
-
- CellCenterPrimaryVariables tmp(0.0);
- tmp[0] = totalMassFlux / avgDensity;
- return tmp;
- };
-
- calculateFluxes(fluxType);
- }
-
std::map<std::string, SurfaceData<surfaceDim ,dim> > surfaces_;
const GridVariables& gridVariables_;
const SolutionVector& sol_;
diff --git a/dumux/freeflow/navierstokes/staggered/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh
index a786cfc1675f5bcf068aae9057c052fff1110666..3ddde3953667bc5e8706645799d5d8afb4948095 100644
--- a/dumux/freeflow/navierstokes/staggered/localresidual.hh
+++ b/dumux/freeflow/navierstokes/staggered/localresidual.hh
@@ -68,6 +68,7 @@ class NavierStokesResidualImpl<TypeTag, DiscretizationMethod::staggered>
using Element = typename GridView::template Codim<0>::Entity;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>;
using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>;
using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>;
@@ -184,7 +185,6 @@ public:
/*!
* \brief Evaluate boundary conditions for a cell center dof
*/
- template<class ElementBoundaryTypes>
CellCenterResidual computeBoundaryFluxForCellCenter(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -232,7 +232,6 @@ public:
/*!
* \brief Evaluate Dirichlet (fixed value) boundary conditions for a face dof
*/
- template<class ElementBoundaryTypes>
void evalDirichletBoundariesForFace(FaceResidual& residual,
const Problem& problem,
const Element& element,
@@ -269,7 +268,6 @@ public:
/*!
* \brief Evaluate boundary boundary fluxes for a face dof
*/
- template<class ElementBoundaryTypes>
FaceResidual computeBoundaryFluxForFace(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc
index dfa5ea5f5c4962f351064b8e0bd246f9b22ea579..bdc19fa31f600a45167bd8c1cd7868203e40a551 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc
@@ -50,6 +50,7 @@
#include <dumux/multidomain/newtonsolver.hh>
#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
+#include <dumux/freeflow/navierstokes/staggered/fluxoversurface.hh>
#include "problem_darcy.hh"
#include "problem_stokes.hh"
@@ -194,6 +195,29 @@ int main(int argc, char** argv) try
couplingManager,
timeLoop);
+ FluxOverSurface<StokesGridVariables,
+ decltype(stokesSol),
+ GetPropType<StokesTypeTag, Properties::ModelTraits>,
+ GetPropType<StokesTypeTag, Properties::LocalResidual>> flux(*stokesGridVariables, stokesSol);
+
+ using StokesGridView = GetPropType<StokesTypeTag, Properties::GridView>;
+ using StokesElement = typename StokesGridView::template Codim<0>::Entity;
+
+ using GlobalPosition = typename StokesElement::Geometry::GlobalCoordinate;
+
+ const Scalar xMin = stokesFvGridGeometry->bBoxMin()[0];
+ const Scalar xMax = stokesFvGridGeometry->bBoxMax()[0];
+ const Scalar yMin = stokesFvGridGeometry->bBoxMin()[1];
+ const Scalar yMax = stokesFvGridGeometry->bBoxMax()[1];
+
+ const auto p0inlet = GlobalPosition{xMin, yMin};
+ const auto p1inlet = GlobalPosition{xMin, yMax};
+ flux.addSurface("inlet", p0inlet, p1inlet);
+
+ const auto p0outlet = GlobalPosition{xMax, yMin};
+ const auto p1outlet = GlobalPosition{xMax, yMax};
+ flux.addSurface("outlet", p0outlet, p1outlet);
+
// the linear solver
using LinearSolver = UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>();
@@ -224,6 +248,16 @@ int main(int argc, char** argv) try
stokesVtkWriter.write(timeLoop->time());
darcyVtkWriter.write(timeLoop->time());
+ flux.calculateMassOrMoleFluxes();
+
+ std::cout << "mole flux at inlet is: " << flux.netFlux("inlet") << std::endl;
+ std::cout << "mole flux at outlet is: " << flux.netFlux("outlet") << std::endl;
+
+ // calculate and print volume fluxes over the planes
+ flux.calculateVolumeFluxes();
+ std::cout << "volume flux at inlet is: " << flux.netFlux("inlet")[0] << std::endl;
+ std::cout << "volume flux at outlet is: " << flux.netFlux("outlet")[0] << std::endl;
+
// report statistics of this time step
timeLoop->reportTimeStep();