diff --git a/dumux/discretization/facecentered/staggered/fvelementgeometry.hh b/dumux/discretization/facecentered/staggered/fvelementgeometry.hh
index 58d008df61c90e894a5d3a32b1080a7a9569912d..357cdd273da04a3cee8bd98277a3ce197007040b 100644
--- a/dumux/discretization/facecentered/staggered/fvelementgeometry.hh
+++ b/dumux/discretization/facecentered/staggered/fvelementgeometry.hh
@@ -197,6 +197,14 @@ public:
const Element& element() const
{ return *elementPtr_; }
+ //! Returns true if the IP of an scvf lies on a concave corner
+ bool scvfIntegrationPointInConcaveCorner(const SubControlVolumeFace& scvf) const
+ { return GG::GeometryHelper::scvfIntegrationPointInConcaveCorner(*this, scvf); }
+
+ //! Returns the the scvf of neighbor element with the same integration point and unit outer normal
+ const SubControlVolumeFace& outsideScvfWithSameIntegrationPoint(const SubControlVolumeFace& scvf) const
+ { return GG::GeometryHelper::outsideScvfWithSameIntegrationPoint(*this, scvf); }
+
private:
const auto& scvfIndices_() const
@@ -249,7 +257,7 @@ public:
, geometryHelper_(gridGeometry.gridView())
{}
- //! Get a sub control volume face with a local scv index
+ //! Get a sub control volume face with a global scvf index
const SubControlVolumeFace& scvf(const GridIndexType scvfIdx) const
{ return scvfs_[findLocalIndex_(scvfIdx, scvfIndices_())]; }
@@ -385,6 +393,14 @@ public:
return *gridGeometry_;
}
+ //! Returns true if the IP of an scvf lies on a concave corner
+ bool scvfIntegrationPointInConcaveCorner(const SubControlVolumeFace& scvf) const
+ { return GG::GeometryHelper::scvfIntegrationPointInConcaveCorner(*this, scvf); }
+
+ //! Returns the the scvf of neighbor element with the same integration point and unit outer normal
+ const SubControlVolumeFace& outsideScvfWithSameIntegrationPoint(const SubControlVolumeFace& scvf) const
+ { return GG::GeometryHelper::outsideScvfWithSameIntegrationPoint(*this, scvf); }
+
private:
//! Binding of an element preparing the geometries of the whole stencil
//! called by the local jacobian to prepare element assembly
diff --git a/dumux/discretization/facecentered/staggered/fvgridgeometry.hh b/dumux/discretization/facecentered/staggered/fvgridgeometry.hh
index 7158921f2f9422e6dd2374b5962dab42fa882074..5830e63d591f15ed4db3dc6f39d0c54ce8924883 100644
--- a/dumux/discretization/facecentered/staggered/fvgridgeometry.hh
+++ b/dumux/discretization/facecentered/staggered/fvgridgeometry.hh
@@ -61,7 +61,7 @@ struct FaceCenteredStaggeredDefaultGridGeometryTraits : public DefaultMapperTrai
using SubControlVolumeFace = FaceCenteredStaggeredSubControlVolumeFace<GridView>;
using IntersectionMapper = ConformingGridIntersectionMapper<GridView>;
using LocalIntersectionMapper = FaceCenteredStaggeredLocalIntersectionIndexMapper<GridView>;
- using GeometryHelper = FaceCenteredStaggeredGeometryHelper<GridView, typename GridView::Grid>;
+ using GeometryHelper = FaceCenteredStaggeredGeometryHelper<GridView>;
template<class GridGeometry>
using ConnectivityMap = FaceCenteredStaggeredConnectivityMap<GridGeometry>;
diff --git a/dumux/discretization/facecentered/staggered/geometryhelper.hh b/dumux/discretization/facecentered/staggered/geometryhelper.hh
index 80bb2433f5c055adcdfe8065dd93267c12aeb127..73013f44f2b392f1b1c2d4bae3e93aec627c8b8d 100644
--- a/dumux/discretization/facecentered/staggered/geometryhelper.hh
+++ b/dumux/discretization/facecentered/staggered/geometryhelper.hh
@@ -24,25 +24,15 @@
#ifndef DUMUX_DISCRETIZATION_FACECENTERED_STAGGERED_GEOMETRY_HELPER_HH
#define DUMUX_DISCRETIZATION_FACECENTERED_STAGGERED_GEOMETRY_HELPER_HH
-#include <dune/common/float_cmp.hh>
+#include <array>
+#include <cassert>
#include <dumux/common/indextraits.hh>
-#include <dune/common/reservedvector.hh>
-#include <dumux/common/math.hh>
-#include <dumux/discretization/basegridgeometry.hh>
-#include <dumux/discretization/checkoverlapsize.hh>
-#include <dumux/discretization/method.hh>
-
-namespace Dune {
-
-template<int dim, class Coordinates>
-class YaspGrid;
-
-}
+#include <dumux/discretization/facecentered/staggered/gridsupportsconcavecorners.hh>
namespace Dumux {
-template<class GridView, class Implementation>
-class FaceCenteredStaggeredGeometryHelperBase
+template<class GridView>
+class FaceCenteredStaggeredGeometryHelper
{
using GridIndexType = typename IndexTraits<GridView>::GridIndex;
using SmallLocalIndexType = typename IndexTraits<GridView>::SmallLocalIndex;
@@ -54,15 +44,15 @@ public:
static constexpr auto numElementFaces = dim * 2;
static constexpr auto numLateralFacesPerScv = 2 * (dim - 1);
- FaceCenteredStaggeredGeometryHelperBase(const GridView& gridView) : gridView_(gridView) {}
+ FaceCenteredStaggeredGeometryHelper(const GridView& gridView) : gridView_(gridView) {}
//! Returns the local index of the opposing face.
static constexpr SmallLocalIndexType localOppositeIdx(const SmallLocalIndexType ownLocalFaceIndex)
{
- return (ownLocalFaceIndex % 2) ? (ownLocalFaceIndex - 1) : (ownLocalFaceIndex + 1);
+ return isOdd_(ownLocalFaceIndex) ? (ownLocalFaceIndex - 1) : (ownLocalFaceIndex + 1);
}
- //! Return the local index of a lateral orthogonal scvf
+ //! Return the local index of a lateral orthogonal scvf
static constexpr int lateralOrthogonalScvfLocalIndex(const SmallLocalIndexType ownLocalScvfIndex)
{
if constexpr(GridView::Grid::dimension == 1)
@@ -128,7 +118,6 @@ public:
}
}
-
//! Returns the local indices of the faces lateral to the own one.
static constexpr auto localLaterFaceIndices(const SmallLocalIndexType ownLocalFaceIndex)
{
@@ -163,144 +152,60 @@ public:
DUNE_THROW(Dune::InvalidStateException, "localFacetIdx " << localFacetIdx << " out of range");
}
- //! Map two local facet indices to a local common entity index (vertex in 2D or edge in 3D)
- SmallLocalIndexType localCommonEntityIndex(SmallLocalIndexType localFacetIdx0, SmallLocalIndexType localFacetIdx1) const
+ //! Returns true if the IP of an scvf lies on a concave corner
+ template<class FVElementGeometry, class SubControlVolumeFace>
+ static bool scvfIntegrationPointInConcaveCorner(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf)
{
- constexpr auto table = []
- {
- using Table = std::array<std::array<SmallLocalIndexType, 2>, dim == 2 ? 4 : 12>;
- if constexpr (dim == 2)
- return Table {{ {0,2}, {1,2}, {0,3}, {1,3} }};
- else
- return Table {{ {0,2}, {1,2}, {0,3}, {1,3}, {0,4}, {1,4}, {2,4}, {3,4}, {0,5}, {1,5}, {2,5}, {3,5} }};
- }();
-
- using std::swap;
- if (localFacetIdx0 > localFacetIdx1)
- swap(localFacetIdx0, localFacetIdx1);
-
- const auto idx = std::find(table.begin(), table.end(), std::array<SmallLocalIndexType, 2>{localFacetIdx0, localFacetIdx1});
+ assert (scvf.isLateral());
- if (idx != std::end(table))
- return std::distance(table.begin(), idx);
+ using Grid = typename FVElementGeometry::GridGeometry::Grid;
+ if constexpr (!GridSupportsConcaveCorners<Grid>::value)
+ return false;
else
- DUNE_THROW(Dune::InvalidStateException, "Faces with local indices " << localFacetIdx0 << " and " << localFacetIdx1 << " do not intersect");
- }
-
- //! Map two facets to a global common entity index (vertex in 2D or edge in 3D)
- auto globalCommonEntityIndex(const Element& element, SmallLocalIndexType localFacetIdx0, SmallLocalIndexType localFacetIdx1) const
- {
- const auto refElement = referenceElement(element);
-
- assert(refElement.size(localFacetIdx0, 1, 2) == refElement.size(localFacetIdx1, 1, 2));
- std::array<GridIndexType, dim == 2 ? 2 : 4> cache = {};
- std::bitset<dim == 2 ? 2 : 4> valueCached = {};
-
- for (int i = 0; i < refElement.size(localFacetIdx0, 1, 2); ++i)
{
- // get local sub index of subentity w.r.t. to the element
- const auto localEntityIdx0 = refElement.subEntity(localFacetIdx0, 1, i, 2);
- // and with the local index we can get the global index of the vertex/edge (2d/3d)
- const auto globalEntityIdx0 = gridView().indexSet().subIndex(element, localEntityIdx0, 2);
- for (int j = 0; j < refElement.size(localFacetIdx1, 1, 2); ++j)
- {
- const auto globalEntityIdx1 = valueCached[j] ? cache[j]
- : gridView().indexSet().subIndex(element, refElement.subEntity(localFacetIdx1, 1, j, 2), 2);
-
- if (globalEntityIdx0 == globalEntityIdx1)
- return globalEntityIdx0;
- else
- {
- cache[j] = globalEntityIdx1;
- valueCached.set(j);
- }
- }
+ if (scvf.boundary())
+ return false;
+
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
+ int onBoundaryCounter = 0;
+ onBoundaryCounter += static_cast<int>(insideScv.boundary());
+ onBoundaryCounter += static_cast<int>(outsideScv.boundary());
+ return onBoundaryCounter == 1;
}
- DUNE_THROW(Dune::InvalidStateException, "No common entity found");
}
- const GridView& gridView() const
- { return gridView_; }
-
-protected:
-
- Implementation &asImp()
- { return *static_cast<Implementation*>(this); }
-
- const Implementation &asImp() const
- { return *static_cast<const Implementation*>(this); }
-
-private:
- GridView gridView_;
-};
+ template<class FVElementGeometry, class SubControlVolumeFace>
+ static const SubControlVolumeFace& outsideScvfWithSameIntegrationPoint(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf)
+ {
+ const auto& lateralOrthogonalScvf = fvGeometry.lateralOrthogonalScvf(scvf);
+ assert(!lateralOrthogonalScvf.boundary());
+ const int offset = (dim == 2) ? 3 : 5;
+ const auto otherLocalIdx = isOdd_(scvf.localIndex()) ? scvf.localIndex() - offset : scvf.localIndex() + offset;
-template<class GridView, class Grid = typename GridView::Grid>
-class FaceCenteredStaggeredGeometryHelper
- : public FaceCenteredStaggeredGeometryHelperBase<GridView, FaceCenteredStaggeredGeometryHelper<GridView, Grid>>
-{
- using ParentType = FaceCenteredStaggeredGeometryHelperBase<GridView, FaceCenteredStaggeredGeometryHelper<GridView, Grid>>;
- using SmallLocalIndexType = typename IndexTraits<GridView>::SmallLocalIndex;
- using Element = typename GridView::template Codim<0>::Entity;
-public:
-
- using ParentType::ParentType;
+ auto outsideFVGeometry = localView(fvGeometry.gridGeometry());
+ const auto outsideElementIdx = fvGeometry.scv(lateralOrthogonalScvf.outsideScvIdx()).elementIndex();
+ outsideFVGeometry.bindElement(fvGeometry.gridGeometry().element(outsideElementIdx));
- //! Update the map for getting the corresponding local face indices in another element.
- void update(const Element& ownElement, const Element& otherElement)
- {
- SmallLocalIndexType ownIdx = 0;
- SmallLocalIndexType otherIdx = 0;
- for (const auto& ownIs : intersections(this->gridView(), ownElement))
+ for (const auto& otherScvf : scvfs(outsideFVGeometry))
{
- // helper lambda to make sure the inner loop stops once the other index is found
- const auto computeOtherIdx = [&] (const auto& ownOuterNormal)
- {
- for (const auto& otherIs : intersections(this->gridView(), otherElement))
- {
- const auto& otherOuterNormal = otherIs.centerUnitOuterNormal();
- if (Dune::FloatCmp::eq<typename GridView::ctype>(ownOuterNormal*otherOuterNormal, 1.0))
- {
- map_[ownIdx++] = otherIdx++;
- return;
- }
- }
- };
-
- computeOtherIdx(ownIs.centerUnitOuterNormal());
+ if (otherScvf.localIndex() == otherLocalIdx)
+ return otherScvf;
}
- if (ownIdx != map_.size())
- DUNE_THROW(Dune::InvalidStateException, "Index could not be mapped");
+
+ DUNE_THROW(Dune::InvalidStateException, "No outside scvf found");
}
- //! Return the local index of the other element's facet with the same position as the own element's facet.
- SmallLocalIndexType localFaceIndexInOtherElement(const SmallLocalIndexType localFaceIndexInOwnElement) const
- { return map_[localFaceIndexInOwnElement]; }
+ const GridView& gridView() const
+ { return gridView_; }
private:
- std::array<SmallLocalIndexType, ParentType::numElementFaces> map_;
-};
-template<class GridView, class YaspCoordinates>
-class FaceCenteredStaggeredGeometryHelper<GridView, Dune::YaspGrid<GridView::Grid::dimension, YaspCoordinates>>
- : public FaceCenteredStaggeredGeometryHelperBase<GridView, FaceCenteredStaggeredGeometryHelper<GridView, typename GridView::Grid>>
-{
- using ParentType = FaceCenteredStaggeredGeometryHelperBase<GridView, FaceCenteredStaggeredGeometryHelper<GridView, typename GridView::Grid>>;
- using SmallLocalIndexType = typename IndexTraits<GridView>::SmallLocalIndex;
- using Element = typename GridView::template Codim<0>::Entity;
-public:
+ static constexpr bool isOdd_(int number)
+ { return number % 2; }
- using ParentType::ParentType;
-
- //! Update the map for getting the corresponding local face indices in another element.
- //! Nothing needs to be done here.
- void update(const Element& ownElement, const Element& otherElement)
- {}
-
- //! Return the local index of the other element's facet with the same position as the own element's facet.
- //! For Yasp grids, this is just the same index.
- SmallLocalIndexType localFaceIndexInOtherElement(const SmallLocalIndexType localFaceIndexInOwnElement) const
- { return localFaceIndexInOwnElement; }
+ GridView gridView_;
};
} // end namespace Dumux
diff --git a/dumux/discretization/facecentered/staggered/gridsupportsconcavecorners.hh b/dumux/discretization/facecentered/staggered/gridsupportsconcavecorners.hh
new file mode 100644
index 0000000000000000000000000000000000000000..60584593dbdee4cfccd2d1df66767ff45d99babf
--- /dev/null
+++ b/dumux/discretization/facecentered/staggered/gridsupportsconcavecorners.hh
@@ -0,0 +1,50 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Typetraits
+ * \copydoc Type trait to determine if a grid supports concave corners (e.g. by cutting out a hole from the domain interior)
+ */
+#ifndef DUMUX_DISCRETIZATION_FACECENTERED_STAGGERED_GRID_SUPPORTS_CONCAVE_CORNERS_HH
+#define DUMUX_DISCRETIZATION_FACECENTERED_STAGGERED_GRID_SUPPORTS_CONCAVE_CORNERS_HH
+
+#include <type_traits>
+
+// forward declare
+namespace Dune {
+
+template<int dim, class Coordinates>
+class YaspGrid;
+
+}
+
+namespace Dumux {
+
+/*!
+ * \brief Type trait to determine if a grid supports concave corners (e.g. by cutting out a hole from the domain interior)
+ */
+template<class T>
+struct GridSupportsConcaveCorners : public std::true_type {};
+
+template<int dim, class Coords>
+struct GridSupportsConcaveCorners<Dune::YaspGrid<dim, Coords>> : public std::false_type {};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/freeflow/navierstokes/momentum/fluxhelper.hh b/dumux/freeflow/navierstokes/momentum/fluxhelper.hh
index 7a023de98d7d18dfddf5a5d54747345108aec2ff..2bbc646ca02d512600aefa1f3e5e736909c0ccc3 100644
--- a/dumux/freeflow/navierstokes/momentum/fluxhelper.hh
+++ b/dumux/freeflow/navierstokes/momentum/fluxhelper.hh
@@ -127,28 +127,39 @@ struct NavierStokesMomentumBoundaryFluxHelper
const auto& orthogonalScvf = fvGeometry.lateralOrthogonalScvf(scvf);
const auto innerTransportingVelocity = elemVolVars[orthogonalScvf.insideScvIdx()].velocity();
- if (scvf.boundary())
+ static const bool useOldScheme = getParam<bool>("FreeFlow.UseOldTransportingVelocity", true); // TODO how to deprecate?
+
+ if (useOldScheme)
{
- if (const auto bcTypes = problem.boundaryTypes(element, scvf); bcTypes.isDirichlet(scvf.normalAxis()))
- return problem.dirichlet(element, scvf)[scvf.normalAxis()];
+ if (scvf.boundary())
+ {
+ if (problem.boundaryTypes(element, scvf).isDirichlet(scvf.normalAxis()))
+ return problem.dirichlet(element, scvf)[scvf.normalAxis()];
+ }
else
- return
- innerTransportingVelocity; // fallback
+ return innerTransportingVelocity;
}
- else
+
+ // use the Dirichlet velocity as transporting velocity if the lateral face is on a Dirichlet boundary
+ if (scvf.boundary())
{
- static const bool useOldScheme = getParam<bool>("FreeFlow.UseOldTransportingVelocity", true); // TODO how to deprecate?
- if (useOldScheme)
- return innerTransportingVelocity;
+ if (problem.boundaryTypes(element, scvf).isDirichlet(scvf.normalAxis()))
+ return 0.5*(problem.dirichlet(element, scvf)[scvf.normalAxis()] + innerTransportingVelocity);
+ }
+
+ if (orthogonalScvf.boundary())
+ {
+ if (problem.boundaryTypes(element, orthogonalScvf).isDirichlet(scvf.normalAxis()))
+ return 0.5*(problem.dirichlet(element, scvf)[scvf.normalAxis()] + innerTransportingVelocity);
else
- {
- // average the transporting velocity by weighting with the scv volumes
- const auto insideVolume = fvGeometry.scv(orthogonalScvf.insideScvIdx()).volume();
- const auto outsideVolume = fvGeometry.scv(orthogonalScvf.outsideScvIdx()).volume();
- const auto outerTransportingVelocity = elemVolVars[orthogonalScvf.outsideScvIdx()].velocity();
- return (insideVolume*innerTransportingVelocity + outsideVolume*outerTransportingVelocity) / (insideVolume + outsideVolume);
- }
+ return innerTransportingVelocity; // fallback value, should actually never be called
}
+
+ // average the transporting velocity by weighting with the scv volumes
+ const auto insideVolume = fvGeometry.scv(orthogonalScvf.insideScvIdx()).volume();
+ const auto outsideVolume = fvGeometry.scv(orthogonalScvf.outsideScvIdx()).volume();
+ const auto outerTransportingVelocity = elemVolVars[orthogonalScvf.outsideScvIdx()].velocity();
+ return (insideVolume*innerTransportingVelocity + outsideVolume*outerTransportingVelocity) / (insideVolume + outsideVolume);
}();
// lateral face normal to boundary (integration point touches boundary)
diff --git a/dumux/freeflow/navierstokes/momentum/fluxvariables.hh b/dumux/freeflow/navierstokes/momentum/fluxvariables.hh
index 765620d2fe55aee919111fe35cfd917e865e48b3..590d1380996331a42d0a28dac66846251acd8ad1 100644
--- a/dumux/freeflow/navierstokes/momentum/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/momentum/fluxvariables.hh
@@ -391,25 +391,33 @@ public:
// get the transporting velocity which is perpendicular to our own (inner) velocity
const Scalar transportingVelocity = [&]()
{
+ const auto& orthogonalScvf = fvGeometry.lateralOrthogonalScvf(scvf);
+ const Scalar innerTransportingVelocity = elemVolVars[orthogonalScvf.insideScvIdx()].velocity();
+
static const bool useOldScheme = getParam<bool>("FreeFlow.UseOldTransportingVelocity", true); // TODO how to deprecate?
+ if (useOldScheme)
+ {
+ if (scvf.boundary() && fvGeometry.scv(scvf.insideScvIdx()).boundary())
+ {
+ if (this->elemBcTypes()[scvf.localIndex()].isDirichlet(scvf.normalAxis()))
+ return problem.dirichlet(this->element(), scvf)[scvf.normalAxis()];
+ }
+ else
+ return innerTransportingVelocity;
+ }
+
// use the Dirichlet velocity as transporting velocity if the lateral face is on a Dirichlet boundary
- if (!useOldScheme && scvf.boundary())
+ if (scvf.boundary())
{
if (this->elemBcTypes()[scvf.localIndex()].isDirichlet(scvf.normalAxis()))
- return problem.dirichlet(this->element(), scvf)[scvf.normalAxis()];
+ return 0.5*(problem.dirichlet(this->element(), scvf)[scvf.normalAxis()] + innerTransportingVelocity);
}
- const auto& orthogonalScvf = fvGeometry.lateralOrthogonalScvf(scvf);
- const Scalar innerTransportingVelocity = elemVolVars[orthogonalScvf.insideScvIdx()].velocity();
-
- if (useOldScheme)
- return innerTransportingVelocity;
-
if (orthogonalScvf.boundary())
{
if (this->elemBcTypes()[orthogonalScvf.localIndex()].isDirichlet(scvf.normalAxis()))
- return problem.dirichlet(this->element(), scvf)[scvf.normalAxis()];
+ return 0.5*(problem.dirichlet(this->element(), scvf)[scvf.normalAxis()] + innerTransportingVelocity);
else
return innerTransportingVelocity; // fallback value, should actually never be called
}
@@ -423,13 +431,32 @@ public:
const Scalar transportedMomentum = [&]()
{
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+
+ auto getDirichletMomentumFlux = [&]()
+ {
+ return problem.dirichlet(this->element(), scvf)[insideScv.dofAxis()] * this->problem().density(this->element(), insideScv);
+ };
+
// use the Dirichlet velocity as for transported momentum if the lateral face is on a Dirichlet boundary
if (scvf.boundary())
{
- if (const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); this->elemBcTypes()[scvf.localIndex()].isDirichlet(scv.dofAxis()))
- return problem.dirichlet(this->element(), scvf)[scv.dofAxis()] * this->problem().density(this->element(), scv);
- else
+ if (!this->elemBcTypes()[scvf.localIndex()].isDirichlet(insideScv.dofAxis()))
DUNE_THROW(Dune::InvalidStateException, "Neither Dirichlet nor Neumann BC set at " << scvf.ipGlobal());
+
+ return getDirichletMomentumFlux();
+ }
+ else
+ {
+ if (fvGeometry.scvfIntegrationPointInConcaveCorner(scvf))
+ {
+ // TODO: we could put this into an assert, as the construction of outsideScvfWithSameIntegrationPoint is quite expensive
+ const auto& outsideScvfWithSameIntegrationPoint = fvGeometry.outsideScvfWithSameIntegrationPoint(scvf);
+ if (!this->problem().boundaryTypes(this->element(), outsideScvfWithSameIntegrationPoint).isDirichlet(insideScv.dofAxis()))
+ DUNE_THROW(Dune::InvalidStateException, "Neither Dirichlet nor Neumann BC set at " << outsideScvfWithSameIntegrationPoint.ipGlobal());
+
+ return getDirichletMomentumFlux();
+ }
}
const bool selfIsUpstream = scvf.directionSign() == sign(transportingVelocity);
diff --git a/test/freeflow/navierstokes/channel/3d/CMakeLists.txt b/test/freeflow/navierstokes/channel/3d/CMakeLists.txt
index 11356d3181361dc964ea2ff5cdb069fedf48ba23..97f0fcd0f6f79580d3913cc0ec83fcbbfb9296c9 100644
--- a/test/freeflow/navierstokes/channel/3d/CMakeLists.txt
+++ b/test/freeflow/navierstokes/channel/3d/CMakeLists.txt
@@ -25,6 +25,19 @@ dumux_add_test(NAME test_ff_stokes_channel_3d_staircase
-Problem.Name test_ff_stokes_channel_3d_staircase -Problem.IsStaircaseGeometry true"
--zeroThreshold {"velocity_liq \(m/s\)":1e-12})
+dumux_add_test(NAME test_ff_stokes_channel_3d_staircase_nocaching
+ SOURCES main.cc
+ LABELS freeflow navierstokes
+ CMAKE_GUARD HAVE_UMFPACK dune-subgrid_FOUND
+ COMPILE_DEFINITIONS ENABLECACHING=0
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_stokes_channel_3d_staircase-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_ff_stokes_channel_3d_staircase_nocaching-00001.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_stokes_channel_3d params.input
+ -Problem.Name test_ff_stokes_channel_3d_staircase_nocaching -Problem.IsStaircaseGeometry true"
+ --zeroThreshold {"velocity_liq \(m/s\)":1e-12})
+
dumux_add_test(NAME test_ff_stokes_channel_pseudo3d
LABELS freeflow navierstokes
SOURCES main.cc
diff --git a/test/freeflow/navierstokes/channel/3d/main.cc b/test/freeflow/navierstokes/channel/3d/main.cc
index 2ba127d34d836c842abf182306e09bfc635bb747..4a3e6912a1fe19f6a4b254071bd1710b24645301 100644
--- a/test/freeflow/navierstokes/channel/3d/main.cc
+++ b/test/freeflow/navierstokes/channel/3d/main.cc
@@ -29,20 +29,22 @@
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
-#include <dune/grid/io/file/vtk.hh>
-#include <dune/istl/io.hh>
-#include <dumux/assembly/staggeredfvassembler.hh>
-#include <dumux/assembly/diffmethod.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>
-#include <dumux/freeflow/navierstokes/staggered/fluxoversurface.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager_sub.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>
-#include <dumux/io/staggeredvtkoutputmodule.hh>
+
#include <dumux/linear/seqsolverbackend.hh>
-#include <dumux/nonlinear/newtonsolver.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/traits.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+
+#include <dumux/freeflow/navierstokes/velocityoutput.hh>
+#include <dumux/freeflow/navierstokes/fluxoveraxisalignedsurface.hh>
#include "properties.hh"
@@ -51,7 +53,8 @@ int main(int argc, char** argv)
using namespace Dumux;
// define the type tag for this problem
- using TypeTag = Properties::TTag::ThreeDChannelTest;
+ using MomentumTypeTag = Properties::TTag::ThreeDChannelTestMomentum;
+ using MassTypeTag = Properties::TTag::ThreeDChannelTestMass;
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
@@ -64,8 +67,8 @@ int main(int argc, char** argv)
Parameters::init(argc, argv);
// create a grid
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using Grid = GetPropType<TypeTag, Properties::Grid>;
+ using Scalar = GetPropType<MassTypeTag, Properties::Scalar>;
+ using Grid = GetPropType<MassTypeTag, Properties::Grid>;
Dumux::GridManager<Grid> gridManager;
#if HAVE_DUNE_SUBGRID && GRID_DIM == 3
@@ -87,6 +90,8 @@ int main(int argc, char** argv)
};
gridManager.init(selector, "Internal");
+#elif HAVE_DUNE_SUBGRID
+ gridManager.init([&](const auto& element) { return true; });
#else
gridManager.init();
#endif
@@ -99,110 +104,110 @@ int main(int argc, char** argv)
const auto& leafGridView = gridManager.grid().leafGridView();
// create the finite volume grid geometry
- using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
- auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
-
- // the problem (boundary conditions)
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- auto problem = std::make_shared<Problem>(gridGeometry);
+ using MomentumGridGeometry = GetPropType<MomentumTypeTag, Properties::GridGeometry>;
+ auto momentumGridGeometry = std::make_shared<MomentumGridGeometry>(leafGridView);
+ using MassGridGeometry = GetPropType<MassTypeTag, Properties::GridGeometry>;
+ auto massGridGeometry = std::make_shared<MassGridGeometry>(leafGridView);
+
+ // the coupling manager
+ using Traits = MultiDomainTraits<MomentumTypeTag, MassTypeTag>;
+ using CouplingManager = StaggeredFreeFlowCouplingManager<Traits>;
+ auto couplingManager = std::make_shared<CouplingManager>();
+
+ // the problems (boundary conditions)
+ using MomentumProblem = GetPropType<MomentumTypeTag, Properties::Problem>;
+ auto momentumProblem = std::make_shared<MomentumProblem>(momentumGridGeometry, couplingManager);
+ using MassProblem = GetPropType<MassTypeTag, Properties::Problem>;
+ auto massProblem = std::make_shared<MassProblem>(massGridGeometry, couplingManager);
// the solution vector
- using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+ constexpr auto momentumIdx = CouplingManager::freeFlowMomentumIndex;
+ constexpr auto massIdx = CouplingManager::freeFlowMassIndex;
+ using SolutionVector = typename Traits::SolutionVector;
SolutionVector x;
- x[GridGeometry::cellCenterIdx()].resize(gridGeometry->numCellCenterDofs());
- x[GridGeometry::faceIdx()].resize(gridGeometry->numFaceDofs());
+ x[momentumIdx].resize(momentumGridGeometry->numDofs());
+ x[massIdx].resize(massGridGeometry->numDofs());
// the grid variables
- using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
- auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
- gridVariables->init(x);
-
- // intialize the vtk output module
- using IOFields = GetPropType<TypeTag, Properties::IOFields>;
- StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
+ using MomentumGridVariables = GetPropType<MomentumTypeTag, Properties::GridVariables>;
+ auto momentumGridVariables = std::make_shared<MomentumGridVariables>(momentumProblem, momentumGridGeometry);
+ using MassGridVariables = GetPropType<MassTypeTag, Properties::GridVariables>;
+ auto massGridVariables = std::make_shared<MassGridVariables>(massProblem, massGridGeometry);
+
+ // compute coupling stencil and afterwards initialize grid variables (need coupling information)
+ couplingManager->init(momentumProblem, massProblem, std::make_tuple(momentumGridVariables, massGridVariables), x);
+ massGridVariables->init(x[massIdx]);
+ momentumGridVariables->init(x[momentumIdx]);
+
+ using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(std::make_tuple(momentumProblem, massProblem),
+ std::make_tuple(momentumGridGeometry, massGridGeometry),
+ std::make_tuple(momentumGridVariables, massGridVariables),
+ couplingManager);
+
+ // initialize the vtk output module
+ using IOFields = GetPropType<MassTypeTag, Properties::IOFields>;
+ VtkOutputModule vtkWriter(*massGridVariables, x[massIdx], massProblem->name());
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
+ vtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<MassGridVariables>>());
vtkWriter.write(0.0);
- // the assembler with time loop for instationary problem
- using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
- auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
-
// the linear solver
using LinearSolver = Dumux::UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
- using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
- NewtonSolver nonLinearSolver(assembler, linearSolver);
-
- // set up two planes over which fluxes are calculated
- FluxOverSurface<GridVariables,
- SolutionVector,
- GetPropType<TypeTag, Properties::ModelTraits>,
- GetPropType<TypeTag, Properties::LocalResidual>> flux(*gridVariables, x);
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+ using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ // set up three planes over which fluxes are calculated
+ FluxOverAxisAlignedSurface flux(*massGridVariables, x[massIdx], assembler->localResidual(massIdx));
+
+ using GridView = typename GetPropType<MassTypeTag, Properties::GridGeometry>::GridView;
using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
- const Scalar xMin = gridGeometry->bBoxMin()[0];
- const Scalar xMax = gridGeometry->bBoxMax()[0];
- const Scalar yMin = gridGeometry->bBoxMin()[1];
- const Scalar yMax = gridGeometry->bBoxMax()[1];
+ const Scalar xMin = massGridGeometry->bBoxMin()[0];
+ const Scalar xMax = massGridGeometry->bBoxMax()[0];
+ const Scalar yMin = massGridGeometry->bBoxMin()[1];
+ const Scalar yMax = massGridGeometry->bBoxMax()[1];
#if GRID_DIM == 3
- const Scalar zMin = gridGeometry->bBoxMin()[2];
- const Scalar zMax = gridGeometry->bBoxMax()[2];
+ const Scalar zMin = massGridGeometry->bBoxMin()[2];
+ const Scalar zMax = massGridGeometry->bBoxMax()[2];
#endif
- // The first plane shall be placed at the middle of the channel.
- // If we have an odd number of cells in x-direction, there would not be any cell faces
- // at the postion of the plane (which is required for the flux calculation).
- // In this case, we add half a cell-width to the x-position in order to make sure that
- // the cell faces lie on the plane. This assumes a regular cartesian grid.
- // The second plane is placed at the outlet of the channel.
+ // the first plane is at the inlet
#if GRID_DIM == 3
- const auto p0inlet = GlobalPosition{xMin, yMin, zMin};
- const auto p1inlet = GlobalPosition{xMin, yMax, zMin};
- const auto p2inlet = GlobalPosition{xMin, yMin, zMax};
- const auto p3inlet = GlobalPosition{xMin, yMax, zMax};
- flux.addSurface("inlet", p0inlet, p1inlet, p2inlet, p3inlet);
+ const auto inletLowerLeft = GlobalPosition{xMin, yMin, zMin};
+ const auto inletUpperRight = GlobalPosition{xMin, yMax, zMax};
+ flux.addAxisAlignedSurface("inlet", inletLowerLeft, inletUpperRight);
#else
- const auto p0inlet = GlobalPosition{xMin, yMin};
- const auto p1inlet = GlobalPosition{xMin, yMax};
- flux.addSurface("inlet", p0inlet, p1inlet);
+ const auto inletLowerLeft = GlobalPosition{xMin, yMin};
+ const auto inletUpperRight = GlobalPosition{xMin, yMax};
+ flux.addAxisAlignedSurface("inlet", inletLowerLeft, inletUpperRight);
#endif
+ // the second plane is at the middle of the channel
const Scalar planePosMiddleX = xMin + 0.5*(xMax - xMin);
- int numCellsX = getParam<std::vector<int>>("Grid.Cells")[0];
-
- const unsigned int refinement = getParam<unsigned int>("Grid.Refinement", 0);
-
- numCellsX *= (1<<refinement);
-
- const Scalar offsetX = (numCellsX % 2 == 0) ? 0.0 : 0.5*((xMax - xMin) / numCellsX);
-
#if GRID_DIM == 3
- const auto p0middle = GlobalPosition{planePosMiddleX + offsetX, yMin, zMin};
- const auto p1middle = GlobalPosition{planePosMiddleX + offsetX, yMax, zMin};
- const auto p2middle = GlobalPosition{planePosMiddleX + offsetX, yMin, zMax};
- const auto p3middle = GlobalPosition{planePosMiddleX + offsetX, yMax, zMax};
- flux.addSurface("middle", p0middle, p1middle, p2middle, p3middle);
+ const auto middleLowerLeft = GlobalPosition{planePosMiddleX, yMin, zMin};
+ const auto middleUpperRight = GlobalPosition{planePosMiddleX, yMax, zMax};
+ flux.addAxisAlignedSurface("middle", middleLowerLeft, middleUpperRight);
#else
- const auto p0middle = GlobalPosition{planePosMiddleX + offsetX, yMin};
- const auto p1middle = GlobalPosition{planePosMiddleX + offsetX, yMax};
-flux.addSurface("middle", p0middle, p1middle);
+ const auto middleLowerLeft = GlobalPosition{planePosMiddleX, yMin};
+ const auto middleUpperRight = GlobalPosition{planePosMiddleX, yMax};
+ flux.addAxisAlignedSurface("middle", middleLowerLeft, middleUpperRight);
#endif
- // The second plane is placed at the outlet of the channel.
+ // The last plane is placed at the outlet of the channel.
#if GRID_DIM == 3
- const auto p0outlet = GlobalPosition{xMax, yMin, zMin};
- const auto p1outlet = GlobalPosition{xMax, yMax, zMin};
- const auto p2outlet = GlobalPosition{xMax, yMin, zMax};
- const auto p3outlet = GlobalPosition{xMax, yMax, zMax};
- flux.addSurface("outlet", p0outlet, p1outlet, p2outlet, p3outlet);
+ const auto outletLowerLeft = GlobalPosition{xMax, yMin, zMin};
+ const auto outletUpperRight = GlobalPosition{xMax, yMax, zMax};
+ flux.addAxisAlignedSurface("outlet", outletLowerLeft, outletUpperRight);
#else
- const auto p0outlet = GlobalPosition{xMax, yMin};
- const auto p1outlet = GlobalPosition{xMax, yMax};
- flux.addSurface("outlet", p0outlet, p1outlet);
+ const auto outletLowerLeft = GlobalPosition{xMax, yMin};
+ const auto outletUpperRight = GlobalPosition{xMax, yMax};
+ flux.addAxisAlignedSurface("outlet", outletLowerLeft, outletUpperRight);
#endif
// linearize & solve
@@ -213,37 +218,19 @@ flux.addSurface("middle", p0middle, p1middle);
vtkWriter.write(1.0);
// calculate and print mass fluxes over the planes
- flux.calculateMassOrMoleFluxes();
- if(GetPropType<TypeTag, Properties::ModelTraits>::enableEnergyBalance())
- {
- std::cout << "mass / energy flux at inlet is: " << flux.netFlux("inlet") << std::endl;
- std::cout << "mass / energy flux at middle is: " << flux.netFlux("middle") << std::endl;
- std::cout << "mass / energy flux at outlet is: " << flux.netFlux("outlet") << std::endl;
- }
- else
- {
- std::cout << "mass flux at inlet is: " << flux.netFlux("inlet") << std::endl;
- std::cout << "mass flux at middle is: " << flux.netFlux("middle") << std::endl;
- std::cout << "mass 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 middle is: " << flux.netFlux("middle")[0] << std::endl;
- std::cout << "volume flux at outlet is: " << flux.netFlux("outlet")[0] << std::endl;
-
+ flux.calculateAllFluxes();
+ std::cout << "mass flux at inlet is: " << flux.flux("inlet") << std::endl;
+ std::cout << "mass flux at middle is: " << flux.flux("middle") << std::endl;
+ std::cout << "mass flux at outlet is: " << flux.flux("outlet") << std::endl;
+ std::cout << "analyticalFlux: " << massProblem->analyticalFlux()*1e3 << std::endl;
timer.stop();
- std::cout << "analyticalFlux: " << problem->analyticalFlux() << std::endl;
-
const auto& comm = Dune::MPIHelper::getCommunication();
std::cout << "Simulation took " << timer.elapsed() << " seconds on "
<< comm.size() << " processes.\n"
<< "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n";
-
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
@@ -256,4 +243,4 @@ flux.addSurface("middle", p0middle, p1middle);
}
return 0;
-} // end main
+}
diff --git a/test/freeflow/navierstokes/channel/3d/params.input b/test/freeflow/navierstokes/channel/3d/params.input
index 9d69b20680936982cff50c6b10a510ab97fca632..f64ce1a23b32cea2955ddbf023f6c31e1276901c 100644
--- a/test/freeflow/navierstokes/channel/3d/params.input
+++ b/test/freeflow/navierstokes/channel/3d/params.input
@@ -1,9 +1,9 @@
[Grid]
UpperRight = 0.005 0.002 0.00025
-Cells = 3 20 10
+Cells = 3 20 10
[Problem]
-Name = test_stokes_channel_3d # name passed to the output routines
+Name = test_stokes_channel_3d
DeltaP = 10
EnableGravity = false
Height = 0.00025
diff --git a/test/freeflow/navierstokes/channel/3d/params_pseudo.input b/test/freeflow/navierstokes/channel/3d/params_pseudo.input
index e450d1a3b7a0259b60763e538f19c8350e47af73..e2f1b9c5717addb68f9982ccc960b600fda822b1 100644
--- a/test/freeflow/navierstokes/channel/3d/params_pseudo.input
+++ b/test/freeflow/navierstokes/channel/3d/params_pseudo.input
@@ -3,7 +3,7 @@ UpperRight = 0.005 0.002
Cells = 3 20
[Problem]
-Name = test_stokes_channel_pseudo3d # name passed to the output routines
+Name = test_stokes_channel_pseudo3d
DeltaP = 10
EnableGravity = false
Height = 0.00025
@@ -12,13 +12,12 @@ Height = 0.00025
LiquidDensity = 1e3
LiquidKinematicViscosity = 1e-6
-[Newton]
-MaxSteps = 10
-MaxRelativeShift = 1e-8
-
[Assembly]
NumericDifference.BaseEpsilon = 1e-8
+[Newton]
+MaxSteps = 10
+MaxRelativeShift = 1e-8
[Vtk]
WriteFaceData = false
diff --git a/test/freeflow/navierstokes/channel/3d/problem.hh b/test/freeflow/navierstokes/channel/3d/problem.hh
index c64755be9ea1dabf014c4f5241f5c801aa1f537a..b342843bd783a7e10978410397665864047049f7 100644
--- a/test/freeflow/navierstokes/channel/3d/problem.hh
+++ b/test/freeflow/navierstokes/channel/3d/problem.hh
@@ -25,17 +25,17 @@
* to mimic the 3D behavior of the flow.
*/
-#ifndef DUMUX_3D_CHANNEL_PROBLEM_HH
-#define DUMUX_3D_CHANNEL_PROBLEM_HH
+#ifndef DUMUX_TEST_FREEFLOW_NAVIERSTOKES_3D_CHANNEL_PROBLEM_HH
+#define DUMUX_TEST_FREEFLOW_NAVIERSTOKES_3D_CHANNEL_PROBLEM_HH
#include <dune/common/float_cmp.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
-#include <dumux/common/numeqvector.hh>
-#include <dumux/freeflow/navierstokes/boundarytypes.hh>
#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/freeflow/navierstokes/momentum/fluxhelper.hh>
+#include <dumux/freeflow/navierstokes/scalarfluxhelper.hh>
#ifndef GRID_DIM
#define GRID_DIM 3
@@ -59,33 +59,31 @@ class ThreeDChannelTestProblem : public NavierStokesProblem<TypeTag>
{
using ParentType = NavierStokesProblem<TypeTag>;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-
- using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
- using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
- using Element = typename GridView::template Codim<0>::Entity;
-
+ using BoundaryTypes = typename ParentType::BoundaryTypes;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ using NumEqVector = typename ParentType::NumEqVector;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+ using PrimaryVariables = typename ParentType::PrimaryVariables;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
- static constexpr int dim = GridView::dimension;
- static constexpr int dimWorld = GridView::dimensionworld;
+ static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
+ using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
- using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
- using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>;
-
- using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
- using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
static constexpr bool enablePseudoThreeDWallFriction = !(GRID_DIM == 3);
+ static constexpr int dim = GridGeometry::GridView::dimension;
public:
- ThreeDChannelTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
- : ParentType(gridGeometry)
+ ThreeDChannelTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(gridGeometry, couplingManager)
{
deltaP_ = getParam<Scalar>("Problem.DeltaP");
height_ = getParam<Scalar>("Problem.Height");
@@ -95,7 +93,7 @@ public:
if(dim == 3 && !Dune::FloatCmp::eq(height_, this->gridGeometry().bBoxMax()[2]))
DUNE_THROW(Dune::InvalidStateException, "z-dimension must equal height");
- if(enablePseudoThreeDWallFriction)
+ if constexpr (enablePseudoThreeDWallFriction)
extrusionFactor_ = 2.0/3.0 * height_;
else
extrusionFactor_ = 1.0;
@@ -117,31 +115,30 @@ public:
/*!
* \brief Evaluates the source term for all phases within a given
- * sub-control volume face.
+ * sub-control volume
*/
- using ParentType::source;
- template<class ElementVolumeVariables, class ElementFaceVariables>
- NumEqVector source(const Element &element,
+ template<class ElementVolumeVariables>
+ NumEqVector source(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
- const ElementFaceVariables& elemFaceVars,
- const SubControlVolumeFace &scvf) const
+ const SubControlVolume& scv) const
{
auto source = NumEqVector(0.0);
-#if GRID_DIM != 3
+ if constexpr (ParentType::isMomentumProblem() && enablePseudoThreeDWallFriction)
+ {
static const Scalar height = getParam<Scalar>("Problem.Height");
static const Scalar factor = getParam<Scalar>("Problem.PseudoWallFractionFactor", 8.0);
- source[scvf.directionIndex()] = this->pseudo3DWallFriction(scvf, elemVolVars, elemFaceVars, height, factor);
-#endif
+ source[scv.dofAxis()] = this->pseudo3DWallFriction(element, fvGeometry, elemVolVars, scv, height, factor);
+ }
return source;
}
+ //! the domain's extrusion factor
Scalar extrusionFactorAtPos(const GlobalPosition& pos) const
{ return extrusionFactor_; }
-
// \}
/*!
* \name Boundary conditions
@@ -158,16 +155,14 @@ public:
{
BoundaryTypes values;
- // set a fixed pressure at the inlet and outlet
- if (isOutlet_(globalPos) || isInlet_(globalPos))
- values.setDirichlet(Indices::pressureIdx);
- else
+ if constexpr (ParentType::isMomentumProblem())
{
- values.setDirichlet(Indices::velocityXIdx);
- values.setDirichlet(Indices::velocityYIdx);
- if(dim == 3)
- values.setDirichlet(Indices::velocityZIdx);
+ values.setAllDirichlet();
+ if (isOutlet_(globalPos) || isInlet_(globalPos))
+ values.setAllNeumann();
}
+ else
+ values.setNeumann(Indices::conti0EqIdx);
return values;
}
@@ -179,30 +174,48 @@ public:
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
- PrimaryVariables values(0.0);
-
- if(isInlet_(globalPos))
- values[Indices::pressureIdx] = 1e5 + deltaP_;
-
- if(isOutlet_(globalPos))
- values[Indices::pressureIdx] = 1e5;
-
- return values;
+ // no-flow/no-slip
+ return PrimaryVariables(0.0);
}
- // \}
-
/*!
- * \brief Evaluates the initial value for a control volume.
+ * \brief Evaluates the boundary conditions for a Neumann control volume.
*
- * \param globalPos The global position
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeometry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param elemFaceVars The element face variables
+ * \param scvf The boundary sub control volume face
*/
- PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ template<class ElementVolumeVariables, class ElementFluxVariablesCache>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
{
- PrimaryVariables values(0.0);
+ NumEqVector values(0.0);
+ const auto& globalPos = scvf.ipGlobal();
+
+ if constexpr (ParentType::isMomentumProblem())
+ {
+ const auto p = isInlet_(globalPos) ? 1e5 + deltaP_ : 1e5;
+ values = NavierStokesMomentumBoundaryFluxHelper::fixedPressureMomentumFlux(
+ *this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache, p
+ );
+ }
+ else
+ {
+ values = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>::scalarOutflowFlux(
+ *this, element, fvGeometry, scvf, elemVolVars
+ );
+ }
+
return values;
}
+ // \}
+
//! Returns the analytical solution for the flux through the rectangular channel
Scalar analyticalFlux() const
{
@@ -218,14 +231,10 @@ public:
private:
bool isInlet_(const GlobalPosition& globalPos) const
- {
- return globalPos[0] < eps_;
- }
+ { return globalPos[0] < eps_; }
bool isOutlet_(const GlobalPosition& globalPos) const
- {
- return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_;
- }
+ { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
static constexpr Scalar eps_=1e-6;
Scalar deltaP_;
@@ -234,6 +243,7 @@ private:
Scalar rho_;
Scalar nu_;
};
+
} // end namespace Dumux
#endif
diff --git a/test/freeflow/navierstokes/channel/3d/properties.hh b/test/freeflow/navierstokes/channel/3d/properties.hh
index 8a0e341762c58163a599c0fe8041c9f0aa22816b..cb23a1d3976de04f9da26d9c860e98621616373e 100644
--- a/test/freeflow/navierstokes/channel/3d/properties.hh
+++ b/test/freeflow/navierstokes/channel/3d/properties.hh
@@ -21,27 +21,38 @@
* \ingroup NavierStokesTests
* \brief The properties for the channel flow test for the staggered grid (Navier-)Stokes model.
*/
-#ifndef DUMUX_3D_CHANNEL_PROPERTIES_HH
-#define DUMUX_3D_CHANNEL_PROPERTIES_HH
+#ifndef DUMUX_TEST_FREEFLOW_NAVIERSTOKES_3D_CHANNEL_PROPERTIES_HH
+#define DUMUX_TEST_FREEFLOW_NAVIERSTOKES_3D_CHANNEL_PROPERTIES_HH
+
+#ifndef ENABLECACHING
+#define ENABLECACHING 1
+#endif
#include <dune/grid/yaspgrid.hh>
-#include <dumux/discretization/staggered/freeflow/properties.hh>
-#include <dumux/freeflow/navierstokes/model.hh>
+#if HAVE_DUNE_SUBGRID
+#include <dune/subgrid/subgrid.hh>
+#endif
+
+#include <dumux/discretization/fcstaggered.hh>
+#include <dumux/discretization/cctpfa.hh>
+#include <dumux/multidomain/staggeredfreeflow/couplingmanager.hh>
+
+#include <dumux/freeflow/navierstokes/momentum/model.hh>
+#include <dumux/freeflow/navierstokes/mass/1p/model.hh>
+
#include <dumux/material/components/constant.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include "problem.hh"
-#if HAVE_DUNE_SUBGRID && GRID_DIM == 3
-#include <dune/subgrid/subgrid.hh>
-#endif
-
namespace Dumux::Properties {
// Create new type tags
namespace TTag {
-struct ThreeDChannelTest { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
+struct ThreeDChannelTest {};
+struct ThreeDChannelTestMomentum { using InheritsFrom = std::tuple<ThreeDChannelTest, NavierStokesMomentum, FaceCenteredStaggeredModel>; };
+struct ThreeDChannelTestMass { using InheritsFrom = std::tuple<ThreeDChannelTest, NavierStokesMassOneP, CCTpfaModel>; };
} // end namespace TTag
// the fluid system
@@ -60,7 +71,7 @@ struct Grid<TypeTag, TTag::ThreeDChannelTest>
using HostGrid = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, dim> >;
-#if HAVE_DUNE_SUBGRID && GRID_DIM == 3
+#if HAVE_DUNE_SUBGRID
using type = Dune::SubGrid<HostGrid::dimension, HostGrid>;
#else
using type = HostGrid;
@@ -72,11 +83,18 @@ template<class TypeTag>
struct Problem<TypeTag, TTag::ThreeDChannelTest> { using type = ThreeDChannelTestProblem<TypeTag> ; };
template<class TypeTag>
-struct EnableGridGeometryCache<TypeTag, TTag::ThreeDChannelTest> { static constexpr bool value = true; };
+struct EnableGridGeometryCache<TypeTag, TTag::ThreeDChannelTest> { static constexpr bool value = ENABLECACHING; };
template<class TypeTag>
-struct EnableGridFluxVariablesCache<TypeTag, TTag::ThreeDChannelTest> { static constexpr bool value = true; };
+struct EnableGridFluxVariablesCache<TypeTag, TTag::ThreeDChannelTest> { static constexpr bool value = ENABLECACHING; };
template<class TypeTag>
-struct EnableGridVolumeVariablesCache<TypeTag, TTag::ThreeDChannelTest> { static constexpr bool value = true; };
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::ThreeDChannelTest> { static constexpr bool value = ENABLECACHING; };
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::ThreeDChannelTest>
+{
+ using Traits = MultiDomainTraits<TTag::ThreeDChannelTestMomentum, TTag::ThreeDChannelTestMass>;
+ using type = StaggeredFreeFlowCouplingManager<Traits>;
+};
} // end namespace Dumux::Properties