From 4cb8ce63aec3ca146209b32ad8000ef35fbba80a Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Fri, 23 Oct 2020 11:15:27 +0200
Subject: [PATCH] [rans][problem] cleanup and separate stored kinematic
 viscosity into viscosity and density

---
 dumux/freeflow/rans/problem.hh                | 21 +++++++++++++------
 dumux/freeflow/rans/twoeq/kepsilon/problem.hh | 11 ++--------
 2 files changed, 17 insertions(+), 15 deletions(-)

diff --git a/dumux/freeflow/rans/problem.hh b/dumux/freeflow/rans/problem.hh
index d4bc6ee852..75675792a3 100644
--- a/dumux/freeflow/rans/problem.hh
+++ b/dumux/freeflow/rans/problem.hh
@@ -126,7 +126,8 @@ public:
         vorticityTensorScalarProduct_.resize(this->gridGeometry().elementMapper().size(), 0.0);
         flowDirectionAxis_.resize(this->gridGeometry().elementMapper().size(), fixedFlowDirectionAxis_);
         wallNormalAxis_.resize(this->gridGeometry().elementMapper().size(), fixedWallNormalAxis_);
-        kinematicViscosity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
+        storedViscosity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
+        storedDensity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
 
         if ( !(hasParamInGroup(this->paramGroup(), "RANS.IsFlatWallBounded")))
         {
@@ -254,7 +255,7 @@ public:
         calculateMaxMinVelocities_();
         calculateStressTensor_();
         calculateVorticityTensor_();
-        storeKinematicViscosity_(curSol);
+        storeViscosities_(curSol);
     }
 
     /*!
@@ -412,8 +413,14 @@ public:
     Scalar vorticityTensorScalarProduct(const int elementIdx) const
     { return vorticityTensorScalarProduct_[elementIdx]; }
 
+    Scalar storedViscosity(const int elementIdx) const
+    { return storedViscosity_[elementIdx]; }
+
+    Scalar storedDensity(const int elementIdx) const
+    { return storedDensity_[elementIdx]; }
+
     Scalar kinematicViscosity(const int elementIdx) const
-    { return kinematicViscosity_[elementIdx]; }
+    { return storedViscosity(elementIdx) / storedDensity(elementIdx); }
 
     bool calledUpdateStaticWallProperties = false;
 
@@ -675,7 +682,7 @@ private:
         }
     }
 
-    void storeKinematicViscosity_(const SolutionVector& curSol)
+    void storeViscosities_(const SolutionVector& curSol)
     {
         // calculate or call all secondary variables
         for (const auto& element : elements(this->gridGeometry().gridView()))
@@ -693,7 +700,8 @@ private:
 
                 VolumeVariables volVars;
                 volVars.update(elemSol, asImp_(), element, scv);
-                kinematicViscosity_[elementIdx] = volVars.viscosity() / volVars.density();
+                storedDensity_[elementIdx] = volVars.density();
+                storedViscosity_[elementIdx] = volVars.viscosity();
             }
         }
     }
@@ -716,7 +724,8 @@ private:
     std::vector<Scalar> stressTensorScalarProduct_;
     std::vector<Scalar> vorticityTensorScalarProduct_;
 
-    std::vector<Scalar> kinematicViscosity_;
+    std::vector<Scalar> storedDensity_;
+    std::vector<Scalar> storedViscosity_;
 
     //! Returns the implementation of the problem (i.e. static polymorphism)
     Implementation &asImp_()
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
index 5ac8179b8d..3f73d9b800 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
@@ -102,7 +102,6 @@ public:
         ParentType::updateStaticWallProperties();
         // update size and initial values of the global vectors
         matchingPointIdx_.resize(this->gridGeometry().elementMapper().size(), 0);
-        storedDensity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
         storedDissipation_.resize(this->gridGeometry().elementMapper().size(), 0.0);
         storedTurbulentKineticEnergy_.resize(this->gridGeometry().elementMapper().size(), 0.0);
         storedDynamicEddyViscosity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
@@ -138,7 +137,6 @@ public:
                 VolumeVariables volVars;
                 volVars.update(elemSol, asImp_(), element, scv);
                 storedDynamicEddyViscosity_[elementIdx] = volVars.calculateEddyViscosity();
-                storedDensity_[elementIdx] = volVars.density();
             }
         }
 
@@ -255,7 +253,7 @@ public:
         unsigned int wallNormalAxis = asImp_().wallNormalAxis(elementIdx);
         unsigned int flowDirectionAxis = asImp_().flowDirectionAxis(elementIdx);
         Scalar velocityGradient = asImp_().velocityGradient(elementIdx, flowDirectionAxis, wallNormalAxis);
-        return mixingLength * mixingLength * abs(velocityGradient) * storedDensity(elementIdx);
+        return mixingLength * mixingLength * abs(velocityGradient) * asImp_().storedDensity(elementIdx);
     }
 
     //! \brief Returns the wall shear stress velocity
@@ -302,8 +300,7 @@ public:
     {
         using std::log;
         Scalar velocityNominal = uStarNominal(elementIdx) * (1.0 / asImp_().karmanConstant() * log(yPlusNominal(elementIdx)) + 5.0);
-        return uStarNominal(elementIdx) * uStarNominal(elementIdx)
-               * velocity / velocityNominal;
+        return uStarNominal(elementIdx) * uStarNominal(elementIdx) * velocity / velocityNominal;
     }
 
     //! \brief Checks whether a wall function should be used
@@ -469,9 +466,6 @@ public:
         return useStoredEddyViscosity;
     }
 
-    Scalar storedDensity(const int elementIdx) const
-    { return storedDensity_[elementIdx]; }
-
     Scalar storedDissipation(const int elementIdx) const
     { return storedDissipation_[elementIdx]; }
 
@@ -489,7 +483,6 @@ public:
 
 private:
     std::vector<unsigned int> matchingPointIdx_;
-    std::vector<Scalar> storedDensity_;
     std::vector<Scalar> storedDissipation_;
     std::vector<Scalar> storedTurbulentKineticEnergy_;
     std::vector<Scalar> storedDynamicEddyViscosity_;
-- 
GitLab