From 02c4522a02bbf51ccb16ad68e006416d003828f8 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Mon, 12 Jun 2017 17:35:41 +0200
Subject: [PATCH] [staggered] Introduce new boundary types

* Allow to set a fixed value for a cell at the boundary directly
  instead of setting a value for the face
---
 .../staggered/globalvolumevariables.hh        |  2 +-
 dumux/freeflow/staggered/boundarytypes.hh     | 98 +++++++++++++++++++
 dumux/freeflow/staggered/localresidual.hh     | 36 +++++--
 dumux/freeflow/staggered/propertydefaults.hh  |  6 ++
 dumux/freeflow/staggerednc/localresidual.hh   | 36 +++++++
 .../staggered/closedsystemtestproblem.hh      |  2 +-
 test/freeflow/staggered/doneatestproblem.hh   |  2 +-
 .../staggerednc/densityflowproblem.hh         |  2 +-
 8 files changed, 173 insertions(+), 11 deletions(-)
 create mode 100644 dumux/freeflow/staggered/boundarytypes.hh

diff --git a/dumux/discretization/staggered/globalvolumevariables.hh b/dumux/discretization/staggered/globalvolumevariables.hh
index 1510817e64..6365d63688 100644
--- a/dumux/discretization/staggered/globalvolumevariables.hh
+++ b/dumux/discretization/staggered/globalvolumevariables.hh
@@ -107,7 +107,7 @@ public:
 
                 for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
                 {
-                    if(bcTypes.isDirichlet(eqIdx))
+                    if(bcTypes.isDirichlet(eqIdx) || bcTypes.isDirichletCell(eqIdx))
                         boundaryPriVars[cellCenterIdx][eqIdx] = problem.dirichlet(element, scvf)[cellCenterIdx][eqIdx];
                     else if(bcTypes.isNeumann(eqIdx) || bcTypes.isOutflow(eqIdx))
                         boundaryPriVars[cellCenterIdx][eqIdx] = sol[cellCenterIdx][scvf.insideScvIdx()][eqIdx];
diff --git a/dumux/freeflow/staggered/boundarytypes.hh b/dumux/freeflow/staggered/boundarytypes.hh
new file mode 100644
index 0000000000..a24d5cc778
--- /dev/null
+++ b/dumux/freeflow/staggered/boundarytypes.hh
@@ -0,0 +1,98 @@
+// -*- 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 2 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
+ * \brief Class to specify the type of a boundary for the staggered Navier-Stokes model.
+ */
+#ifndef STAGGERED_FREEFLOW_BOUNDARY_TYPES_HH
+#define STAGGERED_FREEFLOW_BOUNDARY_TYPES_HH
+
+#include <dumux/common/boundarytypes.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup BC
+ * \brief Class to specify the type of a boundary for the staggered Navier-Stokes model.
+ */
+template <int numEq>
+class StaggeredFreeFlowBoundaryTypes : public Dumux::BoundaryTypes<numEq>
+{
+    using ParentType = Dumux::BoundaryTypes<numEq>;
+
+public:
+    StaggeredFreeFlowBoundaryTypes()
+    {
+
+        for (int eqIdx=0; eqIdx < numEq; ++eqIdx)
+            resetEq(eqIdx);
+    }
+
+
+    /*!
+     * \brief Reset the boundary types for one equation.
+     */
+    void resetEq(int eqIdx)
+    {
+        ParentType::resetEq(eqIdx);
+
+        boundaryInfo_[eqIdx].visited = false;
+        boundaryInfo_[eqIdx].isDirichletCell = false;
+    }
+
+
+    /*!
+     * \brief Sets a fixed Dirichlet value for a cell (such as pressure) at the boundary.
+     *        This is a provisional alternative to setting the Dirichlet value on the boundary directly.
+     *
+     * \param eqIdx The index of the equation which should used to set
+     *              the Dirichlet condition
+     */
+    void setDirichletCell(int eqIdx)
+    {
+        resetEq(eqIdx);
+        boundaryInfo_[eqIdx].visited = true;
+        boundaryInfo_[eqIdx].isDirichletCell = true;
+
+        Valgrind::SetDefined(boundaryInfo_[eqIdx]);
+    }
+
+    /*!
+     * \brief Returns true if an equation is used to specify a
+     *        Dirichlet condition.
+     *
+     * \param eqIdx The index of the equation
+     */
+    bool isDirichletCell(unsigned eqIdx) const
+    { return boundaryInfo_[eqIdx].isDirichletCell; }
+
+
+protected:
+    struct StaggeredFreeFlowBoundaryInfo {
+        bool visited;
+        bool isDirichletCell;
+    };
+
+    std::array<StaggeredFreeFlowBoundaryInfo, numEq> boundaryInfo_;
+};
+
+}
+
+#endif
diff --git a/dumux/freeflow/staggered/localresidual.hh b/dumux/freeflow/staggered/localresidual.hh
index 25d58b08d9..3c3a197194 100644
--- a/dumux/freeflow/staggered/localresidual.hh
+++ b/dumux/freeflow/staggered/localresidual.hh
@@ -274,17 +274,31 @@ protected:
 
                 this->ccResidual_ += boundaryFlux;
 
