From 7ce6c90b69d61bdab49f5c7f0de04427b969c4c4 Mon Sep 17 00:00:00 2001
From: Henry <henry.korb@geo.uu.se>
Date: Tue, 25 Jan 2022 12:15:30 +0100
Subject: [PATCH] performance and consistency update in actuatorLine

---
 .../LBM/ActuatorLine/configActuatorLine.txt   |   4 +-
 .../GPU/CudaMemoryManager.cpp                 |  94 ++++-----
 .../PreCollisionInteractor/ActuatorLine.cu    | 190 ++++++++----------
 .../PreCollisionInteractor/ActuatorLine.h     |  61 +++---
 src/lbm/constants/NumericConstants.h          |   2 +
 5 files changed, 165 insertions(+), 186 deletions(-)

diff --git a/apps/gpu/LBM/ActuatorLine/configActuatorLine.txt b/apps/gpu/LBM/ActuatorLine/configActuatorLine.txt
index 3b590b29f..233994f0d 100644
--- a/apps/gpu/LBM/ActuatorLine/configActuatorLine.txt
+++ b/apps/gpu/LBM/ActuatorLine/configActuatorLine.txt
@@ -1,8 +1,8 @@
 ##################################################
 #informations for Writing
 ##################################################
-Path = "."
+Path = .
 ##################################################
 #informations for reading
 ##################################################
-GridPath="."
+GridPath=.
diff --git a/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp b/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp
index 7074a222b..83fc16cf6 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp
+++ b/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp
@@ -2677,29 +2677,29 @@ void CudaMemoryManager::cudaFreeBladeRadii(ActuatorLine* actuatorLine)
 
 void CudaMemoryManager::cudaAllocBladeCoords(ActuatorLine* actuatorLine)
 {
-    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeCoordsXH, sizeof(real)*actuatorLine->getNumberOfNodes()) );
-    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeCoordsYH, sizeof(real)*actuatorLine->getNumberOfNodes()) );
-    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeCoordsZH, sizeof(real)*actuatorLine->getNumberOfNodes()) );
+    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeCoordsXH, sizeof(real)*actuatorLine->getNNodes()) );
+    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeCoordsYH, sizeof(real)*actuatorLine->getNNodes()) );
+    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeCoordsZH, sizeof(real)*actuatorLine->getNNodes()) );
 
-    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeCoordsXD, sizeof(real)*actuatorLine->getNumberOfNodes()) );
-    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeCoordsYD, sizeof(real)*actuatorLine->getNumberOfNodes()) );
-    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeCoordsZD, sizeof(real)*actuatorLine->getNumberOfNodes()) );
+    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeCoordsXD, sizeof(real)*actuatorLine->getNNodes()) );
+    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeCoordsYD, sizeof(real)*actuatorLine->getNNodes()) );
+    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeCoordsZD, sizeof(real)*actuatorLine->getNNodes()) );
 
-    setMemsizeGPU(3.f*actuatorLine->getNumberOfNodes(), false);
+    setMemsizeGPU(3.f*actuatorLine->getNNodes(), false);
 }
 
 void CudaMemoryManager::cudaCopyBladeCoordsHtoD(ActuatorLine* actuatorLine)
 {
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeCoordsXD, actuatorLine->bladeCoordsXH, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyHostToDevice) );
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeCoordsYD, actuatorLine->bladeCoordsYH, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyHostToDevice) );
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeCoordsZD, actuatorLine->bladeCoordsZH, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyHostToDevice) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeCoordsXD, actuatorLine->bladeCoordsXH, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyHostToDevice) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeCoordsYD, actuatorLine->bladeCoordsYH, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyHostToDevice) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeCoordsZD, actuatorLine->bladeCoordsZH, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyHostToDevice) );
 }
 
 void CudaMemoryManager::cudaCopyBladeCoordsDtoH(ActuatorLine* actuatorLine)
 {
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeCoordsXH, actuatorLine->bladeCoordsXD, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeCoordsYH, actuatorLine->bladeCoordsYD, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeCoordsZH, actuatorLine->bladeCoordsZD, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeCoordsXH, actuatorLine->bladeCoordsXD, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyDeviceToHost) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeCoordsYH, actuatorLine->bladeCoordsYD, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyDeviceToHost) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeCoordsZH, actuatorLine->bladeCoordsZD, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyDeviceToHost) );
 }
 
 void CudaMemoryManager::cudaFreeBladeCoords(ActuatorLine* actuatorLine)
@@ -2715,16 +2715,16 @@ void CudaMemoryManager::cudaFreeBladeCoords(ActuatorLine* actuatorLine)
 
 void CudaMemoryManager::cudaAllocBladeIndices(ActuatorLine* actuatorLine)
 {
-    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeIndicesH, sizeof(uint)*actuatorLine->getNumberOfNodes()) );
+    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeIndicesH, sizeof(uint)*actuatorLine->getNNodes()) );
 
-    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeIndicesD, sizeof(uint)*actuatorLine->getNumberOfNodes()) );
+    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeIndicesD, sizeof(uint)*actuatorLine->getNNodes()) );
 
-    setMemsizeGPU(sizeof(uint)*actuatorLine->getNumberOfNodes(), false);
+    setMemsizeGPU(sizeof(uint)*actuatorLine->getNNodes(), false);
 }
 
 void CudaMemoryManager::cudaCopyBladeIndicesHtoD(ActuatorLine* actuatorLine)
 {
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeIndicesD, actuatorLine->bladeIndicesH, sizeof(uint)*actuatorLine->getNumberOfNodes(), cudaMemcpyHostToDevice) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeIndicesD, actuatorLine->bladeIndicesH, sizeof(uint)*actuatorLine->getNNodes(), cudaMemcpyHostToDevice) );
 }
 
 void CudaMemoryManager::cudaFreeBladeIndices(ActuatorLine* actuatorLine)
