[staggered][localresidual] Eval boundaries for faces directly in flux calculation

dumux/assembly/staggeredlocalresidual.hh

@@ -234,6 +234,8 @@ public:
     {
         if (!scvf.boundary())
             residual += asImp_().computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
+        else
+            residual += asImp_().computeBoundaryFluxForFace(problem, element, fvGeometry, scvf, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache);
     }
 
     //! Evaluate the source terms for a face residual
@@ -298,20 +300,6 @@ public:
         residual += storage;
     }
 
-    //! Evaluate the boundary conditions for a face residual
-    void evalBoundaryForFace(FaceResidualValue& residual,
-                             const Problem& problem,
-                             const Element& element,
-                             const FVElementGeometry& fvGeometry,
-                             const ElementVolumeVariables& elemVolVars,
-                             const ElementFaceVariables& elemFaceVars,
-                             const ElementBoundaryTypes& bcTypes,
-                             const ElementFluxVariablesCache& elemFluxVarsCache,
-                             const SubControlVolumeFace& scvf) const
-    {
-        asImp_().evalBoundaryForFace_(residual, problem, element, fvGeometry, scvf, elemVolVars, elemFaceVars, bcTypes, elemFluxVarsCache);
-    }
-
     //! If no solution has been set, we treat the problem as stationary.
     bool isStationary() const
     { return !timeLoop_; }

dumux/freeflow/navierstokes/staggered/localresidual.hh

@@ -181,8 +181,6 @@ public:
         return fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars);
     }
 
-protected:
-
     /*!
      * \brief Evaluate boundary conditions for a cell center dof
      */
@@ -242,15 +240,15 @@ protected:
      * \brief Evaluate boundary conditions for a face dof
      */
     template<class ElementBoundaryTypes>
-    void evalBoundaryForFace_(FaceResidual& residual,
-                              const Problem& problem,
-                              const Element& element,
-                              const FVElementGeometry& fvGeometry,
-                              const SubControlVolumeFace& scvf,
-                              const ElementVolumeVariables& elemVolVars,
-                              const ElementFaceVariables& elemFaceVars,
-                              const ElementBoundaryTypes& elemBcTypes,
-                              const ElementFluxVariablesCache& elemFluxVarsCache) const
+    void evalBoundaryForFace(FaceResidual& residual,
+                             const Problem& problem,
+                             const Element& element,
+                             const FVElementGeometry& fvGeometry,
+                             const SubControlVolumeFace& scvf,
+                             const ElementVolumeVariables& elemVolVars,
+                             const ElementFaceVariables& elemFaceVars,
+                             const ElementBoundaryTypes& elemBcTypes,
+                             const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         if (scvf.boundary())
         {
@@ -283,16 +281,47 @@ protected:
                 const Scalar fixedValue = 0.0;
                 residual = velocity - fixedValue;
             }
+        }
+    }
+
+    /*!
+     * \brief Evaluate boundary boundary fluxes for a face dof
+     */
+    template<class ElementBoundaryTypes>
+    FaceResidual computeBoundaryFluxForFace(const Problem& problem,
+                                            const Element& element,
+                                            const FVElementGeometry& fvGeometry,
+                                            const SubControlVolumeFace& scvf,
+                                            const ElementVolumeVariables& elemVolVars,
+                                            const ElementFaceVariables& elemFaceVars,
+                                            const ElementBoundaryTypes& elemBcTypes,
+                                            const ElementFluxVariablesCache& elemFluxVarsCache) const
+    {
+        FaceResidual result(0.0);
+
+        if (scvf.boundary())
+        {
+            // handle the actual boundary conditions:
+            const auto bcTypes = problem.boundaryTypes(element, scvf);
+            if (bcTypes.isNeumann(Indices::velocity(scvf.directionIndex())))
+            {
+                // the source term has already been accounted for, here we
+                // add a given Neumann flux for the face on the boundary itself ...
+                const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
+                result = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())]
+                                         * extrusionFactor * scvf.area();
+
+                // ... and treat the fluxes of the remaining (frontal and lateral) faces of the staggered control volume
+                result += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
+            }
             else if(bcTypes.isDirichlet(Indices::pressureIdx))
             {
                 // if none of the above conditions apply, we are at an "fixed pressure" boundary for which the resdiual of the momentum balance needs to be assembled
                 // as if it where inside the domain and not on the boundary (source term has already been acounted for)
-                residual += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
+                result = computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
             }
-            else
-                DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the boundary conditions "
-                           "for the momentum equations at global position " << scvf.center());
         }
+        return result;
     }
 
 private:

dumux/multidomain/subdomainstaggeredlocalassembler.hh

@@ -249,7 +249,9 @@ public:
         if (!this->assembler().isStationaryProblem())
             residual += evalLocalStorageResidualForFace(scvf);
 
-        this->localResidual().evalBoundaryForFace(residual, problem(), this->element(), this->fvGeometry(), elemVolVars, elemFaceVars, this->elemBcTypes(), this->elemFluxVarsCache(), scvf);
+        this->localResidual().evalDirichletBoundariesForFace(residual, this->problem(), this->element(),
+                                                             this->fvGeometry(), scvf, elemVolVars, elemFaceVars,
+                                                             this->elemBcTypes(), this->elemFluxVarsCache());
 
         return residual;
     }