Merge branch 'feature/improve-beavers-joseph' into 'master'

Feature/improve beavers joseph

See merge request !1833

CHANGELOG.md

@@ -5,6 +5,8 @@ Differences Between DuMuX 3.2 and DuMuX 3.1
 
 - __Radially symmetric problems__: We now have support for radially symmetric problems (disc, ball, toroid). The support comes in form of wrappers for sub control volumes and faces that overload the respective `volume()` and `area()` function turning a 1d or 2d problem into a 2d or 3d radially symmetric problem.
 
+- __Improvements of Beavers-Joseph(-Saffman) condition for the free flow model__: The naming for handling BJ(-S) boundary conditions has been adapted from `isBJS()` to `isBeaversJoseph()` / `setBJS()` to `setBeaversJoseph()`. In order to consider the velocity within the porous medium, the old `velocityPorousMedium(element, scvf)` method (returning a Scalar) has been renamed to `porousMediumVelocity(element, scvf)` (returning a velocity vector). The latter defaults to `VelocityVector(0.0)`.
+
 ### Immediate interface changes not allowing/requiring a deprecation period
 
 - Remove `Grid.HeapSize` as dune-ugrid removed the according feature as well.

dumux/discretization/staggered/freeflow/boundarytypes.hh

@@ -31,7 +31,7 @@ namespace Dumux {
 
 /*!
  * \ingroup StaggeredDiscretization
- * \brief Class to specify the type of a boundary for the staggered Navier-Stokes model.
+ * \brief Class to specify the type of a boundary condition for the staggered Navier-Stokes model.
  */
 template <int numEq>
 class StaggeredFreeFlowBoundaryTypes : public Dumux::BoundaryTypes<numEq>
@@ -41,7 +41,6 @@ class StaggeredFreeFlowBoundaryTypes : public Dumux::BoundaryTypes<numEq>
 public:
     StaggeredFreeFlowBoundaryTypes()
     {
-
         for (int eqIdx=0; eqIdx < numEq; ++eqIdx)
             resetEq(eqIdx);
     }
@@ -55,7 +54,7 @@ public:
 
         boundaryInfo_[eqIdx].visited = false;
         boundaryInfo_[eqIdx].isSymmetry = false;
-        boundaryInfo_[eqIdx].isBJS = false;
+        boundaryInfo_[eqIdx].isBeaversJoseph = false;
     }
 
     /*!
@@ -99,11 +98,19 @@ public:
      * \brief Set a boundary condition for a single equation to
      *        Beavers-Joseph-Saffman (special case of Dirichlet b.c.).
      */
+    [[deprecated("Use setBeaversJoseph instead. Will be removed after 3.2")]]
     void setBJS(int eqIdx)
+    { setBeaversJoseph(eqIdx); }
+
+    /*!
+     * \brief Set a boundary condition for a single equation to
+     *        Beavers-Joseph(-Saffmann) (special case of Dirichlet b.c.).
+     */
+    void setBeaversJoseph(unsigned eqIdx)
     {
         resetEq(eqIdx);
         boundaryInfo_[eqIdx].visited = true;
-        boundaryInfo_[eqIdx].isBJS = true;
+        boundaryInfo_[eqIdx].isBeaversJoseph = true;
     }
 
     /*!
@@ -112,28 +119,45 @@ public:
      *
      * \param eqIdx The index of the equation
      */
+    [[deprecated("Use isBeaversJoseph instead. Will be removed after 3.2")]]
     bool isBJS(unsigned eqIdx) const
-    {
-        return boundaryInfo_[eqIdx].isBJS;
-    }
+    { return isBeaversJoseph(eqIdx); }
+
+    /*!
+     * \brief Returns true if an equation is used to specify a
+     *        Beavers-Joseph(-Saffman) boundary condition.
+     *
+     * \param eqIdx The index of the equation
+     */
+    bool isBeaversJoseph(unsigned eqIdx) const
+    { return boundaryInfo_[eqIdx].isBeaversJoseph; }
 
     /*!
      * \brief Returns true if some equation is used to specify a
      *        Beavers-Joseph-Saffman boundary condition.
      */
+    [[deprecated("Use hasBeaversJoseph instead. Will be removed after 3.2")]]
     bool hasBJS() const
