[lowrekepsilon] Various improvements

* add test
* add PDELab solution for comparison
* speed-up by avoid calling getparam
* improve determination of the gradient
* increase verbosity for comparison

dumux/freeflow/rans/problem.hh

@@ -123,7 +123,7 @@ public:
                     wallPositions.push_back(global);
                     for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
                     {
-                        if (abs(intersection.centerUnitOuterNormal()[dimIdx]) > 1e-10)
+                        if (abs(intersection.centerUnitOuterNormal()[dimIdx]) > 1e-8)
                             wallNormalAxisTemp.push_back(dimIdx);
                     }
                 }

dumux/freeflow/rans/twoeq/lowrekepsilon/indices.hh

@@ -24,7 +24,7 @@
 #ifndef DUMUX_LOWREKEPSILON_INDICES_HH
 #define DUMUX_LOWREKEPSILON_INDICES_HH
 
-#include <dumux/freeflow//navierstokes/indices.hh>
+#include <dumux/freeflow/navierstokes/indices.hh>
 
 namespace Dumux {
 

dumux/freeflow/rans/twoeq/lowrekepsilon/model.hh

@@ -113,6 +113,19 @@ public:
     using type = LowReKEpsilonVtkOutputFields<FVGridGeometry>;
 };
 
+//! Set one or different base epsilons for the calculations of the localJacobian's derivatives
+SET_PROP(LowReKEpsilon, BaseEpsilon)
+{
+private:
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+public:
+    static constexpr auto getEps()
+    {
+        return std::array<std::array<Scalar, 2>, 2>{{{1e-6/*dCCdCC*/, 1e-8/*dCCdFace*/},
+                                                     {1e-6/*dFacedCC*/, 1e-8/*dFacedFace*/}}};
+    }
+};
+
 //////////////////////////////////////////////////////////////////
 // default property values for the non-isothermal low-Reynolds k-epsilon model
 //////////////////////////////////////////////////////////////////

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

@@ -73,6 +73,7 @@ public:
         {
             DUNE_THROW(Dune::NotImplemented, "This low-Re k-epsilon model is not implemented: " << lowReKEpsilonModel_);
         }