@@ -2736,29 +2736,29 @@ void CudaMemoryManager::cudaFreeBladeIndices(ActuatorLine* actuatorLine)
 
 void CudaMemoryManager::cudaAllocBladeVelocities(ActuatorLine* actuatorLine)
 {
-    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeVelocitiesXH, sizeof(real)*actuatorLine->getNumberOfNodes()) );
-    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeVelocitiesYH, sizeof(real)*actuatorLine->getNumberOfNodes()) );
-    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeVelocitiesZH, sizeof(real)*actuatorLine->getNumberOfNodes()) );
+    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeVelocitiesXH, sizeof(real)*actuatorLine->getNNodes()) );
+    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeVelocitiesYH, sizeof(real)*actuatorLine->getNNodes()) );
+    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeVelocitiesZH, sizeof(real)*actuatorLine->getNNodes()) );
 
-    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeVelocitiesXD, sizeof(real)*actuatorLine->getNumberOfNodes()) );
-    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeVelocitiesYD, sizeof(real)*actuatorLine->getNumberOfNodes()) );
-    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeVelocitiesZD, sizeof(real)*actuatorLine->getNumberOfNodes()) );
+    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeVelocitiesXD, sizeof(real)*actuatorLine->getNNodes()) );
+    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeVelocitiesYD, sizeof(real)*actuatorLine->getNNodes()) );
+    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeVelocitiesZD, sizeof(real)*actuatorLine->getNNodes()) );
 
-    setMemsizeGPU(3.*sizeof(real)*actuatorLine->getNumberOfNodes(), false);
+    setMemsizeGPU(3.*sizeof(real)*actuatorLine->getNNodes(), false);
 }
 
 void CudaMemoryManager::cudaCopyBladeVelocitiesHtoD(ActuatorLine* actuatorLine)
 {
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeVelocitiesXD, actuatorLine->bladeVelocitiesXH, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyHostToDevice) );
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeVelocitiesYD, actuatorLine->bladeVelocitiesYH, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyHostToDevice) );
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeVelocitiesZD, actuatorLine->bladeVelocitiesZH, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyHostToDevice) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeVelocitiesXD, actuatorLine->bladeVelocitiesXH, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyHostToDevice) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeVelocitiesYD, actuatorLine->bladeVelocitiesYH, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyHostToDevice) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeVelocitiesZD, actuatorLine->bladeVelocitiesZH, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyHostToDevice) );
 }
 
 void CudaMemoryManager::cudaCopyBladeVelocitiesDtoH(ActuatorLine* actuatorLine)
 {
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeVelocitiesXH, actuatorLine->bladeVelocitiesXD, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeVelocitiesYH, actuatorLine->bladeVelocitiesYD, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeVelocitiesZH, actuatorLine->bladeVelocitiesZD, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeVelocitiesXH, actuatorLine->bladeVelocitiesXD, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyDeviceToHost) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeVelocitiesYH, actuatorLine->bladeVelocitiesYD, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyDeviceToHost) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeVelocitiesZH, actuatorLine->bladeVelocitiesZD, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyDeviceToHost) );
 }
 
 void CudaMemoryManager::cudaFreeBladeVelocities(ActuatorLine* actuatorLine)
@@ -2774,29 +2774,29 @@ void CudaMemoryManager::cudaFreeBladeVelocities(ActuatorLine* actuatorLine)
 
 void CudaMemoryManager::cudaAllocBladeForces(ActuatorLine* actuatorLine)
 {   
-    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeForcesXH, sizeof(real)*actuatorLine->getNumberOfNodes()) );
-    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeForcesYH, sizeof(real)*actuatorLine->getNumberOfNodes()) );
-    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeForcesZH, sizeof(real)*actuatorLine->getNumberOfNodes()) );
+    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeForcesXH, sizeof(real)*actuatorLine->getNNodes()) );
+    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeForcesYH, sizeof(real)*actuatorLine->getNNodes()) );
+    checkCudaErrors( cudaMallocHost((void**) &actuatorLine->bladeForcesZH, sizeof(real)*actuatorLine->getNNodes()) );
 
-    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeForcesXD, sizeof(real)*actuatorLine->getNumberOfNodes()) );
-    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeForcesYD, sizeof(real)*actuatorLine->getNumberOfNodes()) );
-    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeForcesZD, sizeof(real)*actuatorLine->getNumberOfNodes()) );
+    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeForcesXD, sizeof(real)*actuatorLine->getNNodes()) );
+    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeForcesYD, sizeof(real)*actuatorLine->getNNodes()) );
+    checkCudaErrors( cudaMalloc((void**) &actuatorLine->bladeForcesZD, sizeof(real)*actuatorLine->getNNodes()) );
 