+    { return hasBeaversJoseph(); }
+
+    /*!
+     * \brief Returns true if some equation is used to specify a
+     *        Beavers-Joseph(-Saffman) boundary condition.
+     */
+    bool hasBeaversJoseph() const
     {
         for (int i = 0; i < numEq; ++i)
-            if (boundaryInfo_[i].isBJS)
+            if (boundaryInfo_[i].isBeaversJoseph)
                 return true;
         return false;
     }
 
 protected:
-    struct StaggeredFreeFlowBoundaryInfo {
+    struct StaggeredFreeFlowBoundaryInfo
+    {
         bool visited;
         bool isSymmetry;
-        bool isBJS;
+        bool isBeaversJoseph;
     };
 
     std::array<StaggeredFreeFlowBoundaryInfo, numEq> boundaryInfo_;

dumux/freeflow/navierstokes/problem.hh

@@ -24,15 +24,11 @@
 #ifndef DUMUX_NAVIERSTOKES_PROBLEM_HH
 #define DUMUX_NAVIERSTOKES_PROBLEM_HH
 
-#include <dune/common/deprecated.hh>
-
 #include <dune/common/exceptions.hh>
 #include <dumux/common/properties.hh>
 #include <dumux/common/staggeredfvproblem.hh>
 #include <dumux/discretization/method.hh>
 
-#include "model.hh"
-
 namespace Dumux {
 
 //! The implementation is specialized for the different discretizations
@@ -88,6 +84,7 @@ class NavierStokesProblem : public NavierStokesParentProblem<TypeTag>
       };
 
     using GlobalPosition = typename SubControlVolumeFace::GlobalPosition;
+    using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
     using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
 
 public:
@@ -150,7 +147,6 @@ public:
         sol[GridGeometry::faceIdx()][scvf.dofIndex()][0] = initSol[Indices::velocity(scvf.directionIndex())];
     }
 
