diff --git a/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp b/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
index 921d1cef3680390c8985b1849ddf80a8d5556f2a..48af64b9ffb8a13c77bc47a5ec1722b383740e7e 100644
--- a/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
+++ b/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
@@ -50,6 +50,7 @@
 #include "VirtualFluids_GPU/PreCollisionInteractor/Probes/PlanarAverageProbe.h"
 #include "VirtualFluids_GPU/PreCollisionInteractor/Probes/WallModelProbe.h"
 #include "VirtualFluids_GPU/BoundaryConditions/BoundaryConditionFactory.h"
+#include "VirtualFluids_GPU/TurbulenceModels/TurbulenceModelFactory.h"
 
 #include "VirtualFluids_GPU/GPU/CudaMemoryManager.h"
 
@@ -161,10 +162,7 @@ void multipleLevel(const std::string& configPath)
     para->setViscosityRatio( dx*dx/dt );
     para->setDensityRatio( 1.0 );
 
-    if(para->getUseAMD())
-        para->setMainKernel("TurbulentViscosityCumulantK17CompChim");
-    else
-        para->setMainKernel("CumulantK17CompChim");
+    para->setMainKernel("TurbulentViscosityCumulantK17CompChim");
 
     para->setIsBodyForce( config.getValue<bool>("bodyForce") );
 
@@ -172,8 +170,13 @@ void multipleLevel(const std::string& configPath)
     para->setTimestepOut( uint(tOut/dt) );
     para->setTimestepEnd( uint(tEnd/dt) );
 
-    // para->setTimestepOut( 100 );
-    // para->setTimestepEnd( 100000 );
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    SPtr<TurbulenceModelFactory> tmFactory = SPtr<TurbulenceModelFactory>( new TurbulenceModelFactory(para) );
+    tmFactory->readConfigFile( config );
+    
+    // tmFactory->setTurbulenceModel(TurbulenceModel::AMD);
+    // tmFactory->setModelConstant(config.getValue<real>("SGSconstant"));
 
     /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
@@ -228,7 +231,7 @@ void multipleLevel(const std::string& configPath)
     auto cudaMemoryManager = std::make_shared<CudaMemoryManager>(para);
     auto gridGenerator = GridProvider::makeGridGenerator(gridBuilder, para, cudaMemoryManager, communicator);
 
-    Simulation sim(para, cudaMemoryManager, communicator, *gridGenerator, &bcFactory);
+    Simulation sim(para, cudaMemoryManager, communicator, *gridGenerator, &bcFactory, tmFactory);
     sim.run();
 }
 
diff --git a/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp b/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
index 0a310e500ad0b5ee262b3d5519591e0b5c4c9c8c..4f08c0cc9b4e5ac4af0d11b5eebcbc275e07cdad 100644
--- a/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
+++ b/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
@@ -4,10 +4,10 @@
 
 #include "Communication/ExchangeData27.h"
 #include "Parameter/CudaStreamManager.h"
-#include "GPU/TurbulentViscosity.h"
 #include "KernelManager/BCKernelManager.h"
 #include "KernelManager/ADKernelManager.h"
 #include "KernelManager/GridScalingKernelManager.h"
+#include "TurbulenceModels/TurbulenceModelFactory.h"
 #include "Kernel/Kernel.h"
 
 #include "CollisionStrategy.h"
@@ -36,11 +36,10 @@ void UpdateGrid27::updateGrid(int level, unsigned int t)
 
     //////////////////////////////////////////////////////////////////////////
 
-    if (para->getUseWale())
+    if (para->getUseWale()) //TODO: make WALE consistent with structure of other turbulence models
         calcMacroscopicQuantities(level);
 
-    if (para->getUseTurbulentViscosity())
-        calcTurbulentViscosity(level);
+    calcTurbulentViscosity(level);
 
     //////////////////////////////////////////////////////////////////////////
 
@@ -364,8 +363,7 @@ void  UpdateGrid27::interactWithProbes(int level, unsigned int t)
 
 void  UpdateGrid27::calcTurbulentViscosity(int level)
 {
-    if(para->getUseAMD())
-        calcTurbulentViscosityAMD(para.get(), level);
+    this->tmFactory->runTurbulenceModelKernel(level);
 }
 
 void UpdateGrid27::exchangeData(int level)
@@ -374,8 +372,8 @@ void UpdateGrid27::exchangeData(int level)
 }
 
 UpdateGrid27::UpdateGrid27(SPtr<Parameter> para, vf::gpu::Communicator &comm, SPtr<CudaMemoryManager> cudaMemoryManager,
-                           std::vector<std::shared_ptr<PorousMedia>> &pm, std::vector<SPtr<Kernel>> &kernels , BoundaryConditionFactory* bcFactory)
-    : para(para), comm(comm), cudaMemoryManager(cudaMemoryManager), pm(pm), kernels(kernels)
+                           std::vector<std::shared_ptr<PorousMedia>> &pm, std::vector<SPtr<Kernel>> &kernels , BoundaryConditionFactory* bcFactory, SPtr<TurbulenceModelFactory>  tmFactory)
+    : para(para), comm(comm), cudaMemoryManager(cudaMemoryManager), pm(pm), kernels(kernels), tmFactory(tmFactory)
 {
     this->collision = getFunctionForCollisionAndExchange(para->getUseStreams(), para->getNumprocs(), para->getKernelNeedsFluidNodeIndicesToRun());
     this->refinement = getFunctionForRefinementAndExchange(para->getUseStreams(), para->getNumprocs(), para->getMaxLevel(), para->useReducedCommunicationAfterFtoC);
diff --git a/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.h b/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.h
index 9001b3db2e7c1b8f9e790a7686090d3af89484a4..2f0779c6ae2b1461f01383cf9fef9d17a4a38b01 100644
--- a/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.h
+++ b/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.h
@@ -10,9 +10,11 @@
 
 class BCKernelManager;
 class ADKernelManager;
+class TurbulenceModelManager;
 class GridScalingKernelManager;
 class Kernel;
 class BoundaryConditionFactory;
+class TurbulenceModelFactory;
 
 class UpdateGrid27;
 using CollisionStrategy = std::function<void (UpdateGrid27* updateGrid, Parameter* para, int level, unsigned int t)>;
@@ -23,7 +25,7 @@ class UpdateGrid27
 {
 public:
     UpdateGrid27(SPtr<Parameter> para, vf::gpu::Communicator &comm, SPtr<CudaMemoryManager> cudaMemoryManager,
-                 std::vector<std::shared_ptr<PorousMedia>> &pm, std::vector<SPtr<Kernel>> &kernels, BoundaryConditionFactory* bcFactory);
+                 std::vector<std::shared_ptr<PorousMedia>> &pm, std::vector<SPtr<Kernel>> &kernels, BoundaryConditionFactory* bcFactory, SPtr<TurbulenceModelFactory> tmFactory);
     void updateGrid(int level, unsigned int t);
     void exchangeData(int level);
 
@@ -79,6 +81,8 @@ private:
     std::shared_ptr<ADKernelManager> adKernelManager;
     //! \property gridScalingKernelManager is a shared pointer to an object of GridScalingKernelManager
     std::shared_ptr<GridScalingKernelManager> gridScalingKernelManager;
+    //! \property tmFactory is a shared pointer to an object of TurbulenceModelFactory
+    std::shared_ptr<TurbulenceModelFactory> tmFactory;
 };
 
 #endif
