Merge branch 'feature/navier_stokes_ctd' into 'master'

Navier slip condition as an option for interface condition

Closes #9

See merge request !14

appl/coupled/stokes_darcy/beaversjoseph/params.input

@@ -29,7 +29,8 @@ NumericDifference.BaseEpsilon = 1e-8
 [Darcy]
 Problem.Name = darcy
 SpatialParams.Permeability = 1e-4
-SpatialParams.AlphaBeaversJoseph = 1.0
+SpatialParams.Porosity = 1.0
+SpatialParams.NavierSlipLength = 0.01
 
 [Vtk]
 WriteFaceData = false

appl/coupled/stokes_darcy/beaversjoseph/problem_stokes.hh

@@ -110,7 +110,7 @@ public:
         {
             values.setCouplingNeumann(Indices::conti0EqIdx);
             values.setCouplingNeumann(scvf.directionIndex());
-            values.setBeaversJoseph(1 - scvf.directionIndex());
+            values.setNavierSlip(1 - scvf.directionIndex());
         }
         else
             values.setAllSymmetry();
@@ -176,6 +176,10 @@ public:
     Scalar alphaBJ(const SubControlVolumeFace& scvf) const
     { return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center()); }
 
+    //! Returns the slip length value required as a parameter for the Navier-Slip boundary condition
+    Scalar navierSlipLength(const SubControlVolumeFace& scvf) const
+    { return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().navierSlipLengthAtPos(scvf.center()); }
+
     //! Set the coupling manager
     void setCouplingManager(std::shared_ptr<CouplingManager> cm)
     { couplingManager_ = cm; }

appl/coupled/stokes_darcy/beaversjoseph/spatialparams.hh

@@ -56,6 +56,8 @@ public:
     {
         permeability_ = getParam<Scalar>("Darcy.SpatialParams.Permeability");
         alphaBJ_ = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
+        porosity_ = getParam<Scalar>("Darcy.SpatialParams.Porosity");
+        navierSlipLength_ = getParam<Scalar>("Darcy.SpatialParams.NavierSlipLength");
     }
 
     /*!
@@ -72,7 +74,7 @@ public:
      * \param globalPos The global position
      */
     Scalar porosityAtPos(const GlobalPosition& globalPos) const
-    { return 0.4; }
+    { return porosity_; }
 
     /*! \brief Defines the Beavers-Joseph coefficient in [-].
      *
@@ -81,10 +83,18 @@ public:
     Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const
     { return alphaBJ_; }
 
+    /*! \brief Defines the Navier-Slip coefficient in [-].
+     *
+     * \param globalPos The global position
+     */
+    Scalar navierSlipLengthAtPos(const GlobalPosition& globalPos) const
+    { return navierSlipLength_; }
 
 private:
     Scalar permeability_;
     Scalar alphaBJ_;
+    Scalar porosity_;
+    Scalar navierSlipLength_;
 };
 
 } // end namespace Dumux

dumux/freeflow/CMakeLists.txt

+#install headers
+install(
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/)

dumux/freeflow/navierstokes/CMakeLists.txt

+#install headers
+install(FILES
+boundarytypes.hh problem.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/navierstokes/)

