diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
index dc73e95de430ca8d8c780709b3100289fde74cbb..ff452ddb752ea92cb95f1577a7052cb87a0cdfde 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
@@ -515,8 +515,9 @@ private:
}
}();
+ const bool selfIsUpstream = lateralFace.directionSign() == sign(transportingVelocity);
StaggeredUpwindHelper<TypeTag, upwindSchemeOrder> upwindHelper(element, fvGeometry, scvf, faceVars, elemVolVars, gridFluxVarsCache.staggeredUpwindMethods());
- return upwindHelper.computeUpwindedLateralMomentum(localSubFaceIdx, currentScvfBoundaryTypes, lateralFaceBoundaryTypes)
+ return upwindHelper.computeUpwindLateralMomentum(selfIsUpstream, lateralFace, localSubFaceIdx, currentScvfBoundaryTypes, lateralFaceBoundaryTypes)
* transportingVelocity * lateralFace.directionSign() * lateralFace.area() * 0.5 * extrusionFactor_(elemVolVars, lateralFace);
}
diff --git a/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh
index ace92a10d35319924e98114f5906196e95777dfc..2b0c14330ae2043f0304c027c90cfade72d9be5e 100644
--- a/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh
@@ -131,25 +131,33 @@ public:
* Then the corresponding set of momenta are collected and the prescribed
* upwinding method is used to calculate the momentum.
*/
- FacePrimaryVariables computeUpwindedLateralMomentum(const int localSubFaceIdx,
- const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
- const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes) const
+ FacePrimaryVariables computeUpwindLateralMomentum(const bool selfIsUpstream,
+ const SubControlVolumeFace& lateralFace,
+ const int localSubFaceIdx,
+ const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
+ const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes) const
{
// Check whether the own or the neighboring element is upstream.
- const auto eIdx = scvf_.insideScvIdx();
- const auto& lateralFace = fvGeometry_.scvf(eIdx, scvf_.pairData(localSubFaceIdx).localLateralFaceIdx);
// Get the transporting velocity, located at the scvf perpendicular to the current scvf where the dof
// of interest is located.
- const Scalar transportingVelocity = faceVars_.velocityLateralInside(localSubFaceIdx);
+ const Scalar insideDensity = elemVolVars_[lateralFace.insideScvIdx()].density();
+ const Scalar outsideDensity = elemVolVars_[lateralFace.outsideScvIdx()].density();
- // Check whether the own or the neighboring element is upstream.
- const bool selfIsUpstream = ( lateralFace.directionSign() == sign(transportingVelocity) );
- const bool canHigherOrder = canLateralSecondOrder_(scvf_, selfIsUpstream, localSubFaceIdx);
- const auto parallelUpwindingMomenta = getLateralUpwindingMomenta_(elemVolVars_.gridVolVars().problem(), fvGeometry_, element_, scvf_, elemVolVars_, faceVars_,
- transportingVelocity, localSubFaceIdx, currentScvfBoundaryTypes,
- lateralFaceBoundaryTypes, canHigherOrder);
- return doLateralMomentumUpwinding_(fvGeometry_, scvf_, parallelUpwindingMomenta, transportingVelocity,
- localSubFaceIdx, upwindScheme_, canHigherOrder);
+ // for higher order schemes do higher order upwind reconstruction
+ if constexpr (useHigherOrder)
+ {
+ // only second order is implemented so far
+ if (canDoLateralSecondOrder_(selfIsUpstream, localSubFaceIdx))
+ {
+ const auto upwindMomenta = getLateralSecondOrderUpwindMomenta_(insideDensity, outsideDensity, selfIsUpstream, localSubFaceIdx, currentScvfBoundaryTypes, lateralFaceBoundaryTypes);
+ const auto distances = getLateralDistances_(localSubFaceIdx, selfIsUpstream);
+ return upwindScheme_.tvd(upwindMomenta, distances, selfIsUpstream, upwindScheme_.tvdApproach());
+ }
+ }
+
+ // otherwise apply first order upwind scheme
+ const auto upwindMomenta = getLateralFirstOrderUpwindMomenta_(insideDensity, outsideDensity, selfIsUpstream, localSubFaceIdx, currentScvfBoundaryTypes, lateralFaceBoundaryTypes);
+ return upwindScheme_.upwind(upwindMomenta[0], upwindMomenta[1]);
}
private:
@@ -225,282 +233,155 @@ private:
*
* This helper function checks if the scvf of interest is not too near to the
* boundary so that a dof upstream with respect to the upstream dof is available.
- *
- * If higher order methods are not prescribed, false is returned.
- *
- * \param ownScvf The SubControlVolumeFace we are considering
- * \param selfIsUpstream @c true if the velocity at ownScvf is upstream wrt the transporting velocity
- * \param localSubFaceIdx The local subface index
*/
- static bool canLateralSecondOrder_(const SubControlVolumeFace& ownScvf,
- [[maybe_unused]] const bool selfIsUpstream,
- [[maybe_unused]] const int localSubFaceIdx)
+ bool canDoLateralSecondOrder_(const bool selfIsUpstream, const int localSubFaceIdx) const
{
- if constexpr (useHigherOrder)
- {
- if (selfIsUpstream)
- {
- // The self velocity is upstream. The downstream velocity can be assigned or retrieved
- // from the boundary, even if there is no parallel neighbor.
- return true;
- }
- else
- {
- // The self velocity is downstream. If there is no parallel neighbor I cannot use a second order approximation.
- return ownScvf.hasParallelNeighbor(localSubFaceIdx, 0);
- }
- }
- else
- return false;
+ // If the self velocity is upstream, the downstream velocity can be assigned or retrieved
+ // from the boundary, even if there is no parallel neighbor.
+ // If the self velocity is downstream and there is no parallel neighbor I cannot use a second order approximation.
+ return selfIsUpstream || scvf_.hasParallelNeighbor(localSubFaceIdx, 0);
}
/*!
* \brief Returns an array of momenta needed for higher order or calls a function to return an array for basic upwinding methods.
*/
- static auto getLateralUpwindingMomenta_(const Problem& problem,
- const FVElementGeometry& fvGeometry,
- const Element& element,
- const SubControlVolumeFace& ownScvf,
- const ElementVolumeVariables& elemVolVars,
- const FaceVariables& faceVars,
- const Scalar transportingVelocity,
- const int localSubFaceIdx,
- const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
- const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
- [[maybe_unused]] const bool canHigherOrder)
+ std::array<Scalar, 3> getLateralSecondOrderUpwindMomenta_(const Scalar insideDensity,
+ const Scalar outsideDensity,
+ const bool selfIsUpstream,
+ const int localSubFaceIdx,
+ const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
+ const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes) const
{
- // Check whether the own or the neighboring element is upstream.
- const SubControlVolumeFace& lateralFace = fvGeometry.scvf(ownScvf.insideScvIdx(), ownScvf.pairData(localSubFaceIdx).localLateralFaceIdx);
+ static_assert(useHigherOrder, "Should only be reached if higher order methods are enabled");
- // Get the volume variables of the own and the neighboring element
- const auto& insideVolVars = elemVolVars[lateralFace.insideScvIdx()];
- const auto& outsideVolVars = elemVolVars[lateralFace.outsideScvIdx()];
+ // If the lateral face lies on a boundary, we assume that the parallel velocity on the boundary is actually known,
+ // thus we always use this value for the computation of the transported momentum.
+ if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0))
+ {
+ const Scalar boundaryVelocity = getParallelVelocityFromBoundary_(element_, scvf_, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
+ const Scalar boundaryMomentum = boundaryVelocity*insideDensity;
+ return {boundaryMomentum, boundaryMomentum, boundaryMomentum};
+ }
- // Check whether the own or the neighboring element is upstream.
- const bool selfIsUpstream = lateralFace.directionSign() == sign(transportingVelocity);
+ if (selfIsUpstream)
+ {
+ std::array<Scalar, 3> momenta;
+ momenta[1] = faceVars_.velocitySelf()*insideDensity;
+ momenta[0] = faceVars_.velocityParallel(localSubFaceIdx, 0)*insideDensity;
- if constexpr (useHigherOrder)
+ // The local index of the faces that is opposite to localSubFaceIdx
+ const int oppositeSubFaceIdx = localSubFaceIdx % 2 ? localSubFaceIdx - 1 : localSubFaceIdx + 1;
+
+ // The "upstream-upstream" velocity is retrieved from the other parallel neighbor or from the boundary.
+ if (scvf_.hasParallelNeighbor(oppositeSubFaceIdx, 0))
+ momenta[2] = faceVars_.velocityParallel(oppositeSubFaceIdx, 0)*insideDensity;
+ else
+ momenta[2] = getParallelVelocityFromOppositeBoundary_(element_, scvf_, currentScvfBoundaryTypes, oppositeSubFaceIdx)*insideDensity;
+ return momenta;
+ }
+ else
{
- if (canHigherOrder)
+ std::array<Scalar, 3> momenta;
+ momenta[0] = faceVars_.velocitySelf()*outsideDensity;
+ momenta[1] = faceVars_.velocityParallel(localSubFaceIdx, 0)*outsideDensity;
+
+ // If there is another parallel neighbor I can assign the "upstream-upstream" velocity, otherwise I retrieve it from the boundary.
+ if (scvf_.hasParallelNeighbor(localSubFaceIdx, 1))
+ momenta[2] = faceVars_.velocityParallel(localSubFaceIdx, 1)*outsideDensity;
+ else
{
- // If the lateral face lies on a boundary, we assume that the parallel velocity on the boundary is actually known,
- // thus we always use this value for the computation of the transported momentum.
- if (!ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
- {
- const Scalar boundaryMomentum = getParallelVelocityFromBoundary_(problem, element, fvGeometry, ownScvf,
- faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes,
- localSubFaceIdx) * insideVolVars.density();
-
- return std::array<Scalar, 3>{boundaryMomentum, boundaryMomentum, boundaryMomentum};
- }
-
- std::array<Scalar, 3> momenta;
- if (selfIsUpstream)
- {
- momenta[1] = faceVars.velocitySelf() * insideVolVars.density();
-
- if (ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
- momenta[0] = faceVars.velocityParallel(localSubFaceIdx, 0) * insideVolVars.density();
- else
- momenta[0] = getParallelVelocityFromBoundary_(problem, element, fvGeometry, ownScvf,
- faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes,
- localSubFaceIdx) * insideVolVars.density();
-
- // The local index of the faces that is opposite to localSubFaceIdx
- const int oppositeSubFaceIdx = localSubFaceIdx % 2 ? localSubFaceIdx - 1 : localSubFaceIdx + 1;
-
- // The "upstream-upstream" velocity is retrieved from the other parallel neighbor or from the boundary.
- if (ownScvf.hasParallelNeighbor(oppositeSubFaceIdx, 0))
- momenta[2] = faceVars.velocityParallel(oppositeSubFaceIdx, 0) * insideVolVars.density();
- else
- momenta[2] = getParallelVelocityFromOppositeBoundary_(problem, element, fvGeometry, ownScvf,
- faceVars, currentScvfBoundaryTypes,
- oppositeSubFaceIdx) * insideVolVars.density();
- }
- else
- {
- momenta[0] = faceVars.velocitySelf() * outsideVolVars.density();
- momenta[1] = faceVars.velocityParallel(localSubFaceIdx, 0) * outsideVolVars.density();
-
- // If there is another parallel neighbor I can assign the "upstream-upstream" velocity, otherwise I retrieve it from the boundary.
- if (ownScvf.hasParallelNeighbor(localSubFaceIdx, 1))
- momenta[2] = faceVars.velocityParallel(localSubFaceIdx, 1) * outsideVolVars.density();
- else
- {
- const Element& elementParallel = fvGeometry.gridGeometry().element(lateralFace.outsideScvIdx());
- const SubControlVolumeFace& firstParallelScvf = fvGeometry.scvf(lateralFace.outsideScvIdx(), ownScvf.localFaceIdx());
-
- momenta[2] = getParallelVelocityFromOppositeBoundary_(problem, elementParallel, fvGeometry, firstParallelScvf,
- faceVars, problem.boundaryTypes(elementParallel, firstParallelScvf),
- localSubFaceIdx) * outsideVolVars.density();
- }
- }
- return momenta;
+ const auto& lateralFace = fvGeometry_.scvf(scvf_.insideScvIdx(), scvf_.pairData(localSubFaceIdx).localLateralFaceIdx);
+ const auto& elementParallel = fvGeometry_.gridGeometry().element(lateralFace.outsideScvIdx());
+ const auto& firstParallelScvf = fvGeometry_.scvf(lateralFace.outsideScvIdx(), scvf_.localFaceIdx());
+ const auto& problem = elemVolVars_.gridVolVars().problem();
+ const auto& boundaryTypes = problem.boundaryTypes(elementParallel, firstParallelScvf);
+ momenta[2] = getParallelVelocityFromOppositeBoundary_(elementParallel, firstParallelScvf, boundaryTypes, localSubFaceIdx)*outsideDensity;
}
- else
- return getFirstOrderLateralUpwindingMomenta_(problem, element, fvGeometry, ownScvf, faceVars,
- currentScvfBoundaryTypes, lateralFaceBoundaryTypes,
- localSubFaceIdx, selfIsUpstream, insideVolVars, outsideVolVars);
+ return momenta;
}
- else
- return getFirstOrderLateralUpwindingMomenta_(problem, element, fvGeometry, ownScvf, faceVars,
- currentScvfBoundaryTypes, lateralFaceBoundaryTypes,
- localSubFaceIdx, selfIsUpstream, insideVolVars, outsideVolVars);
}
/*!
* \brief Returns an array of momenta needed for basic upwinding methods.
*/
- static auto getFirstOrderLateralUpwindingMomenta_(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const SubControlVolumeFace& ownScvf,
- const FaceVariables& faceVars,
- const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
- const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
- const int localSubFaceIdx,
- const bool selfIsUpstream,
- const VolumeVariables& insideVolVars,
- const VolumeVariables& outsideVolVars)
+ std::array<Scalar, 2> getLateralFirstOrderUpwindMomenta_(const Scalar insideDensity,
+ const Scalar outsideDensity,
+ const bool selfIsUpstream,
+ const int localSubFaceIdx,
+ const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
+ const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes) const
{
- const Scalar momentumParallel = ownScvf.hasParallelNeighbor(localSubFaceIdx, 0)
- ? faceVars.velocityParallel(localSubFaceIdx, 0) * outsideVolVars.density()
- : (getParallelVelocityFromBoundary_(problem, element, fvGeometry, ownScvf,
- faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes,
- localSubFaceIdx) * insideVolVars.density());
// If the lateral face lies on a boundary, we assume that the parallel velocity on the boundary is actually known,
// thus we always use this value for the computation of the transported momentum.
- if (!ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
+ if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0))
{
- if constexpr (useHigherOrder)
- return std::array<Scalar, 3>{momentumParallel, momentumParallel};
- else
- return std::array<Scalar, 2>{momentumParallel, momentumParallel};
+ const Scalar boundaryVelocity = getParallelVelocityFromBoundary_(element_, scvf_, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
+ const Scalar boundaryMomentum = boundaryVelocity*insideDensity;
+ return {boundaryMomentum, boundaryMomentum};
}
- const Scalar momentumSelf = faceVars.velocitySelf() * insideVolVars.density();
- if constexpr (useHigherOrder)
- return selfIsUpstream ? std::array<Scalar, 3>{momentumParallel, momentumSelf}
- : std::array<Scalar, 3>{momentumSelf, momentumParallel};
- else
- return selfIsUpstream ? std::array<Scalar, 2>{momentumParallel, momentumSelf}
- : std::array<Scalar, 2>{momentumSelf, momentumParallel};
-
- }
-
- /*!
- * \brief Returns the upwinded momentum for higher order or basic upwind schemes
- *
- * If higher order methods are not enabled, this is fowarded to the Frontal Momentum Upwinding method
- *
- * \param scvf The sub control volume face
- * \param momenta The momenta to be upwinded
- * \param transportingVelocity The average of the self and opposite velocities.
- * \param localSubFaceIdx The local index of the subface
- * \param gridFluxVarsCache The grid flux variables cache
- */
- template<class MomentaArray>
- static Scalar doLateralMomentumUpwinding_([[maybe_unused]] const FVElementGeometry& fvGeometry,
- const SubControlVolumeFace& scvf,
- const MomentaArray& momenta,
- [[maybe_unused]] const Scalar transportingVelocity,
- [[maybe_unused]] const int localSubFaceIdx,
- const UpwindScheme& upwindScheme,
- [[maybe_unused]] const bool canHigherOrder)
- {
- if constexpr (useHigherOrder)
- {
- if (canHigherOrder)
- {
- const auto eIdx = scvf.insideScvIdx();
- const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
- const bool selfIsUpstream = ( lateralFace.directionSign() == sign(transportingVelocity) );
- const std::array<Scalar,3> distances = getLateralDistances_(scvf, localSubFaceIdx, selfIsUpstream);
- return upwindScheme.tvd(momenta, distances, selfIsUpstream, upwindScheme.tvdApproach());
- }
- else
- return upwindScheme.upwind(momenta[0], momenta[1]);
- }
+ const Scalar momentumParallel = faceVars_.velocityParallel(localSubFaceIdx, 0)*outsideDensity;
+ const Scalar momentumSelf = faceVars_.velocitySelf()*insideDensity;
+ if (selfIsUpstream)
+ return {momentumParallel, momentumSelf};
else
- return upwindScheme.upwind(momenta[0], momenta[1]);
-
+ return {momentumSelf, momentumParallel};
}
/*!
* \brief Returns an array of distances needed for non-uniform higher order upwind schemes
- *
- * \param ownScvf The sub control volume face
- * \param localSubFaceIdx The local index of the subface
- * \param selfIsUpstream bool describing upstream face.
+ * computes lateral distances {upstream to downstream distance, up-upstream to upstream distance, downstream staggered cell size}
*/
- static std::array<Scalar, 3> getLateralDistances_(const SubControlVolumeFace& ownScvf,
- const int localSubFaceIdx,
- const bool selfIsUpstream)
+ std::array<Scalar, 3> getLateralDistances_(const int localSubFaceIdx, const bool selfIsUpstream) const
{
static_assert(useHigherOrder, "Should only be reached if higher order methods are enabled");
- // distances {upstream to downstream distance, up-upstream to upstream distance, downstream staggered cell size}
- std::array<Scalar, 3> distances;
-
if (selfIsUpstream)
{
// The local index of the faces that is opposite to localSubFaceIdx
const int oppositeSubFaceIdx = localSubFaceIdx % 2 ? localSubFaceIdx - 1 : localSubFaceIdx + 1;
- distances[0] = ownScvf.parallelDofsDistance(localSubFaceIdx, 0);
- distances[1] = ownScvf.parallelDofsDistance(oppositeSubFaceIdx, 0);
- if (ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
- distances[2] = ownScvf.pairData(localSubFaceIdx).parallelCellWidths[0];
+ std::array<Scalar, 3> distances;
+ distances[0] = scvf_.parallelDofsDistance(localSubFaceIdx, 0);
+ distances[1] = scvf_.parallelDofsDistance(oppositeSubFaceIdx, 0);
+ if (scvf_.hasParallelNeighbor(localSubFaceIdx, 0))
+ distances[2] = scvf_.pairData(localSubFaceIdx).parallelCellWidths[0];
else
- distances[2] = ownScvf.area() / 2.0;
+ distances[2] = scvf_.area() / 2.0;
+ return distances;
}
else
{
- distances[0] = ownScvf.parallelDofsDistance(localSubFaceIdx, 0);
- distances[1] = ownScvf.parallelDofsDistance(localSubFaceIdx, 1);
- distances[2] = ownScvf.pairData(localSubFaceIdx).parallelCellWidths[0];
+ std::array<Scalar, 3> distances;
+ distances[0] = scvf_.parallelDofsDistance(localSubFaceIdx, 0);
+ distances[1] = scvf_.parallelDofsDistance(localSubFaceIdx, 1);
+ distances[2] = scvf_.pairData(localSubFaceIdx).parallelCellWidths[0];
+ return distances;
}
-
- return distances;
}
/*!
* \brief Return the outer parallel velocity for normal faces that are on the boundary and therefore have no neighbor.
- *
* Calls the problem to retrieve a fixed value set on the boundary.
- *
- * \param problem The problem
- * \param scvf The SubControlVolumeFace that is normal to the boundary
- * \param velocitySelf the velocity at scvf
- * \param localSubFaceIdx The local index of the face that is on the boundary
- * \param element The element that is on the boundary
- * \param lateralFaceBoundaryTypes stores the type of boundary at the lateral face
*/
- static Scalar getParallelVelocityFromBoundary_(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const SubControlVolumeFace& scvf,
- const FaceVariables& faceVars,
- const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
- const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
- const int localSubFaceIdx)
+ Scalar getParallelVelocityFromBoundary_(const Element& element,
+ const SubControlVolumeFace& scvf,
+ const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
+ const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
+ const int localSubFaceIdx) const
{
- // Find out what boundary type is set on the lateral face
- const bool useZeroGradient = lateralFaceBoundaryTypes && (lateralFaceBoundaryTypes->isSymmetry()
- || lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx));
- const bool lateralFaceHasBJS = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex()));
- const bool lateralFaceHasDirichletVelocity = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex()));
- const Scalar velocitySelf = faceVars.velocitySelf();
-
// If there is a Dirichlet condition for the pressure we assume zero gradient for the velocity,
// so the velocity at the boundary equal to that on the scvf.
+ const bool useZeroGradient = lateralFaceBoundaryTypes && (lateralFaceBoundaryTypes->isSymmetry() || lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx));
if (useZeroGradient)
- return velocitySelf;
+ return faceVars_.velocitySelf();
+ const bool lateralFaceHasBJS = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex()));
if (lateralFaceHasBJS)
- return VelocityGradients::beaversJosephVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars,
+ return VelocityGradients::beaversJosephVelocityAtLateralScvf(elemVolVars_.gridVolVars().problem(), element, fvGeometry_, scvf, faceVars_,
currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
- else if (lateralFaceHasDirichletVelocity)
+ const bool lateralFaceHasDirichletVelocity = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex()));
+ if (lateralFaceHasDirichletVelocity)
{
// ________________
// ---------------o || frontal face of staggered half-control-volume
@@ -511,38 +392,27 @@ private:
// | || # o position at which the boundary conditions will be evaluated
// ---------------#
- const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(element);
- const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
-
+ const auto& lateralFace = fvGeometry_.scvf(scvf.insideScvIdx(), scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
const auto ghostFace = lateralFace.makeBoundaryFace(scvf.pairData(localSubFaceIdx).lateralStaggeredFaceCenter);
+ const auto& problem = elemVolVars_.gridVolVars().problem();
return problem.dirichlet(element, ghostFace)[Indices::velocity(scvf.directionIndex())];
}
- else
- {
- // Neumann conditions are not well implemented
- DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the boundary conditions for the momentum equations at global position " << fvGeometry.scvf(scvf.insideScvIdx(), scvf.pairData(localSubFaceIdx).localLateralFaceIdx).center());
- }
+
+ // Neumann conditions are not well implemented
+ DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the boundary conditions for the momentum equations at global position "
+ << fvGeometry_.scvf(scvf.insideScvIdx(), scvf.pairData(localSubFaceIdx).localLateralFaceIdx).center());
}
/*!
* \brief Return a velocity value from a boundary for which the boundary conditions have to be checked.
- *
- * \param problem The problem
- * \param scvf The SubControlVolumeFace that is normal to the boundary
- * \param localIdx The local index of the face that is on the boundary
- * \param boundaryElement The element that is on the boundary
- * \param parallelVelocity The velocity over scvf
*/
- static Scalar getParallelVelocityFromOppositeBoundary_(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const SubControlVolumeFace& scvf,
- const FaceVariables& faceVars,
- const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
- const int localOppositeSubFaceIdx)
+ Scalar getParallelVelocityFromOppositeBoundary_(const Element& element,
+ const SubControlVolumeFace& scvf,
+ const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
+ const int localOppositeSubFaceIdx) const
{
// A ghost subface at the boundary is created, featuring the location of the sub face's center
- const SubControlVolumeFace& lateralOppositeScvf = fvGeometry.scvf(scvf.insideScvIdx(), scvf.pairData(localOppositeSubFaceIdx).localLateralFaceIdx);
+ const auto& lateralOppositeScvf = fvGeometry_.scvf(scvf.insideScvIdx(), scvf.pairData(localOppositeSubFaceIdx).localLateralFaceIdx);
GlobalPosition center = scvf.pairData(localOppositeSubFaceIdx).lateralStaggeredFaceCenter + lateralOppositeScvf.center();
center *= 0.5;
@@ -558,11 +428,9 @@ private:
// Get the boundary types of the lateral opposite boundary face
+ const auto& problem = elemVolVars_.gridVolVars().problem();
const auto lateralOppositeFaceBoundaryTypes = problem.boundaryTypes(element, lateralOppositeScvf.makeBoundaryFace(center));
-
- return getParallelVelocityFromBoundary_(problem, element, fvGeometry, scvf,
- faceVars, currentScvfBoundaryTypes, lateralOppositeFaceBoundaryTypes,
- localOppositeSubFaceIdx);
+ return getParallelVelocityFromBoundary_(element, scvf, currentScvfBoundaryTypes, lateralOppositeFaceBoundaryTypes, localOppositeSubFaceIdx);
}
const Element& element_;