diff --git a/src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosity.h b/src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosity.h
deleted file mode 100644
index c0f104f408469ebb65bca5ada4c53ba9a6b13829..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosity.h
+++ /dev/null
@@ -1,10 +0,0 @@
-#ifndef TURBULENT_VISCOSITY_H_
-#define TURBULENT_VISCOSITY_H_
-
-#include "Core/DataTypes.h"
-
-class Parameter;
-
-void calcTurbulentViscosityAMD(Parameter* para, int level);
-
-#endif //TURBULENT_VISCOSITY_H_
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosityInlines.cuh b/src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosityInlines.cuh
new file mode 100644
index 0000000000000000000000000000000000000000..eb301515527a9e8a3056676b0d4dffe8197c7dbe
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosityInlines.cuh
@@ -0,0 +1,61 @@
+//=======================================================================================
+// ____          ____    __    ______     __________   __      __       __        __         
+// \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |        
+//  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |        
+//   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |        
+//    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____    
+//     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|   
+//      \    \  |    |   ________________________________________________________________    
+//       \    \ |    |  |  ______________________________________________________________|   
+//        \    \|    |  |  |         __          __     __     __     ______      _______    
+//         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)   
+//          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______    
+//           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+//            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/   
+//
+//  This file is part of VirtualFluids. VirtualFluids is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file TurbulentViscosity.h
+//! \ingroup GPU
+//! \author Henry Korb, Henrik Asmuth
+//======================================================================================
+
+#ifndef TURBULENT_VISCOSITY_INLINES_CUH_
+#define TURBULENT_VISCOSITY_INLINES_CUH_
+
+#include <cuda.h>
+#include <cuda_runtime.h>
+
+#include "LBM/LB.h" 
+#include <lbm/constants/NumericConstants.h>
+
+using namespace vf::lbm::constant;
+
+__inline__ __device__ real calcTurbulentViscositySmagorinsky(real Cs, real dxux, real dyuy, real dzuz, real Dxy, real Dxz , real Dyz)
+{
+    return Cs*Cs * sqrt( c2o1 * ( dxux*dxux + dyuy*dyuy + dzuz*dzuz ) + Dxy*Dxy + Dxz*Dxz + Dyz*Dyz );
+}
+
+__inline__ __device__ real calcTurbulentViscosityQR(real C, real dxux, real dyuy, real dzuz, real Dxy, real Dxz , real Dyz)
+{
+        // ! Verstappen's QR model
+        //! Second invariant of the strain-rate tensor
+        real Q = c1o2*( dxux*dxux + dyuy*dyuy + dzuz*dzuz ) + c1o4*( Dxy*Dxy + Dxz*Dxz + Dyz*Dyz);
+        //! Third invariant of the strain-rate tensor (determinant)
+        real R = - dxux*dyuy*dzuz - c1o4*( Dxy*Dxz*Dyz + dxux*Dyz*Dyz + dyuy*Dxz*Dxz + dzuz*Dxy*Dxy );
+        
+        return C * max(R, c0o1) / Q;
+}
+
+#endif //TURBULENT_VISCOSITY_H_e
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosity.cu b/src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosityKernels.cu
similarity index 56%
rename from src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosity.cu
rename to src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosityKernels.cu
index ed7a199e414709c6a3beff69374989fef2884dc2..df1dc571cf9627d84e940b2e0f53d55216ca6532 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosity.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosityKernels.cu
@@ -1,4 +1,37 @@
-#include "TurbulentViscosity.h"
+//=======================================================================================
+// ____          ____    __    ______     __________   __      __       __        __         
+// \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |        
+//  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |        
+//   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |        
+//    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____    
+//     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|   
+//      \    \  |    |   ________________________________________________________________    
+//       \    \ |    |  |  ______________________________________________________________|   
+//        \    \|    |  |  |         __          __     __     __     ______      _______    
+//         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)   
+//          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______    
+//           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+//            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/   
+//
+//  This file is part of VirtualFluids. VirtualFluids is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file TurbulentViscosityKernels.cu
+//! \ingroup GPU
+//! \author Henry Korb, Henrik Asmuth
+//======================================================================================
+
+#include "TurbulentViscosityKernels.h"
 #include "lbm/constants/NumericConstants.h"
 #include "Parameter/Parameter.h"
 #include "cuda/CudaGrid.h"
@@ -90,4 +123,13 @@ void calcTurbulentViscosityAMD(Parameter* para, int level)
     );
     getLastCudaError("calcAMD execution failed");
 }