-
     /*!
      * \brief An additional drag term can be included as source term for the momentum balance
      *        to mimic 3D flow behavior in 2D:
@@ -208,9 +204,7 @@ public:
     }
 
     /*!
-     * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
-     *
-     * This member function must be overloaded in the problem implementation, if the BJS boundary condition is used.
+     * \brief Returns the beta value, or the alpha value divided by the square root of the intrinsic permeability.
      */
     Scalar betaBJ(const Element& element, const SubControlVolumeFace& scvf) const
     {
@@ -220,22 +214,25 @@ public:
 
     /*!
      * \brief Returns the velocity in the porous medium (which is 0 by default according to Saffmann).
+     * \note This method is deprecated. Use porousMediumVelocity(element, scvf) instead, returning a velocity vector. Will be removed after 3.2
      */
     Scalar velocityPorousMedium(const Element& element, const SubControlVolumeFace& scvf) const
-    { return 0.0; }
+    {
+        // Redirect to helper method to avoid spurious deprecation warnings. TODO: Remove after 3.2!
+        return deprecatedVelocityPorousMedium_();
+    }
 
-    //! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph-Saffman condition is used
-    DUNE_DEPRECATED_MSG("Use beaversJosephVelocity(element, scv, faceOnPorousBoundary, velocitySelf, tangentialVelocityGradient) instead")
-    const Scalar bjsVelocity(const Element& element,
-                             const SubControlVolume& scv,
-                             const SubControlVolumeFace& faceOnPorousBoundary,
-                             const Scalar velocitySelf) const
+    /*!
+     * \brief Returns the velocity in the porous medium (which is 0 by default according to Saffmann).
+     */
+    VelocityVector porousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const
     {
-        // assume tangential velocity gradient of zero and
-        return beaversJosephVelocity(element, scv, faceOnPorousBoundary, velocitySelf, 0.0);
+        // TODO: return VelocityVector(0.0) after 3.2!
+        return VelocityVector(getVelPM_(element, scvf));
     }
 
     //! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph condition is used
+    [[deprecated("Use beaversJosephVelocity(element, scv, ownScvf, faceOnPorousBoundary, velocitySelf, tangentialVelocityGradient) instead. Will be removed after 3.2")]]
     const Scalar beaversJosephVelocity(const Element& element,
                                        const SubControlVolume& scv,
                                        const SubControlVolumeFace& faceOnPorousBoundary,
@@ -249,8 +246,55 @@ public:
         return (tangentialVelocityGradient*distance + asImp_().velocityPorousMedium(element,faceOnPorousBoundary)*betaBJ*distance + velocitySelf) / (betaBJ*distance + 1.0);
     }
 
+    //! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph condition is used
+    const Scalar beaversJosephVelocity(const Element& element,
+                                       const SubControlVolume& scv,
+                                       const SubControlVolumeFace& ownScvf,
+                                       const SubControlVolumeFace& faceOnPorousBoundary,
+                                       const Scalar velocitySelf,
+                                       const Scalar tangentialVelocityGradient) const
+    {
+        // du/dy + dv/dx = alpha/sqrt(K) * (u_boundary-uPM)
+        // beta = alpha/sqrt(K)
+        const Scalar betaBJ = asImp_().betaBJ(element, faceOnPorousBoundary);
+        const Scalar distanceNormalToBoundary = (faceOnPorousBoundary.center() - scv.center()).two_norm();
+
+        // create a unit normal vector oriented in positive coordinate direction
+        GlobalPosition orientation = ownScvf.unitOuterNormal();
+        orientation[ownScvf.directionIndex()] = 1.0;
+
+        return (tangentialVelocityGradient*distanceNormalToBoundary
+              + asImp_().porousMediumVelocity(element, faceOnPorousBoundary) * orientation * betaBJ * distanceNormalToBoundary
+              + velocitySelf) / (betaBJ*distanceNormalToBoundary + 1.0);
+    }
+
 private:
 
+    // Auxiliary method handling deprecation warnings. TODO: Remove after 3.2!
+    Scalar getVelPM_(const Element& element, const SubControlVolumeFace& scvf) const
+    {
+        // Check if the user problem implements the deprecated velocityPorousMedium method
+        static constexpr bool implHasVelocityPorousMedium = !std::is_same<decltype(&Implementation::velocityPorousMedium), decltype(&NavierStokesProblem::velocityPorousMedium)>::value;
+        // This check would always trigger a spurious deprecation warning if the base class' (NavierStokesProblem) velocityPorousMedium method was equipped with a deprecation warning.
+        // This is why we need another level of redirection there.
+
+        // Forward either to user impl (thereby raising a deprecation warning) or return 0.0 by default
+        return deprecationHelper_(element, scvf, std::integral_constant<bool, implHasVelocityPorousMedium>{});
+    }
+
+    [[deprecated("\nvelocityPorousMedium(element, scvf) is deprecated. Use porousMediumVelocity(element, scvf) instead, returning a velocity vector. Will be removed after 3.2")]]
+    Scalar deprecationHelper_(const Element& element, const SubControlVolumeFace& scvf, std::true_type) const
+    { return asImp_().velocityPorousMedium(element, scvf); }
+
+    // Return 0.0 by default. TODO: Remove this after 3.2
+    Scalar deprecationHelper_(const Element& element, const SubControlVolumeFace& scvf, std::false_type) const
+    { return 0.0; }
+
+    // Auxiliary method to trigger a deprecation warning.
+    [[deprecated("\nvelocityPorousMedium(element, scvf) is deprecated. Use porousMediumVelocity(element, scvf) instead, returning a velocity vector. Will be removed after 3.2")]]
+    Scalar deprecatedVelocityPorousMedium_() const
+    { return 0.0; }
+
     //! Returns the implementation of the problem (i.e. static polymorphism)
     Implementation &asImp_()
     { return *static_cast<Implementation *>(this); }

dumux/freeflow/navierstokes/staggered/fluxvariables.hh

@@ -331,7 +331,7 @@ public:
                 // It is not clear how to evaluate the BJ condition here.
                 // For symmetry reasons, our own scvf should then have the same Neumann flux as the lateral face.
                 // TODO: We should clarify if this is the correct approach.
-                if (currentScvfBoundaryTypes->isBJS(Indices::velocity(lateralScvf.directionIndex())) && lateralFaceBoundaryTypes &&
+                if (currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex())) && lateralFaceBoundaryTypes &&
                     lateralFaceBoundaryTypes->isNeumann(Indices::velocity(scvf.directionIndex())))
                 {
                     const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
@@ -375,7 +375,7 @@ public:
                 std::bitset<3> admittableBcTypes;
                 admittableBcTypes.set(0, lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx));
                 admittableBcTypes.set(1, lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex())));