dumux/freeflow/navierstokes/boundarytypes.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup NavierStokesModel
+ * \copydoc Dumux::NavierStokesBoundaryTypes
+ */
+#ifndef FREEFLOW_NAVIERSTOKES_BOUNDARY_TYPES_HH
+#define FREEFLOW_NAVIERSTOKES_BOUNDARY_TYPES_HH
+
+#include <dumux/common/boundarytypes.hh>
+#include <dumux/common/typetraits/typetraits.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief Class to specify the type of a boundary condition for the Navier-Stokes model.
+ */
+template <int numEq>
+class NavierStokesBoundaryTypes : public BoundaryTypes<numEq>
+{
+    using ParentType = BoundaryTypes<numEq>;
+
+public:
+    NavierStokesBoundaryTypes()
+    {
+        for (int eqIdx=0; eqIdx < numEq; ++eqIdx)
+            resetEq(eqIdx);
+    }
+
+    /*!
+     * \brief Reset the boundary types for one equation.
+     */
+    void resetEq(const int eqIdx)
+    {
+        ParentType::resetEq(eqIdx);
+
+        boundaryInfo_[eqIdx].isSymmetry = false;
+        boundaryInfo_[eqIdx].isBeaversJoseph = false;
+        boundaryInfo_[eqIdx].isOutflow = false;
+        boundaryInfo_[eqIdx].isNavierSlip = false;
+    }
+
+    /*!
+     * \brief Sets a symmetry boundary condition for all equations
+     */
+    void setAllSymmetry()
+    {
+        for (int eqIdx=0; eqIdx < numEq; ++eqIdx)
+        {
+            resetEq(eqIdx);
+            boundaryInfo_[eqIdx].isSymmetry = true;
+        }
+    }
+
+    /*!
+     * \brief Returns true if the there is a symmetry boundary condition
+     */
+    bool isSymmetry() const
+    { return boundaryInfo_[0].isSymmetry; }
+
+    /*!
+     * \brief  Prevent setting all boundary conditions to Dirichlet.
+     */
+    template<class T = void>
+    void setAllDirichlet()
+    {
+        static_assert(AlwaysFalse<T>::value, "Setting all boundary types to Dirichlet not permitted!");
+    }
+
+    /*!
+     * \brief  Prevent setting all boundary conditions to Neumann.
+     */
+    template<class T = void>
+    void setAllNeumann()
+    {
+        static_assert(AlwaysFalse<T>::value, "Setting all boundary types to Neumann not permitted!");
+    }
+
+    /*!
+     * \brief Set a boundary condition for a single equation to
+     *        Beavers-Joseph(-Saffmann) (special case of Dirichlet b.c.).
+     */
+    void setBeaversJoseph(const int eqIdx)
+    {
+        resetEq(eqIdx);
+        boundaryInfo_[eqIdx].isBeaversJoseph = true;
+    }
+
+    /*!
+     * \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(const int eqIdx) const
+    { return boundaryInfo_[eqIdx].isBeaversJoseph; }
+
+    /*!
+     * \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].isBeaversJoseph)
+                return true;
+        return false;
+    }
+
+    /*!
+     * \brief Set a boundary condition for a single equation to
+     *        Beavers-Joseph(-Saffmann) (special case of Dirichlet b.c.).
+     */
+    void setNavierSlip(const int eqIdx)
+    {
+        resetEq(eqIdx);
+        boundaryInfo_[eqIdx].isNavierSlip = true;
+    }
+
+    /*!
+     * \brief Returns true if an equation is used to specify a
+     *        Beavers-Joseph(-Saffman) boundary condition.
+     *
+     * \param eqIdx The index of the equation
+     */
+    bool isNavierSlip(const int eqIdx) const
+    { return boundaryInfo_[eqIdx].isNavierSlip; }
+
+    /*!
+     * \brief Returns true if some equation is used to specify a
+     *        Beavers-Joseph(-Saffman) boundary condition.
+     */
+    bool hasNavierSlip() const
+    {
+        for (int i = 0; i < numEq; ++i)
+            if (boundaryInfo_[i].isNavierSlip)
+                return true;
+        return false;
+    }
+
+
+    /*!
+     * \brief Set an outflow boundary condition
+     */
+    void setOutflow(const int eqIdx)
+    {
+        resetEq(eqIdx);
+        boundaryInfo_[eqIdx].isOutflow = true;
+    }
+
+    /*!
+     * \brief Returns true if an outflow boundary condition was set
+     * \param eqIdx The index of the equation
+     */
+    bool isOutflow(const int eqIdx) const
+    { return boundaryInfo_[eqIdx].isOutflow; }
+
+    /*!
+     * \brief Returns true if some equation has an outflow boundary condition
+     */
+    bool hasOutflow() const
+    {
+        for (int i = 0; i < numEq; ++i)
+            if (boundaryInfo_[i].isOutflow)
+                return true;
+        return false;
+    }
+
+protected:
+    //! use bitfields to minimize the size
+    struct NavierStokesBoundaryInfo
+    {
+        bool isSymmetry : 1;
+        bool isOutflow : 1;
+        bool isBeaversJoseph : 1;
+        bool isNavierSlip : 1;
+    };
+
+    std::array<NavierStokesBoundaryInfo, numEq> boundaryInfo_;
+};
+
+} // end namespace Dumux
+
+#endif

dumux/freeflow/navierstokes/staggered/CMakeLists.txt

+#install headers
+install(FILES
+fluxvariables.hh velocitygradients.hh staggeredupwindhelper.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/navierstokes/staggered/)