From fd2171cbdc73885f147f36bc5d42629fdeeda753 Mon Sep 17 00:00:00 2001
From: Henry <henry.korb@geo.uu.se>
Date: Thu, 12 Oct 2023 15:49:50 +0200
Subject: [PATCH] move calculateOmega to lbm library

---
 .../GPU/GridScaling/scaleCF_compressible.cu   | 33 ++++++++++---------
 .../GPU/GridScaling/scaleFC_compressible.cu   |  2 +-
 .../LBM/GPUHelperFunctions/ScalingUtilities.h | 17 +++++-----
 3 files changed, 27 insertions(+), 25 deletions(-)

diff --git a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleCF_compressible.cu b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleCF_compressible.cu
index 4df1d28b1..05ec84179 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 2492cb2ac..c1202e5a7 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/LBM/GPUHelperFunctions/ScalingUtilities.h b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h
index 2d71eaf42..e24aed949 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::dir;
@@ -94,7 +95,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];
 
@@ -114,7 +115,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);
     momentsSet.calculateMMP(f_fine, omega_);
@@ -132,7 +133,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);
     momentsSet.calculatePMP(f_fine, omega_);
@@ -150,7 +151,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);
     momentsSet.calculatePMM(f_fine, omega_);
@@ -178,7 +179,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);
     momentsSet.calculateMPM(f_fine, omega_);
@@ -196,7 +197,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);
     momentsSet.calculateMPP(f_fine, omega_);
@@ -214,7 +215,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);
     momentsSet.calculatePPP(f_fine, omega_);
@@ -232,7 +233,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);
     momentsSet.calculatePPM(f_fine, omega_);
-- 
GitLab