-                admittableBcTypes.set(2, lateralFaceBoundaryTypes->isBJS(Indices::velocity(scvf.directionIndex())));
+                admittableBcTypes.set(2, lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())));
                 if (admittableBcTypes.count() != 1)
                 {
                     DUNE_THROW(Dune::InvalidStateException, "Invalid boundary conditions for lateral scvf "
@@ -502,7 +502,7 @@ private:
                     const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
                     return problem.dirichlet(element, scvf.makeBoundaryFace(lateralBoundaryFacePos))[Indices::velocity(lateralFace.directionIndex())];
                 }
-                else if (bcTypes.isBJS(Indices::velocity(lateralFace.directionIndex())))
+                else if (bcTypes.isBeaversJoseph(Indices::velocity(lateralFace.directionIndex())))
                 {
                     return VelocityGradients::beaversJosephVelocityAtCurrentScvf(problem, element, fvGeometry, scvf,  faceVars,
                                                                                  currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
@@ -573,7 +573,7 @@ private:
         {
             if (!scvf.boundary() ||
                 currentScvfBoundaryTypes->isDirichlet(Indices::velocity(lateralFace.directionIndex())) ||
-                currentScvfBoundaryTypes->isBJS(Indices::velocity(lateralFace.directionIndex())))
+                currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralFace.directionIndex())))
             {
                 const Scalar velocityGrad_ji = VelocityGradients::velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
                 // Account for the orientation of the staggered normal face's outer normal vector.

dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh

@@ -582,7 +582,7 @@ private:
         // 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->isBJS(Indices::velocity(scvf.directionIndex()));
+        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();
 

dumux/freeflow/navierstokes/staggered/velocitygradients.hh

@@ -121,7 +121,7 @@ public:
                 const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
                 return problem.dirichlet(element, lateralScvf.makeBoundaryFace(lateralBoundaryFacePos))[Indices::velocity(scvf.directionIndex())];
             }
-            else if (lateralFaceBoundaryTypes->isBJS(Indices::velocity(scvf.directionIndex())))
+            else if (lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())))
             {
                 return beaversJosephVelocityAtLateralScvf(problem, element, fvGeometry, scvf,  faceVars,
                                                           currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
@@ -200,7 +200,7 @@ public:
                 const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
                 return problem.dirichlet(element, scvf.makeBoundaryFace(lateralBoundaryFacePos))[Indices::velocity(lateralScvf.directionIndex())];
             }
-            else if (currentScvfBoundaryTypes->isBJS(Indices::velocity(lateralScvf.directionIndex())))
+            else if (currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex())))
             {
                 return beaversJosephVelocityAtCurrentScvf(problem, element, fvGeometry, scvf,  faceVars,
                                                           currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
@@ -265,15 +265,18 @@ public:
             if (lateralScvf.boundary())
             {
                 if (lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
-                    lateralFaceBoundaryTypes->isBJS(Indices::velocity(scvf.directionIndex())))
+                    lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())))
                     return 0.0;
             }
 
             return velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
         }();
 
-        return problem.beaversJosephVelocity(element, fvGeometry.scv(scvf.insideScvIdx()),
-                                             scvf, innerLateralVelocity,
+        return problem.beaversJosephVelocity(element,
+                                             fvGeometry.scv(scvf.insideScvIdx()),
+                                             lateralScvf,
+                                             scvf, /*on boundary*/
+                                             innerLateralVelocity,
                                              tangentialVelocityGradient);
     }
 
@@ -322,16 +325,18 @@ public:
             if (scvf.boundary())
             {
                 if (currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
-                    currentScvfBoundaryTypes->isBJS(Indices::velocity(lateralScvf.directionIndex())))
+                    currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex())))
                     return 0.0;
             }
 
             return velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
         }();
 
-
-        return problem.beaversJosephVelocity(element, fvGeometry.scv(scvf.insideScvIdx()),
-                                             lateralScvf, innerParallelVelocity,
+        return problem.beaversJosephVelocity(element,
+                                             fvGeometry.scv(scvf.insideScvIdx()),
+                                             scvf,
+                                             lateralScvf, /*on boundary*/
+                                             innerParallelVelocity,
                                              tangentialVelocityGradient);
     }
 

dumux/freeflow/rans/problem.hh

@@ -337,14 +337,14 @@ public:
 
                     // adapt calculations for Beavers-Joseph-Saffman condition
                     unsigned int normalNormDim = lateralFace.directionIndex();
