[freeflow] Implement Neumann BC for momentum

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.
      *

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)
     {}
 

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);

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;