-    setMemsizeGPU(3.*sizeof(real)*actuatorLine->getNumberOfNodes(), false);
+    setMemsizeGPU(3.*sizeof(real)*actuatorLine->getNNodes(), false);
 }
 
 void CudaMemoryManager::cudaCopyBladeForcesHtoD(ActuatorLine* actuatorLine)
 {
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeForcesXD, actuatorLine->bladeForcesXH, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyHostToDevice) );
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeForcesYD, actuatorLine->bladeForcesYH, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyHostToDevice) );
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeForcesZD, actuatorLine->bladeForcesZH, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyHostToDevice) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeForcesXD, actuatorLine->bladeForcesXH, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyHostToDevice) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeForcesYD, actuatorLine->bladeForcesYH, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyHostToDevice) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeForcesZD, actuatorLine->bladeForcesZH, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyHostToDevice) );
 }
 
 void CudaMemoryManager::cudaCopyBladeForcesDtoH(ActuatorLine* actuatorLine)
 {
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeForcesXH, actuatorLine->bladeForcesXD, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeForcesYH, actuatorLine->bladeForcesYD, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
-    checkCudaErrors( cudaMemcpy(actuatorLine->bladeForcesZH, actuatorLine->bladeForcesZD, sizeof(real)*actuatorLine->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeForcesXH, actuatorLine->bladeForcesXD, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyDeviceToHost) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeForcesYH, actuatorLine->bladeForcesYD, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyDeviceToHost) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->bladeForcesZH, actuatorLine->bladeForcesZD, sizeof(real)*actuatorLine->getNNodes(), cudaMemcpyDeviceToHost) );
 }
 
 void CudaMemoryManager::cudaFreeBladeForces(ActuatorLine* actuatorLine)
@@ -2812,14 +2812,14 @@ void CudaMemoryManager::cudaFreeBladeForces(ActuatorLine* actuatorLine)
 
 void CudaMemoryManager::cudaAllocSphereIndices(ActuatorLine* actuatorLine)
 {    
-    checkCudaErrors( cudaMallocHost((void**) &(actuatorLine->boundingSphereIndicesH), sizeof(int)*actuatorLine->getNumberOfIndices()));
-    checkCudaErrors( cudaMalloc((void**) &(actuatorLine->boundingSphereIndicesD), sizeof(int)*actuatorLine->getNumberOfIndices()));
-    setMemsizeGPU(sizeof(int)*actuatorLine->getNumberOfIndices(), false);
+    checkCudaErrors( cudaMallocHost((void**) &(actuatorLine->boundingSphereIndicesH), sizeof(int)*actuatorLine->getNIndices()));
+    checkCudaErrors( cudaMalloc((void**) &(actuatorLine->boundingSphereIndicesD), sizeof(int)*actuatorLine->getNIndices()));
+    setMemsizeGPU(sizeof(int)*actuatorLine->getNIndices(), false);
 }
 
 void CudaMemoryManager::cudaCopySphereIndicesHtoD(ActuatorLine* actuatorLine)
 {
-    checkCudaErrors( cudaMemcpy(actuatorLine->boundingSphereIndicesD, actuatorLine->boundingSphereIndicesH, sizeof(int)*actuatorLine->getNumberOfIndices(), cudaMemcpyHostToDevice) );
+    checkCudaErrors( cudaMemcpy(actuatorLine->boundingSphereIndicesD, actuatorLine->boundingSphereIndicesH, sizeof(int)*actuatorLine->getNIndices(), cudaMemcpyHostToDevice) );
 }
 
 void CudaMemoryManager::cudaFreeSphereIndices(ActuatorLine* actuatorLine)
diff --git a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.cu b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.cu
index 2d1ee506c..4f88c8895 100644
--- a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.cu
+++ b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.cu
@@ -5,7 +5,6 @@
 #include <helper_cuda.h>
 
 #include <cuda/CudaGrid.h>
-#include "lbm/constants/NumericConstants.h"
 #include "VirtualFluids_GPU/GPU/GeometryUtils.h"
 
 #include "Parameter/Parameter.h"
