From bfc30f5e559035c2190de68d316e48e8712957ee Mon Sep 17 00:00:00 2001
From: Henry <henry.korb@geo.uu.se>
Date: Fri, 7 Oct 2022 17:57:00 +0200
Subject: [PATCH] adds actuator farm
---
apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp | 52 +-
.../GPU/CudaMemoryManager.cpp | 248 ++++++++
.../VirtualFluids_GPU/GPU/CudaMemoryManager.h | 35 ++
.../PreCollisionInteractor/ActuatorFarm.cu | 551 ++++++++++++++++++
.../PreCollisionInteractor/ActuatorFarm.h | 168 ++++++
5 files changed, 1031 insertions(+), 23 deletions(-)
create mode 100644 src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.cu
create mode 100644 src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.h
diff --git a/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp b/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
index 58e5aede1..a6583214d 100644
--- a/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
+++ b/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
@@ -45,6 +45,7 @@
#include "VirtualFluids_GPU/Parameter/Parameter.h"
#include "VirtualFluids_GPU/Output/FileWriter.h"
#include "VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.h"
+#include "VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.h"
#include "VirtualFluids_GPU/PreCollisionInteractor/Probes/PointProbe.h"
#include "VirtualFluids_GPU/PreCollisionInteractor/Probes/PlaneProbe.h"
#include "VirtualFluids_GPU/Factories/BoundaryConditionFactory.h"
@@ -65,9 +66,9 @@ LbmOrGks lbmOrGks = LBM;
const real reference_diameter = 126.0; // diameter in m
-const real L_x = 10*reference_diameter;
-const real L_y = 6*reference_diameter;
-const real L_z = 6*reference_diameter;
+const real L_x = 30*reference_diameter;
+const real L_y = 10*reference_diameter;
+const real L_z = 10*reference_diameter;
const real viscosity = 1.56e-5;
@@ -75,7 +76,7 @@ const real velocity = 9.0;
const real mach = 0.1;
-const uint nodes_per_diameter = 16;
+const uint nodes_per_diameter = 32;
std::string path(".");
@@ -182,25 +183,30 @@ void multipleLevel(const std::string& configPath)
int level = 0;
uint nBlades = 3;
uint nBladeNodes = 32;
-
- SPtr<ActuatorLine> actuator_line =SPtr<ActuatorLine>( new ActuatorLine(nBlades, density, nBladeNodes, epsilon, turbPos[0], turbPos[1], turbPos[2], reference_diameter, level, dt, dx) );
- para->addActuator( actuator_line );
-
- SPtr<PointProbe> pointProbe = SPtr<PointProbe>( new PointProbe("pointProbe", para->getOutputPath(), 100, 1, 500, 100) );
- std::vector<real> probeCoordsX = {reference_diameter,2*reference_diameter,5*reference_diameter};
- std::vector<real> probeCoordsY = {3*reference_diameter,3*reference_diameter,3*reference_diameter};
- std::vector<real> probeCoordsZ = {3*reference_diameter,3*reference_diameter,3*reference_diameter};
- pointProbe->addProbePointsFromList(probeCoordsX, probeCoordsY, probeCoordsZ);
- // pointProbe->addProbePointsFromXNormalPlane(2*D, 0.0, 0.0, L_y, L_z, (uint)L_y/dx, (uint)L_z/dx);
-
- pointProbe->addStatistic(Statistic::Means);
- pointProbe->addStatistic(Statistic::Variances);
- para->addProbe( pointProbe );
-
- SPtr<PlaneProbe> planeProbe = SPtr<PlaneProbe>( new PlaneProbe("planeProbe", para->getOutputPath(), 100, 500, 100, 100) );
- planeProbe->setProbePlane(5*reference_diameter, 0, 0, dx, L_y, L_z);
- planeProbe->addStatistic(Statistic::Means);
- para->addProbe( planeProbe );
+ real omega = 1.0f;
+ // SPtr<ActuatorLine> actuator_line =SPtr<ActuatorLine>( new ActuatorLine(nBlades, density, nBladeNodes, epsilon, turbPos[0], turbPos[1], turbPos[2], reference_diameter, level, dt, dx) );
+ SPtr<ActuatorFarm> actuator_farm = std::make_shared<ActuatorFarm>(nBlades, density, nBladeNodes, epsilon, level, dt, dx, false);
+ std::vector<real> bladeRadii;
+ real dr = reference_diameter/nBladeNodes;
+ for(uint node=0; node<nBladeNodes; node++){bladeRadii.emplace_back(dr*(node+1));}
+ actuator_farm->addTurbine(turbPos[0], turbPos[1], turbPos[2], reference_diameter, omega, 0, 0, bladeRadii);
+ para->addActuator( actuator_farm );
+
+ // SPtr<PointProbe> pointProbe = SPtr<PointProbe>( new PointProbe("pointProbe", para->getOutputPath(), 100, 1, 500, 100) );
+ // std::vector<real> probeCoordsX = {reference_diameter,2*reference_diameter,5*reference_diameter};
+ // std::vector<real> probeCoordsY = {3*reference_diameter,3*reference_diameter,3*reference_diameter};
+ // std::vector<real> probeCoordsZ = {3*reference_diameter,3*reference_diameter,3*reference_diameter};
+ // pointProbe->addProbePointsFromList(probeCoordsX, probeCoordsY, probeCoordsZ);
+ // // pointProbe->addProbePointsFromXNormalPlane(2*D, 0.0, 0.0, L_y, L_z, (uint)L_y/dx, (uint)L_z/dx);
+
+ // pointProbe->addStatistic(Statistic::Means);
+ // pointProbe->addStatistic(Statistic::Variances);
+ // para->addProbe( pointProbe );
+
+ // SPtr<PlaneProbe> planeProbe = SPtr<PlaneProbe>( new PlaneProbe("planeProbe", para->getOutputPath(), 100, 500, 100, 100) );
+ // planeProbe->setProbePlane(5*reference_diameter, 0, 0, dx, L_y, L_z);
+ // planeProbe->addStatistic(Statistic::Means);
+ // para->addProbe( planeProbe );
auto cudaMemoryManager = std::make_shared<CudaMemoryManager>(para);
diff --git a/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp b/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp
index 221922169..f3d46757a 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp
+++ b/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp
@@ -7,6 +7,7 @@
#include <Parameter/Parameter.h>
#include "Parameter/CudaStreamManager.h"
#include "PreCollisionInteractor/ActuatorLine.h"
+#include "PreCollisionInteractor/ActuatorFarm.h"
#include "PreCollisionInteractor/Probes/Probe.h"
#include "Calculation/PorousMedia.h"
@@ -3147,6 +3148,253 @@ void CudaMemoryManager::cudaFreeSphereIndices(ActuatorLine* actuatorLine)
checkCudaErrors( cudaFreeHost(actuatorLine->boundingSphereIndicesH) );
checkCudaErrors( cudaFree(actuatorLine->boundingSphereIndicesD) );
}
+////////////////////////////////////////////////////////////////////////////////////
+// ActuatorFarm
+///////////////////////////////////////////////////////////////////////////////
+void CudaMemoryManager::cudaAllocBladeGeometries(ActuatorFarm* actuatorFarm)
+{
+ uint sizeRealTurbine = sizeof(real)*actuatorFarm->getNTurbines();
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->bladeRadiiH, sizeRealTurbine*actuatorFarm->getNBladeNodes()) );
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->diametersH, sizeRealTurbine) );
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->turbinePosXH, sizeRealTurbine) );
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->turbinePosYH, sizeRealTurbine) );
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->turbinePosZH, sizeRealTurbine) );
+
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->bladeRadiiD, sizeRealTurbine*actuatorFarm->getNBladeNodes()) );
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->diametersD, sizeRealTurbine) );
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->turbinePosXD, sizeRealTurbine) );
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->turbinePosYD, sizeRealTurbine) );
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->turbinePosZD, sizeRealTurbine) );
+ setMemsizeGPU(sizeof(real)*(actuatorFarm->getNBladeNodes()+4)*actuatorFarm->getNTurbines(), false);
+
+}
+void CudaMemoryManager::cudaCopyBladeGeometriesHtoD(ActuatorFarm* actuatorFarm)
+{
+ uint sizeRealTurbine = sizeof(real)*actuatorFarm->getNTurbines();
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeRadiiD, actuatorFarm->bladeRadiiH, sizeRealTurbine*actuatorFarm->getNBladeNodes(), cudaMemcpyHostToDevice) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->diametersD, actuatorFarm->diametersH, sizeRealTurbine, cudaMemcpyHostToDevice) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->turbinePosXD, actuatorFarm->turbinePosXH, sizeRealTurbine, cudaMemcpyHostToDevice) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->turbinePosYD, actuatorFarm->turbinePosYH, sizeRealTurbine, cudaMemcpyHostToDevice) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->turbinePosZD, actuatorFarm->turbinePosZH, sizeRealTurbine, cudaMemcpyHostToDevice) );
+
+}
+void CudaMemoryManager::cudaCopyBladeGeometriesDtoH(ActuatorFarm* actuatorFarm)
+{
+ uint sizeRealTurbine = sizeof(real)*actuatorFarm->getNTurbines();
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeRadiiH, actuatorFarm->bladeRadiiD, sizeRealTurbine*actuatorFarm->getNBladeNodes(), cudaMemcpyDeviceToHost) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->diametersH, actuatorFarm->diametersD, sizeRealTurbine, cudaMemcpyDeviceToHost) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->turbinePosXH, actuatorFarm->turbinePosXD, sizeRealTurbine, cudaMemcpyDeviceToHost) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->turbinePosYH, actuatorFarm->turbinePosYD, sizeRealTurbine, cudaMemcpyDeviceToHost) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->turbinePosZH, actuatorFarm->turbinePosZD, sizeRealTurbine, cudaMemcpyDeviceToHost) );
+
+}
+void CudaMemoryManager::cudaFreeBladeGeometries(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaFree(actuatorFarm->bladeRadiiD) );
+ checkCudaErrors( cudaFree(actuatorFarm->diametersD) );
+ checkCudaErrors( cudaFree(actuatorFarm->turbinePosXD) );
+ checkCudaErrors( cudaFree(actuatorFarm->turbinePosYD) );
+ checkCudaErrors( cudaFree(actuatorFarm->turbinePosZD) );
+
+ checkCudaErrors( cudaFreeHost(actuatorFarm->bladeRadiiH) );
+ checkCudaErrors( cudaFreeHost(actuatorFarm->diametersH) );
+ checkCudaErrors( cudaFreeHost(actuatorFarm->turbinePosXH) );
+ checkCudaErrors( cudaFreeHost(actuatorFarm->turbinePosYH) );
+ checkCudaErrors( cudaFreeHost(actuatorFarm->turbinePosZH) );
+}
+
+void CudaMemoryManager::cudaAllocBladeOrientations(ActuatorFarm* actuatorFarm)
+{
+ uint sizeRealTurbine = sizeof(real)*actuatorFarm->getNTurbines();
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->omegasH, sizeRealTurbine) );
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->azimuthsH, sizeRealTurbine) );
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->yawsH, sizeRealTurbine) );
+
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->omegasD, sizeRealTurbine) );
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->azimuthsD, sizeRealTurbine) );
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->yawsD, sizeRealTurbine) );
+
+ setMemsizeGPU(3*sizeRealTurbine, false);
+
+}
+void CudaMemoryManager::cudaCopyBladeOrientationsHtoD(ActuatorFarm* actuatorFarm)
+{
+ uint sizeRealTurbine = sizeof(real)*actuatorFarm->getNTurbines();
+ checkCudaErrors( cudaMemcpy(actuatorFarm->omegasD, actuatorFarm->omegasH, sizeRealTurbine, cudaMemcpyHostToDevice) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->azimuthsD, actuatorFarm->azimuthsH, sizeRealTurbine, cudaMemcpyHostToDevice) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->yawsD, actuatorFarm->yawsH, sizeRealTurbine, cudaMemcpyHostToDevice) );
+
+}
+void CudaMemoryManager::cudaCopyBladeOrientationsDtoH(ActuatorFarm* actuatorFarm)
+{
+ uint sizeRealTurbine = sizeof(real)*actuatorFarm->getNTurbines();
+ checkCudaErrors( cudaMemcpy(actuatorFarm->omegasH, actuatorFarm->omegasD, sizeRealTurbine, cudaMemcpyDeviceToHost) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->azimuthsH, actuatorFarm->azimuthsD, sizeRealTurbine, cudaMemcpyDeviceToHost) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->yawsH, actuatorFarm->yawsD, sizeRealTurbine, cudaMemcpyDeviceToHost) );
+}
+void CudaMemoryManager::cudaFreeBladeOrientations(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaFree((void**) &actuatorFarm->omegasD) );
+ checkCudaErrors( cudaFree((void**) &actuatorFarm->azimuthsD) );
+ checkCudaErrors( cudaFree((void**) &actuatorFarm->yawsD) );
+
+ checkCudaErrors( cudaFreeHost((void**) &actuatorFarm->omegasH) );
+ checkCudaErrors( cudaFreeHost((void**) &actuatorFarm->azimuthsH) );
+ checkCudaErrors( cudaFreeHost((void**) &actuatorFarm->yawsH) );
+}
+
+void CudaMemoryManager::cudaAllocBladeCoords(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->bladeCoordsXH, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->bladeCoordsYH, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->bladeCoordsZH, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->bladeCoordsXD, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->bladeCoordsYD, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->bladeCoordsZD, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+
+ setMemsizeGPU(3.f*actuatorFarm->getNumberOfNodes(), false);
+}
+
+void CudaMemoryManager::cudaCopyBladeCoordsHtoD(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeCoordsXD, actuatorFarm->bladeCoordsXH, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyHostToDevice) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeCoordsYD, actuatorFarm->bladeCoordsYH, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyHostToDevice) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeCoordsZD, actuatorFarm->bladeCoordsZH, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyHostToDevice) );
+}
+
+void CudaMemoryManager::cudaCopyBladeCoordsDtoH(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeCoordsXH, actuatorFarm->bladeCoordsXD, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeCoordsYH, actuatorFarm->bladeCoordsYD, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeCoordsZH, actuatorFarm->bladeCoordsZD, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
+}
+
+void CudaMemoryManager::cudaFreeBladeCoords(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaFree(actuatorFarm->bladeCoordsXD) );
+ checkCudaErrors( cudaFree(actuatorFarm->bladeCoordsYD) );
+ checkCudaErrors( cudaFree(actuatorFarm->bladeCoordsZD) );
+
+ checkCudaErrors( cudaFreeHost(actuatorFarm->bladeCoordsXH) );
+ checkCudaErrors( cudaFreeHost(actuatorFarm->bladeCoordsYH) );
+ checkCudaErrors( cudaFreeHost(actuatorFarm->bladeCoordsZH) );
+}
+
+void CudaMemoryManager::cudaAllocBladeIndices(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->bladeIndicesH, sizeof(uint)*actuatorFarm->getNumberOfNodes()) );
+
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->bladeIndicesD, sizeof(uint)*actuatorFarm->getNumberOfNodes()) );
+
+ setMemsizeGPU(sizeof(uint)*actuatorFarm->getNumberOfNodes(), false);
+}
+
+void CudaMemoryManager::cudaCopyBladeIndicesHtoD(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeIndicesD, actuatorFarm->bladeIndicesH, sizeof(uint)*actuatorFarm->getNumberOfNodes(), cudaMemcpyHostToDevice) );
+}
+
+void CudaMemoryManager::cudaFreeBladeIndices(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaFree(actuatorFarm->bladeIndicesD) );
+
+ checkCudaErrors( cudaFreeHost(actuatorFarm->bladeIndicesH) );
+}
+
+void CudaMemoryManager::cudaAllocBladeVelocities(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->bladeVelocitiesXH, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->bladeVelocitiesYH, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->bladeVelocitiesZH, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->bladeVelocitiesXD, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->bladeVelocitiesYD, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->bladeVelocitiesZD, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+
+ setMemsizeGPU(3.*sizeof(real)*actuatorFarm->getNumberOfNodes(), false);
+}
+
+void CudaMemoryManager::cudaCopyBladeVelocitiesHtoD(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeVelocitiesXD, actuatorFarm->bladeVelocitiesXH, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyHostToDevice) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeVelocitiesYD, actuatorFarm->bladeVelocitiesYH, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyHostToDevice) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeVelocitiesZD, actuatorFarm->bladeVelocitiesZH, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyHostToDevice) );
+}
+
+void CudaMemoryManager::cudaCopyBladeVelocitiesDtoH(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeVelocitiesXH, actuatorFarm->bladeVelocitiesXD, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeVelocitiesYH, actuatorFarm->bladeVelocitiesYD, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeVelocitiesZH, actuatorFarm->bladeVelocitiesZD, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
+}
+
+void CudaMemoryManager::cudaFreeBladeVelocities(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaFree(actuatorFarm->bladeVelocitiesXD) );
+ checkCudaErrors( cudaFree(actuatorFarm->bladeVelocitiesYD) );
+ checkCudaErrors( cudaFree(actuatorFarm->bladeVelocitiesZD) );
+
+ checkCudaErrors( cudaFreeHost(actuatorFarm->bladeVelocitiesXH) );
+ checkCudaErrors( cudaFreeHost(actuatorFarm->bladeVelocitiesYH) );
+ checkCudaErrors( cudaFreeHost(actuatorFarm->bladeVelocitiesZH) );
+}
+
+void CudaMemoryManager::cudaAllocBladeForces(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->bladeForcesXH, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->bladeForcesYH, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+ checkCudaErrors( cudaMallocHost((void**) &actuatorFarm->bladeForcesZH, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->bladeForcesXD, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->bladeForcesYD, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+ checkCudaErrors( cudaMalloc((void**) &actuatorFarm->bladeForcesZD, sizeof(real)*actuatorFarm->getNumberOfNodes()) );
+
+ setMemsizeGPU(3.*sizeof(real)*actuatorFarm->getNumberOfNodes(), false);
+}
+
+void CudaMemoryManager::cudaCopyBladeForcesHtoD(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeForcesXD, actuatorFarm->bladeForcesXH, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyHostToDevice) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeForcesYD, actuatorFarm->bladeForcesYH, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyHostToDevice) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeForcesZD, actuatorFarm->bladeForcesZH, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyHostToDevice) );
+}
+
+void CudaMemoryManager::cudaCopyBladeForcesDtoH(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeForcesXH, actuatorFarm->bladeForcesXD, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeForcesYH, actuatorFarm->bladeForcesYD, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
+ checkCudaErrors( cudaMemcpy(actuatorFarm->bladeForcesZH, actuatorFarm->bladeForcesZD, sizeof(real)*actuatorFarm->getNumberOfNodes(), cudaMemcpyDeviceToHost) );
+}
+
+void CudaMemoryManager::cudaFreeBladeForces(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaFree(actuatorFarm->bladeForcesXD) );
+ checkCudaErrors( cudaFree(actuatorFarm->bladeForcesYD) );
+ checkCudaErrors( cudaFree(actuatorFarm->bladeForcesZD) );
+
+ checkCudaErrors( cudaFreeHost(actuatorFarm->bladeForcesXH) );
+ checkCudaErrors( cudaFreeHost(actuatorFarm->bladeForcesYH) );
+ checkCudaErrors( cudaFreeHost(actuatorFarm->bladeForcesZH) );
+}
+
+void CudaMemoryManager::cudaAllocSphereIndices(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaMallocHost((void**) &(actuatorFarm->boundingSphereIndicesH), sizeof(int)*actuatorFarm->getNumberOfIndices()));
+ checkCudaErrors( cudaMalloc((void**) &(actuatorFarm->boundingSphereIndicesD), sizeof(int)*actuatorFarm->getNumberOfIndices()));
+ setMemsizeGPU(sizeof(int)*actuatorFarm->getNumberOfIndices(), false);
+}
+
+void CudaMemoryManager::cudaCopySphereIndicesHtoD(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaMemcpy(actuatorFarm->boundingSphereIndicesD, actuatorFarm->boundingSphereIndicesH, sizeof(int)*actuatorFarm->getNumberOfIndices(), cudaMemcpyHostToDevice) );
+}
+
+void CudaMemoryManager::cudaFreeSphereIndices(ActuatorFarm* actuatorFarm)
+{
+ checkCudaErrors( cudaFreeHost(actuatorFarm->boundingSphereIndicesH) );
+ checkCudaErrors( cudaFree(actuatorFarm->boundingSphereIndicesD) );
+}
////////////////////////////////////////////////////////////////////////////////////
// Probe
diff --git a/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.h b/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.h
index 8df1b412c..beb87ba63 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.h
+++ b/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.h
@@ -19,6 +19,7 @@
class Parameter;
class PorousMedia;
class ActuatorLine;
+class ActuatorFarm;
class Probe;
class VIRTUALFLUIDS_GPU_EXPORT CudaMemoryManager
@@ -381,6 +382,40 @@ public:
void cudaCopySphereIndicesHtoD(ActuatorLine* actuatorLine);
void cudaFreeSphereIndices(ActuatorLine* actuatorLine);
+ // ActuatorFarm
+ void cudaAllocBladeGeometries(ActuatorFarm* actuatorFarm);
+ void cudaCopyBladeGeometriesHtoD(ActuatorFarm* actuatorFarm);
+ void cudaCopyBladeGeometriesDtoH(ActuatorFarm* actuatorFarm);
+ void cudaFreeBladeGeometries(ActuatorFarm* actuatorFarm);
+
+ void cudaAllocBladeOrientations(ActuatorFarm* actuatorFarm);
+ void cudaCopyBladeOrientationsHtoD(ActuatorFarm* actuatorFarm);
+ void cudaCopyBladeOrientationsDtoH(ActuatorFarm* actuatorFarm);
+ void cudaFreeBladeOrientations(ActuatorFarm* actuatorFarm);
+
+ void cudaAllocBladeCoords(ActuatorFarm* actuatorFarm);
+ void cudaCopyBladeCoordsHtoD(ActuatorFarm* actuatorFarm);
+ void cudaCopyBladeCoordsDtoH(ActuatorFarm* actuatorFarm);
+ void cudaFreeBladeCoords(ActuatorFarm* actuatorFarm);
+
+ void cudaAllocBladeIndices(ActuatorFarm* actuatorFarm);
+ void cudaCopyBladeIndicesHtoD(ActuatorFarm* actuatorFarm);
+ void cudaFreeBladeIndices(ActuatorFarm* actuatorFarm);
+
+ void cudaAllocBladeVelocities(ActuatorFarm* actuatorFarm);
+ void cudaCopyBladeVelocitiesHtoD(ActuatorFarm* actuatorFarm);
+ void cudaCopyBladeVelocitiesDtoH(ActuatorFarm* actuatorFarm);
+ void cudaFreeBladeVelocities(ActuatorFarm* actuatorFarm);
+
+ void cudaAllocBladeForces(ActuatorFarm* actuatorFarm);
+ void cudaCopyBladeForcesHtoD(ActuatorFarm* actuatorFarm);
+ void cudaCopyBladeForcesDtoH(ActuatorFarm* actuatorFarm);
+ void cudaFreeBladeForces(ActuatorFarm* actuatorFarm);
+
+ void cudaAllocSphereIndices(ActuatorFarm* actuatorFarm);
+ void cudaCopySphereIndicesHtoD(ActuatorFarm* actuatorFarm);
+ void cudaFreeSphereIndices(ActuatorFarm* actuatorFarm);
+ // Probes
void cudaAllocProbeDistances(Probe* probe, int level);
void cudaCopyProbeDistancesHtoD(Probe* probe, int level);
void cudaCopyProbeDistancesDtoH(Probe* probe, int level);
diff --git a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.cu b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.cu
new file mode 100644
index 000000000..c5e1926d0
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.cu
@@ -0,0 +1,551 @@
+#include "ActuatorFarm.h"
+
+#include <cuda.h>
+#include <cuda_runtime.h>
+#include <helper_cuda.h>
+
+#include <cuda/CudaGrid.h>
+#include "VirtualFluids_GPU/GPU/GeometryUtils.h"
+#include "VirtualFluids_GPU/Kernel/Utilities/DistributionHelper.cuh"
+
+#include "Parameter/Parameter.h"
+#include "DataStructureInitializer/GridProvider.h"
+#include "GPU/CudaMemoryManager.h"
+#include "lbm/constants/NumericConstants.h"
+
+using namespace vf::lbm::constant;
+
+
+__host__ __device__ __inline__ uint calcNode(uint bladeNode, uint nBladeNodes, uint blade, uint nBlades, uint turbine, uint nTurbines)
+{
+
+ return bladeNode+nBladeNodes*(blade+turbine*nBlades);
+}
+
+__host__ __device__ __inline__ void calcTurbineBladeAndBladeNode(uint node, uint& bladeNode, uint nBladeNodes, uint& blade, uint nBlades, uint& turbine, uint nTurbines)
+{
+ turbine = node/(nBladeNodes*nBlades);
+ uint x_off = turbine*nBladeNodes*nBlades;
+ blade = (node - x_off)/nBlades;
+ uint y_off = nBladeNodes*blade+x_off;
+ bladeNode = (node - y_off)/nBladeNodes;
+}
+
+__host__ __device__ __forceinline__ real distSqrd(real distX, real distY, real distZ)
+{
+ return distX*distX+distY*distY+distZ*distZ;
+}
+
+__host__ __device__ __inline__ void rotateFromBladeToGlobal(
+ real& bladeCoordX_BF, real& bladeCoordY_BF, real& bladeCoordZ_BF,
+ real& bladeCoordX_GF, real& bladeCoordY_GF, real& bladeCoordZ_GF,
+ real& azimuth, real& yaw)
+{
+ real tmpX, tmpY, tmpZ;
+
+ 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(
+ real& bladeCoordX_BF, real& bladeCoordY_BF, real& bladeCoordZ_BF,
+ real& bladeCoordX_GF, real& bladeCoordY_GF, real& bladeCoordZ_GF,
+ real& azimuth, real& yaw)
+{
+ real tmpX, tmpY, tmpZ;
+
+ invRotateAboutZ3D(yaw, bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF, tmpX, tmpY, tmpZ);
+ invRotateAboutX3D(azimuth, tmpX, tmpY, tmpZ, bladeCoordX_BF, bladeCoordY_BF, bladeCoordZ_BF);
+}
+
+__global__ void interpolateVelocities(real* gridCoordsX, real* gridCoordsY, real* gridCoordsZ,
+ uint* neighborsX, uint* neighborsY, uint* neighborsZ, uint* neighborsWSB,
+ real* vx, real* vy, real* vz,
+ real* bladeCoordsX, real* bladeCoordsY, real* bladeCoordsZ,
+ real* bladeVelocitiesX, real* bladeVelocitiesY, real* bladeVelocitiesZ,
+ uint nTurbines, uint nBlades, uint nBladeNodes,
+ real* azimuths, real* yaws, real* omegas,
+ real* turbPosX, real* turbPosY, real* turbPosZ,
+ uint* bladeIndices, real velocityRatio, real invDeltaX)
+{
+
+ const uint node = vf::gpu::getNodeIndex();
+
+ if(node>=nBladeNodes*nBlades*nTurbines) return;
+
+ uint turbine, bladeNode, blade;
+
+ calcTurbineBladeAndBladeNode(node, bladeNode, nBladeNodes, blade, nBlades, turbine, nTurbines);
+
+ real bladeCoordX_BF = bladeCoordsX[node];
+ real bladeCoordY_BF = bladeCoordsY[node];
+ real bladeCoordZ_BF = bladeCoordsZ[node];
+
+ real bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF;
+
+ real localAzimuth = azimuths[turbine]+blade*c2Pi/nBlades;
+ 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];
+
+ uint k, ke, kn, kt;
+ uint kne, kte, ktn, ktne;
+
+ k = findNearestCellBSW(bladeIndices[node],
+ gridCoordsX, gridCoordsY, gridCoordsZ,
+ bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF,
+ neighborsX, neighborsY, neighborsZ, neighborsWSB);
+
+ bladeIndices[node] = k;
+
+ getNeighborIndicesOfBSW(k, ke, kn, kt, kne, kte, ktn, ktne, neighborsX, neighborsY, neighborsZ);
+
+ real dW, dE, dN, dS, dT, dB;
+
+ 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);
+
+ real bladeVelX_GF = trilinearInterpolation(dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, vx)*velocityRatio;
+ real bladeVelY_GF = trilinearInterpolation(dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, vy)*velocityRatio;
+ real bladeVelZ_GF = trilinearInterpolation(dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, vz)*velocityRatio;
+
+ real bladeVelX_BF, bladeVelY_BF, bladeVelZ_BF;
+
+
+ rotateFromGlobalToBlade(bladeVelX_BF, bladeVelY_BF, bladeVelZ_BF,
+ bladeVelX_GF, bladeVelY_GF, bladeVelZ_GF,
+ localAzimuth, yaw);
+
+ bladeVelocitiesX[node] = bladeVelX_BF;
+ bladeVelocitiesY[node] = bladeVelY_BF+omegas[turbine]*bladeCoordZ_BF;
+ bladeVelocitiesZ[node] = bladeVelZ_BF;
+}
+
+
+__global__ void applyBodyForces(real* gridCoordsX, real* gridCoordsY, real* gridCoordsZ,
+ real* gridForcesX, real* gridForcesY, real* gridForcesZ,
+ real* bladeCoordsX, real* bladeCoordsY, real* bladeCoordsZ,
+ real* bladeForcesX, real* bladeForcesY,real* bladeForcesZ,
+ uint nTurbines, uint nBlades, uint nBladeNodes,
+ real* azimuths, real* yaws, real* diameters,
+ real* turbPosX, real* turbPosY, real* turbPosZ,
+ uint* gridIndices, uint nIndices,
+ real invEpsilonSqrd, real factorGaussian)
+{
+ const uint x = threadIdx.x;
+ const uint y = blockIdx.x;
+ const uint z = blockIdx.y;
+
+ const uint nx = blockDim.x;
+ const uint ny = gridDim.x;
+
+ const uint index = nx*(ny*z + y) + x;
+
+ if(index>=nIndices) return;
+
+
+ uint gridIndex = gridIndices[index];
+
+ real gridCoordX_GF = gridCoordsX[gridIndex];
+ real gridCoordY_GF = gridCoordsY[gridIndex];
+ real gridCoordZ_GF = gridCoordsZ[gridIndex];
+
+ real gridForceX_RF = c0o1;
+ real gridForceY_RF = c0o1;
+ real gridForceZ_RF = c0o1;
+
+ real dAzimuth = c2Pi/nBlades;
+
+ for(uint turbine = 0; turbine<nTurbines; turbine++)
+ {
+ real radius = c1o2*diameters[turbine];
+ real gridCoordX_RF = gridCoordX_GF - turbPosX[turbine];
+ real gridCoordY_RF = gridCoordY_GF - turbPosY[turbine];
+ real gridCoordZ_RF = gridCoordZ_GF - turbPosZ[turbine];
+
+ if(distSqrd(gridCoordX_RF, gridCoordY_RF, gridCoordZ_RF)*invEpsilonSqrd< radius*radius*invEpsilonSqrd+c7o1)
+ continue;
+
+ real azimuth = azimuths[turbine];
+ real yaw = yaws[turbine];
+
+ for( uint blade=0; blade<nBlades; blade++)
+ {
+ real localAzimuth = azimuth+blade*dAzimuth;
+
+
+ real gridCoordX_BF, gridCoordY_BF, gridCoordZ_BF;
+
+ rotateFromGlobalToBlade(gridCoordX_BF, gridCoordY_BF, gridCoordZ_BF,
+ gridCoordX_RF, gridCoordY_RF, gridCoordZ_RF,
+ localAzimuth, yaw);
+
+ //Handle first and last node separately
+
+ uint node = calcNode(0, nBladeNodes, blade, nBlades, turbine, nTurbines);
+ uint nextNode = calcNode(1, nBladeNodes, blade, nBlades, turbine, nTurbines);
+
+ real last_z = c0o1;
+ real current_z = c0o1;
+ real next_z = bladeCoordsZ[node];
+
+ real x, y, dz;
+
+ x = bladeCoordsX[node];
+ y = bladeCoordsY[node];
+ last_z = current_z;
+ current_z = next_z;
+ next_z = bladeCoordsZ[nextNode];
+
+ dz = c1o2*(next_z+current_z);
+
+ real eta = dz*factorGaussian*exp(-distSqrd(x-gridCoordX_BF, y-gridCoordY_BF, current_z-gridCoordZ_BF)*invEpsilonSqrd);
+
+ real forceX_RF, forceY_RF, forceZ_RF;
+
+ 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;
+
+ for( uint bladeNode=1; bladeNode<nBladeNodes-1; bladeNode++)
+ {
+ node = nextNode;
+ nextNode = calcNode(bladeNode+1, nBladeNodes, blade, nBlades, turbine, nTurbines);
+
+ x = bladeCoordsX[node];
+ y = bladeCoordsY[node];
+ last_z = current_z;
+ current_z = next_z;
+ next_z = bladeCoordsZ[node+1];
+
+ dz = c1o2*(next_z-last_z);
+
+ eta = dz*factorGaussian*exp(-distSqrd(x-gridCoordX_BF, y-gridCoordY_BF, current_z-gridCoordZ_BF)*invEpsilonSqrd);
+
+ 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;
+ }
+
+ node = calcNode(nBladeNodes-1, nBladeNodes, blade, nBlades, turbine, nTurbines);
+
+ x = bladeCoordsX[node];
+ y = bladeCoordsY[node];
+ last_z = current_z;
+ current_z = next_z;
+
+ dz = radius-c1o2*(current_z+last_z);
+
+ eta = dz*factorGaussian*exp(-distSqrd(x-gridCoordX_BF, y-gridCoordY_BF, current_z-gridCoordZ_BF)*invEpsilonSqrd);
+
+ 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;
+ }
+
+ gridForcesX[gridIndex] = gridForceX_RF;
+ gridForcesY[gridIndex] = gridForceY_RF;
+ gridForcesZ[gridIndex] = gridForceZ_RF;
+ }
+}
+
+void ActuatorFarm::addTurbine(real posX, real posY, real posZ, real diameter, real omega, real azimuth, real yaw, std::vector<real> bladeRadii)
+{
+ preInitPosX.push_back(posX);
+ preInitPosY.push_back(posY);
+ preInitPosZ.push_back(posZ);
+ preInitOmegas.push_back(omega);
+ preInitAzimuths.push_back(azimuth);
+ preInitYaws.push_back(yaw);
+ preInitDiameters.push_back(diameter);
+ preInitBladeRadii.push_back(bladeRadii);
+}
+
+void ActuatorFarm::init(Parameter* para, GridProvider* gridProvider, CudaMemoryManager* cudaMemoryManager)
+{
+ if(!para->getIsBodyForce()) throw std::runtime_error("try to allocate ActuatorFarm but BodyForce is not set in Parameter.");
+ this->initTurbineGeometries(cudaMemoryManager);
+ this->initBladeCoords(cudaMemoryManager);
+ this->initBladeIndices(para, cudaMemoryManager);
+ this->initBladeVelocities(cudaMemoryManager);
+ this->initBladeForces(cudaMemoryManager);
+ this->initBoundingSphere(para, cudaMemoryManager);
+}
+
+
+void ActuatorFarm::interact(Parameter* para, CudaMemoryManager* cudaMemoryManager, int level, unsigned int t)
+{
+ if (level != this->level) return;
+
+ if(useHostArrays) cudaMemoryManager->cudaCopyBladeCoordsHtoD(this);
+
+ vf::cuda::CudaGrid bladeGrid = vf::cuda::CudaGrid(para->getParH(level)->numberofthreads, this->numberOfNodes);
+
+ interpolateVelocities<<< bladeGrid.grid, bladeGrid.threads >>>(
+ para->getParD(this->level)->coordinateX, para->getParD(this->level)->coordinateY, para->getParD(this->level)->coordinateZ,
+ para->getParD(this->level)->neighborX, para->getParD(this->level)->neighborY, para->getParD(this->level)->neighborZ, para->getParD(this->level)->neighborInverse,
+ para->getParD(this->level)->velocityX, para->getParD(this->level)->velocityY, para->getParD(this->level)->velocityZ,
+ this->bladeCoordsXD, this->bladeCoordsYD, this->bladeCoordsZD,
+ this->bladeVelocitiesXD, this->bladeVelocitiesYD, this->bladeVelocitiesZD,
+ this->nTurbines, this->nBlades, this->nBladeNodes,
+ this->azimuthsD, this->yawsD, this->omegasD,
+ this->turbinePosXD, this->turbinePosYD, this->turbinePosZD,
+ this->bladeIndicesD, para->getVelocityRatio(), this->invDeltaX);
+
+ if(useHostArrays) cudaMemoryManager->cudaCopyBladeVelocitiesDtoH(this);
+
+ this->calcBladeForces();
+
+ if(useHostArrays) cudaMemoryManager->cudaCopyBladeForcesHtoD(this);
+
+ vf::cuda::CudaGrid sphereGrid = vf::cuda::CudaGrid(para->getParH(level)->numberofthreads, this->numberOfIndices);
+
+ applyBodyForces<<<sphereGrid.grid, sphereGrid.threads>>>(
+ para->getParD(this->level)->coordinateX, para->getParD(this->level)->coordinateY, para->getParD(this->level)->coordinateZ,
+ para->getParD(this->level)->forceX_SP, para->getParD(this->level)->forceY_SP, para->getParD(this->level)->forceZ_SP,
+ this->bladeCoordsXD, this->bladeCoordsYD, this->bladeCoordsZD,
+ this->bladeForcesXD, this->bladeForcesYD, this->bladeForcesZD,
+ this->nTurbines, this->nBlades, this->nBladeNodes,
+ this->azimuthsD, this->yawsD, this->diametersD,
+ this->turbinePosXD, this->turbinePosYD, this->turbinePosZD,
+ this->boundingSphereIndicesD, this->numberOfIndices,
+ this->invEpsilonSqrd, this->factorGaussian);
+
+ for(uint turbine=0; turbine<this->nTurbines; turbine++)
+ this->azimuthsH[turbine] = fmod(this->azimuthsH[turbine]+this->omegasH[turbine]*this->deltaT, c2Pi);
+
+ cudaMemoryManager->cudaCopyBladeOrientationsHtoD(this);
+}
+
+
+void ActuatorFarm::free(Parameter* para, CudaMemoryManager* cudaMemoryManager)
+{
+ cudaMemoryManager->cudaFreeBladeGeometries(this);
+ cudaMemoryManager->cudaFreeBladeOrientations(this);
+ cudaMemoryManager->cudaFreeBladeCoords(this);
+ cudaMemoryManager->cudaFreeBladeVelocities(this);
+ cudaMemoryManager->cudaFreeBladeForces(this);
+ cudaMemoryManager->cudaFreeBladeIndices(this);
+ cudaMemoryManager->cudaFreeSphereIndices(this);
+}
+
+
+void ActuatorFarm::calcForcesEllipticWing()
+{
+ real u_rel, v_rel, u_rel_sq;
+ real phi;
+ real Cl = c1o1;
+ real Cd = c0o1;
+ real c0 = c1o1;
+
+ real c, Cn, Ct;
+ for(uint turbine=0; turbine<this->nTurbines; turbine++)
+ {
+ real diameter = this->diametersH[turbine];
+ for( uint blade=0; blade<this->nBlades; blade++)
+ {
+ for( uint bladeNode=0; bladeNode<this->nBladeNodes; bladeNode++)
+ {
+ uint node = calcNode(bladeNode, this->nBladeNodes, blade, this->nBlades, turbine, this->nTurbines);
+
+ u_rel = this->bladeVelocitiesXH[node];
+ v_rel = this->bladeVelocitiesYH[node];
+ u_rel_sq = u_rel*u_rel+v_rel*v_rel;
+ phi = atan2(u_rel, v_rel);
+
+ real tmp = c4o1*this->bladeRadiiH[bladeNode]/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] = -c1o2*u_rel_sq*c*this->density*Cn;
+ this->bladeForcesYH[node] = -c1o2*u_rel_sq*c*this->density*Ct;
+ this->bladeForcesZH[node] = c0o1;
+ }
+ }
+ }
+}
+
+void ActuatorFarm::calcBladeForces()
+{
+ // this->calcForcesEllipticWing();
+}
+void ActuatorFarm::initTurbineGeometries(CudaMemoryManager* cudaMemoryManager)
+{
+ this->nTurbines = this->preInitDiameters.size();
+ this->numberOfNodes = nTurbines*nBladeNodes*nBlades;
+
+ cudaMemoryManager->cudaAllocBladeGeometries(this);
+ cudaMemoryManager->cudaAllocBladeOrientations(this);
+
+ for(uint turbine=0; turbine<this->nTurbines; turbine++)
+ {
+ for(uint node=0; node<this->nBladeNodes; node++)
+ {
+ this->bladeRadiiH[calcNode(node, nBladeNodes, 0, 1, turbine, nTurbines)] = this->preInitBladeRadii[turbine][node];
+ }
+
+ }
+ std::copy(preInitPosX.begin(), preInitPosX.end(), turbinePosXH);
+ std::copy(preInitPosY.begin(), preInitPosY.end(), turbinePosYH);
+ std::copy(preInitPosZ.begin(), preInitPosZ.end(), turbinePosZH);
+ std::copy(preInitDiameters.begin(), preInitDiameters.end(), diametersH);
+
+ cudaMemoryManager->cudaCopyBladeGeometriesHtoD(this);
+ std::copy(preInitAzimuths.begin(), preInitAzimuths.end(), this->azimuthsH);
+ std::copy(preInitOmegas.begin(), preInitOmegas.end(), this->omegasH);
+ std::copy(preInitYaws.begin(), preInitYaws.end(), this->yawsH);
+
+ cudaMemoryManager->cudaCopyBladeOrientationsHtoD(this);
+ real dxOPiSqrtEps = pow(this->deltaX/(this->epsilon*sqrt(cPi)),c3o1);
+ this->factorGaussian = dxOPiSqrtEps/this->forceRatio;
+}
+
+void ActuatorFarm::initBladeCoords(CudaMemoryManager* cudaMemoryManager)
+{
+ cudaMemoryManager->cudaAllocBladeCoords(this);
+
+ for(uint turbine=0; turbine<nTurbines; turbine++)
+ {
+ for(uint blade=0; blade<this->nBlades; blade++)
+ {
+ for(uint bladeNode=0; bladeNode<this->nBladeNodes; bladeNode++)
+ {
+ uint node = calcNode(bladeNode, this->nBladeNodes, blade, this->nBlades, turbine, this->nTurbines);
+
+ this->bladeCoordsXH[node] = c0o1;
+ this->bladeCoordsYH[node] = c0o1;
+ this->bladeCoordsZH[node] = this->bladeRadiiH[calcNode(bladeNode, nBladeNodes, 0, 1, turbine, nTurbines)];
+ }
+ }
+ }
+ cudaMemoryManager->cudaCopyBladeCoordsHtoD(this);
+}
+
+void ActuatorFarm::initBladeVelocities(CudaMemoryManager* cudaMemoryManager)
+{
+ cudaMemoryManager->cudaAllocBladeVelocities(this);
+
+ std::fill_n(this->bladeVelocitiesXH, this->numberOfNodes, c0o1);
+ std::fill_n(this->bladeVelocitiesYH, this->numberOfNodes, c0o1);
+ std::fill_n(this->bladeVelocitiesZH, this->numberOfNodes, c0o1);
+
+ cudaMemoryManager->cudaCopyBladeVelocitiesHtoD(this);
+}
+
+void ActuatorFarm::initBladeForces(CudaMemoryManager* cudaMemoryManager)
+{
+ cudaMemoryManager->cudaAllocBladeForces(this);
+
+ std::fill_n(this->bladeForcesXH, this->numberOfNodes, c0o1);
+ std::fill_n(this->bladeForcesYH, this->numberOfNodes, c0o1);
+ std::fill_n(this->bladeForcesZH, this->numberOfNodes, c0o1);
+
+ cudaMemoryManager->cudaCopyBladeForcesHtoD(this);
+}
+
+void ActuatorFarm::initBladeIndices(Parameter* para, CudaMemoryManager* cudaMemoryManager)
+{
+ cudaMemoryManager->cudaAllocBladeIndices(this);
+
+ std::fill_n(this->bladeIndicesH, this->numberOfNodes, 1);
+
+ cudaMemoryManager->cudaCopyBladeIndicesHtoD(this);
+}
+
+void ActuatorFarm::initBoundingSphere(Parameter* para, CudaMemoryManager* cudaMemoryManager)
+{
+ // Actuator line exists only on 1 level
+ std::vector<int> nodesInSphere;
+
+ for(uint turbine=0; turbine<this->nTurbines; turbine++)
+ {
+ real sphereRadius = c1o2*this->diametersH[turbine]+c4o1*this->epsilon;
+
+ real posX = this->turbinePosXH[level];
+ real posY = this->turbinePosYH[level];
+ real posZ = this->turbinePosZH[level];
+
+ real sphereRadiusSqrd = sphereRadius*sphereRadius;
+
+ for (uint j = 1; j <= para->getParH(this->level)->numberOfNodes; j++)
+ {
+ const real distX = para->getParH(this->level)->coordinateX[j]-posX;
+ const real distY = para->getParH(this->level)->coordinateY[j]-posY;
+ const real distZ = para->getParH(this->level)->coordinateZ[j]-posZ;
+ if(distSqrd(distX,distY,distZ) < sphereRadiusSqrd) nodesInSphere.push_back(j);
+ }
+ }
+
+ this->numberOfIndices = uint(nodesInSphere.size());
+
+ cudaMemoryManager->cudaAllocSphereIndices(this);
+ std::copy(nodesInSphere.begin(), nodesInSphere.end(), this->boundingSphereIndicesH);
+ cudaMemoryManager->cudaCopySphereIndicesHtoD(this);
+}
+
+void ActuatorFarm::setAllBladeCoords(real* _bladeCoordsX, real* _bladeCoordsY, real* _bladeCoordsZ)
+{
+ std::copy_n(this->bladeCoordsXH, this->numberOfNodes, _bladeCoordsX);
+ std::copy_n(this->bladeCoordsYH, this->numberOfNodes, _bladeCoordsY);
+ std::copy_n(this->bladeCoordsZH, this->numberOfNodes, _bladeCoordsZ);
+
+}
+
+void ActuatorFarm::setAllBladeVelocities(real* _bladeVelocitiesX, real* _bladeVelocitiesY, real* _bladeVelocitiesZ)
+{
+ std::copy_n(this->bladeVelocitiesXH, this->numberOfNodes, _bladeVelocitiesX);
+ std::copy_n(this->bladeVelocitiesYH, this->numberOfNodes, _bladeVelocitiesY);
+ std::copy_n(this->bladeVelocitiesZH, this->numberOfNodes, _bladeVelocitiesZ);
+
+}
+
+void ActuatorFarm::setAllBladeForces(real* _bladeForcesX, real* _bladeForcesY, real* _bladeForcesZ)
+{
+ std::copy_n(this->bladeForcesXH, this->numberOfNodes, _bladeForcesX);
+ std::copy_n(this->bladeForcesYH, this->numberOfNodes, _bladeForcesY);
+ std::copy_n(this->bladeForcesZH, this->numberOfNodes, _bladeForcesZ);
+
+}void ActuatorFarm::setTurbineBladeCoords(uint turbine, real* _bladeCoordsX, real* _bladeCoordsY, real* _bladeCoordsZ)
+{
+ std::copy_n(&this->bladeCoordsXH[turbine*nBladeNodes*nBlades], nBladeNodes*nBlades, _bladeCoordsX);
+ std::copy_n(&this->bladeCoordsYH[turbine*nBladeNodes*nBlades], nBladeNodes*nBlades, _bladeCoordsY);
+ std::copy_n(&this->bladeCoordsZH[turbine*nBladeNodes*nBlades], nBladeNodes*nBlades, _bladeCoordsZ);
+}
+
+void ActuatorFarm::setTurbineBladeVelocities(uint turbine, real* _bladeVelocitiesX, real* _bladeVelocitiesY, real* _bladeVelocitiesZ)
+{
+ std::copy_n(&this->bladeVelocitiesXH[turbine*nBladeNodes*nBlades], nBladeNodes*nBlades, _bladeVelocitiesX);
+ std::copy_n(&this->bladeVelocitiesYH[turbine*nBladeNodes*nBlades], nBladeNodes*nBlades, _bladeVelocitiesY);
+ std::copy_n(&this->bladeVelocitiesZH[turbine*nBladeNodes*nBlades], nBladeNodes*nBlades, _bladeVelocitiesZ);
+}
+
+void ActuatorFarm::setTurbineBladeForces(uint turbine, real* _bladeForcesX, real* _bladeForcesY, real* _bladeForcesZ)
+{
+ std::copy_n(&this->bladeForcesXH[turbine*nBladeNodes*nBlades], nBladeNodes*nBlades, _bladeForcesX);
+ std::copy_n(&this->bladeForcesYH[turbine*nBladeNodes*nBlades], nBladeNodes*nBlades, _bladeForcesY);
+ std::copy_n(&this->bladeForcesZH[turbine*nBladeNodes*nBlades], nBladeNodes*nBlades, _bladeForcesZ);
+}
diff --git a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.h b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.h
new file mode 100644
index 000000000..4e7fc220e
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.h
@@ -0,0 +1,168 @@
+#ifndef ActuatorFarm_H
+#define ActuatorFarm_H
+
+#include "PreCollisionInteractor.h"
+#include "PointerDefinitions.h"
+
+class Parameter;
+class GridProvider;
+
+class ActuatorFarm : public PreCollisionInteractor
+{
+public:
+ ActuatorFarm(
+ const uint _nBlades,
+ const real _density,
+ const uint _nBladeNodes,
+ const real _epsilon,
+ int _level,
+ const real _deltaT,
+ const real _deltaX,
+ const bool _useHostArrays
+ ) :
+ nBlades(_nBlades),
+ density(_density),
+ nBladeNodes(_nBladeNodes),
+ epsilon(_epsilon),
+ level(_level),
+ useHostArrays(_useHostArrays),
+ nTurbines(0),
+ numberOfNodes(0),
+ PreCollisionInteractor()
+ {
+ this->deltaT = _deltaX/pow(2,this->level);
+ this->deltaX = _deltaX/pow(2,this->level);
+ }
+
+ virtual ~ActuatorFarm()
+ {
+
+ }
+ void addTurbine(real turbinePosX, real turbinePosY, real turbinePosZ, real diameter, real omega, real azimuth, real yaw, std::vector<real> bladeRadii);
+ void init(Parameter* para, GridProvider* gridProvider, CudaMemoryManager* cudaManager);
+ void interact(Parameter* para, CudaMemoryManager* cudaManager, int level, uint t);
+ void free(Parameter* para, CudaMemoryManager* cudaManager);
+ void write(uint t);
+
+ uint getNTurbines(){return this->nTurbines; };
+ uint getNBladeNodes(){return this->nBladeNodes;};
+ uint getNBlades(){return this->nBlades;};
+
+ real* getAllAzimuths(){ return azimuthsH; };
+ real* getAllOmegas(){ return omegasH; };
+ real* getAllYaws(){ return yawsH; };
+
+ real* getAllTurbinePosX(){ return turbinePosXH; };
+ real* getAllTurbinePosY(){ return turbinePosYH; };
+ real* getAllTurbinePosZ(){ return turbinePosZH; };
+
+ real getTurbineAzimuth(uint turbine){ return azimuthsH[turbine]; };
+ real getTurbineOmega (uint turbine){ return omegasH[turbine]; };
+ real getTurbineYaw (uint turbine){ return yawsH[turbine]; };
+
+ real getTurbinePosX(uint turbine){ return turbinePosXH[turbine]; };
+ real getTurbinePosY(uint turbine){ return turbinePosYH[turbine]; };
+ real getTurbinePosZ(uint turbine){ return turbinePosZH[turbine]; };
+
+ uint getNumberOfIndices(){return this->numberOfIndices;};
+ uint getNumberOfNodes(){return this->numberOfNodes;};
+ real* getAllBladeCoordsX(){return this->bladeCoordsXH;};
+ real* getAllBladeCoordsY(){return this->bladeCoordsYH;};
+ real* getAllBladeCoordsZ(){return this->bladeCoordsZH;};
+ real* getAllBladeVelocitiesX(){return this->bladeVelocitiesXH;};
+ real* getAllBladeVelocitiesY(){return this->bladeVelocitiesYH;};
+ real* getAllBladeVelocitiesZ(){return this->bladeVelocitiesZH;};
+ real* getAllBladeForcesX(){return this->bladeForcesXH;};
+ real* getAllBladeForcesY(){return this->bladeForcesYH;};
+ real* getAllBladeForcesZ(){return this->bladeForcesZH;};
+
+ real* getTurbineBladeCoordsX(uint turbine){return &this->bladeCoordsXH[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeCoordsY(uint turbine){return &this->bladeCoordsYH[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeCoordsZ(uint turbine){return &this->bladeCoordsZH[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeVelocitiesX(uint turbine){return &this->bladeVelocitiesXH[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeVelocitiesY(uint turbine){return &this->bladeVelocitiesYH[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeVelocitiesZ(uint turbine){return &this->bladeVelocitiesZH[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeForcesX(uint turbine){return &this->bladeForcesXH[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeForcesY(uint turbine){return &this->bladeForcesYH[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeForcesZ(uint turbine){return &this->bladeForcesZH[turbine*nBladeNodes*nBlades];};
+
+ real* getAllBladeCoordsXDevice(){return this->bladeCoordsXD;};
+ real* getAllBladeCoordsYDevice(){return this->bladeCoordsYD;};
+ real* getAllBladeCoordsZDevice(){return this->bladeCoordsZD;};
+ real* getAllBladeVelocitiesXDevice(){return this->bladeVelocitiesXD;};
+ real* getAllBladeVelocitiesYDevice(){return this->bladeVelocitiesYD;};
+ real* getAllBladeVelocitiesZDevice(){return this->bladeVelocitiesZD;};
+ real* getAllBladeForcesXDevice(){return this->bladeForcesXD;};
+ real* getAllBladeForcesYDevice(){return this->bladeForcesYD;};
+ real* getAllBladeForcesZDevice(){return this->bladeForcesZD;};
+
+ real* getTurbineBladeCoordsXDevice(uint turbine){return &this->bladeCoordsXD[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeCoordsYDevice(uint turbine){return &this->bladeCoordsYD[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeCoordsZDevice(uint turbine){return &this->bladeCoordsZD[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeVelocitiesXDevice(uint turbine){return &this->bladeVelocitiesXD[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeVelocitiesYDevice(uint turbine){return &this->bladeVelocitiesYD[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeVelocitiesZDevice(uint turbine){return &this->bladeVelocitiesZD[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeForcesXDevice(uint turbine){return &this->bladeForcesXD[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeForcesYDevice(uint turbine){return &this->bladeForcesYD[turbine*nBladeNodes*nBlades];};
+ real* getTurbineBladeForcesZDevice(uint turbine){return &this->bladeForcesZD[turbine*nBladeNodes*nBlades];};
+
+ void setAllBladeCoords(real* _bladeCoordsX, real* _bladeCoordsY, real* _bladeCoordsZ);
+ void setAllBladeVelocities(real* _bladeVelocitiesX, real* _bladeVelocitiesY, real* _bladeVelocitiesZ);
+ void setAllBladeForces(real* _bladeForcesX, real* _bladeForcesY, real* _bladeForcesZ);
+
+ void setTurbineBladeCoords(uint turbine, real* _bladeCoordsX, real* _bladeCoordsY, real* _bladeCoordsZ);
+ void setTurbineBladeVelocities(uint turbine, real* _bladeVelocitiesX, real* _bladeVelocitiesY, real* _bladeVelocitiesZ);
+ void setTurbineBladeForces(uint turbine, real* _bladeForcesX, real* _bladeForcesY, real* _bladeForcesZ);
+
+
+private:
+ void initTurbineGeometries(CudaMemoryManager* cudaManager);
+ void initBoundingSphere(Parameter* para, CudaMemoryManager* cudaManager);
+ void initBladeCoords(CudaMemoryManager* cudaManager);
+ void initBladeVelocities(CudaMemoryManager* cudaManager);
+ void initBladeForces(CudaMemoryManager* cudaManager);
+ void initBladeIndices(Parameter* para, CudaMemoryManager* cudaManager);
+
+ void calcForcesEllipticWing();
+ void calcBladeForces();
+ void rotateBlades(real angle, uint turbineID);
+
+ void writeBladeCoords(uint t);
+ void writeBladeForces(uint t);
+ void writeBladeVelocities(uint t);
+
+
+
+
+public:
+ real* bladeRadiiH;
+ real* bladeRadiiD;
+ real* bladeCoordsXH, * bladeCoordsYH, * bladeCoordsZH;
+ real* bladeCoordsXD, * bladeCoordsYD, * bladeCoordsZD;
+ real* bladeVelocitiesXH, * bladeVelocitiesYH, * bladeVelocitiesZH;
+ real* bladeVelocitiesXD, * bladeVelocitiesYD, * bladeVelocitiesZD;
+ real* bladeForcesXH, * bladeForcesYH, * bladeForcesZH;
+ real* bladeForcesXD, * bladeForcesYD, * bladeForcesZD;
+ uint* bladeIndicesH;
+ uint* bladeIndicesD;
+ uint* boundingSphereIndicesH;
+ uint* boundingSphereIndicesD;
+ real* turbinePosXH, *turbinePosYH, *turbinePosZH, *omegasH, *azimuthsH, *yawsH, *diametersH;
+ real* turbinePosXD, *turbinePosYD, *turbinePosZD, *omegasD, *azimuthsD, *yawsD, *diametersD;
+
+private:
+ std::vector<real> preInitPosX, preInitPosY, preInitPosZ, preInitDiameters, preInitOmegas, preInitAzimuths, preInitYaws;
+ std::vector<std::vector<real>> preInitBladeRadii;
+ const bool useHostArrays;
+ const real density;
+ real deltaT, deltaX;
+ const uint nBladeNodes, nBlades;
+ uint nTurbines;
+ const real epsilon; // in m
+ const int level;
+ uint numberOfIndices;
+ uint numberOfNodes;
+ real forceRatio, factorGaussian, invEpsilonSqrd, invDeltaX;
+};
+
+#endif
\ No newline at end of file
--
GitLab