-                // set a fixed pressure for cells adjacent to a wall
-                if(bcTypes.isDirichlet(massBalanceIdx) && bcTypes.isDirichlet(momentumBalanceIdx))
-                {
-                    const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
-                    const auto& insideVolVars = elemVolVars[insideScv];
-                    this->ccResidual_[pressureIdx] = insideVolVars.pressure() - this->problem().dirichletAtPos(insideScv.center())[cellCenterIdx][pressureIdx];
-                }
+                asImp_().setFixedCell_(fvGeometry.scv(scvf.insideScvIdx()), elemVolVars, bcTypes);
             }
         }
     }
 
+    /*!
+     * \brief Sets a fixed Dirichlet value for a cell (such as pressure) at the boundary.
+     *        This is a provisional alternative to setting the Dirichlet value on the boundary directly.
+     *
+     * \param insideScv The sub control volume
+     * \param elemVolVars The current or previous element volVars
+     * \param bcTypes The boundary types
+     */
+    void setFixedCell_(const SubControlVolume& insideScv,
+                       const ElementVolumeVariables& elemVolVars,
+                       const BoundaryTypes& bcTypes)
+    {
+        // set a fixed pressure for cells adjacent to a wall
+        if(bcTypes.isDirichletCell(massBalanceIdx))
+        {
+            const auto& insideVolVars = elemVolVars[insideScv];
+            this->ccResidual_[pressureIdx] = insideVolVars.pressure() - this->problem().dirichletAtPos(insideScv.center())[cellCenterIdx][pressureIdx];
+        }
+    }
+
      /*!
      * \brief Evaluate boundary conditions for a face dof
      */
@@ -353,6 +367,14 @@ private:
         }
         return result;
     }
+
+    //! Returns the implementation of the problem (i.e. static polymorphism)
+    Implementation &asImp_()
+    { return *static_cast<Implementation *>(this); }
+
+    //! \copydoc asImp_()
+    const Implementation &asImp_() const
+    { return *static_cast<const Implementation *>(this); }
 };
 }
 
diff --git a/dumux/freeflow/staggered/propertydefaults.hh b/dumux/freeflow/staggered/propertydefaults.hh
index 8259358bf9..4206e1fa9c 100644
--- a/dumux/freeflow/staggered/propertydefaults.hh
+++ b/dumux/freeflow/staggered/propertydefaults.hh
@@ -38,6 +38,7 @@
 #include "fluxvariablescache.hh"
 #include "velocityoutput.hh"
 #include "vtkoutputmodule.hh"
+#include "boundarytypes.hh"
 
 #include <dumux/implicit/staggered/localresidual.hh>
 #include <dumux/freeflow/staggeredni/localresidual.hh>
@@ -180,6 +181,11 @@ public:
     using type = StaggeredPrimaryVariables<TypeTag, CellCenterBoundaryValues, FaceBoundaryValues>;
 };
 
+//! Boundary types at a single degree of freedom
+SET_TYPE_PROP(NavierStokes,
+              BoundaryTypes,
+              StaggeredFreeFlowBoundaryTypes<GET_PROP_VALUE(TypeTag, NumEq)>);
+
 SET_TYPE_PROP(NavierStokes, VtkOutputModule, FreeFlowStaggeredVtkOutputModule<TypeTag>);
 
 SET_TYPE_PROP(NavierStokes, VelocityOutput, StaggeredFreeFlowVelocityOutput<TypeTag>);