-                    if (lateralFace.boundary() && (asImp_().boundaryTypes(element, lateralFace).isBJS(Indices::velocity(velIdx))))
+                    if (lateralFace.boundary() && (asImp_().boundaryTypes(element, lateralFace).isBeaversJoseph(Indices::velocity(velIdx))))
                     {
                         unsigned int neighborIdx = neighborIdx_[elementIdx][normalNormDim][0];
                         if (lateralFace.center()[normalNormDim] < cellCenter_[elementIdx][normalNormDim])
                             neighborIdx = neighborIdx_[elementIdx][normalNormDim][1];
 
                         const SubControlVolume& scv = fvGeometry.scv(scvf.insideScvIdx());
-                        bjsVelocityAverage[normalNormDim] += ParentType::beaversJosephVelocity(element, scv, lateralFace, velocity_[elementIdx][velIdx], 0.0);
+                        bjsVelocityAverage[normalNormDim] += ParentType::beaversJosephVelocity(element, scv, scvf, lateralFace, velocity_[elementIdx][velIdx], 0.0);
                         if (bjsNumFaces[normalNormDim] > 0 && neighborIdx != bjsNeighbor[normalNormDim])
                             DUNE_THROW(Dune::InvalidStateException, "Two different neighborIdx should not occur");
                         bjsNeighbor[normalNormDim] = neighborIdx;

test/multidomain/boundary/stokesdarcy/1p2c_1p2c/diffusionlawcomparison/problem_stokes.hh

@@ -198,7 +198,7 @@ public:
             values.setCouplingNeumann(Indices::conti0EqIdx);
             values.setCouplingNeumann(Indices::conti0EqIdx + 1);
             values.setCouplingNeumann(Indices::momentumYBalanceIdx);
-            values.setBJS(Indices::momentumXBalanceIdx);
+            values.setBeaversJoseph(Indices::momentumXBalanceIdx);
         }
 
         return values;

test/multidomain/boundary/stokesdarcy/1p2c_1p2c/problem_stokes.hh

@@ -209,7 +209,7 @@ public:
             values.setCouplingNeumann(Indices::conti0EqIdx);
             values.setCouplingNeumann(Indices::conti0EqIdx + 1);
             values.setCouplingNeumann(Indices::momentumYBalanceIdx);
-            values.setBJS(Indices::momentumXBalanceIdx);
+            values.setBeaversJoseph(Indices::momentumXBalanceIdx);
         }
 
         return values;

test/multidomain/boundary/stokesdarcy/1p2c_2p2c/problem_stokes.hh

@@ -207,7 +207,7 @@ public:
             values.setCouplingNeumann(Indices::conti0EqIdx);
             values.setCouplingNeumann(Indices::conti0EqIdx + 1);
             values.setCouplingNeumann(Indices::momentumYBalanceIdx);
-            values.setBJS(Indices::momentumXBalanceIdx);
+            values.setBeaversJoseph(Indices::momentumXBalanceIdx);
         }
         return values;
     }

test/multidomain/boundary/stokesdarcy/1p3c_1p3c/problem_stokes.hh

@@ -189,7 +189,7 @@ public:
             values.setNeumann(Indices::conti0EqIdx+1);
             values.setNeumann(Indices::conti0EqIdx+2);
             values.setNeumann(Indices::momentumYBalanceIdx);
-            values.setBJS(Indices::velocityXIdx);
+            values.setBeaversJoseph(Indices::velocityXIdx);
         }
 
         return values;

test/multidomain/boundary/stokesdarcy/1p_1p/problem_stokes.hh

@@ -190,7 +190,7 @@ public:
         {
             values.setCouplingNeumann(Indices::conti0EqIdx);
             values.setCouplingNeumann(Indices::momentumYBalanceIdx);
-            values.setBJS(Indices::momentumXBalanceIdx);
+            values.setBeaversJoseph(Indices::momentumXBalanceIdx);
         }
 
         return values;

test/multidomain/boundary/stokesdarcy/1p_2p/problem_stokes.hh

@@ -153,7 +153,7 @@ public:
         {
             values.setCouplingNeumann(Indices::conti0EqIdx);
             values.setCouplingNeumann(Indices::momentumYBalanceIdx);
-            values.setBJS(Indices::momentumXBalanceIdx);
+            values.setBeaversJoseph(Indices::momentumXBalanceIdx);
         }
         else
         {