diff --git a/src/gpu/VirtualFluids_GPU/CMakeLists.txt b/src/gpu/VirtualFluids_GPU/CMakeLists.txt
index d7e4a776271c5bc243a4ed663ffecd23dec5cebb..8228d25394c5bc4b3a55d7f978438d2ad187e285 100644
--- a/src/gpu/VirtualFluids_GPU/CMakeLists.txt
+++ b/src/gpu/VirtualFluids_GPU/CMakeLists.txt
@@ -26,4 +26,5 @@ if(BUILD_VF_UNIT_TESTS)
set_source_files_properties(Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17Test.cpp PROPERTIES LANGUAGE CUDA)
set_source_files_properties(Kernel/Utilities/DistributionHelperTests.cpp PROPERTIES LANGUAGE CUDA)
set_source_files_properties(Parameter/ParameterTest.cpp PROPERTIES LANGUAGE CUDA)
+ set_source_files_properties(PreCollisionInteractor/ActuatorFarmInlinesTest.cpp PROPERTIES LANGUAGE CUDA)
endif()
diff --git a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.cu b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.cu
index 22f5c4791f4207ceda7c78d6b437ffffdd21ad80..2a8e96c95682aa4ee8f7e47ae2a09ca75f752600 100644
--- a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.cu
+++ b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.cu
@@ -31,6 +31,7 @@
//! \author Henrik Asmuth, Henry Korb
//======================================================================================
#include "ActuatorFarm.h"
+#include "ActuatorFarmInlines.h"
#include <cuda.h>
#include <cuda_runtime.h>
@@ -50,22 +51,6 @@
using namespace vf::basics::constant;
-
-__host__ __device__ __inline__ uint calcNode(uint bladeNode, uint numberOfBladeNodes, uint blade, uint numberOfBlades, uint turbine, uint numberOfTurbines)
-{
-
- return bladeNode+numberOfBladeNodes*(blade+numberOfBlades*turbine);
-}
-
-__host__ __device__ __inline__ void calcTurbineBladeAndBladeNode(uint node, uint& bladeNode, uint numberOfBladeNodes, uint& blade, uint numberOfBlades, uint& turbine, uint numberOfTurbines)
-{
- turbine = node/(numberOfBladeNodes*numberOfBlades);
- uint x_off = turbine*numberOfBladeNodes*numberOfBlades;
- blade = (node - x_off)/numberOfBlades;
- uint y_off = numberOfBladeNodes*blade+x_off;
- bladeNode = (node - y_off);
-}
-
__host__ __device__ __forceinline__ real distSqrd(real distX, real distY, real distZ)
{
return distX*distX+distY*distY+distZ*distZ;
@@ -101,6 +86,24 @@ __host__ __device__ __inline__ void rotateFromGlobalToBlade(
invRotateAboutX3D(azimuth, tmpX, tmpY, tmpZ, bladeCoordX_BF, bladeCoordY_BF, bladeCoordZ_BF);
}
+__host__ __device__ __inline__ void calcCoordinateOfBladeNodeInGridFrameOfReference(uint bladeNodeIndex, uint turbine, real localAzimuth, real yaw,
+ const real* bladeCoordsX, const real* bladeCoordsY, const real* bladeCoordsZ,
+ const real* turbPosX, const real* turbPosY, const real* turbPosZ,
+ real& bladeCoordX_GF, real& bladeCoordY_GF, real& bladeCoordZ_GF)
+{
+
+ const real bladeCoordX_BF = bladeCoordsX[bladeNodeIndex];
+ const real bladeCoordY_BF = bladeCoordsY[bladeNodeIndex];
+ const real bladeCoordZ_BF = bladeCoordsZ[bladeNodeIndex];
+
+ rotateFromBladeToGlobal(bladeCoordX_BF, bladeCoordY_BF, bladeCoordZ_BF, bladeCoordX_GF, bladeCoordY_GF,
+ bladeCoordZ_GF, localAzimuth, yaw);
+
+ bladeCoordX_GF += turbPosX[turbine];
+ bladeCoordY_GF += turbPosY[turbine];
+ bladeCoordZ_GF += turbPosZ[turbine];
+}
+
__global__ void interpolateVelocities(real* gridCoordsX, real* gridCoordsY, real* gridCoordsZ,
uint* neighborsX, uint* neighborsY, uint* neighborsZ, uint* neighborsWSB,
real* vx, real* vy, real* vz,
@@ -117,29 +120,19 @@ __global__ void interpolateVelocities(real* gridCoordsX, real* gridCoordsY, real
//!
const unsigned nodeIndex = vf::gpu::getNodeIndex();
- if(nodeIndex>=numberOfBladeNodes*numberOfBlades*numberOfTurbines) return;
-
- uint turbine, bladeNode, blade;
+ if (nodeIndex >= numberOfBladeNodes * numberOfBlades * numberOfTurbines) return;
- calcTurbineBladeAndBladeNode(nodeIndex, bladeNode, numberOfBladeNodes, blade, numberOfBlades, turbine, numberOfTurbines);
+ uint turbine, bladeNodeIndexOnBlade, blade;
+ calcTurbineBladeAndBladeNode(nodeIndex, bladeNodeIndexOnBlade, numberOfBladeNodes, blade, numberOfBlades, turbine, numberOfTurbines);
- real bladeCoordX_BF = bladeCoordsX[nodeIndex];
- real bladeCoordY_BF = bladeCoordsY[nodeIndex];
- real bladeCoordZ_BF = bladeCoordsZ[nodeIndex];
+ const real localAzimuth = azimuths[turbine] + blade * c2Pi / numberOfBlades;
+ const real yaw = yaws[turbine];
real bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF;
-
- real localAzimuth = azimuths[turbine]+blade*c2Pi/numberOfBlades;
- real yaw = yaws[turbine];
-
-
- rotateFromBladeToGlobal(bladeCoordX_BF, bladeCoordY_BF, bladeCoordZ_BF,
- bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF,
- localAzimuth, yaw);
-
- bladeCoordX_GF += turbPosX[turbine];
- bladeCoordY_GF += turbPosY[turbine];
- bladeCoordZ_GF += turbPosZ[turbine];
+ calcCoordinateOfBladeNodeInGridFrameOfReference(nodeIndex, turbine, localAzimuth, yaw,
+ bladeCoordsX, bladeCoordsY, bladeCoordsZ,
+ turbPosX, turbPosY, turbPosZ,
+ bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF);
uint k, ke, kn, kt;
uint kne, kte, ktn, ktne;
@@ -148,7 +141,7 @@ __global__ void interpolateVelocities(real* gridCoordsX, real* gridCoordsY, real
gridCoordsX, gridCoordsY, gridCoordsZ,
bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF,
neighborsX, neighborsY, neighborsZ, neighborsWSB);
-
+
bladeIndices[nodeIndex] = k;
getNeighborIndicesOfBSW(k, ke, kn, kt, kne, kte, ktn, ktne, neighborsX, neighborsY, neighborsZ);
@@ -172,7 +165,7 @@ __global__ void interpolateVelocities(real* gridCoordsX, real* gridCoordsY, real
localAzimuth, yaw);
bladeVelocitiesX[nodeIndex] = bladeVelX_BF;
- bladeVelocitiesY[nodeIndex] = bladeVelY_BF+omegas[turbine]*bladeCoordZ_BF;
+ bladeVelocitiesY[nodeIndex] = bladeVelY_BF + omegas[turbine] * bladeCoordsZ[nodeIndex];
bladeVelocitiesZ[nodeIndex] = bladeVelZ_BF;
}
@@ -192,7 +185,6 @@ __global__ void applyBodyForces(real* gridCoordsX, real* gridCoordsY, real* grid
if(index>=nIndices) return;
-
uint gridIndex = gridIndices[index];
real gridCoordX_GF = gridCoordsX[gridIndex];
@@ -230,7 +222,7 @@ __global__ void applyBodyForces(real* gridCoordsX, real* gridCoordsY, real* grid
localAzimuth, yaw);
uint node;
- uint nextNode = calcNode(0, numberOfBladeNodes, blade, numberOfBlades, turbine, numberOfTurbines);
+ uint nextNode = calcNodeIndexInBladeArrays(0, numberOfBladeNodes, blade, numberOfBlades, turbine, numberOfTurbines);
real last_z = c0o1;
real current_z = c0o1;
@@ -241,7 +233,7 @@ __global__ void applyBodyForces(real* gridCoordsX, real* gridCoordsY, real* grid
for( uint bladeNode=0; bladeNode<numberOfBladeNodes-1; bladeNode++)
{
node = nextNode;
- nextNode = calcNode(bladeNode+1, numberOfBladeNodes, blade, numberOfBlades, turbine, numberOfTurbines);
+ nextNode = calcNodeIndexInBladeArrays(bladeNode+1, numberOfBladeNodes, blade, numberOfBlades, turbine, numberOfTurbines);
x = bladeCoordsX[node];
y = bladeCoordsY[node];
@@ -255,7 +247,7 @@ __global__ void applyBodyForces(real* gridCoordsX, real* gridCoordsY, real* grid
rotateFromBladeToGlobal(bladeForcesX[node], bladeForcesY[node], bladeForcesZ[node],
forceX_RF, forceY_RF, forceZ_RF,
localAzimuth, yaw);
-
+
gridForceX_RF += forceX_RF*eta;
gridForceY_RF += forceY_RF*eta;
gridForceZ_RF += forceZ_RF*eta;
@@ -277,7 +269,7 @@ __global__ void applyBodyForces(real* gridCoordsX, real* gridCoordsY, real* grid
rotateFromBladeToGlobal(bladeForcesX[node], bladeForcesY[node], bladeForcesZ[node],
forceX_RF, forceY_RF, forceZ_RF,
localAzimuth, yaw);
-
+
gridForceX_RF += forceX_RF*eta;
gridForceY_RF += forceY_RF*eta;
gridForceZ_RF += forceZ_RF*eta;
@@ -320,7 +312,7 @@ void ActuatorFarm::interact(Parameter* para, CudaMemoryManager* cudaMemoryManage
cudaStream_t stream = para->getStreamManager()->getStream(CudaStreamIndex::ActuatorFarm, this->streamIndex);
- if(useHostArrays) cudaMemoryManager->cudaCopyBladeCoordsHtoD(this);
+ if (useHostArrays) cudaMemoryManager->cudaCopyBladeCoordsHtoD(this);
vf::cuda::CudaGrid bladeGrid = vf::cuda::CudaGrid(para->getParH(level)->numberofthreads, this->numberOfGridNodes);
@@ -389,7 +381,7 @@ void ActuatorFarm::calcForcesEllipticWing()
{
for( uint bladeNode=0; bladeNode<this->numberOfBladeNodes; bladeNode++)
{
- uint node = calcNode(bladeNode, this->numberOfBladeNodes, blade, this->numberOfBlades, turbine, this->numberOfTurbines);
+ uint node = calcNodeIndexInBladeArrays(bladeNode, this->numberOfBladeNodes, blade, this->numberOfBlades, turbine, this->numberOfTurbines);
u_rel = this->bladeVelocitiesXH[node];
v_rel = this->bladeVelocitiesYH[node];
@@ -436,7 +428,7 @@ void ActuatorFarm::initTurbineGeometries(CudaMemoryManager* cudaMemoryManager)
{
for(uint node=0; node<this->numberOfBladeNodes; node++)
{
- this->bladeRadiiH[calcNode(node, numberOfBladeNodes, 0, 1, turbine, numberOfTurbines)] = this->preInitBladeRadii[turbine][node];
+ this->bladeRadiiH[calcNodeIndexInBladeArrays(node, numberOfBladeNodes, 0, 1, turbine, numberOfTurbines)] = this->preInitBladeRadii[turbine][node];
}
}
@@ -464,11 +456,11 @@ void ActuatorFarm::initBladeCoords(CudaMemoryManager* cudaMemoryManager)
{
for(uint bladeNode=0; bladeNode<this->numberOfBladeNodes; bladeNode++)
{
- uint node = calcNode(bladeNode, this->numberOfBladeNodes, blade, this->numberOfBlades, turbine, this->numberOfTurbines);
+ uint node = calcNodeIndexInBladeArrays(bladeNode, this->numberOfBladeNodes, blade, this->numberOfBlades, turbine, this->numberOfTurbines);
this->bladeCoordsXH[node] = c0o1;
this->bladeCoordsYH[node] = c0o1;
- this->bladeCoordsZH[node] = this->bladeRadiiH[calcNode(bladeNode, numberOfBladeNodes, 0, 1, turbine, numberOfTurbines)];
+ this->bladeCoordsZH[node] = this->bladeRadiiH[calcNodeIndexInBladeArrays(bladeNode, numberOfBladeNodes, 0, 1, turbine, numberOfTurbines)];
}
}
}
@@ -596,7 +588,8 @@ void ActuatorFarm::setAllBladeForces(real* _bladeForcesX, real* _bladeForcesY, r
std::copy_n(_bladeForcesY, this->numberOfGridNodes, this->bladeForcesYH);
std::copy_n(_bladeForcesZ, this->numberOfGridNodes, this->bladeForcesZH);
-}void ActuatorFarm::setTurbineBladeCoords(uint turbine, real* _bladeCoordsX, real* _bladeCoordsY, real* _bladeCoordsZ)
+}
+void ActuatorFarm::setTurbineBladeCoords(uint turbine, real* _bladeCoordsX, real* _bladeCoordsY, real* _bladeCoordsZ)
{
std::copy_n(_bladeCoordsX, numberOfBladeNodes*numberOfBlades, &this->bladeCoordsXH[turbine*numberOfBladeNodes*numberOfBlades]);
std::copy_n(_bladeCoordsY, numberOfBladeNodes*numberOfBlades, &this->bladeCoordsYH[turbine*numberOfBladeNodes*numberOfBlades]);
diff --git a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarmInlines.h b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarmInlines.h
new file mode 100644
index 0000000000000000000000000000000000000000..ee19c4b55183eb610a10980951fc0a6621655b8b
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarmInlines.h
@@ -0,0 +1,20 @@
+#ifndef ACTUATOR_FARM_INLINES
+#define ACTUATOR_FARM_INLINES
+
+#include "basics/DataTypes.h"
+
+__host__ __device__ __inline__ uint calcNodeIndexInBladeArrays(uint bladeNode, uint numberOfNodesPerBlade, uint blade, uint numberOfBlades, uint turbine, uint numberOfTurbines)
+{
+ return bladeNode + numberOfNodesPerBlade * (blade + numberOfBlades * turbine);
+}
+
+__host__ __device__ __inline__ void calcTurbineBladeAndBladeNode(uint node, uint &bladeNode, uint numberOfNodesPerBlade, uint &blade, uint numberOfBlades, uint &turbine, uint numberOfTurbines)
+{
+ turbine = node / (numberOfNodesPerBlade * numberOfBlades);
+ uint x_off = turbine * numberOfNodesPerBlade * numberOfBlades;
+ blade = (node - x_off) / numberOfNodesPerBlade;
+ uint y_off = numberOfNodesPerBlade * blade + x_off;
+ bladeNode = node - y_off;
+}
+
+#endif
diff --git a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarmInlinesTest.cpp b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarmInlinesTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..720ff3c3c3217456ace0284e73d4953d08d7baa8
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarmInlinesTest.cpp
@@ -0,0 +1,99 @@
+#include <gmock/gmock.h>
+
+#include "ActuatorFarmInlines.h"
+
+TEST(ActuatorFarmInlinesTest, calcNodeIndexInBladeArrays)
+{
+ const uint numberOfNodesPerBlade = 4;
+ const uint numberOfBlades = 3;
+ const uint numberOfTurbines = 2;
+
+ // first node on first blade
+ uint bladeNode = 0;
+ uint blade = 0;
+ uint turbine = 0;
+ auto nodeIndexInBladeArrays = calcNodeIndexInBladeArrays(bladeNode, numberOfNodesPerBlade, blade, numberOfBlades, turbine, numberOfTurbines);
+ EXPECT_THAT(nodeIndexInBladeArrays, testing::Eq(0));
+
+ // last node on first blade
+ bladeNode = 3;
+ blade = 0;
+ turbine = 0;
+ nodeIndexInBladeArrays = calcNodeIndexInBladeArrays(bladeNode, numberOfNodesPerBlade, blade, numberOfBlades, turbine, numberOfTurbines);
+ EXPECT_THAT(nodeIndexInBladeArrays, testing::Eq(3));
+
+ // first node on third blade
+ bladeNode = 0;
+ blade = 2;
+ turbine = 0;
+ nodeIndexInBladeArrays = calcNodeIndexInBladeArrays(bladeNode, numberOfNodesPerBlade, blade, numberOfBlades, turbine, numberOfTurbines);
+ EXPECT_THAT(nodeIndexInBladeArrays, testing::Eq(8));
+
+ // last node on third blade, also last node on first turbine
+ bladeNode = 3;
+ blade = 2;
+ turbine = 0;
+ nodeIndexInBladeArrays = calcNodeIndexInBladeArrays(bladeNode, numberOfNodesPerBlade, blade, numberOfBlades, turbine, numberOfTurbines);
+ EXPECT_THAT(nodeIndexInBladeArrays, testing::Eq(11));
+
+ // first node on second turbine
+ bladeNode = 0;
+ blade = 0;
+ turbine = 1;
+ nodeIndexInBladeArrays = calcNodeIndexInBladeArrays(bladeNode, numberOfNodesPerBlade, blade, numberOfBlades, turbine, numberOfTurbines);
+ EXPECT_THAT(nodeIndexInBladeArrays, testing::Eq(12));
+
+ // last node on second turbine
+ bladeNode = 3;
+ blade = 2;
+ turbine = 1;
+ nodeIndexInBladeArrays = calcNodeIndexInBladeArrays(bladeNode, numberOfNodesPerBlade, blade, numberOfBlades, turbine, numberOfTurbines);
+ EXPECT_THAT(nodeIndexInBladeArrays, testing::Eq(23));
+}
+
+TEST(ActuatorFarmInlinesTest, calcTurbineBladeAndBladeNode)
+{
+ const uint numberOfNodesPerBlade = 4;
+ const uint numberOfBlades = 3;
+ const uint numberOfTurbines = 2;
+
+ uint bladeNode;
+ uint blade;
+ uint turbine;
+
+ uint node = 0; // first node on first blade
+ calcTurbineBladeAndBladeNode(node, bladeNode, numberOfNodesPerBlade, blade, numberOfBlades, turbine, numberOfTurbines);
+ EXPECT_THAT(bladeNode, testing::Eq(0));
+ EXPECT_THAT(blade, testing::Eq(0));
+ EXPECT_THAT(turbine, testing::Eq(0));
+
+ node = 3; // last node on first blade
+ calcTurbineBladeAndBladeNode(node, bladeNode, numberOfNodesPerBlade, blade, numberOfBlades, turbine, numberOfTurbines);
+ EXPECT_THAT(bladeNode, testing::Eq(3));
+ EXPECT_THAT(blade, testing::Eq(0));
+ EXPECT_THAT(turbine, testing::Eq(0));
+
+ node = 8; // first node on third blade
+ calcTurbineBladeAndBladeNode(node, bladeNode, numberOfNodesPerBlade, blade, numberOfBlades, turbine, numberOfTurbines);
+ EXPECT_THAT(bladeNode, testing::Eq(0));
+ EXPECT_THAT(blade, testing::Eq(2));
+ EXPECT_THAT(turbine, testing::Eq(0));
+
+ node = 11; // last node on third blade, also last node on first turbine
+ calcTurbineBladeAndBladeNode(node, bladeNode, numberOfNodesPerBlade, blade, numberOfBlades, turbine, numberOfTurbines);
+ EXPECT_THAT(bladeNode, testing::Eq(3));
+ EXPECT_THAT(blade, testing::Eq(2));
+ EXPECT_THAT(turbine, testing::Eq(0));
+
+ node = 12; // first node on second turbine
+ calcTurbineBladeAndBladeNode(node, bladeNode, numberOfNodesPerBlade, blade, numberOfBlades, turbine, numberOfTurbines);
+ EXPECT_THAT(bladeNode, testing::Eq(0));
+ EXPECT_THAT(blade, testing::Eq(0));
+ EXPECT_THAT(turbine, testing::Eq(1));
+
+ node = 23; // last node on second turbine
+ calcTurbineBladeAndBladeNode(node, bladeNode, numberOfNodesPerBlade, blade, numberOfBlades, turbine, numberOfTurbines);
+ EXPECT_THAT(bladeNode, testing::Eq(3));
+ EXPECT_THAT(blade, testing::Eq(2));
+ EXPECT_THAT(turbine, testing::Eq(1));
+}