diff --git a/apps/gpu/ActuatorLine/configActuatorLine.txt b/apps/gpu/ActuatorLine/configActuatorLine.txt
index 5799f24716777295b2f835ab00561ff767ba87b9..fbba0b89c689957831ce6ce7625da1441d8909a6 100644
--- a/apps/gpu/ActuatorLine/configActuatorLine.txt
+++ b/apps/gpu/ActuatorLine/configActuatorLine.txt
@@ -24,11 +24,11 @@ tStartOutProbe=100
 tOutProbe=100
 
 ##################################################
-#TurbulenceModel = QR
-#SGSconstant = 0.3333333
-#
-#QuadricLimiterP = 100000.0
-#QuadricLimiterM = 100000.0
-#QuadricLimiterD = 100000.0
+TurbulenceModel = QR
+SGSconstant = 0.3333333
+
+QuadricLimiterP = 10000.0
+QuadricLimiterM = 10000.0
+QuadricLimiterD = 10000.0
 ##################################################
 
diff --git a/src/cpu/VirtualFluidsCore/LBM/K17CompressibleNavierStokes.cpp b/src/cpu/VirtualFluidsCore/LBM/K17CompressibleNavierStokes.cpp
index 186574807cab3cea5dc3d46edb8111a392cf19ed..6fdc2b96ca744906a26266e0101a0145cc3928dd 100644
--- a/src/cpu/VirtualFluidsCore/LBM/K17CompressibleNavierStokes.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/K17CompressibleNavierStokes.cpp
@@ -153,11 +153,9 @@ void K17CompressibleNavierStokes::calculate(int step)
                 real quadricLimiter[3] = { 0.01, 0.01, 0.01 }; // TODO: Where do we configure the quadricLimiter?
                 parameter.quadricLimiter = quadricLimiter;
 
-                const bool writeMacroscopicVariables = false;
                 vf::lbm::MacroscopicValues mv;  // not used
                 vf::lbm::TurbulentViscosity tv; // not used
-                vf::lbm::runK17CompressibleNavierStokes<vf::lbm::TurbulenceModel::None, writeMacroscopicVariables>(parameter,
-                                                                                                                   mv, tv);
+                vf::lbm::runK17CompressibleNavierStokes<vf::lbm::TurbulenceModel::None>(parameter, mv, tv);
                 dataSet->getFdistributions()->setDistribution(parameter.distribution, x1, x2, x3);
             }
         }
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleCF_compressible.cu b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleCF_compressible.cu
index 4df1d28b148fa00f5efb81a5683d1df33b71b73d..05ec8417976bd1b71d0a9ab08b1a44ab06ffd7b0 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleCF_compressible.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleCF_compressible.cu
@@ -38,6 +38,7 @@
 
 #include <lbm/refinement/InterpolationCF.h>
 #include <lbm/interpolation/InterpolationCoefficients.h>
+#include <lbm/collision/TurbulentViscosity.h>
 
 
 
