diff --git a/dumux/common/staggeredfvproblem.hh b/dumux/common/staggeredfvproblem.hh
index c66515c0f607e1e2f2b138d778f4c9cd99722301..a517417d81c0234901f11337bf1668effe789188 100644
--- a/dumux/common/staggeredfvproblem.hh
+++ b/dumux/common/staggeredfvproblem.hh
@@ -53,6 +53,13 @@ class StaggeredFVProblem : public FVProblem<TypeTag>
using Implementation = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
+
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
+ using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
+ using GridFaceVariables = typename GridVariables::GridFaceVariables;
+ using ElementFaceVariables = typename GridFaceVariables::LocalView;
+
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
@@ -106,6 +113,28 @@ public:
return asImp_().sourceAtPos(e.center());
}
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * This is the method for the case where the Neumann condition is
+ * potentially solution dependent
+ * \param elemFaceVars All face variables for the element
+ * \param scvf The sub control volume face
+ *
+ * Negative values mean influx.
+ * E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
+ */
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ // forward it to the interface with only the global position
+ return asImp_().neumannAtPos(scvf.ipGlobal());
+ }
+
/*!
* \brief Evaluate the initial value for a control volume.
*
diff --git a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
index 426657563cfc9ce92f8fac15456551a3391028f9..bc1c61fddd6a26444f168552b33434f9a7b3d035 100644
--- a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
@@ -100,11 +100,13 @@ public:
//! Constructor for a ghost face outside of the domain. Only needed to retrieve the center and scvIndices
FreeFlowStaggeredSubControlVolumeFace(const GlobalPosition& dofPosition,
const std::vector<GridIndexType>& scvIndices,
- const int dofIdx = -1)
+ const int dofIdx = -1,
+ const int scvfIndex = -1)
: center_(dofPosition),
- scvfIndex_(-1),
+ scvfIndex_(scvfIndex),
scvIndices_(scvIndices),
dofIdx_(dofIdx),
+ selfToOppositeDistance_(0.0),
isGhostFace_(true)
{}
diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
index 996bf049b2e75d090287fc194c2d2cfb5f007f45..5780abd862d1d06b90981f8321dd506dc0cfe53c 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
@@ -272,7 +272,7 @@ public:
* scvf
* ---------####### || and # staggered half-control-volume
* | || | current scvf
- * | || | # normal staggered faces over wich fluxes are calculated
+ * | || | # normal staggered sub faces over wich fluxes are calculated
* | || x~~~~> vel.Self
* | || | x dof position
* scvf | || |
@@ -291,22 +291,43 @@ public:
auto& faceVars = elemFaceVars[scvf];
const int numSubFaces = scvf.pairData().size();
- // Account for all sub-faces.
+ // Account for all sub faces.
for(int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx)
{
const auto eIdx = scvf.insideScvIdx();
+ // Get the face normal to the face the dof lives on. The staggered sub face conincides with half of this normal face.
const auto& normalFace = fvGeometry.scvf(eIdx, scvf.pairData()[localSubFaceIdx].localNormalFaceIdx);
- // Check if we have a symmetry boundary condition. If yes, the tangental part of the momentum flux can be neglected.
+ // Check if there is face/element parallel to our face of interest where the dof lives on. If there is no parallel neighbor,
+ // we are on a boundary where we have to check for boundary conditions.
if(!scvf.hasParallelNeighbor(localSubFaceIdx))
{
- // Use the ghost face to check if there is a symmetry boundary condition and skip any further steps if yes.
- const auto bcTypes = problem.boundaryTypes(element, makeParallelGhostFace_(scvf, localSubFaceIdx));
+ // Construct a temporary scvf which corresponds to the staggered sub face, featuring the location
+ // the sub faces's center.
+ auto localSubFaceCenter = scvf.pairData(localSubFaceIdx).virtualOuterParallelFaceDofPos - normalFace.center();
+ localSubFaceCenter *= 0.5;
+ localSubFaceCenter += normalFace.center();
+ const auto localSubFace = makeGhostFace_(normalFace, localSubFaceCenter);
+
+ // Retrieve the boundary types that correspond to the sub face.
+ const auto bcTypes = problem.boundaryTypes(element, localSubFace);
+
+ // Check if we have a symmetry boundary condition. If yes, the tangental part of the momentum flux can be neglected
+ // and we may skip any further calculations for the given sub face.
if(bcTypes.isSymmetry())
continue;
+
+ // Handle Neumann boundary conditions. No further calculations are then required for the given sub face.
+ if(bcTypes.isNeumann(Indices::velocity(scvf.directionIndex())))
+ {
+ const auto extrusionFactor = elemVolVars[normalFace.insideScvIdx()].extrusionFactor();
+ normalFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, localSubFace)[Indices::velocity(scvf.directionIndex())]
+ * extrusionFactor * normalFace.area() * 0.5;
+ continue;
+ }
}
- // If there is no symmetry boundary condition, proceed to calculate the tangential momentum flux.
+ // If there is no symmetry or Neumann boundary condition for the given sub face, proceed to calculate the tangential momentum flux.
if(enableInertiaTerms)
normalFlux += computeAdvectivePartOfNormalMomentumFlux_(problem, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx);
@@ -541,19 +562,19 @@ private:
private:
- //! helper functiuob to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
+ //! helper function to conveniently create a ghost face used to retrieve boundary values from the problem
SubControlVolumeFace makeGhostFace_(const SubControlVolumeFace& ownScvf, const GlobalPosition& pos) const
{
- return SubControlVolumeFace(pos, std::vector<unsigned int>{ownScvf.insideScvIdx(), ownScvf.outsideScvIdx()}, ownScvf.dofIndex());
+ return SubControlVolumeFace(pos, std::vector<unsigned int>{ownScvf.insideScvIdx(), ownScvf.outsideScvIdx()}, ownScvf.dofIndex(), ownScvf.index());
};
- //! helper functiuob to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
+ //! helper function to conveniently create a ghost face which is outside the domain, normal to the scvf of interest
SubControlVolumeFace makeNormalGhostFace_(const SubControlVolumeFace& ownScvf, const int localSubFaceIdx) const
{
return makeGhostFace_(ownScvf, ownScvf.pairData(localSubFaceIdx).virtualOuterNormalFaceDofPos);
};
- //! helper functiuob to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
+ //! helper function to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
SubControlVolumeFace makeParallelGhostFace_(const SubControlVolumeFace& ownScvf, const int localSubFaceIdx) const
{
return makeGhostFace_(ownScvf, ownScvf.pairData(localSubFaceIdx).virtualOuterParallelFaceDofPos);
diff --git a/dumux/freeflow/navierstokes/staggered/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh
index 9c1c37a7fe86a5da12d47d30b4ac21c19eb29daf..7ef4336a84df321e2a90f91a02b915c4c68cec91 100644
--- a/dumux/freeflow/navierstokes/staggered/localresidual.hh
+++ b/dumux/freeflow/navierstokes/staggered/localresidual.hh
@@ -233,8 +233,8 @@ protected:
{
if(bcTypes.isNeumann(eqIdx + cellCenterOffset))
{
- const auto extrusionFactor = 1.0; //TODO: get correct extrusion factor
- boundaryFlux[eqIdx] = problem.neumann(element, fvGeometry, elemVolVars, scvf)[eqIdx + cellCenterOffset]
+ const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
+ boundaryFlux[eqIdx] = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[eqIdx + cellCenterOffset]
* extrusionFactor * scvf.area();
}
}
@@ -293,11 +293,23 @@ protected:
residual = velocity - dirichletValue;
}
+ // handle Neumann boundary conditions
+ if(bcTypes.isNeumann(Indices::velocity(scvf.directionIndex())))
+ {
+ // set the given Neumann flux for the face on the boundary itself
+ const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
+ residual = problem.neumann(element, fvGeometry, elemVolVars, elementFaceVars, scvf)[Indices::velocity(scvf.directionIndex())]
+ * extrusionFactor * scvf.area();
+
+ // treat the remaining (normal) faces of the staggered control volume
+ FluxVariables fluxVars;
+ residual += fluxVars.computeNormalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elementFaceVars);
+ }
+
// For symmetry boundary conditions, there is no flow accross the boundary and
// we therefore treat it like a Dirichlet boundary conditions with zero velocity
if(bcTypes.isSymmetry())
{
- // const Scalar velocity = faceVars.faceVars(scvf.dofIndex()).velocity();
const Scalar velocity = elementFaceVars[scvf].velocitySelf();
const Scalar fixedValue = 0.0;
residual = velocity - fixedValue;