From c72379e12ceabc403f35b402f781bfece4569987 Mon Sep 17 00:00:00 2001
From: Timo Koch <timokoch@math.uio.no>
Date: Fri, 17 Mar 2023 19:52:47 +0100
Subject: [PATCH] [fcdiamond] Treat fcdiamond on the boundary like a CVFE
scheme
---
dumux/common/fvproblem.hh | 19 +++++----
.../cvfe/elementboundarytypes.hh | 40 ++++---------------
.../freeflow/navierstokes/momentum/problem.hh | 22 +++++-----
.../navierstokes/analyticalsolutionvectors.hh | 4 +-
.../navierstokes/unstructured/main.cc | 33 ++++-----------
5 files changed, 35 insertions(+), 83 deletions(-)
diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh
index 5fd5ff9fcf..bfdc48a2c5 100644
--- a/dumux/common/fvproblem.hh
+++ b/dumux/common/fvproblem.hh
@@ -71,8 +71,7 @@ class FVProblem
using PointSourceMap = std::map< std::pair<std::size_t, std::size_t>,
std::vector<PointSource> >;
- static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethods::box;
- static constexpr bool isPQ1Bubble = GridGeometry::discMethod == DiscretizationMethods::pq1bubble;
+ static constexpr bool isCVFE = DiscretizationMethods::isCVFE<typename GridGeometry::DiscretizationMethod>;
static constexpr bool isStaggered = GridGeometry::discMethod == DiscretizationMethods::staggered;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@@ -143,9 +142,9 @@ public:
auto boundaryTypes(const Element &element,
const SubControlVolume &scv) const
{
- if (!isBox && !isPQ1Bubble)
+ if (!isCVFE)
DUNE_THROW(Dune::InvalidStateException,
- "boundaryTypes(..., scv) called for cell-centered method.");
+ "boundaryTypes(..., scv) called for non-CVFE method.");
// forward it to the method which only takes the global coordinate
return asImp_().boundaryTypesAtPos(scv.dofPosition());
@@ -161,9 +160,9 @@ public:
auto boundaryTypes(const Element &element,
const SubControlVolumeFace &scvf) const
{
- if (isBox || isPQ1Bubble)
+ if (isCVFE)
DUNE_THROW(Dune::InvalidStateException,
- "boundaryTypes(..., scvf) called for box method.");
+ "boundaryTypes(..., scvf) called for CVFE method.");
// forward it to the method which only takes the global coordinate
return asImp_().boundaryTypesAtPos(scvf.ipGlobal());
@@ -195,9 +194,9 @@ public:
PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
{
// forward it to the method which only takes the global coordinate
- if (isBox || isPQ1Bubble)
+ if (isCVFE)
{
- DUNE_THROW(Dune::InvalidStateException, "dirichlet(scvf) called for box method.");
+ DUNE_THROW(Dune::InvalidStateException, "dirichlet(scvf) called for CVFE method.");
}
else
return asImp_().dirichletAtPos(scvf.ipGlobal());
@@ -214,9 +213,9 @@ public:
PrimaryVariables dirichlet(const Element &element, const SubControlVolume &scv) const
{
// forward it to the method which only takes the global coordinate
- if (!isBox && !isStaggered && !isPQ1Bubble)
+ if (!isCVFE && !isStaggered)
{
- DUNE_THROW(Dune::InvalidStateException, "dirichlet(scv) called for other than box or staggered method.");
+ DUNE_THROW(Dune::InvalidStateException, "dirichlet(scv) called for other than CVFE or staggered method.");
}
else
return asImp_().dirichletAtPos(scv.dofPosition());
diff --git a/dumux/discretization/cvfe/elementboundarytypes.hh b/dumux/discretization/cvfe/elementboundarytypes.hh
index ab58b07004..d3abb4eb46 100644
--- a/dumux/discretization/cvfe/elementboundarytypes.hh
+++ b/dumux/discretization/cvfe/elementboundarytypes.hh
@@ -54,47 +54,21 @@ public:
const typename FVElementGeometry::Element& element,
const FVElementGeometry& fvGeometry)
{
- using GridGeometry = typename FVElementGeometry::GridGeometry;
-
- if constexpr (GridGeometry::discMethod == DiscretizationMethods::fcdiamond)
- if (!fvGeometry.hasBoundaryScvf())
- return;
-
bcTypes_.resize(fvGeometry.numScv());
hasDirichlet_ = false;
hasNeumann_ = false;
- // boundary dofs always have boundary corresponding faces
- if constexpr (GridGeometry::discMethod == DiscretizationMethods::fcdiamond)
+ for (const auto& scv : scvs(fvGeometry))
{
- for (const auto& scvf : scvfs(fvGeometry))
- {
- if (scvf.boundary())
- {
- const auto localIndex = fvGeometry.scv(scvf.insideScvIdx()).localDofIndex();
- bcTypes_[localIndex] = problem.boundaryTypes(element, scvf);
- hasDirichlet_ = hasDirichlet_ || bcTypes_[localIndex].hasDirichlet();
- hasNeumann_ = hasNeumann_ || bcTypes_[localIndex].hasNeumann();
- }
- }
- }
+ const auto scvIdxLocal = scv.localDofIndex();
+ bcTypes_[scvIdxLocal].reset();
- // generally, single vertices maybe on the boundary even if the face is not
- // in this case we need a different mechanism
- else
- {
- for (const auto& scv : scvs(fvGeometry))
+ if (fvGeometry.gridGeometry().dofOnBoundary(scv.dofIndex()))
{
- const auto scvIdxLocal = scv.localDofIndex();
- bcTypes_[scvIdxLocal].reset();
-
- if (fvGeometry.gridGeometry().dofOnBoundary(scv.dofIndex()))
- {
- bcTypes_[scvIdxLocal] = problem.boundaryTypes(element, scv);
- hasDirichlet_ = hasDirichlet_ || bcTypes_[scvIdxLocal].hasDirichlet();
- hasNeumann_ = hasNeumann_ || bcTypes_[scvIdxLocal].hasNeumann();
- }
+ bcTypes_[scvIdxLocal] = problem.boundaryTypes(element, scv);
+ hasDirichlet_ = hasDirichlet_ || bcTypes_[scvIdxLocal].hasDirichlet();
+ hasNeumann_ = hasNeumann_ || bcTypes_[scvIdxLocal].hasNeumann();
}
}
}
diff --git a/dumux/freeflow/navierstokes/momentum/problem.hh b/dumux/freeflow/navierstokes/momentum/problem.hh
index 4000f78c7f..a8b1ca719b 100644
--- a/dumux/freeflow/navierstokes/momentum/problem.hh
+++ b/dumux/freeflow/navierstokes/momentum/problem.hh
@@ -568,9 +568,9 @@ class NavierStokesMomentumProblemImpl<TypeTag, DiscretizationMethods::CVFE<DM>>
using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
- static constexpr bool isVertexCentered
- = GridGeometry::discMethod == DiscretizationMethods::box
- || GridGeometry::discMethod == DiscretizationMethods::pq1bubble;
+ static constexpr bool isCVFE = DiscretizationMethods::isCVFE<
+ typename GridGeometry::DiscretizationMethod
+ >;
public:
//! These types are used in place of the typical NumEqVector type.
@@ -674,9 +674,9 @@ public:
BoundaryTypes boundaryTypes(const Element& element,
const SubControlVolume& scv) const
{
- if (!isVertexCentered)
+ if (!isCVFE)
DUNE_THROW(Dune::InvalidStateException,
- "boundaryTypes(..., scv) called for for a non vertex-centered method.");
+ "boundaryTypes(..., scv) called for for a non CVFE method.");
// forward it to the method which only takes the global coordinate
return asImp_().boundaryTypesAtPos(scv.dofPosition());
@@ -692,9 +692,9 @@ public:
BoundaryTypes boundaryTypes(const Element& element,
const SubControlVolumeFace& scvf) const
{
- if (isVertexCentered)
+ if (isCVFE)
DUNE_THROW(Dune::InvalidStateException,
- "boundaryTypes(..., scvf) called for a vertex-centered method.");
+ "boundaryTypes(..., scvf) called for a CVFE method.");
// forward it to the method which only takes the global coordinate
return asImp_().boundaryTypesAtPos(scvf.ipGlobal());
@@ -710,8 +710,8 @@ public:
DirichletValues dirichlet(const Element& element, const SubControlVolume& scv) const
{
// forward it to the method which only takes the global coordinate
- if (!isVertexCentered)
- DUNE_THROW(Dune::InvalidStateException, "dirichlet(scv) called for a non vertex-centered method.");
+ if (!isCVFE)
+ DUNE_THROW(Dune::InvalidStateException, "dirichlet(scv) called for a non CVFE method.");
else
return asImp_().dirichletAtPos(scv.dofPosition());
}
@@ -726,8 +726,8 @@ public:
DirichletValues dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
{
// forward it to the method which only takes the global coordinate
- if (isVertexCentered)
- DUNE_THROW(Dune::InvalidStateException, "dirichlet(scvf) called for vertex-centered method.");
+ if (isCVFE)
+ DUNE_THROW(Dune::InvalidStateException, "dirichlet(scvf) called for CVFE method.");
else
return asImp_().dirichletAtPos(scvf.ipGlobal());
}
diff --git a/test/freeflow/navierstokes/analyticalsolutionvectors.hh b/test/freeflow/navierstokes/analyticalsolutionvectors.hh
index d67b86d059..f1c4046757 100644
--- a/test/freeflow/navierstokes/analyticalsolutionvectors.hh
+++ b/test/freeflow/navierstokes/analyticalsolutionvectors.hh
@@ -185,9 +185,7 @@ public:
analyticalVelocityAtDofs_[scv.dofIndex()][scv.dofAxis()]
= momentumProblem_->analyticalSolution(scv.center(), time)[MomIndices::velocity(scv.dofAxis())];
- else if constexpr (MomentumGridGeometry::discMethod == DiscretizationMethods::fcdiamond
- || MomentumGridGeometry::discMethod == DiscretizationMethods::pq1bubble
- || MomentumGridGeometry::discMethod == DiscretizationMethods::box)
+ else if constexpr (DiscretizationMethods::isCVFE<typename MomentumGridGeometry::DiscretizationMethod>)
for (const auto& scv : scvs(fvGeometry))
for (int dirIdx = 0; dirIdx < dimWorld; ++dirIdx)
analyticalVelocityAtDofs_[scv.dofIndex()][dirIdx]
diff --git a/test/freeflow/navierstokes/unstructured/main.cc b/test/freeflow/navierstokes/unstructured/main.cc
index cc39cf559f..563f829921 100644
--- a/test/freeflow/navierstokes/unstructured/main.cc
+++ b/test/freeflow/navierstokes/unstructured/main.cc
@@ -58,36 +58,17 @@ auto dirichletDofs(std::shared_ptr<MomGG> momentumGridGeometry,
auto fvGeometry = localView(*momentumGridGeometry);
for (const auto& element : elements(momentumGridGeometry->gridView()))
{
- if constexpr (MomGG::discMethod == Dumux::DiscretizationMethods::pq1bubble)
+ fvGeometry.bind(element);
+ for (const auto& scv : scvs(fvGeometry))
{
- fvGeometry.bind(element);
- for (const auto& scv : scvs(fvGeometry))
+ if (momentumGridGeometry->dofOnBoundary(scv.dofIndex()))
{
- if (momentumGridGeometry->dofOnBoundary(scv.dofIndex()))
- {
- const auto bcTypes = momentumProblem->boundaryTypes(element, scv);
- for (int i = 0; i < bcTypes.size(); ++i)
- if (bcTypes.isDirichlet(i))
- dirichletDofs[momentumIdx][scv.dofIndex()][i] = 1.0;
- }
+ const auto bcTypes = momentumProblem->boundaryTypes(element, scv);
+ for (int i = 0; i < bcTypes.size(); ++i)
+ if (bcTypes.isDirichlet(i))
+ dirichletDofs[momentumIdx][scv.dofIndex()][i] = 1.0;
}
}
- else if constexpr (MomGG::discMethod == Dumux::DiscretizationMethods::fcdiamond)
- {
- fvGeometry.bind(element);
- for (const auto& scvf : scvfs(fvGeometry))
- {
- if (scvf.boundary())
- {
- const auto bcTypes = momentumProblem->boundaryTypes(element, scvf);
- for (int i = 0; i < bcTypes.size(); ++i)
- if (bcTypes.isDirichlet(i))
- dirichletDofs[momentumIdx][fvGeometry.scv(scvf.insideScvIdx()).dofIndex()][i] = 1.0;
- }
- }
- }
- else
- DUNE_THROW(Dune::NotImplemented, "dirichletDof help for discretization " << MomGG::discMethod);
}
return dirichletDofs;
--
GitLab