+        useStoredEddyViscosity_ = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "RANS.UseStoredEddyViscosity", true);
     }
 
     /*!
@@ -131,6 +132,7 @@ public:
     std::vector<Scalar> storedKinematicEddyViscosity_;
     std::vector<Scalar> storedTurbulentKineticEnergy_;
     int lowReKEpsilonModel_;
+    bool useStoredEddyViscosity_;
 
 private:
     //! Returns the implementation of the problem (i.e. static polymorphism)

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

@@ -95,8 +95,8 @@ public:
 
         // calculate advective flux
         const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
-        const bool isOutflowK = bcTypes.isOutflow(turbulentKineticEnergyEqIdx);
-        const bool isOutflowEpsilon = bcTypes.isOutflow(dissipationEqIdx);
+        const bool isOutflowK = scvf.boundary() && bcTypes.isOutflow(turbulentKineticEnergyEqIdx);
+        const bool isOutflowEpsilon = scvf.boundary() && bcTypes.isOutflow(dissipationEqIdx);
         auto upwindTermK = [](const auto& volVars)
         {
             return volVars.turbulentKineticEnergy();

dumux/freeflow/rans/twoeq/lowrekepsilon/volumevariables.hh

@@ -110,9 +110,8 @@ public:
         storedDissipationTilde_ = problem.storedDissipationTilde_[ParentType::elementID()];
         storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[ParentType::elementID()];
         stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[ParentType::elementID()];
-        dofPosition_ = scv.dofPosition();
-        if (getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "RANS.UseStoredEddyViscosity", true))
-            ParentType::setDynamicEddyViscosity(problem.storedKinematicEddyViscosity_[ParentType::elementID()] * ParentType::density());
+        if (problem.useStoredEddyViscosity_)
+            ParentType::setKinematicEddyViscosity(problem.storedKinematicEddyViscosity_[ParentType::elementID()]);
         else
             calculateEddyViscosity();
     }
@@ -124,8 +123,18 @@ public:
     {
         Scalar kinematicEddyViscosity
             = cMu() * fMu() * turbulentKineticEnergy() * turbulentKineticEnergy()
-              / std::max(dissipationTilde(), getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "RANS.MinDissipation", -1));
-        ParentType::setDynamicEddyViscosity(kinematicEddyViscosity * ParentType::density());
+              / dissipationTilde();
+//         if (std::isnan(kinematicEddyViscosity) || dissipationTilde() < 1e-8)
+//         {
+        Dune::dinfo /*<< " scv.dofPosition() " << dofPosition_
+                    */<< " cMu() " << cMu()
+                    << " fMu() " << fMu()
+                    << " turbulentKineticEnergy() " << turbulentKineticEnergy()
+                    << " dissipationTilde() " << dissipationTilde()
+                    << " kinematicEddyViscosity " << kinematicEddyViscosity
+                    << std::endl;
+//         }
+        ParentType::setKinematicEddyViscosity(kinematicEddyViscosity);
     }
 
     /*!
@@ -303,7 +312,6 @@ public:
     }
 
 protected:
-    Dune::FieldVector<Scalar, 2> dofPosition_;
     int lowReKEpsilonModel_;
     Scalar turbulentKineticEnergy_;
     Scalar dissipationTilde_;

dumux/freeflow/rans/twoeq/lowrekepsilon/vtkoutputfields.hh

@@ -54,6 +54,19 @@ public:
         vtk.addVolumeVariable([](const auto& v){ return v.turbulentKineticEnergy(); }, "turbulentKineticEnergy");
         vtk.addVolumeVariable([](const auto& v){ return v.dissipationTilde(); }, "dissipation");
         vtk.addVolumeVariable([](const auto& v){ return v.stressTensorScalarProduct(); }, "stressTensorScalarProduct");
+        vtk.addVolumeVariable([](const auto& v){ return v.kinematicEddyViscosity(); }, "eddyViscosity");
+        vtk.addVolumeVariable([](const auto& v){ return 2.0 * v.kinematicEddyViscosity()
+                                               * v.stressTensorScalarProduct(); }, "production_k");
+        vtk.addVolumeVariable([](const auto& volVars){ return volVars.cOneEpsilon() * volVars.fOne()
+                                    * volVars.dissipationTilde() / volVars.turbulentKineticEnergy()
+                                    * 2.0 * volVars.kinematicEddyViscosity()
+                                    * volVars.stressTensorScalarProduct(); }, "production_eps");
+        vtk.addVolumeVariable([](const auto& volVars){ return volVars.dissipationTilde(); }, "destruction_k");
+        vtk.addVolumeVariable([](const auto& volVars){ return volVars.cTwoEpsilon() * volVars.fTwo()
+                                    * volVars.dissipationTilde() * volVars.dissipationTilde()
+                                    / volVars.turbulentKineticEnergy(); }, "destruction_eps");
+        vtk.addVolumeVariable([](const auto& v){ return v.dValue(); }, "dValue");
+        vtk.addVolumeVariable([](const auto& v){ return v.eValue(); }, "eValue");
     }
 };
 

dumux/freeflow/rans/volumevariables.hh

@@ -149,11 +149,11 @@ public:
     { return uPlus_; }
 
     /*!
-     * \brief Return the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$ of the flow within the
+     * \brief Return the kinematic eddy viscosity \f$\mathrm{[m^2/s]}\f$ of the flow within the
      *        control volume.
      */
-    Scalar dynamicEddyViscosity() const
-    { return dynamicEddyViscosity_; }
+    Scalar kinematicEddyViscosity() const
+    { return dynamicEddyViscosity() / asImp_().density(); }
 
     /*!
      * \brief Return the kinematic viscosity \f$\mathrm{[m^2/s]}\f$ of the fluid within the
@@ -163,11 +163,11 @@ public:
     { return asImp_().viscosity() / asImp_().density(); }
 
     /*!
-     * \brief Return the kinematic eddy viscosity \f$\mathrm{[Pa s]}\f$ of the flow within the
+     * \brief Return the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$ of the flow within the
      *        control volume.
      */
