From 20525641d0fee44eeca1b6e519f2f83616987311 Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Mon, 19 Oct 2020 17:06:35 +0200
Subject: [PATCH] [rans][problem] add geometric check for isFlatWallBounded

---
 dumux/common/parameters.hh     |  2 --
 dumux/freeflow/rans/problem.hh | 42 ++++++++++++++++++++++++++++++++--
 2 files changed, 40 insertions(+), 4 deletions(-)

diff --git a/dumux/common/parameters.hh b/dumux/common/parameters.hh
index 801b2b126c..c54fe99fc8 100644
--- a/dumux/common/parameters.hh
+++ b/dumux/common/parameters.hh
@@ -333,8 +333,6 @@ private:
         // parameters in the mpfa group
         defaultParams["MPFA.Q"] = "0.0";
 
-        // parameters in the RANS group
-        defaultParams["RANS.IsFlatWallBounded"] = "false";
 
         // merge the global default tree but do not overwrite if the parameter already exists
         mergeTree_(params, defaultParams, false);
diff --git a/dumux/freeflow/rans/problem.hh b/dumux/freeflow/rans/problem.hh
index 997528b0c0..d4bc6ee852 100644
--- a/dumux/freeflow/rans/problem.hh
+++ b/dumux/freeflow/rans/problem.hh
@@ -24,6 +24,8 @@
 #ifndef DUMUX_RANS_PROBLEM_HH
 #define DUMUX_RANS_PROBLEM_HH
 
+#include <algorithm>
+
 #include <dune/common/fmatrix.hh>
 #include <dumux/common/properties.hh>
 #include <dumux/common/staggeredfvproblem.hh>
@@ -126,6 +128,13 @@ public:
         wallNormalAxis_.resize(this->gridGeometry().elementMapper().size(), fixedWallNormalAxis_);
         kinematicViscosity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
 
+        if ( !(hasParamInGroup(this->paramGroup(), "RANS.IsFlatWallBounded")))
+        {
+            std::cout << "The parameter \"Rans.IsFlatWallBounded\" is not specified. \n"
+                    << " -- Based on the grid and the isOnWallAtPos function specified by the user,"
+                    << " this parameter is set to be "<< std::boolalpha << isFlatWallBounded() << "\n";
+        }
+
         std::vector<WallElementInformation> wallElements;
 
         const auto gridView = this->gridGeometry().gridView();
@@ -307,8 +316,8 @@ public:
 
     bool isFlatWallBounded() const
     {
-        static const bool isFlatWallBounded = getParamFromGroup<bool>(this->paramGroup(), "RANS.IsFlatWallBounded");
-        return isFlatWallBounded;
+        static const bool hasAlignedWalls = hasAlignedWalls_();
+        return hasAlignedWalls;
     }
 
     /*!
@@ -409,6 +418,35 @@ public:
     bool calledUpdateStaticWallProperties = false;
 
 private:
+
+    bool hasAlignedWalls_() const
+    {
+        if ( hasParamInGroup(this->paramGroup(), "RANS.IsFlatWallBounded"))
+        {
+            static const bool isFlatWallBounded = getParamFromGroup<bool>(this->paramGroup(), "RANS.IsFlatWallBounded");
+            return isFlatWallBounded;
+        }
+
+        std::vector<int> wallFaceAxis;
+        wallFaceAxis.reserve(this->gridGeometry().numBoundaryScvf());
+
+        const auto gridView = this->gridGeometry().gridView();
+        auto fvGeometry = localView(this->gridGeometry());
+        for (const auto& element : elements(gridView))
+        {
+            fvGeometry.bindElement(element);
+            for (const auto& scvf : scvfs(fvGeometry))
+            {
+                // only search for walls at a global boundary
+                if (!scvf.boundary() && asImp_().isOnWall(scvf))
+                    wallFaceAxis.push_back(scvf.directionIndex());
+            }
+        }
+
+        // Returns if all wall directions are the same
+        return std::all_of(wallFaceAxis.begin(), wallFaceAxis.end(), [firstDir=wallFaceAxis[0]](auto dir){ return (dir == firstDir);} ) ;
+    }
+
     void calculateCCVelocities_(const SolutionVector& curSol)
     {
         // calculate cell-center-averaged velocities
-- 
GitLab