@@ -20,16 +19,12 @@ __host__ __device__ __inline__ uint calcNode(uint bladeNode, uint nBladeNodes, u
 __host__ __device__ __inline__ void calcBladeAndBladeNode(uint node, uint& bladeNode, uint nBladeNodes, uint& blade, uint nBlades)
 {
     blade = node/nBladeNodes;
-    bladeNode = node % nBladeNodes;
+    bladeNode = node - blade*nBladeNodes;
 }
 
-__host__ __device__ __inline__ real calcGaussian3D(real posX, real posY, real posZ, real destX, real destY, real destZ, real epsilon)
+__host__ __device__ __forceinline__ real distSqrd(real distX, real distY, real distZ)
 {
-    real distX = destX-posX;
-    real distY = destY-posY;
-    real distZ = destZ-posZ;
-    real dist = sqrt(distX*distX+distY*distY+distZ*distZ);
-    return pow(epsilon,-3)*pow(vf::lbm::constant::cPi,-1.5f)*exp(-pow(dist/epsilon,2));
+    return distX*distX+distY*distY+distZ*distZ;
 }
 
 __host__ __device__ __inline__ void rotateFromBladeToGlobal(
@@ -41,7 +36,6 @@ __host__ __device__ __inline__ void rotateFromBladeToGlobal(
 
     rotateAboutX3D(azimuth, bladeCoordX_BF, bladeCoordY_BF, bladeCoordZ_BF, tmpX, tmpY, tmpZ);
     rotateAboutZ3D(yaw, tmpX, tmpY, tmpZ, bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF);
-
 }
 
 __host__ __device__ __inline__ void rotateFromGlobalToBlade(
@@ -62,8 +56,8 @@ __global__ void interpolateVelocities(real* gridCoordsX, real* gridCoordsY, real
                                       real* bladeVelocitiesX, real* bladeVelocitiesY, real* bladeVelocitiesZ, 
                                       uint nBlades, uint nBladeNodes, 
                                       real azimuth, real yaw, real omega, 
-                                      real turbPosX, real turbPosZ, real turbPosY,
-                                      uint* bladeIndices, real velocityRatio)
+                                      real turbPosX, real turbPosY, real turbPosZ,
+                                      uint* bladeIndices, real velocityRatio, real invDeltaX)
 {
     const uint x = threadIdx.x; 
     const uint y = blockIdx.x;
@@ -80,21 +74,21 @@ __global__ void interpolateVelocities(real* gridCoordsX, real* gridCoordsY, real
 
     if(node>=nBladeNodes*nBlades) return;
 
-    real bladePosX_BF = bladeCoordsX[node];
-    real bladePosY_BF = bladeCoordsY[node];
-    real bladePosZ_BF = bladeCoordsZ[node];
+    real bladeCoordX_BF = bladeCoordsX[node];
+    real bladeCoordY_BF = bladeCoordsY[node];
+    real bladeCoordZ_BF = bladeCoordsZ[node];
 
-    real bladePosX_GF, bladePosY_GF, bladePosZ_GF;
+    real bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF;
 
-    real localAzimuth = azimuth+blade*2*vf::lbm::constant::cPi/nBlades;
+    real localAzimuth = azimuth+blade*c2Pi/nBlades;
 
-    rotateFromBladeToGlobal(bladePosX_BF, bladePosY_BF, bladePosZ_BF, 
-                            bladePosX_GF, bladePosY_GF, bladePosZ_GF,
+    rotateFromBladeToGlobal(bladeCoordX_BF, bladeCoordY_BF, bladeCoordZ_BF, 
+                            bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF,
                             localAzimuth, yaw);
 
-    bladePosX_GF += turbPosX;
-    bladePosY_GF += turbPosY;
-    bladePosZ_GF += turbPosZ;
+    bladeCoordX_GF += turbPosX;
+    bladeCoordY_GF += turbPosY;
+    bladeCoordZ_GF += turbPosZ;
 
     uint old_index = bladeIndices[node];
 
@@ -103,7 +97,7 @@ __global__ void interpolateVelocities(real* gridCoordsX, real* gridCoordsY, real
 
     k = findNearestCellBSW(old_index, 
                            gridCoordsX, gridCoordsY, gridCoordsZ, 
-                           bladePosX_GF, bladePosY_GF, bladePosZ_GF, 
+                           bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF, 
                            neighborsX, neighborsY, neighborsZ, neighborsWSB);
         
     bladeIndices[node] = k;
@@ -112,10 +106,9 @@ __global__ void interpolateVelocities(real* gridCoordsX, real* gridCoordsY, real
 
     real dW, dE, dN, dS, dT, dB;
 
-    real invDeltaX = 1.f/(gridCoordsX[ktne]-gridCoordsX[k]);
-    real distX = invDeltaX*(bladePosX_GF-gridCoordsX[k]);
-    real distY = invDeltaX*(bladePosY_GF-gridCoordsY[k]);
-    real distZ = invDeltaX*(bladePosZ_GF-gridCoordsZ[k]);
+    real distX = invDeltaX*(bladeCoordX_GF-gridCoordsX[k]);
+    real distY = invDeltaX*(bladeCoordY_GF-gridCoordsY[k]);
+    real distZ = invDeltaX*(bladeCoordZ_GF-gridCoordsZ[k]);
 
     getInterpolationWeights(dW, dE, dN, dS, dT, dB, distX, distY, distZ);
 
@@ -128,7 +121,7 @@ __global__ void interpolateVelocities(real* gridCoordsX, real* gridCoordsY, real
     rotateFromGlobalToBlade(bladeVelX_BF, bladeVelY_BF, bladeVelZ_BF, bladeVelX_GF, bladeVelY_GF, bladeVelZ_GF, localAzimuth, yaw);
 
     bladeVelocitiesX[node] = bladeVelX_BF;
-    bladeVelocitiesY[node] = bladeVelY_BF+omega*bladePosZ_BF;
+    bladeVelocitiesY[node] = bladeVelY_BF+omega*bladeCoordZ_BF;
     bladeVelocitiesZ[node] = bladeVelZ_BF;
 }
 
@@ -139,10 +132,9 @@ __global__ void applyBodyForces(real* gridCoordsX, real* gridCoordsY, real* grid
                                 real* bladeForcesX, real* bladeForcesY,real* bladeForcesZ,
                                 uint nBlades, uint nBladeNodes,
                                 real azimuth, real yaw, real omega, 
-                                real turbPosX, real turbPosZ, real turbPosY,
-                                real* bladeRadii, real forceRatio,
-                                uint* gridIndices, uint numberOfIndices, 
-                                real radius, real epsilon, real delta_x)
+                                real turbPosX, real turbPosY, real turbPosZ,
+                                uint* gridIndices, uint nIndices, 
+                                real invEpsilonSqrd, real factorGaussian)
 {
     const uint x = threadIdx.x; 
     const uint y = blockIdx.x;
@@ -153,63 +145,49 @@ __global__ void applyBodyForces(real* gridCoordsX, real* gridCoordsY, real* grid
 
     const uint index = nx*(ny*z + y) + x;
 
-    if(index>=numberOfIndices) return;
+    if(index>=nIndices) return;
 
     int gridIndex = gridIndices[index];
 
-    real posX = gridCoordsX[gridIndex];
-    real posY = gridCoordsY[gridIndex];
-    real posZ = gridCoordsZ[gridIndex];
-
-    real gridForceX = 0.0f;
-    real gridForceY = 0.0f;
-    real gridForceZ = 0.0f;
+    real gridCoordX_RF = gridCoordsX[gridIndex] - turbPosX;
+    real gridCoordY_RF = gridCoordsY[gridIndex] - turbPosY;
+    real gridCoordZ_RF = gridCoordsZ[gridIndex] - turbPosZ;
 
-    real deltaXCubed = pow(delta_x,3);
+    real gridForceX_RF = c0o1;
+    real gridForceY_RF = c0o1;
+    real gridForceZ_RF = c0o1;
 
-    real dAzimuth = 2*vf::lbm::constant::cPi/nBlades;
+    real dAzimuth = c2Pi/nBlades;
 
     for( uint blade=0; blade<nBlades; blade++)
     { 
-        real last_r = 0.0f;
-        real r = 0.0f;
         real localAzimuth = azimuth+blade*dAzimuth;
         
         for( uint bladeNode=0; bladeNode<nBladeNodes; bladeNode++)
         {
-            r = bladeRadii[bladeNode];
-
             uint node = calcNode(bladeNode, nBladeNodes, blade, nBlades);
 
-            real bladePosX_GF, bladePosY_GF, bladePosZ_GF;
+            real bladeCoordX_RF, bladeCoordY_RF, bladeCoordZ_RF;
 
             rotateFromBladeToGlobal(bladeCoordsX[node], bladeCoordsY[node], bladeCoordsZ[node], 
-                                    bladePosX_GF, bladePosY_GF, bladePosZ_GF,
+                                    bladeCoordX_RF, bladeCoordY_RF, bladeCoordZ_RF,
                                     localAzimuth, yaw);
 
-            bladePosX_GF += turbPosX;
-            bladePosY_GF += turbPosY;
-            bladePosZ_GF += turbPosZ;
+            real eta = factorGaussian*exp(-distSqrd(bladeCoordX_RF-gridCoordX_RF, bladeCoordY_RF-gridCoordY_RF, bladeCoordZ_RF-gridCoordZ_RF)*invEpsilonSqrd);
 
-            real eta = (r-last_r)*calcGaussian3D(posX, posY, posZ, bladePosX_GF, bladePosY_GF, bladePosZ_GF, epsilon)*deltaXCubed;
+            real forceX_RF, forceY_RF, forceZ_RF;
 
-            real forceX_GF, forceY_GF, forceZ_GF;
-
-            rotateFromBladeToGlobal(bladeForcesX[node], bladeForcesY[node], bladeForcesZ[node], forceX_GF, forceY_GF, forceZ_GF, localAzimuth, yaw);
+            rotateFromBladeToGlobal(bladeForcesX[node], bladeForcesY[node], bladeForcesZ[node], forceX_RF, forceY_RF, forceZ_RF, localAzimuth, yaw);
             
-            gridForceX += forceX_GF*eta;
-            gridForceY += forceY_GF*eta;
-            gridForceZ += forceZ_GF*eta;
-
-            last_r = r;
-        }         
+            gridForceX_RF += forceX_RF*eta;
+            gridForceY_RF += forceY_RF*eta;
+            gridForceZ_RF += forceZ_RF*eta;
+        }
     }
 
-    real invForceRatio = 1.f/forceRatio;
-
-    gridForcesX[gridIndex] = gridForceX*invForceRatio;
-    gridForcesY[gridIndex] = gridForceY*invForceRatio;
-    gridForcesZ[gridIndex] = gridForceZ*invForceRatio;
+    atomicAdd(&gridForcesX[gridIndex],gridForceX_RF);
+    atomicAdd(&gridForcesY[gridIndex],gridForceY_RF);
+    atomicAdd(&gridForcesZ[gridIndex],gridForceZ_RF);
 }
 
 
@@ -231,8 +209,7 @@ void ActuatorLine::interact(Parameter* para, CudaMemoryManager* cudaManager, int
 
     cudaManager->cudaCopyBladeCoordsHtoD(this);
 
-    uint numberOfThreads = para->getParH(level)->numberofthreads;
-    vf::cuda::CudaGrid bladeGrid = vf::cuda::CudaGrid(numberOfThreads, this->numberOfNodes);
+    vf::cuda::CudaGrid bladeGrid = vf::cuda::CudaGrid(para->getParH(level)->numberofthreads, this->nNodes);
 
     interpolateVelocities<<< bladeGrid.grid, bladeGrid.threads >>>(
         para->getParD(this->level)->coordX_SP, para->getParD(this->level)->coordY_SP, para->getParD(this->level)->coordZ_SP,        
@@ -242,8 +219,8 @@ void ActuatorLine::interact(Parameter* para, CudaMemoryManager* cudaManager, int
         this->bladeVelocitiesXD, this->bladeVelocitiesYD, this->bladeVelocitiesZD,  
         this->nBlades, this->nBladeNodes,
         this->azimuth, this->yaw, this->omega, 
-        this->turbinePosX, this->turbinePosZ, this->turbinePosY,
-        this->bladeIndicesD, para->getVelocityRatio());
+        this->turbinePosX, this->turbinePosY, this->turbinePosZ,
+        this->bladeIndicesD, para->getVelocityRatio(), this->invDeltaX);
 
     cudaManager->cudaCopyBladeVelocitiesDtoH(this);
 
@@ -251,7 +228,7 @@ void ActuatorLine::interact(Parameter* para, CudaMemoryManager* cudaManager, int
 
     cudaManager->cudaCopyBladeForcesHtoD(this);
 
-    vf::cuda::CudaGrid sphereGrid = vf::cuda::CudaGrid(numberOfThreads, this->numberOfIndices);
+    vf::cuda::CudaGrid sphereGrid = vf::cuda::CudaGrid(para->getParH(level)->numberofthreads, this->nIndices);
 
     applyBodyForces<<<sphereGrid.grid, sphereGrid.threads>>>(
         para->getParD(this->level)->coordX_SP, para->getParD(this->level)->coordY_SP, para->getParD(this->level)->coordZ_SP,        
@@ -260,12 +237,11 @@ void ActuatorLine::interact(Parameter* para, CudaMemoryManager* cudaManager, int
         this->bladeForcesXD, this->bladeForcesYD, this->bladeForcesZD,
         this->nBlades, this->nBladeNodes,
         this->azimuth, this->yaw, this->omega, 
-        this->turbinePosX, this->turbinePosZ, this->turbinePosY,
-        this->bladeRadiiD, this->density*pow(para->getViscosityRatio(),2),
-        this->boundingSphereIndicesD, this->numberOfIndices,
-        this->diameter*0.5f, this->epsilon, this->delta_x);
+        this->turbinePosX, this->turbinePosY, this->turbinePosZ,
+        this->boundingSphereIndicesD, this->nIndices,
+        this->invEpsilonSqrd, this->factorGaussian);
 
-    this->azimuth = fmod(this->azimuth+this->omega*this->delta_t,2*vf::lbm::constant::cPi);
+    this->azimuth = fmod(this->azimuth+this->omega*this->deltaT,c2Pi);
 }
 
 
@@ -285,9 +261,9 @@ void ActuatorLine::calcForcesEllipticWing()
     uint node;
     real u_rel, v_rel, u_rel_sq;
     real phi;
-    real Cl = 1.f;
-    real Cd = 0.f;
-    real c0 = 1.f;
+    real Cl = c1o1;
+    real Cd = c0o1;
+    real c0 = c1o1;
 
     real c, Cn, Ct;
 
@@ -301,14 +277,15 @@ void ActuatorLine::calcForcesEllipticWing()
             v_rel = this->bladeVelocitiesYH[node];
             u_rel_sq = u_rel*u_rel+v_rel*v_rel;
             phi = atan2(u_rel, v_rel);
-
-            c = c0 * sqrt( 1.f- pow(4.f*this->bladeRadiiH[bladeNode]/this->diameter-1.f, 2.f) );
-            Cn =   Cl*cos(phi)+Cd*sin(phi);
-            Ct =  -Cl*sin(phi)+Cd*cos(phi);
+            
+            real tmp = c4o1*this->bladeRadiiH[bladeNode]/this->diameter-c1o1;
+            c = c0 * sqrt( c1o1- tmp*tmp );
+            Cn = Cl*cos(phi)+Cd*sin(phi);
+            Ct = Cl*sin(phi)-Cd*cos(phi);
         
-            this->bladeForcesXH[node] = -0.5f*u_rel_sq*c*this->density*Cn;
-            this->bladeForcesYH[node] = -0.5f*u_rel_sq*c*this->density*Ct;
-            this->bladeForcesZH[node] = 0.0f;
+            this->bladeForcesXH[node] = -c1o2*u_rel_sq*c*this->density*Cn;
+            this->bladeForcesYH[node] = -c1o2*u_rel_sq*c*this->density*Ct;
+            this->bladeForcesZH[node] = c0o1;
         }
     }
 }
@@ -322,13 +299,16 @@ void ActuatorLine::initBladeRadii(CudaMemoryManager* cudaManager)
 {   
     cudaManager->cudaAllocBladeRadii(this);
 
-    real dx = 0.5f*this->diameter/this->nBladeNodes;  
+    real dr = c1o2*this->diameter/this->nBladeNodes;  
 
     for(uint node=0; node<this->nBladeNodes; node++)
     {
-        this->bladeRadiiH[node] = dx*(node+1);
+        this->bladeRadiiH[node] = dr*(node+1);
     }
     cudaManager->cudaCopyBladeRadiiHtoD(this);
+
+    real dxOPiSqrtEps = this->deltaX/(this->epsilon*sqrt(cPi));
+    this->factorGaussian = dr*dxOPiSqrtEps*dxOPiSqrtEps*dxOPiSqrtEps/this->forceRatio;
 }
 
 void ActuatorLine::initBladeCoords(CudaMemoryManager* cudaManager)
@@ -341,8 +321,8 @@ void ActuatorLine::initBladeCoords(CudaMemoryManager* cudaManager)
         {
             uint node = calcNode(bladeNode, this->nBladeNodes, blade, this->nBlades);
 
-            this->bladeCoordsXH[node] = 0.f;
-            this->bladeCoordsYH[node] = 0.f;
+            this->bladeCoordsXH[node] = c0o1;
+            this->bladeCoordsYH[node] = c0o1;
             this->bladeCoordsZH[node] = this->bladeRadiiH[bladeNode];
         }
     }
@@ -353,11 +333,11 @@ void ActuatorLine::initBladeVelocities(CudaMemoryManager* cudaManager)
 {   
     cudaManager->cudaAllocBladeVelocities(this);
 
-    for(uint node=0; node<this->numberOfNodes; node++)
+    for(uint node=0; node<this->nNodes; node++)
     {
-        this->bladeVelocitiesXH[node] = 0.f;
-        this->bladeVelocitiesYH[node] = 0.f;
-        this->bladeVelocitiesZH[node] = 0.f;
+        this->bladeVelocitiesXH[node] = c0o1;
+        this->bladeVelocitiesYH[node] = c0o1;
+        this->bladeVelocitiesZH[node] = c0o1;
     }
     cudaManager->cudaCopyBladeVelocitiesHtoD(this);
 }
@@ -366,11 +346,11 @@ void ActuatorLine::initBladeForces(CudaMemoryManager* cudaManager)
 {   
     cudaManager->cudaAllocBladeForces(this);
 
-    for(uint node=0; node<this->numberOfNodes; node++)
+    for(uint node=0; node<this->nNodes; node++)
     {
-        this->bladeForcesXH[node] = 0.f;
-        this->bladeForcesYH[node] = 0.f;
-        this->bladeForcesZH[node] = 0.f;
+        this->bladeForcesXH[node] = c0o1;
+        this->bladeForcesYH[node] = c0o1;
+        this->bladeForcesZH[node] = c0o1;
     }
     cudaManager->cudaCopyBladeForcesHtoD(this);
 }
@@ -379,7 +359,7 @@ void ActuatorLine::initBladeIndices(Parameter* para, CudaMemoryManager* cudaMana
 {   
     cudaManager->cudaAllocBladeIndices(this);
 
-    for(uint node=0; node<this->numberOfNodes; node++)
+    for(uint node=0; node<this->nNodes; node++)
     {
         this->bladeIndicesH[node] = 1;
     }
@@ -390,18 +370,18 @@ void ActuatorLine::initBoundingSphere(Parameter* para, CudaMemoryManager* cudaMa
 {
     // Actuator line exists only on 1 level
     std::vector<int> nodesInSphere;
-    real sphereRadius = 0.5*this->diameter+4.f*this->epsilon;
+    real sphereRadius = c1o2*this->diameter+c4o1*this->epsilon;
+    real sphereRadiusSqrd = sphereRadius*sphereRadius;
 
     for (uint j = 1; j <= para->getParH(this->level)->size_Mat_SP; j++)
     {
         const real distX = para->getParH(this->level)->coordX_SP[j]-this->turbinePosX;
         const real distY = para->getParH(this->level)->coordY_SP[j]-this->turbinePosY;
         const real distZ = para->getParH(this->level)->coordZ_SP[j]-this->turbinePosZ;
-        const real dist = sqrt(pow(distX,2)+pow(distY,2)+pow(distZ,2));
-        if(dist < sphereRadius) nodesInSphere.push_back(j);
+        if(distSqrd(distX,distY,distZ) < sphereRadiusSqrd) nodesInSphere.push_back(j);
     }
 
-    this->numberOfIndices = uint(nodesInSphere.size());
+    this->nIndices = uint(nodesInSphere.size());
     cudaManager->cudaAllocSphereIndices(this);
     std::copy(nodesInSphere.begin(), nodesInSphere.end(), this->boundingSphereIndicesH);
     cudaManager->cudaCopySphereIndicesHtoD(this);
@@ -410,7 +390,7 @@ void ActuatorLine::initBoundingSphere(Parameter* para, CudaMemoryManager* cudaMa
 void ActuatorLine::setBladeCoords(real* _bladeCoordsX, real* _bladeCoordsY, real* _bladeCoordsZ)
 { 
 
-    for(uint node=0; node<this->numberOfNodes; node++)
+    for(uint node=0; node<this->nNodes; node++)
     {
         this->bladeCoordsXH[node] = _bladeCoordsX[node];
         this->bladeCoordsYH[node] = _bladeCoordsY[node];
@@ -420,7 +400,7 @@ void ActuatorLine::setBladeCoords(real* _bladeCoordsX, real* _bladeCoordsY, real
 
 void ActuatorLine::setBladeVelocities(real* _bladeVelocitiesX, real* _bladeVelocitiesY, real* _bladeVelocitiesZ)
 { 
-    for(uint node=0; node<this->numberOfNodes; node++)
+    for(uint node=0; node<this->nNodes; node++)
     {
         this->bladeVelocitiesXH[node] = _bladeVelocitiesX[node];
         this->bladeVelocitiesYH[node] = _bladeVelocitiesY[node];
@@ -430,7 +410,7 @@ void ActuatorLine::setBladeVelocities(real* _bladeVelocitiesX, real* _bladeVeloc
 
 void ActuatorLine::setBladeForces(real* _bladeForcesX, real* _bladeForcesY, real* _bladeForcesZ)
 { 
-    for(uint node=0; node<this->numberOfNodes; node++)
+    for(uint node=0; node<this->nNodes; node++)
     {
         this->bladeForcesXH[node] = _bladeForcesX[node];
         this->bladeForcesYH[node] = _bladeForcesY[node];
diff --git a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.h b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.h
index 6a61cb7bd..af08c6efc 100644
--- a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.h
+++ b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.h
@@ -4,10 +4,12 @@
 #include "PreCollisionInteractor.h"
 #include "PointerDefinitions.h"
 #include "VirtualFluids_GPU_export.h"
+#include "lbm/constants/NumericConstants.h"
 
 class Parameter;
 class GridProvider;
 
+using namespace vf::lbm::constant;
 class VIRTUALFLUIDS_GPU_EXPORT ActuatorLine : public PreCollisionInteractor
 {
 public:
@@ -19,8 +21,8 @@ public:
         real _turbinePosX, real _turbinePosY, real _turbinePosZ,
         const real _diameter,
         int _level,
-        const real _delta_t,
-        const real _delta_x
+        const real _deltaT,
+        const real _deltaX
     ) : nBlades(_nBlades),
         density(_density),
         nBladeNodes(_nBladeNodes), 
@@ -28,15 +30,17 @@ public:
         turbinePosX(_turbinePosX), turbinePosY(_turbinePosY), turbinePosZ(_turbinePosZ),
         diameter(_diameter),
         level(_level),
-        delta_x(_delta_x),
         PreCollisionInteractor()
     {
-        this->delta_t = _delta_t/pow(2,this->level);
-        this->delta_x = _delta_x/pow(2,this->level);
-        this->numberOfNodes = this->nBladeNodes*this->nBlades;
-        this->omega = 1.0f;
-        this->azimuth = 0.0f;
-        this->yaw = 0.0f;
+        this->deltaT = _deltaT/pow(2,this->level);
+        this->deltaX = _deltaX/pow(2,this->level);
+        this->invDeltaX = c1o1/this->deltaX;
+        this->forceRatio = this->density*pow(this->deltaX,4)*pow(this->deltaT,-2);
+        this->invEpsilonSqrd = c1o1/(this->epsilon*this->epsilon);
+        this->nNodes = this->nBladeNodes*this->nBlades;
+        this->omega = c1o1;
+        this->azimuth = c0o1;
+        this->yaw = c0o1;
     };
 
     virtual ~ActuatorLine(){};
@@ -46,26 +50,26 @@ public:
     void free(Parameter* para, CudaMemoryManager* cudaManager) override;
     void write(uint t);
 
-    uint getNBladeNodes(){return this->nBladeNodes;};
-    uint getNBlades(){return this->nBlades;};
-    uint getNumberOfIndices(){return this->numberOfIndices;};
-    uint getNumberOfNodes(){return this->numberOfNodes;};
+    uint getNBladeNodes(){ return this->nBladeNodes; };
+    uint getNBlades(){ return this->nBlades;};
+    uint getNIndices(){ return this->nIndices; };
+    uint getNNodes(){ return this->nNodes; };
     real getOmega(){ return this->omega; };
     real getAzimuth(){ return this->azimuth; };
     real getDensity(){ return this->density; };
     real getPositionX(){ return this->turbinePosX; };
     real getPositionY(){ return this->turbinePosY; };
     real getPositionZ(){ return this->turbinePosZ; };
-    real* getBladeRadii(){return this->bladeRadiiH;};
-    real* getBladeCoordsX(){return this->bladeCoordsXH;};
-    real* getBladeCoordsY(){return this->bladeCoordsYH;};
-    real* getBladeCoordsZ(){return this->bladeCoordsZH;};
-    real* getBladeVelocitiesX(){return this->bladeVelocitiesXH;};
-    real* getBladeVelocitiesY(){return this->bladeVelocitiesYH;};
-    real* getBladeVelocitiesZ(){return this->bladeVelocitiesZH;};
-    real* getBladeForcesX(){return this->bladeForcesXH;};
-    real* getBladeForcesY(){return this->bladeForcesYH;};
-    real* getBladeForcesZ(){return this->bladeForcesZH;};
+    real* getBladeRadii(){ return this->bladeRadiiH; };
+    real* getBladeCoordsX(){ return this->bladeCoordsXH; };
+    real* getBladeCoordsY(){ return this->bladeCoordsYH; };
+    real* getBladeCoordsZ(){ return this->bladeCoordsZH; };
+    real* getBladeVelocitiesX(){ return this->bladeVelocitiesXH; };
+    real* getBladeVelocitiesY(){ return this->bladeVelocitiesYH; };
+    real* getBladeVelocitiesZ(){ return this->bladeVelocitiesZH; };
+    real* getBladeForcesX(){ return this->bladeForcesXH; };
+    real* getBladeForcesY(){ return this->bladeForcesYH; };
+    real* getBladeForcesZ(){ return this->bladeForcesZH; };
 
     void setOmega(real _omega){ this->omega = _omega; };
     void setAzimuth(real _azimuth){ this->azimuth = _azimuth; };
@@ -85,12 +89,6 @@ private:
     void initBladeIndices(Parameter* para, CudaMemoryManager* cudaManager);
 
     void calcForcesEllipticWing();
-    void rotateBlades(real angle);
-
-    void writeBladeCoords(uint t){};
-    void writeBladeForces(uint t){};
-    void writeBladeVelocities(uint t){};
-    
 
 public:
     real* bladeRadiiH;
@@ -109,14 +107,13 @@ public:
 private:
     const real density;
     real turbinePosX, turbinePosY, turbinePosZ;
-    real omega, azimuth, yaw, delta_t, delta_x;
+    real omega, azimuth, yaw, deltaT, deltaX, invDeltaX, forceRatio, factorGaussian, invEpsilonSqrd;
     const real diameter;
     const uint nBladeNodes;
     const uint nBlades;
     const real epsilon; // in m
     const int level;
-    uint numberOfIndices;
-    uint numberOfNodes;
+    uint nIndices, nNodes;
 };
 
 #endif
\ No newline at end of file
diff --git a/src/lbm/constants/NumericConstants.h b/src/lbm/constants/NumericConstants.h
index 5c7556ceb..5f023276e 100644
--- a/src/lbm/constants/NumericConstants.h
+++ b/src/lbm/constants/NumericConstants.h
@@ -119,6 +119,7 @@ static constexpr double c10eM10 = 1e-10;
 static constexpr double smallSingle = 0.0000000002;
 
 static constexpr double cPi = 3.1415926535;
+static constexpr double c2Pi = 6.28318530717;
 static constexpr double cPio180 = 1.74532925199e-2;
 static constexpr double c180oPi = 57.2957795131;
 
@@ -232,6 +233,7 @@ static constexpr float c10eM10 = 1e-10f;
 static constexpr float smallSingle = 0.0000000002f;
 
 static constexpr float cPi = 3.1415926535f;
+static constexpr double c2Pi = 6.2831853071f;
 static constexpr float cPio180 = 1.74532925199e-2f;
 static constexpr float c180oPi = 57.2957795131f;
 
-- 
GitLab