diff --git a/dumux/freeflow/staggerednc/localresidual.hh b/dumux/freeflow/staggerednc/localresidual.hh
index 5481018e17..0fb4538ede 100644
--- a/dumux/freeflow/staggerednc/localresidual.hh
+++ b/dumux/freeflow/staggerednc/localresidual.hh
@@ -60,13 +60,17 @@ class StaggeredNavierStokesResidualImpl;
 template<class TypeTag>
 class StaggeredNavierStokesResidualImpl<TypeTag, true> : public StaggeredNavierStokesResidualImpl<TypeTag, false>
 {
+    using ParentType = StaggeredNavierStokesResidualImpl<TypeTag, false>;
     friend class StaggeredLocalResidual<TypeTag>;
+    friend ParentType;
 
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
     using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
     using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+    using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
 
 
     using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
@@ -126,6 +130,38 @@ class StaggeredNavierStokesResidualImpl<TypeTag, true> : public StaggeredNavierS
 
         return storage;
     }
+
+protected:
+
+    /*!
+     * \brief Sets a fixed Dirichlet value for a cell (such as pressure) at the boundary.
+     *        This is a provisional alternative to setting the Dirichlet value on the boundary directly.
+     *
+     * \param insideScv The sub control volume
+     * \param elemVolVars The current or previous element volVars
+     * \param bcTypes The boundary types
+     */
+    void setFixedCell_(const SubControlVolume& insideScv,
+                       const ElementVolumeVariables& elemVolVars,
+                       const BoundaryTypes& bcTypes)
+    {
+        ParentType::setFixedCell_(insideScv, elemVolVars, bcTypes);
+
+        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+        {
+            // get equation index
+            const auto eqIdx = conti0EqIdx + compIdx;
+
+            // set a fixed mole fraction for cells
+            if(eqIdx != conti0EqIdx && bcTypes.isDirichletCell(eqIdx))
+            {
+                const auto& insideVolVars = elemVolVars[insideScv];
+                const Scalar massOrMoleFraction = useMoles ? insideVolVars.moleFraction(phaseIdx, compIdx) : insideVolVars.massFraction(phaseIdx, compIdx);
+                this->ccResidual_[eqIdx] = massOrMoleFraction - this->problem().dirichletAtPos(insideScv.center())[cellCenterIdx][eqIdx];
+            }
+        }
+
+    }
 };
 }
 
diff --git a/test/freeflow/staggered/closedsystemtestproblem.hh b/test/freeflow/staggered/closedsystemtestproblem.hh
index 36c4c1ab3a..0dd74e48bd 100644
--- a/test/freeflow/staggered/closedsystemtestproblem.hh
+++ b/test/freeflow/staggered/closedsystemtestproblem.hh
@@ -201,7 +201,7 @@ public:
 
         // set a fixed pressure in one cell
         if (isLowerLeftCell_(globalPos))
-            values.setDirichlet(massBalanceIdx);
+            values.setDirichletCell(massBalanceIdx);
         else
             values.setNeumann(massBalanceIdx);
 
diff --git a/test/freeflow/staggered/doneatestproblem.hh b/test/freeflow/staggered/doneatestproblem.hh
index 4541d8c559..49fe61b776 100644
--- a/test/freeflow/staggered/doneatestproblem.hh
+++ b/test/freeflow/staggered/doneatestproblem.hh
@@ -229,7 +229,7 @@ public:
 
         // set Dirichlet values for the velocity and pressure everywhere
         values.setDirichlet(momentumBalanceIdx);
-        values.setDirichlet(massBalanceIdx);
+        values.setDirichletCell(massBalanceIdx);
 
         return values;
     }
diff --git a/test/freeflow/staggerednc/densityflowproblem.hh b/test/freeflow/staggerednc/densityflowproblem.hh
index 685c2b0f94..8143c2891a 100644
--- a/test/freeflow/staggerednc/densityflowproblem.hh
+++ b/test/freeflow/staggerednc/densityflowproblem.hh
@@ -212,7 +212,7 @@ public:
         values.setOutflow(massBalanceIdx);
 
         if(globalPos[1] <  eps_)
-            values.setDirichlet(massBalanceIdx);
+            values.setDirichletCell(massBalanceIdx);
 
         if(globalPos[1] > this->bBoxMax()[1] - eps_)
         {
-- 
GitLab