-    Scalar kinematicEddyViscosity() const
-    { return dynamicEddyViscosity() / asImp_().density(); }
+    Scalar dynamicEddyViscosity() const
+    { return dynamicEddyViscosity_; }
 
 protected:
     DimVector velocity_;

test/freeflow/rans/CMakeLists.txt

 add_input_file_links()
-dune_symlink_to_source_files(FILES laufer_re50000.csv laufer_re50000_u+y+.csv)
+dune_symlink_to_source_files(FILES laufer_re50000.csv laufer_re50000_u+y+.csv pdelab-lowrekepsilon.csv pdelab-zeroeq.csv)
 
 dune_add_test(NAME test_pipe_laufer_zeroeq
               SOURCES test_pipe_laufer.cc
@@ -7,7 +7,7 @@ dune_add_test(NAME test_pipe_laufer_zeroeq
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS  --script fuzzy
                         --files ${CMAKE_SOURCE_DIR}/test/references/pipe_laufer_zeroeq.vtu
-                                ${CMAKE_CURRENT_BINARY_DIR}/pipe_laufer_zeroeq_reference-00016.vtu
+                                ${CMAKE_CURRENT_BINARY_DIR}/pipe_laufer_reference-00044.vtu
                         --command "${CMAKE_CURRENT_BINARY_DIR}/test_pipe_laufer_zeroeq test_pipe_laufer_reference.input")
 
 dune_add_test(NAME test_pipe_laufer_lowrekepsilon
@@ -15,9 +15,9 @@ dune_add_test(NAME test_pipe_laufer_lowrekepsilon
               CMAKE_GUARD HAVE_UMFPACK
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS  --script fuzzy
-                        --files ${CMAKE_SOURCE_DIR}/test/references/pipe_laufer_zeroeq.vtu
-                                ${CMAKE_CURRENT_BINARY_DIR}/pipe_laufer_zeroeq_reference-00020.vtu
-                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_pipe_laufer_zeroeq test_pipe_laufer_reference.input")
+                        --files ${CMAKE_SOURCE_DIR}/test/references/pipe_laufer_lowrekepsilon.vtu
+                                ${CMAKE_CURRENT_BINARY_DIR}/pipe_laufer_reference-00074.vtu
+                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_pipe_laufer_lowrekepsilon test_pipe_laufer_reference.input")
 target_compile_definitions(test_pipe_laufer_lowrekepsilon PUBLIC "LOWREKEPSILON=1")
 
 dune_add_test(NAME test_pipe_zeroeqni

test/freeflow/rans/pipelauferproblem.hh

@@ -192,6 +192,12 @@ public:
         values.setOutflow(Indices::conti0EqIdx);
         values.setDirichlet(Indices::momentumXBalanceIdx);
         values.setDirichlet(Indices::momentumYBalanceIdx);
+        if (isOutlet(globalPos))
+        {
+            values.setDirichlet(Indices::conti0EqIdx);
+            values.setOutflow(Indices::momentumXBalanceIdx);
+            values.setOutflow(Indices::momentumYBalanceIdx);
+        }
 
 #if NONISOTHERMAL
         values.setDirichlet(Indices::energyBalanceIdx);
@@ -211,14 +217,6 @@ public:
         }
 #endif
 
-        // set a fixed pressure in one cell
-        if (isOutlet(globalPos))
-        {
-            values.setDirichlet(Indices::conti0EqIdx);
-            values.setOutflow(Indices::momentumXBalanceIdx);
-            values.setOutflow(Indices::momentumYBalanceIdx);
-        }
-
         return values;
     }
 

test/freeflow/rans/test_pipe_laufer.cc

@@ -24,16 +24,16 @@
  * This test simulates is based on pipe flow experiments by
  * John Laufers experiments in 1954 \cite Laufer1954a.
  */
- #include <config.h>
+#include <config.h>
 
- #include <ctime>
- #include <iostream>
+#include <ctime>
+#include <iostream>
 
- #include <dune/common/parallel/mpihelper.hh>
- #include <dune/common/timer.hh>
- #include <dune/grid/io/file/dgfparser/dgfexception.hh>
- #include <dune/grid/io/file/vtk.hh>
- #include <dune/istl/io.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+#include <dune/istl/io.hh>
 
 #include "pipelauferproblem.hh"
 
@@ -204,8 +204,8 @@ int main(int argc, char** argv) try
     ////////////////////////////////////////////////////////////
 
 #if HAVE_PVPYTHON
-    bool plotLawOfTheWall = getParam<bool>("Output.PlotLawOfTheWall", false);
-    bool plotVelocityProfile = getParam<bool>("Output.PlotVelocityProfile", false);
+    static const bool plotLawOfTheWall = getParam<bool>("Output.PlotLawOfTheWall", false);
+    static const bool plotVelocityProfile = getParam<bool>("Output.PlotVelocityProfile", false);
     if (plotLawOfTheWall || plotVelocityProfile)
     {
         char fileName[255];

test/freeflow/rans/test_pipe_laufer.input

 [TimeLoop]
-DtInitial = 1e-3 # [s]
+DtInitial = 1e-5 # [s]
 TEnd = 100 # [s]
-# MaxTimeStepSize = 1e-2
 
 [Grid]
 Verbosity = true
 Positions0 = 0.0 10.0
 Positions1 = 0.0 0.12345 0.2469
-Cells0 = 5
-# Cells1 = 25 25
-# Grading1 = 1.2 -1.2
+Cells0 = 25
+Cells1 = 25 25
+Grading1 = 1.2 -1.2
+# Cells0 = 5
 # Cells1 = 15 15
 # Grading1 = 1.5 -1.5
-Cells1 = 5 5
-Grading1 = 4.0 -4.0
+# Cells0 = 2
+# Cells1 = 1 1
+# Grading1 = 4.0 -4.0
 
 [Output]
 PlotLawOfTheWall = true
 PlotVelocityProfile = true
 
 [Problem]
-Name = pipe_laufer_zeroeq
+Name = pipe_laufer
 InletVelocity = 2.5 # [m/s]
 EnableGravity = false
 SandGrainRoughness = 0.0 # [m] # not implemented for EddyViscosityModel = 3
 
 [RANS]
 EddyViscosityModel = 3
-MinDissipation = 1e-8
 UseStoredEddyViscosity = false
-# LowReKEpsilonModel = 1
-# StartWithZeroVelocity = true
 
 [Assembly]
 NumericDifferenceMethod = 0
 NumericDifference.BaseEpsilon = 1e-8
 
 [Newton]
-MaxSteps = 8
-TargetSteps = 5
-MaxRelativeShift = 1e-5
-# MaxTimeStepDivisions = 2
+MaxSteps = 7
+TargetSteps = 6
+MaxRelativeShift = 1e-8
+UseLineSearch = true
 
 [Vtk]
 AddVelocity = true

test/freeflow/rans/test_pipe_laufer_reference.input

 [TimeLoop]
-DtInitial = 1e-1 # [s]
+DtInitial = 1e-3 # [s]
 TEnd = 100 # [s]
 
 [Grid]
@@ -11,7 +11,7 @@ Cells1 = 8 8
 Grading1 = 1.4 -1.4
 
 [Problem]
-Name = pipe_laufer_zeroeq_reference
+Name = pipe_laufer_reference
 InletVelocity = 2.5 # [m/s]
 EnableGravity = false
 
@@ -24,7 +24,8 @@ NumericDifference.BaseEpsilon = 1e-8
 
 [Newton]
 MaxSteps = 10
-MaxRelativeShift = 1e-5
+TargetSteps = 5
+MaxRelativeShift = 1e-6
 
 [Vtk]
 AddVelocity = true