-    
\ No newline at end of file
+    
+__inline__ __device__ real calcTurbulentViscosityQR(real C, real dxux, real dyuy, real dzuz, real Dxy, real Dxz , real Dyz)
+{
+        // ! Verstappen's QR model
+        //! Second invariant of the strain-rate tensor
+        real Q = c1o2*( dxux*dxux + dyuy*dyuy + dzuz*dzuz ) + c1o4*( Dxy*Dxy + Dxz*Dxz + Dyz*Dyz);
+        //! Third invariant of the strain-rate tensor (determinant)
+        real R = - dxux*dyuy*dzuz - c1o4*( Dxy*Dxz*Dyz + dxux*Dyz*Dyz + dyuy*Dxz*Dxz + dzuz*Dxy*Dxy );
+        return C * max(R, c0o1) / Q;
+}
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosityKernels.h b/src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosityKernels.h
new file mode 100644
index 0000000000000000000000000000000000000000..b227e680301cd4639d48a5cf3ce74f08eb7e1b9f
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosityKernels.h
@@ -0,0 +1,52 @@
+//=======================================================================================
+// ____          ____    __    ______     __________   __      __       __        __         
+// \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |        
+//  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |        
+//   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |        
+//    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____    
+//     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|   
+//      \    \  |    |   ________________________________________________________________    
+//       \    \ |    |  |  ______________________________________________________________|   
+//        \    \|    |  |  |         __          __     __     __     ______      _______    
+//         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)   
+//          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______    
+//           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+//            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/   
+//
+//  This file is part of VirtualFluids. VirtualFluids is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file TurbulentViscosityKernels.h
+//! \ingroup GPU
+//! \author Henry Korb, Henrik Asmuth
+//======================================================================================
+
+#ifndef TURBULENT_VISCOSITY_KERNELS_H_
+#define TURBULENT_VISCOSITY_KERNELS_H_
+
+#include <cuda.h>
+#include <cuda_runtime.h>
+
+#include "LBM/LB.h" 
+#include "Core/DataTypes.h"
+#include <lbm/constants/NumericConstants.h>
+
+using namespace vf::lbm::constant;
+
+class Parameter;
+
+void calcTurbulentViscosityAMD(Parameter* para, int level);
+
+
+
+#endif //TURBULENT_VISCOSITY_H_e
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim.cu b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim.cu
index c1be283f7fb0e5585c28bdaf5e4519f468c32c8f..fb39d57db2101e1edecf1a7707b5bc9c306c6b1c 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim.cu
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim.cu
@@ -12,27 +12,37 @@ void TurbulentViscosityCumulantK17CompChim::run()
 {
 	vf::cuda::CudaGrid grid = vf::cuda::CudaGrid(para->getParH(level)->numberofthreads, para->getParH(level)->numberOfNodes);
 
-	LB_Kernel_TurbulentViscosityCumulantK17CompChim <<< grid.grid, grid.threads >>>(
-		para->getParD(level)->omega,
-		para->getParD(level)->typeOfGridNode,
-		para->getParD(level)->neighborX,
-		para->getParD(level)->neighborY,
-		para->getParD(level)->neighborZ,
-		para->getParD(level)->distributions.f[0],
-		para->getParD(level)->rho,
-		para->getParD(level)->velocityX,
-		para->getParD(level)->velocityY,
-		para->getParD(level)->velocityZ,
-		para->getParD(level)->turbViscosity,
-		(unsigned long)para->getParD(level)->numberOfNodes,
-		level,
-		para->getIsBodyForce(),
-		para->getForcesDev(),
-		para->getParD(level)->forceX_SP,
-		para->getParD(level)->forceY_SP,
-		para->getParD(level)->forceZ_SP,
-        para->getQuadricLimitersDev(),
-		para->getParD(level)->isEvenTimestep);
+	TurbulenceModel turbulenceModel = para->getTurbulenceModel();
+	switch(para->getTurbulenceModel())
+	{
+		case TurbulenceModel::None: 		
+			LB_Kernel_TurbulentViscosityCumulantK17CompChim < TurbulenceModel::None > <<< grid.grid, grid.threads >>>(  para->getParD(level)->omega, 	para->getParD(level)->typeOfGridNode, 	para->getParD(level)->neighborX,	para->getParD(level)->neighborY,	para->getParD(level)->neighborZ,	para->getParD(level)->distributions.f[0],	
+																														para->getParD(level)->rho,		para->getParD(level)->velocityX,		para->getParD(level)->velocityY,	para->getParD(level)->velocityZ,	para->getParD(level)->turbViscosity,para->getSGSConstant(),
+																														(unsigned long)para->getParD(level)->numberOfNodes,	level,				para->getIsBodyForce(),				para->getForcesDev(),				para->getParD(level)->forceX_SP,	para->getParD(level)->forceY_SP,
+																														para->getParD(level)->forceZ_SP,para->getQuadricLimitersDev(),			para->getParD(level)->isEvenTimestep);
+			break;
+		case TurbulenceModel::AMD: 		
+			LB_Kernel_TurbulentViscosityCumulantK17CompChim < TurbulenceModel::AMD  > <<< grid.grid, grid.threads >>>(  para->getParD(level)->omega, 	para->getParD(level)->typeOfGridNode, 	para->getParD(level)->neighborX,	para->getParD(level)->neighborY,	para->getParD(level)->neighborZ,	para->getParD(level)->distributions.f[0],	
+																														para->getParD(level)->rho,		para->getParD(level)->velocityX,		para->getParD(level)->velocityY,	para->getParD(level)->velocityZ,	para->getParD(level)->turbViscosity,para->getSGSConstant(),
+																														(unsigned long)para->getParD(level)->numberOfNodes,	level,				para->getIsBodyForce(),				para->getForcesDev(),				para->getParD(level)->forceX_SP,	para->getParD(level)->forceY_SP,
+																														para->getParD(level)->forceZ_SP,para->getQuadricLimitersDev(),			para->getParD(level)->isEvenTimestep);
+			break;						
+		case TurbulenceModel::Smagorinsky: 	
+			LB_Kernel_TurbulentViscosityCumulantK17CompChim < TurbulenceModel::AMD  > <<< grid.grid, grid.threads >>>(  para->getParD(level)->omega, 	para->getParD(level)->typeOfGridNode, 	para->getParD(level)->neighborX,	para->getParD(level)->neighborY,	para->getParD(level)->neighborZ,	para->getParD(level)->distributions.f[0],	
+																														para->getParD(level)->rho,		para->getParD(level)->velocityX,		para->getParD(level)->velocityY,	para->getParD(level)->velocityZ,	para->getParD(level)->turbViscosity,para->getSGSConstant(),
+																														(unsigned long)para->getParD(level)->numberOfNodes,	level,				para->getIsBodyForce(),				para->getForcesDev(),				para->getParD(level)->forceX_SP,	para->getParD(level)->forceY_SP,
+																														para->getParD(level)->forceZ_SP,para->getQuadricLimitersDev(),			para->getParD(level)->isEvenTimestep);
+			break;
+		case TurbulenceModel::QR: 	
+			LB_Kernel_TurbulentViscosityCumulantK17CompChim < TurbulenceModel::QR  > <<< grid.grid, grid.threads >>>(  para->getParD(level)->omega, 	para->getParD(level)->typeOfGridNode, 	para->getParD(level)->neighborX,	para->getParD(level)->neighborY,	para->getParD(level)->neighborZ,	para->getParD(level)->distributions.f[0],	
+																														para->getParD(level)->rho,		para->getParD(level)->velocityX,		para->getParD(level)->velocityY,	para->getParD(level)->velocityZ,	para->getParD(level)->turbViscosity,para->getSGSConstant(),
+																														(unsigned long)para->getParD(level)->numberOfNodes,	level,				para->getIsBodyForce(),				para->getForcesDev(),				para->getParD(level)->forceX_SP,	para->getParD(level)->forceY_SP,
+																														para->getParD(level)->forceZ_SP,para->getQuadricLimitersDev(),			para->getParD(level)->isEvenTimestep);
+			break;
+		default:
+			throw std::runtime_error("TurbulentViscosityCumulantK17CompChim: Invalid turbulence mpdel ");
+			break;
+	}
 	getLastCudaError("LB_Kernel_TurbulentViscosityCumulantK17CompChim execution failed");
 }
 
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim_Device.cu b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim_Device.cu
index 685d5742c5bbd3a25632b8a7ad8c6f5be0857643..af7a67ce3f0b0cd9c4c6279d0728ab838c76fd4a 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim_Device.cu
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim_Device.cu
@@ -42,6 +42,7 @@
 #include "LBM/LB.h" 
 #include "lbm/constants/D3Q27.h"
 #include <lbm/constants/NumericConstants.h>