@@ -89,7 +90,7 @@ template <bool hasTurbulentViscosity> __device__ void interpolate(
     //! - Set moments (zeroth to sixth order) on destination node
     //!
     real omegaF  = omegaFine;
-    omegaF = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
+    omegaF = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
 
     const real epsilon_new = c1o2; // ratio of grid resolutions
     real f[27];
@@ -120,7 +121,7 @@ template <bool hasTurbulentViscosity> __device__ void interpolate(
     ////////////////////////////////////////////////////////////////////////////////
     // Set moments (zeroth to sixth orders) on destination node
 
-    omegaF = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
+    omegaF = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
 
     vf::lbm::interpolateCF(f, omegaF, epsilon_new, coefficients, x, y, z);
 
@@ -148,7 +149,7 @@ template <bool hasTurbulentViscosity> __device__ void interpolate(
     ////////////////////////////////////////////////////////////////////////////////
     // Set moments (zeroth to sixth orders) on destination node
 
-    omegaF = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
+    omegaF = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
 
     vf::lbm::interpolateCF(f, omegaF, epsilon_new, coefficients, x, y, z);
 
@@ -176,7 +177,7 @@ template <bool hasTurbulentViscosity> __device__ void interpolate(
     ////////////////////////////////////////////////////////////////////////////////
     // Set moments (zeroth to sixth orders) on destination node
 
-    omegaF = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
+    omegaF = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
 
     vf::lbm::interpolateCF(f, omegaF, epsilon_new, coefficients, x, y, z);
 
@@ -215,7 +216,7 @@ template <bool hasTurbulentViscosity> __device__ void interpolate(
     ////////////////////////////////////////////////////////////////////////////////
     // Set moments (zeroth to sixth orders) on destination node
 
-    omegaF = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
+    omegaF = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
 
     vf::lbm::interpolateCF(f, omegaF, epsilon_new, coefficients, x, y, z);
 
@@ -239,7 +240,7 @@ template <bool hasTurbulentViscosity> __device__ void interpolate(
     indices.k_0MM = neighborZfine[indices.k_0MM];
     indices.k_MMM = neighborZfine[indices.k_MMM];
 
-    omegaF = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
+    omegaF = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
 
     vf::lbm::interpolateCF(f, omegaF, epsilon_new, coefficients, x, y, z);
 
@@ -263,7 +264,7 @@ template <bool hasTurbulentViscosity> __device__ void interpolate(
     indices.k_0MM = indices.k_MMM;
     indices.k_MMM = neighborXfine[indices.k_MMM];
 
-    omegaF = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
+    omegaF = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
 
     vf::lbm::interpolateCF(f, omegaF, epsilon_new, coefficients, x, y, z);
 
@@ -287,7 +288,7 @@ template <bool hasTurbulentViscosity> __device__ void interpolate(
     indices.k_0M0 = k_base_MM0;
     indices.k_MM0 = neighborXfine[k_base_MM0];
 
-    omegaF = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
+    omegaF = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaFine, turbulentViscosityFine[indices.k_000]) : omegaFine;
 
     vf::lbm::interpolateCF(f, omegaF, epsilon_new, coefficients, x, y, z);
 
@@ -365,7 +366,7 @@ template<bool hasTurbulentViscosity> __global__ void scaleCF_compressible(
     indices.k_0MM = k_base_0MM;
     indices.k_MMM = k_base_MMM;
 
-    omegaC = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
+    omegaC = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
 
     real f_coarse[27];
 
@@ -385,7 +386,7 @@ template<bool hasTurbulentViscosity> __global__ void scaleCF_compressible(
     indices.k_0MM = neighborZcoarse[indices.k_0MM];
     indices.k_MMM = neighborZcoarse[indices.k_MMM];
 
-    omegaC = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
+    omegaC = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
 
     read(f_coarse, distCoarse, indices);
     momentsSet.calculateMMP(f_coarse, omegaC);
@@ -403,7 +404,7 @@ template<bool hasTurbulentViscosity> __global__ void scaleCF_compressible(
     indices.k_0MM = indices.k_MMM;
     indices.k_MMM = neighborXcoarse[indices.k_MMM];
 
-    omegaC = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
+    omegaC = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
 
     read(f_coarse, distCoarse, indices);
     momentsSet.calculatePMP(f_coarse, omegaC);
@@ -421,7 +422,7 @@ template<bool hasTurbulentViscosity> __global__ void scaleCF_compressible(
     indices.k_0M0 = k_base_MM0;
     indices.k_MM0 = neighborXcoarse[k_base_MM0];
 
-    omegaC = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
+    omegaC = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
 
     read(f_coarse, distCoarse, indices);
     momentsSet.calculatePMM(f_coarse, omegaC);
@@ -449,7 +450,7 @@ template<bool hasTurbulentViscosity> __global__ void scaleCF_compressible(
     indices.k_0MM = k_base_0MM;
     indices.k_MMM = k_base_MMM;
 
-    omegaC = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
+    omegaC = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
 
     read(f_coarse, distCoarse, indices);
     momentsSet.calculateMPM(f_coarse, omegaC);
@@ -467,7 +468,7 @@ template<bool hasTurbulentViscosity> __global__ void scaleCF_compressible(
     indices.k_0MM = neighborZcoarse[indices.k_0MM];
     indices.k_MMM = neighborZcoarse[indices.k_MMM];
 
-    omegaC = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
+    omegaC = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
 
     read(f_coarse, distCoarse, indices);
     momentsSet.calculateMPP(f_coarse, omegaC);
@@ -484,7 +485,7 @@ template<bool hasTurbulentViscosity> __global__ void scaleCF_compressible(
     indices.k_0MM = indices.k_MMM;
     indices.k_MMM = neighborXcoarse[indices.k_MMM];
 
-    omegaC = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
+    omegaC = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
 
     read(f_coarse, distCoarse, indices);
     momentsSet.calculatePPP(f_coarse, omegaC);
@@ -501,7 +502,7 @@ template<bool hasTurbulentViscosity> __global__ void scaleCF_compressible(
     indices.k_0M0 = k_base_MM0;
     indices.k_MM0 = neighborXcoarse[k_base_MM0];
 
-    omegaC = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
+    omegaC = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
 
     read(f_coarse, distCoarse, indices);
     momentsSet.calculatePPM(f_coarse, omegaC);
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleFC_compressible.cu b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleFC_compressible.cu
index c0f6abe6bc83826f9ade5860419432017f7a42c2..9decb44aab7d1f0944d98cd83fcfb081eabdabee 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleFC_compressible.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleFC_compressible.cu
@@ -60,7 +60,7 @@ template <bool hasTurbulentViscosity> __device__ void interpolate(
     vf::gpu::ListIndices indices(indicesCoarse000[nodeIndex], neighborXcoarse, neighborYcoarse, neighborZcoarse);
 
     const real epsilonNew = c2o1; // ratio of grid resolutions
-    const real omegaCoarseNew = hasTurbulentViscosity ? vf::gpu::calculateOmega(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
+    const real omegaCoarseNew = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omegaCoarse, turbulentViscosityCoarse[indices.k_000]) : omegaCoarse;
     real fCoarse[27];
     vf::lbm::interpolateFC(fCoarse, epsilonNew, omegaCoarseNew, coefficients);
 
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes.cu b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes.cu
index 02df4e149018668f830b161f86639f102ca7b12d..60d4a031451b73e1c4b8de1b317d6488d32e69e1 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes.cu
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes.cu
@@ -60,7 +60,7 @@ void K17CompressibleNavierStokes<turbulenceModel>::run()
 
     auto collision = [] __device__(vf::lbm::CollisionParameter & parameter, vf::lbm::MacroscopicValues & macroscopicValues,
                                    vf::lbm::TurbulentViscosity & turbulentViscosity) {
-        return vf::lbm::runK17CompressibleNavierStokes<turbulenceModel, false>(parameter, macroscopicValues, turbulentViscosity);
+        return vf::lbm::runK17CompressibleNavierStokes<turbulenceModel>(parameter, macroscopicValues, turbulentViscosity);
     };
 
     vf::gpu::runCollision<decltype(collision), turbulenceModel, false, false><<<cudaGrid.grid, cudaGrid.threads>>>(collision, kernelParameter);
@@ -75,24 +75,19 @@ void K17CompressibleNavierStokes<turbulenceModel>::runOnIndices(const unsigned i
 {
     cudaStream_t stream = para->getStreamManager()->getStream(streamIndex);
 
+    auto collision = [] __device__(vf::lbm::CollisionParameter& parameter, vf::lbm::MacroscopicValues& macroscopicValues, vf::lbm::TurbulentViscosity& turbulentViscosity) {
+        return vf::lbm::runK17CompressibleNavierStokes<turbulenceModel>(parameter, macroscopicValues, turbulentViscosity);
+    };
+
     switch (collisionTemplate) {
         case CollisionTemplate::Default: {
             vf::gpu::GPUCollisionParameter kernelParameter = getCollisionParameter(para, level, indices, size_indices);
-
-            auto collision = [] __device__(vf::lbm::CollisionParameter& parameter, vf::lbm::MacroscopicValues& macroscopicValues, vf::lbm::TurbulentViscosity& turbulentViscosity) {
-                return vf::lbm::runK17CompressibleNavierStokes<turbulenceModel, false>(parameter, macroscopicValues, turbulentViscosity);
-            };
-
             vf::gpu::runCollision<decltype(collision), turbulenceModel, false, false><<<cudaGrid.grid, cudaGrid.threads, 0, stream>>>(collision, kernelParameter);
 
             break;
         }
         case CollisionTemplate::WriteMacroVars: {
             vf::gpu::GPUCollisionParameter kernelParameter = getCollisionParameter(para, level, indices, size_indices);
-
-            auto collision = [] __device__(vf::lbm::CollisionParameter& parameter, vf::lbm::MacroscopicValues& macroscopicValues, vf::lbm::TurbulentViscosity& turbulentViscosity) {
-                return vf::lbm::runK17CompressibleNavierStokes<turbulenceModel, true>(parameter, macroscopicValues, turbulentViscosity);
-            };
             vf::gpu::runCollision<decltype(collision), turbulenceModel, true, false><<<cudaGrid.grid, cudaGrid.threads, 0, stream>>>(collision, kernelParameter);
 
             break;
@@ -100,20 +95,12 @@ void K17CompressibleNavierStokes<turbulenceModel>::runOnIndices(const unsigned i
         case CollisionTemplate::SubDomainBorder:
         case CollisionTemplate::AllFeatures: {
             vf::gpu::GPUCollisionParameter kernelParameter = getCollisionParameter(para, level, indices, size_indices);
-
-            auto collision = [] __device__(vf::lbm::CollisionParameter& parameter, vf::lbm::MacroscopicValues& macroscopicValues, vf::lbm::TurbulentViscosity& turbulentViscosity) {
-                return vf::lbm::runK17CompressibleNavierStokes<turbulenceModel, true>(parameter, macroscopicValues, turbulentViscosity);
-            };
             vf::gpu::runCollision<decltype(collision), turbulenceModel, true, true><<<cudaGrid.grid, cudaGrid.threads, 0, stream>>>(collision, kernelParameter);
 
             break;
         }
         case CollisionTemplate::ApplyBodyForce: {
             vf::gpu::GPUCollisionParameter kernelParameter = getCollisionParameter(para, level, indices, size_indices);
-
-            auto collision = [] __device__(vf::lbm::CollisionParameter& parameter, vf::lbm::MacroscopicValues& macroscopicValues, vf::lbm::TurbulentViscosity& turbulentViscosity) {
-                return vf::lbm::runK17CompressibleNavierStokes<turbulenceModel, false>(parameter, macroscopicValues, turbulentViscosity);
-            };
             vf::gpu::runCollision<decltype(collision), turbulenceModel, false, true><<<cudaGrid.grid, cudaGrid.threads, 0, stream>>>(collision, kernelParameter);
 
             break;
diff --git a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/KernelUtilities.h b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/KernelUtilities.h
index f67a2d752741f840f9a0189ce67299cc41e7a4d2..2403a7acd9153bf2a1b9aad4d19c5f8d5fa60766 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/KernelUtilities.h
+++ b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/KernelUtilities.h
@@ -391,11 +391,6 @@ __inline__ __device__ void writeInverse(Distributions27& destination, const List
     (destination.f[DIR_MMM])[indices.k_MMM] = f[DIR_PPP];
 }
 
-__inline__ __device__ real calculateOmega(const real omega_old, real turbulenceViscosity)
-{
-    return omega_old / (c1o1 + c3o1 * omega_old * turbulenceViscosity);
-}
-
 }
 
 #endif
diff --git a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/RunCollision.cuh b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/RunCollision.cuh
index ba1034ac59b52522c8afd881d9f91f90e701d5ae..9c5b17adf7b5c41fa774511a0797f311a6f1296d 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/RunCollision.cuh
+++ b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/RunCollision.cuh
@@ -58,9 +58,9 @@ __global__ void runCollision(CollisionFunctor collision, GPUCollisionParameter c
         para.forceZ = (collisionParameter.forces[2] + collisionParameter.bodyForceZ[k_000]) * c1o2 * collisionParameter.forceFactor;
 
         // Reset body force. To be used when not using round-off correction.
-        collisionParameter.bodyForceX[k_000] = 0.0f;
-        collisionParameter.bodyForceX[k_000] = 0.0f;
-        collisionParameter.bodyForceX[k_000] = 0.0f;
+        collisionParameter.bodyForceX[k_000] = c0o1;
+        collisionParameter.bodyForceY[k_000] = c0o1;
+        collisionParameter.bodyForceZ[k_000] = c0o1;
 
         ////////////////////////////////////////////////////////////////////////////////////
         //!> Round-off correction
@@ -97,7 +97,7 @@ __global__ void runCollision(CollisionFunctor collision, GPUCollisionParameter c
     vf::lbm::MacroscopicValues macroscopicValues;
     collision(para, macroscopicValues, turbulentViscosity);
 
-    if (writeMacroscopicVariables) {
+    if (writeMacroscopicVariables || turbulenceModel == vf::lbm::TurbulenceModel::AMD) {
         collisionParameter.vx[k_000] = macroscopicValues.vx;
         collisionParameter.vy[k_000] = macroscopicValues.vy;
         collisionParameter.vz[k_000] = macroscopicValues.vz;
diff --git a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h
index 850f56ead98c03b01d80dbc0c184f7eaf757af39..d337d7330c0fccf412d4ab2710c35ee6df39f509 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h
+++ b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h
@@ -41,6 +41,7 @@
 #include <basics/DataTypes.h>
 
 #include <lbm/interpolation/InterpolationCoefficients.h>
+#include <lbm/collision/TurbulentViscosity.h>
 
 using namespace vf::basics::constant;
 using namespace vf::lbm;
@@ -97,7 +98,7 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
     indices.k_0MM = k_base_0MM;
     indices.k_MMM = k_base_MMM;
 
-    omega_ = hasTurbulentViscosity ? calculateOmega(omega, turbulentViscosity[indices.k_000]) : omega;
+    omega_ = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omega, turbulentViscosity[indices.k_000]) : omega;
 
     real f_fine[27];
 
@@ -123,7 +124,7 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
     indices.k_0MM = neighborZ[indices.k_0MM];
     indices.k_MMM = neighborZ[indices.k_MMM];
 
-    omega_ = hasTurbulentViscosity ? calculateOmega(omega, turbulentViscosity[indices.k_000]) : omega;
+    omega_ = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omega, turbulentViscosity[indices.k_000]) : omega;
 
     read(f_fine, distFine, indices);
 
@@ -146,7 +147,7 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
     indices.k_0MM = indices.k_MMM;
     indices.k_MMM = neighborX[indices.k_MMM];
 
-    omega_ = hasTurbulentViscosity ? calculateOmega(omega, turbulentViscosity[indices.k_000]) : omega;
+    omega_ = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omega, turbulentViscosity[indices.k_000]) : omega;
 
     read(f_fine, distFine, indices);
     if (forcesForAllNodesX)
@@ -168,7 +169,7 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
     indices.k_0M0 = k_base_MM0;
     indices.k_MM0 = neighborX[k_base_MM0];
 
-    omega_ = hasTurbulentViscosity ? calculateOmega(omega, turbulentViscosity[indices.k_000]) : omega;
+    omega_ = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omega, turbulentViscosity[indices.k_000]) : omega;
 
     read(f_fine, distFine, indices);
     if (forcesForAllNodesX)
@@ -201,7 +202,7 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
     indices.k_0MM = k_base_0MM;
     indices.k_MMM = k_base_MMM;
 
-    omega_ = hasTurbulentViscosity ? calculateOmega(omega, turbulentViscosity[indices.k_000]) : omega;
+    omega_ = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omega, turbulentViscosity[indices.k_000]) : omega;
 
     read(f_fine, distFine, indices);
     if (forcesForAllNodesX)
@@ -223,7 +224,7 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
     indices.k_0MM = neighborZ[indices.k_0MM];
     indices.k_MMM = neighborZ[indices.k_MMM];
 
-    omega_ = hasTurbulentViscosity ? calculateOmega(omega, turbulentViscosity[indices.k_000]) : omega;
+    omega_ = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omega, turbulentViscosity[indices.k_000]) : omega;
 
     read(f_fine, distFine, indices);
     if (forcesForAllNodesX)
@@ -245,7 +246,7 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
     indices.k_0MM = indices.k_MMM;
     indices.k_MMM = neighborX[indices.k_MMM];
 
-    omega_ = hasTurbulentViscosity ? calculateOmega(omega, turbulentViscosity[indices.k_000]) : omega;
+    omega_ = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omega, turbulentViscosity[indices.k_000]) : omega;
 
     read(f_fine, distFine, indices);
     if (forcesForAllNodesX)
@@ -266,8 +267,8 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
     indices.k_M00 = neighborX[k_base_M00];
     indices.k_0M0 = k_base_MM0;
     indices.k_MM0 = neighborX[k_base_MM0];
-
-    omega_ = hasTurbulentViscosity ? calculateOmega(omega, turbulentViscosity[indices.k_000]) : omega;
+    
+    omega_ = hasTurbulentViscosity ? vf::lbm::calculateOmegaWithturbulentViscosity(omega, turbulentViscosity[indices.k_000]) : omega;
 
     read(f_fine, distFine, indices);
     if (forcesForAllNodesX)
diff --git a/src/lbm/collision/K17CompressibleNavierStokes.h b/src/lbm/collision/K17CompressibleNavierStokes.h
index af0d55e419212265cd4760eb6d82198f80e61499..33fd0afcef01f580427b2cc940ff06e13216bc6f 100644
--- a/src/lbm/collision/K17CompressibleNavierStokes.h
+++ b/src/lbm/collision/K17CompressibleNavierStokes.h
@@ -79,7 +79,7 @@ namespace vf::lbm
 //! ]</b></a> and \ref <a href="https://doi.org/10.1016/j.jcp.2017.07.004"><b>[ M. Geier et al. (2017),
 //! DOI:10.1016/j.jcp.2017.07.004 ]</b></a>
 ////////////////////////////////////////////////////////////////////////////////
-template <TurbulenceModel turbulenceModel, bool writeMacroscopicVariables>
+template <TurbulenceModel turbulenceModel>
 __host__ __device__ void runK17CompressibleNavierStokes(CollisionParameter& parameter, MacroscopicValues& macroscopicValues, TurbulentViscosity& turbulentViscosity)
 {
     auto& distribution = parameter.distribution;
@@ -237,10 +237,7 @@ __host__ __device__ void runK17CompressibleNavierStokes(CollisionParameter& para
     ////////////////////////////////////////////////////////////////////////////////////
     //! - Calculate modified omega with turbulent viscosity
     //!
-    real omega = parameter.omega;
-    if (turbulenceModel != TurbulenceModel::None) {
-        omega /= (c1o1 + c3o1 * parameter.omega * turbulentViscosity.value);
-    }
+    const real omega = turbulenceModel == TurbulenceModel::None ? parameter.omega : vf::lbm::calculateOmegaWithturbulentViscosity(parameter.omega, turbulentViscosity.value);
     ////////////////////////////////////////////////////////////
     // 2.
     real OxxPyyPzz = c1o1;
@@ -510,14 +507,6 @@ __host__ __device__ void runK17CompressibleNavierStokes(CollisionParameter& para
     m_010 = -m_010;
     m_001 = -m_001;
 
-    // Write to array here to distribute read/write
-    if (writeMacroscopicVariables || turbulenceModel == TurbulenceModel::AMD) {
-        macroscopicValues.rho = drho;
-        macroscopicValues.vx = vvx;
-        macroscopicValues.vy = vvy;
-        macroscopicValues.vz = vvz;
-    }
-
     ////////////////////////////////////////////////////////////////////////////////////
     //! - Chimera transform from central moments to well conditioned distributions as defined in Appendix J in
     //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015),
@@ -588,6 +577,11 @@ __host__ __device__ void runK17CompressibleNavierStokes(CollisionParameter& para
     distribution[DIR_MMP] = f_MMP;
     distribution[DIR_MPM] = f_MPM;
     distribution[DIR_MMM] = f_MMM;
+
+    macroscopicValues.rho = drho;
+    macroscopicValues.vx = vvx;
+    macroscopicValues.vy = vvy;
+    macroscopicValues.vz = vvz;
 }
 
 } // namespace vf::lbm
diff --git a/src/lbm/collision/TurbulentViscosity.h b/src/lbm/collision/TurbulentViscosity.h
index 0122dda7710f0832511afb83103c97ebcda9a384..7284638f1c21f894ba84cc121eb7ebd48a64704c 100644
--- a/src/lbm/collision/TurbulentViscosity.h
+++ b/src/lbm/collision/TurbulentViscosity.h
@@ -71,7 +71,7 @@ inline __host__ __device__ real calcTurbulentViscositySmagorinsky(real Cs, real
 }
 
 template <typename T>
-__host__ __device__ int max( T a, T b )
+__host__ __device__ T max( T a, T b )
 {
     return ( a > b ) ? a : b;
 }
@@ -87,6 +87,11 @@ inline __host__ __device__ real calcTurbulentViscosityQR(real C, real dxux, real
     return C * max(R, c0o1) / Q;
 }
 
+inline __host__ __device__ real calculateOmegaWithturbulentViscosity(real omega, real turbulenceViscosity)
+{
+    return omega / (c1o1 + c3o1 * omega * turbulenceViscosity);
+}
+
 } // namespace vf::lbm
 
 #endif //TURBULENT_VISCOSITY_H
diff --git a/utilities/ci-regression-tests/regression-tests-ci.yml.j2 b/utilities/ci-regression-tests/regression-tests-ci.yml.j2
index f162ff0c855130ffb8f2cc1c391a62360577b82a..98ddaccd48204f652d4768fb8cff3fa31eb2c7de 100644
--- a/utilities/ci-regression-tests/regression-tests-ci.yml.j2
+++ b/utilities/ci-regression-tests/regression-tests-ci.yml.j2
@@ -16,6 +16,11 @@ stages:
     - chmod +x ./regression-tests/*
     - pip install fieldcompare
 
+  artifacts:
+    expire_in: 1 hrs
+    paths:
+      - output/
+
 {% for regression_test in regression_tests %}
 run-regression-test-{{ regression_test }}:
   extends: .regression-test