Merge branch 'fix/ransVelGradForSingleElement' into 'master'

[freeflow][rans] Set velocity gradient to zero, if it is a single element in on direction

See merge request !1130

dumux/freeflow/rans/problem.hh

@@ -110,25 +110,24 @@ public:
         std::vector<unsigned int> wallElements;
         std::vector<GlobalPosition> wallPositions;
         std::vector<unsigned int> wallNormalAxisTemp;
-        auto& gridView(this->fvGridGeometry().gridView());
+
+        const auto gridView = this->fvGridGeometry().gridView();
+        auto fvGeometry = localView(this->fvGridGeometry());
+
         for (const auto& element : elements(gridView))
         {
-            for (const auto& intersection : intersections(gridView, element))
+            fvGeometry.bindElement(element);
+            for (const auto& scvf : scvfs(fvGeometry))
             {
                 // only search for walls at a global boundary
-                if (!intersection.boundary())
+                if (!scvf.boundary())
                     continue;
 
-                GlobalPosition global = intersection.geometry().center();
-                if (asImp_().isOnWall(global))
+                if (asImp_().isOnWall(scvf))
                 {
                     wallElements.push_back(this->fvGridGeometry().elementMapper().index(element));
-                    wallPositions.push_back(global);
-                    for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
-                    {
-                        if (abs(intersection.centerUnitOuterNormal()[dimIdx]) > 1e-8)
-                            wallNormalAxisTemp.push_back(dimIdx);
-                    }
+                    wallPositions.push_back(scvf.center());
+                    wallNormalAxisTemp.push_back(scvf.directionIndex());
                 }
             }
         }
@@ -265,6 +264,9 @@ public:
                               - velocity_[neighborIdx_[elementIdx][dimIdx][0]][velIdx])
                           / (cellCenter_[neighborIdx_[elementIdx][dimIdx][1]][dimIdx]
                               - cellCenter_[neighborIdx_[elementIdx][dimIdx][0]][dimIdx]);
+                    if (abs(cellCenter_[neighborIdx_[elementIdx][dimIdx][1]][dimIdx]
+                            - cellCenter_[neighborIdx_[elementIdx][dimIdx][0]][dimIdx]) < 1e-8)
+                        velocityGradients_[elementIdx][velIdx][dimIdx] = 0.0;
                 }
 
                 if (abs(velocity_[elementIdx][dimIdx]) > abs(velocityMaximum_[wallElementIdx][dimIdx]))
@@ -414,12 +416,22 @@ public:
         }
     }
 
+    /*!
+     * \brief Returns whether a given sub control volume face is on a wall
+     *
+     * \param scvf The sub control volume face.
+     */
+    bool isOnWall(const SubControlVolumeFace& scvf) const
+    {
+        return asImp_().isOnWallAtPos(scvf.center());
+    }
+
     /*!
      * \brief Returns whether a given point is on a wall
      *
      * \param globalPos The position in global coordinates.
      */