+#include "GPU/TurbulentViscosityInlines.cuh"
 
 using namespace vf::lbm::constant;
 using namespace vf::lbm::dir;
@@ -49,6 +50,7 @@ using namespace vf::lbm::dir;
 
 
 ////////////////////////////////////////////////////////////////////////////////
+template<TurbulenceModel turbulenceModel>
 __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
 	real omega_in,
 	uint* typeOfGridNode,
@@ -61,6 +63,7 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
     real* vy,
     real* vz,
     real* turbulentViscosity,
+    real SGSconstant,
 	unsigned long size_Mat,
 	int level,
     bool bodyForce,
@@ -466,23 +469,21 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
         real dyuy = dxux + omega * c3o2 * mxxMyy;
         real dzuz = dxux + omega * c3o2 * mxxMzz;
 
-        //Smagorinsky for debugging
-        // if(true)
-        // {   
-            // if(false && k==99976)
-            // {
-            //     printf("dudz+dwdu: \t %1.14f \n", Dxz );
-            //     printf("dvdz+dudy: \t %1.14f \n", Dxy );  
-            //     printf("dwdy+dvdz: \t %1.14f \n", Dyz );  
-            //     printf("nu_t * dudz+dwdu: \t %1.14f \n", turbulentViscosity[k]*Dxz );
-            //     printf("nu_t * dvdz+dudy: \t %1.14f \n", turbulentViscosity[k]*Dxy );  
-            //     printf("nu_t * dwdy+dvdz: \t %1.14f \n", turbulentViscosity[k]*Dyz );      
-            // } 
-        //     real Sbar = sqrt(c2o1*(dxux*dxux+dyuy*dyuy+dzuz*dzuz)+Dxy*Dxy+Dxz*Dxz+Dyz*Dyz);
-        //     real Cs = 0.08f;
-        //     turbulentViscosity[k] = Cs*Cs*Sbar;
-        // }
-
+        ////////////////////////////////////////////////////////////////////////////////////
+        switch (turbulenceModel)
+        {
+        case TurbulenceModel::None:
+        case TurbulenceModel::AMD:  //AMD is computed in separate kernel
+            break;
+        case TurbulenceModel::Smagorinsky:
+            turbulentViscosity[k] = calcTurbulentViscositySmagorinsky(SGSconstant, dxux, dyuy, dzuz, Dxy, Dxz , Dyz);
+            break;
+        case TurbulenceModel::QR:
+            turbulentViscosity[k] = calcTurbulentViscosityQR(SGSconstant, dxux, dyuy, dzuz, Dxy, Dxz , Dyz);
+            break;
+        default:
+            break;
+        }
         ////////////////////////////////////////////////////////////
         //! - Relaxation of second order cumulants with correction terms according to Eq. (33)-(35) in
         //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017),
@@ -724,4 +725,12 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
 
 
     }
-}
\ No newline at end of file
+}
+
+template __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim < TurbulenceModel::AMD > ( real omega_in, uint* typeOfGridNode, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long size_Mat, int level, bool bodyForce, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep);
+
+template __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim < TurbulenceModel::Smagorinsky > ( real omega_in, uint* typeOfGridNode, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long size_Mat, int level, bool bodyForce, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep);
+
+template __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim < TurbulenceModel::QR > ( real omega_in, uint* typeOfGridNode, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long size_Mat, int level, bool bodyForce, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep);
+
+template __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim < TurbulenceModel::None > ( real omega_in, uint* typeOfGridNode, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long size_Mat, int level, bool bodyForce, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep);
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim_Device.cuh b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim_Device.cuh
index 3633dfa4e18d057d60acf68faf348d3c5f5855d8..5ef37557399f263d25edf03b02b00f6a03c6e1cb 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim_Device.cuh
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim_Device.cuh
@@ -4,7 +4,7 @@
 #include <DataTypes.h>
 #include <curand.h>
 
