Merge branch 'feature/symmetry-bc-for-staggered' into 'next'

[staggered][freeflow] Implement symmetry boundary condition

See merge request !550

dumux/discretization/staggered/globalvolumevariables.hh

@@ -109,7 +109,7 @@ public:
                 {
                     if(bcTypes.isDirichlet(eqIdx) || bcTypes.isDirichletCell(eqIdx))
                         boundaryPriVars[cellCenterIdx][eqIdx] = problem.dirichlet(element, scvf)[cellCenterIdx][eqIdx];
-                    else if(bcTypes.isNeumann(eqIdx) || bcTypes.isOutflow(eqIdx))
+                    else if(bcTypes.isNeumann(eqIdx) || bcTypes.isOutflow(eqIdx) || bcTypes.isSymmetry())
                         boundaryPriVars[cellCenterIdx][eqIdx] = sol[cellCenterIdx][scvf.insideScvIdx()][eqIdx];
                     //TODO: this assumes a zero-gradient for e.g. the pressure on the boundary
                     // could be made more general by allowing a non-zero-gradient, provided in problem file

dumux/freeflow/staggered/boundarytypes.hh

@@ -55,6 +55,7 @@ public:
 
         boundaryInfo_[eqIdx].visited = false;
         boundaryInfo_[eqIdx].isDirichletCell = false;
+        boundaryInfo_[eqIdx].isSymmetry = false;
     }
 
 
@@ -83,11 +84,32 @@ public:
     bool isDirichletCell(unsigned eqIdx) const
     { return boundaryInfo_[eqIdx].isDirichletCell; }
 
+    /*!
+     * \brief Sets a symmetry boundary condition for all equations
+     */
+    void setSymmetry()
+    {
+        for (int eqIdx=0; eqIdx < numEq; ++eqIdx)
+        {
+            resetEq(eqIdx);
+            boundaryInfo_[eqIdx].visited = true;
+            boundaryInfo_[eqIdx].isSymmetry = true;
+            Valgrind::SetDefined(boundaryInfo_[eqIdx]);
+        }
+    }
+
+    /*!
+     * \brief Returns true if the there is a symmetry boundary condition
+     */
+    bool isSymmetry() const
+    { return boundaryInfo_[0].isSymmetry; }
+
 
 protected:
     struct StaggeredFreeFlowBoundaryInfo {
         bool visited;
         bool isDirichletCell;
+        bool isSymmetry;
     };
 
     std::array<StaggeredFreeFlowBoundaryInfo, numEq> boundaryInfo_;

dumux/freeflow/staggered/fluxvariables.hh

@@ -264,11 +264,27 @@ public:
       };
 
       // account for all sub-faces
-      for(auto subFaceData : scvf.pairData())
+      for(const auto& subFaceData : scvf.pairData())
       {
-          const int eIdx = scvf.insideScvIdx();
+          const auto eIdx = scvf.insideScvIdx();
           const auto& normalFace = fvGeometry.scvf(eIdx, subFaceData.localNormalFaceIdx);
 
+          // Check if we have a symmetry boundary condition. If yes, the tangental part of the momentum flux can be neglected.
+          if(subFaceData.outerParallelFaceDofIdx < 0)
+          {
+              // lambda to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
+              auto makeGhostFace = [eIdx] (const GlobalPosition& pos)
+              {
+                  return SubControlVolumeFace(pos, std::vector<unsigned int>{eIdx,eIdx});
+              };
+
+              // 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, makeGhostFace(subFaceData.virtualOuterParallelFaceDofPos));
+              if(bcTypes.isSymmetry())
+                continue;
+          }
+
+          // if there is no symmetry boundary condition, proceed to calculate the tangential momentum flux
           if(navierStokes)
               tangentialFlux += computeAdvectivePartOfTangentialMomentumFlux_(problem, element, scvf, normalFace, subFaceData, elemVolVars, velocity);
 
@@ -292,8 +308,8 @@ private:
       const auto insideScvIdx = normalFace.insideScvIdx();
       const auto outsideScvIdx = normalFace.outsideScvIdx();
 
-      // convenience function to create a ghost face, outside of the domain
-    //   auto makeGhostFace = [insideScvIdx, outsideScvIdx] (const GlobalPosition& pos)
+      // lambda to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
+    //   auto makeGhostFace = [insideScvIdx] (const GlobalPosition& pos)
     //   {
     //       return SubControlVolumeFace(pos, std::vector<unsigned int>{insideScvIdx,insideScvIdx});
     //   };
@@ -340,8 +356,8 @@ private:
       const auto& insideVolVars = elemVolVars[insideScvIdx];
       const auto& outsideVolVars = elemVolVars[outsideScvIdx];
 
-      // convenience function to create a ghost face, outside of the domain
-      auto makeGhostFace = [insideScvIdx, outsideScvIdx] (const GlobalPosition& pos)
+      // lambda to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
+      auto makeGhostFace = [insideScvIdx] (const GlobalPosition& pos)
       {
           return SubControlVolumeFace(pos, std::vector<unsigned int>{insideScvIdx,insideScvIdx});
       };

dumux/freeflow/staggered/localresidual.hh

@@ -315,7 +315,7 @@ protected:
             // handle the actual boundary conditions:
             const auto bcTypes = this->problem().boundaryTypes(element, scvf);
 
-            // set a fixed value for the velocity
+            // set a fixed value for the velocity for Dirichlet boundary conditions
             if(bcTypes.isDirichlet(momentumBalanceIdx))
             {
                 const Scalar velocity = faceVars.faceVars(scvf.dofIndex()).velocity();
@@ -323,6 +323,15 @@ protected:
                 this->faceResiduals_[scvf.localFaceIdx()] = velocity - dirichletValue;
             }
 
+            // 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 fixedValue = 0.0;
+                this->faceResiduals_[scvf.localFaceIdx()] = velocity - fixedValue;
+            }
+
             // outflow condition for the momentum balance equation
             if(bcTypes.isOutflow(momentumBalanceIdx))
             {