-    bool isOnWall(const GlobalPosition &globalPos) const
+    bool isOnWallAtPos(const GlobalPosition &globalPos) const
     {
         // Throw an exception if no walls are implemented
         DUNE_THROW(Dune::InvalidStateException,

dumux/freeflow/rans/twoeq/kepsilon/problem.hh

@@ -312,7 +312,7 @@ public:
     {
         unsigned int elementIdx = asImp_().fvGridGeometry().elementMapper().index(element);
         auto bcTypes = asImp_().boundaryTypes(element, localSubFace);
-        return asImp_().isOnWall(localSubFace.center())
+        return asImp_().isOnWall(localSubFace)
                && bcTypes.isDirichlet(eqIdx)
                && isMatchingPoint(elementIdx);
     }

dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh

@@ -154,7 +154,7 @@ public:
         const auto bcTypes = problem.boundaryTypes(element, scvf);
         if (!(scvf.boundary() && (bcTypes.isOutflow(Indices::turbulentKineticEnergyEqIdx)
                                   || bcTypes.isSymmetry()
-                                  || problem.isOnWall(scvf.center()))))
+                                  || problem.isOnWall(scvf))))
         {
             if (!(insideVolVars.isMatchingPoint() && outsideVolVars.isMatchingPoint())
                 || !(insideVolVars.isMatchingPoint() && outsideVolVars.inNearWallRegion())

test/freeflow/rans/pipelauferproblem.hh

@@ -185,7 +185,7 @@ public:
      */
     // \{
 
-    bool isOnWall(const GlobalPosition &globalPos) const
+    bool isOnWallAtPos(const GlobalPosition &globalPos) const
     {
         return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_
                || globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_;
@@ -260,7 +260,7 @@ public:
 #endif
 #if KOMEGA
             // set a fixed dissipation (omega) in one cell
-            if (isOnWall(globalPos))
+            if (isOnWallAtPos(globalPos))
                 values.setDirichletCell(Indices::dissipationIdx);
 #endif
 #if ONEEQ
@@ -285,7 +285,7 @@ public:
         if (time() > 10.0)
         {
             values[Indices::temperatureIdx] = inletTemperature_;
-            if (isOnWall(globalPos))
+            if (isOnWallAtPos(globalPos))
             {
                 values[Indices::temperatureIdx] = wallTemperature_;
             }
@@ -327,14 +327,14 @@ public:
         values[Indices::velocityXIdx] = time() > initializationTime_
                                         ? inletVelocity_
                                         : time() / initializationTime_ * inletVelocity_;
-        if (isOnWall(globalPos))
+        if (isOnWallAtPos(globalPos))
         {
             values[Indices::velocityXIdx] = 0.0;
         }
 
 #if NONISOTHERMAL
         values[Indices::temperatureIdx] = inletTemperature_;
-        if (isOnWall(globalPos))
+        if (isOnWallAtPos(globalPos))
         {
             values[Indices::temperatureIdx] = wallTemperature_;
         }
@@ -343,7 +343,7 @@ public:
 #if LOWREKEPSILON || KEPSILON || KOMEGA
         values[Indices::turbulentKineticEnergyIdx] = turbulentKineticEnergy_;
         values[Indices::dissipationIdx] = dissipation_;
-        if (isOnWall(globalPos))
+        if (isOnWallAtPos(globalPos))
         {
             values[Indices::turbulentKineticEnergyIdx] = 0.0;
             values[Indices::dissipationIdx] = 0.0;
@@ -352,7 +352,7 @@ public:
 
 #if ONEEQ
         values[Indices::viscosityTildeIdx] = viscosityTilde_;
-        if (isOnWall(globalPos))
+        if (isOnWallAtPos(globalPos))
         {
             values[Indices::viscosityTildeIdx] = 0.0;
         }

test/freeflow/ransnc/flatplatetestproblem.hh

@@ -266,7 +266,7 @@ public:
             values.setOutflow(Indices::viscosityTildeIdx);
 #endif
         }
-        else if(isOnWall(globalPos))
+        else if(isOnWallAtPos(globalPos))
         {
             values.setDirichlet(Indices::velocityXIdx);
             values.setDirichlet(Indices::velocityYIdx);
@@ -316,7 +316,7 @@ public:
                 values[transportCompIdx] = 1e-3;
             }
 #if NONISOTHERMAL
-            if (isOnWall(globalPos))
+            if (isOnWallAtPos(globalPos))
             {
                 values[Indices::temperatureIdx] += 30.;
             }
@@ -370,14 +370,14 @@ public:
 
         // block velocity profile
         values[Indices::velocityXIdx] = 0.0;
-        if (!isOnWall(globalPos))
+        if (!isOnWallAtPos(globalPos))
             values[Indices::velocityXIdx] =  inletVelocity_;
         values[Indices::velocityYIdx] = 0.0;
 
 #if KEPSILON || KOMEGA || LOWREKEPSILON
         values[Indices::turbulentKineticEnergyEqIdx] = turbulentKineticEnergy_;
         values[Indices::dissipationEqIdx] = dissipation_;
-        if (isOnWall(globalPos))
+        if (isOnWallAtPos(globalPos))
         {
             values[Indices::turbulentKineticEnergyEqIdx] = 0.0;
             values[Indices::dissipationEqIdx] = 0.0;
@@ -386,7 +386,7 @@ public:
 
 #if ONEEQ
         values[Indices::viscosityTildeIdx] = viscosityTilde_;
-        if (isOnWall(globalPos))
+        if (isOnWallAtPos(globalPos))
         {
             values[Indices::viscosityTildeIdx] = 0.0;
         }
@@ -407,7 +407,7 @@ public:
         return timeLoop_->time();
     }
 
-    bool isOnWall(const GlobalPosition& globalPos) const
+    bool isOnWallAtPos(const GlobalPosition& globalPos) const
     {
         return globalPos[1] < eps_;
     }