-__global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
+template< TurbulenceModel turbulenceModel > __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
 	real omega_in,
 	uint* typeOfGridNode,
 	uint* neighborX,
@@ -16,6 +16,7 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
     real* vy,
     real* vz,
 	real* turbulentViscosity,
+	real SGSconstant,
 	unsigned long size_Mat,
 	int level,
 	bool bodyForce,
diff --git a/src/gpu/VirtualFluids_GPU/LBM/LB.h b/src/gpu/VirtualFluids_GPU/LBM/LB.h
index 10bfa5977f7cccb421d05f94c876a185cd667fcb..eea4adfda3c1ef0862f39ef58fc6e065af7bab1b 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/LB.h
+++ b/src/gpu/VirtualFluids_GPU/LBM/LB.h
@@ -50,6 +50,19 @@
 #include <string>
 #include <vector>
 
+//! \brief An enumeration for selecting a turbulence model
+enum class TurbulenceModel {
+   //! - Smagorinsky
+    Smagorinsky,
+    //! - AMD (Anisotropic Minimum Dissipation) model, see e.g. Rozema et al., Phys. Fluids 27, 085107 (2015), https://doi.org/10.1063/1.4928700
+    AMD,
+    //! - QR model by Verstappen 
+    QR,
+    //! - TODO: move the WALE model here from the old kernels
+    //WALE
+    //! - No turbulence model
+    None
+};
 
 struct InitCondition
 {
@@ -122,8 +135,8 @@ struct InitCondition
    bool calcMedian {false};
    bool isConc {false};
    bool isWale {false};
+   TurbulenceModel turbulenceModel {TurbulenceModel::None};
    bool isTurbulentViscosity {false};
-   bool isAMD {false};
    real SGSConstant {0.0};
    bool isMeasurePoints {false};
    bool isInitNeq {false};
@@ -285,6 +298,8 @@ typedef struct PLP{
 	uint memSizeID, memSizeTimestep, memSizerealAll, memSizereal, memSizeBool, memSizeBoolBC;
 }PathLineParticles;
 
+
+
 //////////////////////////////////////////////////////////////////////////
 inline int vectorPosition(int i, int j, int k, int Lx, int Ly )
 {
diff --git a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
index 92e6883bfe5f4fc2fcab4dbbd98c164a38bf7192..811dd3796b849a9ba0f10b54308969a8e7038f32 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
+++ b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
@@ -55,7 +55,7 @@
 #include "PreProcessor/PreProcessorFactory/PreProcessorFactoryImp.h"
 #include "Kernel/Utilities/KernelFactory/KernelFactoryImp.h"
 #include "Kernel/Kernel.h"
-
+#include "TurbulenceModels/TurbulenceModelFactory.h"
 #include <cuda/DeviceInfo.h>
 
 #include <logger/Logger.h>
@@ -72,10 +72,19 @@ Simulation::Simulation(std::shared_ptr<Parameter> para, std::shared_ptr<CudaMemo
     : para(para), cudaMemoryManager(memoryManager), communicator(communicator), kernelFactory(std::make_unique<KernelFactoryImp>()),
       preProcessorFactory(std::make_shared<PreProcessorFactoryImp>()), dataWriter(std::make_unique<FileWriter>())
 {
-	init(gridProvider, bcFactory);
+	this->tmFactory = SPtr<TurbulenceModelFactory>( new TurbulenceModelFactory(para) );
+	init(gridProvider, bcFactory, tmFactory);
+}
+
+Simulation::Simulation(std::shared_ptr<Parameter> para, std::shared_ptr<CudaMemoryManager> memoryManager,
+                       vf::gpu::Communicator &communicator, GridProvider &gridProvider, BoundaryConditionFactory* bcFactory, SPtr<TurbulenceModelFactory> tmFactory)
+    : para(para), cudaMemoryManager(memoryManager), communicator(communicator), kernelFactory(std::make_unique<KernelFactoryImp>()),
+      preProcessorFactory(std::make_shared<PreProcessorFactoryImp>()), dataWriter(std::make_unique<FileWriter>())
+{
+	init(gridProvider, bcFactory, tmFactory);
 }
 
-void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFactory)
+void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFactory, SPtr<TurbulenceModelFactory> tmFactory)
 {
     gridProvider.initalGridInformations();
 
@@ -138,8 +147,6 @@ void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFa
     output << "delta_rho:  " << para->getDensityRatio() << "\n";
     output << "QuadricLimiters:  " << para->getQuadricLimitersHost()[0] << "\t" << para->getQuadricLimitersHost()[1]
            << "\t" << para->getQuadricLimitersHost()[2] << "\n";
-    if (para->getUseAMD())
-        output << "AMD SGS model:  " << para->getSGSConstant() << "\n";
     //////////////////////////////////////////////////////////////////////////
 
     /////////////////////////////////////////////////////////////////////////
@@ -382,7 +389,7 @@ void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFa
     //////////////////////////////////////////////////////////////////////////
     // Init UpdateGrid
     //////////////////////////////////////////////////////////////////////////
-    this->updateGrid27 = std::make_unique<UpdateGrid27>(para, communicator, cudaMemoryManager, pm, kernels, bcFactory);
+    this->updateGrid27 = std::make_unique<UpdateGrid27>(para, communicator, cudaMemoryManager, pm, kernels, bcFactory, tmFactory);
 
     //////////////////////////////////////////////////////////////////////////
     // Write Initialized Files
diff --git a/src/gpu/VirtualFluids_GPU/LBM/Simulation.h b/src/gpu/VirtualFluids_GPU/LBM/Simulation.h
index 98dca98747d0eb6603812f75428242e05edd23e4..9220b92b26ac7f3f665293c017aed400aa29d976 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/Simulation.h
+++ b/src/gpu/VirtualFluids_GPU/LBM/Simulation.h
@@ -30,12 +30,16 @@ class UpdateGrid27;
 class KineticEnergyAnalyzer;
 class EnstrophyAnalyzer;
 class BoundaryConditionFactory;
+class TurbulenceModelFactory;
 
 class Simulation
 {
 public:
     Simulation(std::shared_ptr<Parameter> para, std::shared_ptr<CudaMemoryManager> memoryManager,
                vf::gpu::Communicator &communicator, GridProvider &gridProvider, BoundaryConditionFactory* bcFactory);
+	
+	Simulation(std::shared_ptr<Parameter> para, std::shared_ptr<CudaMemoryManager> memoryManager,
+               vf::gpu::Communicator &communicator, GridProvider &gridProvider, BoundaryConditionFactory* bcFactory, SPtr<TurbulenceModelFactory> tmFactory);
 
     ~Simulation();
     void run();
@@ -47,7 +51,7 @@ public:
     void addEnstrophyAnalyzer(uint tAnalyse);
 
 private:
-	void init(GridProvider &gridProvider, BoundaryConditionFactory *bcFactory);
+	void init(GridProvider &gridProvider, BoundaryConditionFactory *bcFactory, SPtr<TurbulenceModelFactory> tmFactory);
     void allocNeighborsOffsetsScalesAndBoundaries(GridProvider& gridProvider);
     void porousMedia();
     void definePMarea(std::shared_ptr<PorousMedia>& pm);
@@ -75,6 +79,7 @@ private:
 	std::vector < SPtr< Kernel>> kernels;
 	std::vector < SPtr< ADKernel>> adKernels;
 	std::shared_ptr<PreProcessor> preProcessor;
+	SPtr<TurbulenceModelFactory> tmFactory;
 
     SPtr<RestartObject> restart_object;
 
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
index 31fe3bef2a2916f0920aef2321126915fe0fe003..3118d17801500e63fa7bff9bd98f55b2f02dbd86 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
@@ -125,12 +125,6 @@ void Parameter::readConfigData(const vf::basics::ConfigurationFile &configData)
     if (configData.contains("UseWale"))
         this->setUseWale(configData.getValue<bool>("UseWale"));
     //////////////////////////////////////////////////////////////////////////
-    if (configData.contains("UseAMD"))
-        this->setUseAMD(configData.getValue<bool>("UseAMD"));
-    //////////////////////////////////////////////////////////////////////////
-    if (configData.contains("SGSconstant"))
-        this->setSGSConstant(configData.getValue<real>("SGSconstant"));
-    //////////////////////////////////////////////////////////////////////////
     if (configData.contains("UseInitNeq"))
         this->setUseInitNeq(configData.getValue<bool>("UseInitNeq"));
     //////////////////////////////////////////////////////////////////////////
@@ -935,11 +929,9 @@ void Parameter::setUseWale(bool useWale)
     if (useWale)
         setUseTurbulentViscosity(true);
 }
-void Parameter::setUseAMD(bool useAMD)
+void Parameter::setTurbulenceModel(TurbulenceModel turbulenceModel)
 {
-    ic.isAMD = useAMD;
-    if (useAMD)
-        setUseTurbulentViscosity(true);
+    ic.turbulenceModel = turbulenceModel;
 }
 void Parameter::setSGSConstant(real SGSConstant)
 {
@@ -2326,9 +2318,9 @@ bool Parameter::getUseWale()
 {
     return ic.isWale;
 }
-bool Parameter::getUseAMD()
+TurbulenceModel Parameter::getTurbulenceModel()
 {
-    return ic.isAMD;
+    return ic.turbulenceModel;
 }
 bool Parameter::getUseTurbulentViscosity()
 {
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
index 941e30413ff83eb18b2cdf2a511848afdc886b70..e5935336b5d800c8b132c9a54f7dc8f94e1c7c7d 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
@@ -42,6 +42,7 @@
 #include "lbm/constants/D3Q27.h"
 #include "LBM/LB.h"
 #include "PreCollisionInteractor/PreCollisionInteractor.h"
+#include "TurbulenceModels/TurbulenceModelFactory.h"
 
 #include "VirtualFluids_GPU_export.h"
 
@@ -520,6 +521,7 @@ public:
     void setStreetVelocityFile(bool streetVelocityFile);
     void setUseMeasurePoints(bool useMeasurePoints);
     void setUseWale(bool useWale);
+    void setTurbulenceModel(TurbulenceModel turbulenceModel);
     void setUseTurbulentViscosity(bool useTurbulentViscosity);
     void setUseAMD(bool useAMD);
     void setSGSConstant(real SGSConstant);
@@ -800,8 +802,8 @@ public:
     bool isStreetVelocityFile();
     bool getUseMeasurePoints();
     bool getUseWale();
+    TurbulenceModel getTurbulenceModel();
     bool getUseTurbulentViscosity();
-    bool getUseAMD();
     real getSGSConstant();
     bool getHasWallModelMonitor();
     bool getUseInitNeq();
diff --git a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/PlanarAverageProbe.cu b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/PlanarAverageProbe.cu
index d762717a49d122484bf07970e32df0797be17e7b..60932417ef9a1939fdaf935550f202e4a45ec8d2 100644
--- a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/PlanarAverageProbe.cu
+++ b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/PlanarAverageProbe.cu
@@ -142,14 +142,16 @@ std::vector<PostProcessingVariable> PlanarAverageProbe::getPostProcessingVariabl
     switch (statistic)
     {
     case Statistic::SpatialMeans:
-        postProcessingVariables.push_back( PostProcessingVariable("vx_spatMean",  velocityRatio) );
+        postProcessingVariables.push_back( PostProcessingVariable("vx_spatMean",  this->velocityRatio) );
         postProcessingVariables.push_back( PostProcessingVariable("vy_spatMean",  this->velocityRatio) );
         postProcessingVariables.push_back( PostProcessingVariable("vz_spatMean",  this->velocityRatio) );
+        postProcessingVariables.push_back( PostProcessingVariable("nut_spatMean", this->viscosityRatio) );
         break;
     case Statistic::SpatioTemporalMeans:
         postProcessingVariables.push_back( PostProcessingVariable("vx_spatTmpMean",  this->velocityRatio) );
         postProcessingVariables.push_back( PostProcessingVariable("vy_spatTmpMean",  this->velocityRatio) );
         postProcessingVariables.push_back( PostProcessingVariable("vz_spatTmpMean",  this->velocityRatio) );
+        postProcessingVariables.push_back( PostProcessingVariable("nut_spatTmpMean", this->viscosityRatio) );
         break;
     case Statistic::SpatialCovariances:
         postProcessingVariables.push_back( PostProcessingVariable("vxvx_spatMean",  pow(this->velocityRatio, 2.0)) );
@@ -294,6 +296,7 @@ void PlanarAverageProbe::calculateQuantities(SPtr<ProbeStruct> probeStruct, Para
     thrust::device_ptr<real> vx_thrust = thrust::device_pointer_cast(para->getParD(level)->velocityX);
     thrust::device_ptr<real> vy_thrust = thrust::device_pointer_cast(para->getParD(level)->velocityY);
     thrust::device_ptr<real> vz_thrust = thrust::device_pointer_cast(para->getParD(level)->velocityZ);
+    thrust::device_ptr<real> nut_thrust = thrust::device_pointer_cast(para->getParD(level)->turbViscosity);
 
     real N = (real)probeStruct->nIndices;
     real n = (real)probeStruct->vals;
@@ -307,6 +310,8 @@ void PlanarAverageProbe::calculateQuantities(SPtr<ProbeStruct> probeStruct, Para
     thrust::permutation_iterator<valIterator, indIterator> vy_iter_end  (vy_thrust, indices_thrust+probeStruct->nIndices);
     thrust::permutation_iterator<valIterator, indIterator> vz_iter_begin(vz_thrust, indices_thrust);
     thrust::permutation_iterator<valIterator, indIterator> vz_iter_end  (vz_thrust, indices_thrust+probeStruct->nIndices);
+    thrust::permutation_iterator<valIterator, indIterator> nut_iter_begin(nut_thrust, indices_thrust);
+    thrust::permutation_iterator<valIterator, indIterator> nut_iter_end  (nut_thrust, indices_thrust+probeStruct->nIndices);
 
     for( uint i=0; i<nPoints; i++ )
     {
@@ -319,10 +324,14 @@ void PlanarAverageProbe::calculateQuantities(SPtr<ProbeStruct> probeStruct, Para
             real spatMean_vy = thrust::reduce(vy_iter_begin, vy_iter_end)/N;
             real spatMean_vz = thrust::reduce(vz_iter_begin, vz_iter_end)/N;
 
+            real spatMean_nut;
+            if(para->getUseTurbulentViscosity()) spatMean_nut = thrust::reduce(nut_iter_begin, nut_iter_end)/N;
+
             uint arrOff = probeStruct->arrayOffsetsH[int(Statistic::SpatialMeans)];
             probeStruct->quantitiesArrayH[(arrOff+0)*nPoints+node] = spatMean_vx;
             probeStruct->quantitiesArrayH[(arrOff+1)*nPoints+node] = spatMean_vy;
             probeStruct->quantitiesArrayH[(arrOff+2)*nPoints+node] = spatMean_vz;
+            if(para->getUseTurbulentViscosity()) probeStruct->quantitiesArrayH[(arrOff+3)*nPoints+node] = spatMean_nut;
 
             if(probeStruct->quantitiesH[int(Statistic::SpatioTemporalMeans)] && doTmpAveraging)
             {
@@ -330,10 +339,14 @@ void PlanarAverageProbe::calculateQuantities(SPtr<ProbeStruct> probeStruct, Para
             real spatTmpMean_vx_old = probeStruct->quantitiesArrayH[(arrOff+0)*nPoints+node];
             real spatTmpMean_vy_old = probeStruct->quantitiesArrayH[(arrOff+1)*nPoints+node];
             real spatTmpMean_vz_old = probeStruct->quantitiesArrayH[(arrOff+2)*nPoints+node];
+            real spatTmpMean_nut_old;
+            if(para->getUseTurbulentViscosity()) spatTmpMean_nut_old = probeStruct->quantitiesArrayH[(arrOff+3)*nPoints+node];;
 
             probeStruct->quantitiesArrayH[(arrOff+0)*nPoints+node] += (spatMean_vx-spatTmpMean_vx_old)/n;
             probeStruct->quantitiesArrayH[(arrOff+1)*nPoints+node] += (spatMean_vy-spatTmpMean_vy_old)/n;
             probeStruct->quantitiesArrayH[(arrOff+2)*nPoints+node] += (spatMean_vz-spatTmpMean_vz_old)/n;
+            if(para->getUseTurbulentViscosity()) probeStruct->quantitiesArrayH[(arrOff+3)*nPoints+node] += (spatMean_nut-spatTmpMean_nut_old)/n;
+
             }
         
             if(probeStruct->quantitiesH[int(Statistic::SpatialCovariances)])
diff --git a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/Probe.cu b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/Probe.cu
index 5c0d76cf892d581baade960b7836ced10ca7b591..e13f411273a824cfee4d52c712dd81245301396d 100644
--- a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/Probe.cu
+++ b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/Probe.cu
@@ -184,6 +184,7 @@ void Probe::init(Parameter* para, GridProvider* gridProvider, CudaMemoryManager*
     this->forceRatio         = para->getForceRatio();
     this->stressRatio        = para->getDensityRatio()*pow(para->getVelocityRatio(), 2.0);
     this->accelerationRatio  = para->getVelocityRatio()/para->getTimeRatio();
+    this->viscosityRatio     = para->getViscosityRatio();
 
     probeParams.resize(para->getMaxLevel()+1);
 
diff --git a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/Probe.h b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/Probe.h
index a83f4d2b614f02fc8e013c6fe3d11f40eadc9924..07ffd19c11f7c929a2c353e1bc16cece519a1db7 100644
--- a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/Probe.h
+++ b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/Probe.h
@@ -213,11 +213,13 @@ protected:
 
     uint tProbe = 0; //!> counter for number of probe evaluations. Only used when outputting timeseries
 
+
     real velocityRatio;
     real densityRatio;
     real forceRatio;
     real stressRatio;
     real accelerationRatio;
+    real viscosityRatio;
 };
 
 #endif
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/TurbulenceModels/TurbulenceModelFactory.cpp b/src/gpu/VirtualFluids_GPU/TurbulenceModels/TurbulenceModelFactory.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..980a199b03710ac91310d26424c20c7371f1097d
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/TurbulenceModels/TurbulenceModelFactory.cpp
@@ -0,0 +1,88 @@
+//=======================================================================================
+// ____          ____    __    ______     __________   __      __       __        __
+// \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
+//  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
+//   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
+//    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
+//     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
+//      \    \  |    |   ________________________________________________________________
+//       \    \ |    |  |  ______________________________________________________________|
+//        \    \|    |  |  |         __          __     __     __     ______      _______
+//         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
+//          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
+//           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+//            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
+//
+//  This file is part of VirtualFluids. VirtualFluids is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file TurbulentViscosityFactory.cpp
+//! \ingroup TurbulentViscosity
+//! \author Henrik Asmuth
+//=======================================================================================
+#include "LBM/LB.h"
+#include "TurbulenceModelFactory.h"
+#include "GPU/TurbulentViscosityKernels.h"
+#include "Parameter/Parameter.h"
+#include <logger/Logger.h>
+
+#include <variant>
+
+void TurbulenceModelFactory::setTurbulenceModel(TurbulenceModel _turbulenceModel)
+{
+    this->turbulenceModel = _turbulenceModel;
+    para->setTurbulenceModel(_turbulenceModel);
+    if(this->turbulenceModel != TurbulenceModel::None) para->setUseTurbulentViscosity(true);
+
+    switch (this->turbulenceModel) {
+        case TurbulenceModel::AMD:
+            this->turbulenceModelKernel = calcTurbulentViscosityAMD;
+            break;
+        default:
+            this->turbulenceModelKernel = nullptr;
+    }
+}
+
+void TurbulenceModelFactory::setModelConstant(real modelConstant)
+{
+    para->setSGSConstant(modelConstant);
+}
+
+void TurbulenceModelFactory::readConfigFile(const vf::basics::ConfigurationFile &configData)
+{
+    if (configData.contains("TurbulenceModel"))
+    {
+        std::string config = configData.getValue<std::string>("TurbulenceModel");
+        
+        if      (config == "Smagorinsky") this->setTurbulenceModel( TurbulenceModel::Smagorinsky ); 
+        else if (config == "AMD")         this->setTurbulenceModel( TurbulenceModel::AMD );               
+        else if (config == "QR" )         this->setTurbulenceModel( TurbulenceModel::QR );             
+        else if (config == "None")        this->setTurbulenceModel( TurbulenceModel::None );           
+        else    std::runtime_error("TurbulenceModelFactory: Invalid turbulence model!");           
+
+        VF_LOG_INFO("Turbulence model: {}", config);
+        
+    }
+
+    if (configData.contains("SGSconstant"))
+    {
+        para->setSGSConstant(configData.getValue<real>("SGSconstant"));
+
+        VF_LOG_INFO("SGS constant: {}", para->getSGSConstant() );
+    }
+}
+
+void TurbulenceModelFactory::runTurbulenceModelKernel(const int level) const
+{
+    if(this->turbulenceModelKernel) this->turbulenceModelKernel(para.get(), level);
+}
diff --git a/src/gpu/VirtualFluids_GPU/TurbulenceModels/TurbulenceModelFactory.h b/src/gpu/VirtualFluids_GPU/TurbulenceModels/TurbulenceModelFactory.h
new file mode 100644
index 0000000000000000000000000000000000000000..e71c8ed5f7be016a4a800b83cfb3252ee6b8246e
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/TurbulenceModels/TurbulenceModelFactory.h
@@ -0,0 +1,71 @@
+//=======================================================================================
+// ____          ____    __    ______     __________   __      __       __        __
+// \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
+//  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
+//   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
+//    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
+//     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
+//      \    \  |    |   ________________________________________________________________
+//       \    \ |    |  |  ______________________________________________________________|
+//        \    \|    |  |  |         __          __     __     __     ______      _______
+//         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
+//          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
+//           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+//            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
+//
+//  This file is part of VirtualFluids. VirtualFluids is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file TurbulentViscosityFactory.h
+//! \ingroup TurbulentViscosity
+//! \author Henrik Asmuth
+//=======================================================================================
+#ifndef TurbulenceModelFactory_H
+#define TurbulenceModelFactory_H
+
+#include <functional>
+#include <map>
+#include <string>
+#include <variant>
+
+#include "LBM/LB.h"
+#include "Parameter/Parameter.h"
+
+#include <basics/config/ConfigurationFile.h>
+
+class Parameter;
+
+using TurbulenceModelKernel = std::function<void(Parameter *, int )>;
+
+class TurbulenceModelFactory
+{
+public:
+    
+    TurbulenceModelFactory(SPtr<Parameter> parameter): para(parameter) {}
+
+    void setTurbulenceModel(TurbulenceModel _turbulenceModel);
+
+    void setModelConstant(real modelConstant);
+
+    void readConfigFile(const vf::basics::ConfigurationFile &configData);
+
+    void runTurbulenceModelKernel(const int level) const;
+
+private:
+    TurbulenceModel turbulenceModel = TurbulenceModel::None;
+    TurbulenceModelKernel turbulenceModelKernel = nullptr;
+    SPtr<Parameter> para;
+
+};
+
+#endif