diff --git a/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarm.cu b/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarm.cu index 1676a883671ae7a410415e569a49fda74d9b4e0b..776ebc6dcfdfe6281dcd2b558552d5a6965a58dd 100644 --- a/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarm.cu +++ b/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarm.cu @@ -37,17 +37,17 @@ #include <cuda_runtime.h> #include <helper_cuda.h> -#include <logger/Logger.h> -#include <cuda_helper/CudaGrid.h> #include <basics/constants/NumericConstants.h> #include <basics/writer/WbWriterVtkXmlBinary.h> +#include <cuda_helper/CudaGrid.h> +#include <logger/Logger.h> +#include "DataStructureInitializer/GridProvider.h" +#include "GPU/CudaMemoryManager.h" #include "GPU/GeometryUtils.h" #include "LBM/GPUHelperFunctions/KernelUtilities.h" -#include "Parameter/Parameter.h" #include "Parameter/CudaStreamManager.h" -#include "DataStructureInitializer/GridProvider.h" -#include "GPU/CudaMemoryManager.h" +#include "Parameter/Parameter.h" using namespace vf::basics::constant; @@ -55,16 +55,16 @@ struct GridData { const uint* indices; const uint nIndices; - const real* coordsX, * coordsY, *coordsZ; - const uint* neighborsX, * neighborsY, * neighborsZ, * neighborsWSB; - const real* vx, *vy, *vz; - real* fx, *fy, *fz; + const real *coordsX, *coordsY, *coordsZ; + const uint *neighborsX, *neighborsY, *neighborsZ, *neighborsWSB; + const real *vx, *vy, *vz; + real *fx, *fy, *fz; const real inverseDeltaX, velocityRatio; }; struct TurbineData { - const real* posX, *posY, *posZ; + const real *posX, *posY, *posZ; const uint numberOfTurbines; const real smearingWidth, factorGaussian; }; @@ -73,18 +73,18 @@ struct ComponentData { const real referenceLength; const uint numberOfNodesPerTurbine; - const real* coordsX, * coordsY, * coordsZ; - real* velocitiesX, * velocitiesY, * velocitiesZ; - const real* forcesX, * forcesY, * forcesZ; + const real *coordsX, *coordsY, *coordsZ; + real *velocitiesX, *velocitiesY, *velocitiesZ; + const real *forcesX, *forcesY, *forcesZ; uint* gridIndices; }; - __global__ void interpolateVelocities(const GridData gridData, const TurbineData turbineData, ComponentData componentData) { const unsigned nodeIndex = vf::gpu::getNodeIndex(); - if (nodeIndex >= componentData.numberOfNodesPerTurbine * turbineData.numberOfTurbines) return; + if (nodeIndex >= componentData.numberOfNodesPerTurbine * turbineData.numberOfTurbines) + return; const real coordX = componentData.coordsX[nodeIndex]; const real coordY = componentData.coordsY[nodeIndex]; @@ -100,7 +100,8 @@ __global__ void interpolateVelocities(const GridData gridData, const TurbineData componentData.gridIndices[nodeIndex] = k; - getNeighborIndicesOfBSW(k, ke, kn, kt, kne, kte, ktn, ktne, gridData.neighborsX, gridData.neighborsY, gridData.neighborsZ); + getNeighborIndicesOfBSW(k, ke, kn, kt, kne, kte, ktn, ktne, gridData.neighborsX, gridData.neighborsY, + gridData.neighborsZ); real dW, dE, dN, dS, dT, dB; @@ -110,9 +111,15 @@ __global__ void interpolateVelocities(const GridData gridData, const TurbineData getInterpolationWeights(dW, dE, dN, dS, dT, dB, distX, distY, distZ); - componentData.velocitiesX[nodeIndex] = trilinearInterpolation(dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, gridData.vx) * gridData.velocityRatio; - componentData.velocitiesY[nodeIndex] = trilinearInterpolation(dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, gridData.vy) * gridData.velocityRatio; - componentData.velocitiesZ[nodeIndex] = trilinearInterpolation(dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, gridData.vz) * gridData.velocityRatio; + componentData.velocitiesX[nodeIndex] = + trilinearInterpolation(dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, gridData.vx) * + gridData.velocityRatio; + componentData.velocitiesY[nodeIndex] = + trilinearInterpolation(dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, gridData.vy) * + gridData.velocityRatio; + componentData.velocitiesZ[nodeIndex] = + trilinearInterpolation(dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, gridData.vz) * + gridData.velocityRatio; } __global__ void applyBodyForces(GridData gridData, const TurbineData turbineData, const ComponentData componentData) @@ -120,7 +127,8 @@ __global__ void applyBodyForces(GridData gridData, const TurbineData turbineData const uint index = vf::gpu::getNodeIndex(); - if (index >= gridData.nIndices) return; + if (index >= gridData.nIndices) + return; const uint gridIndex = gridData.indices[index]; const real gridCoordX = gridData.coordsX[gridIndex]; @@ -136,10 +144,11 @@ __global__ void applyBodyForces(GridData gridData, const TurbineData turbineData const real distToHubY = gridCoordY - turbineData.posY[turbine]; const real distToHubZ = gridCoordZ - turbineData.posZ[turbine]; - if (!inBoundingSphere(distToHubX, distToHubY, distToHubZ, componentData.referenceLength, turbineData.smearingWidth)) continue; + if (!inBoundingSphere(distToHubX, distToHubY, distToHubZ, componentData.referenceLength, turbineData.smearingWidth)) + continue; for (uint turbineNode = 0; turbineNode < componentData.numberOfNodesPerTurbine; turbineNode++) { - const uint node = turbine*componentData.numberOfNodesPerTurbine + turbineNode; + const uint node = turbine * componentData.numberOfNodesPerTurbine + turbineNode; const real distX = componentData.coordsX[node] - gridCoordX; const real distY = componentData.coordsY[node] - gridCoordY; @@ -156,9 +165,10 @@ __global__ void applyBodyForces(GridData gridData, const TurbineData turbineData gridData.fz[gridIndex] += gridForceZ; } -void ActuatorFarm::init(Parameter *para, GridProvider *gridProvider, CudaMemoryManager *cudaManager) +void ActuatorFarm::init(Parameter* para, GridProvider* gridProvider, CudaMemoryManager* cudaManager) { - if (!para->getIsBodyForce()) throw std::runtime_error("try to allocate ActuatorFarm but BodyForce is not set in Parameter."); + if (!para->getIsBodyForce()) + throw std::runtime_error("try to allocate ActuatorFarm but BodyForce is not set in Parameter."); this->forceRatio = para->getForceRatio(); this->initTurbineGeometries(cudaManager); this->initBladeCoords(cudaManager); @@ -169,17 +179,18 @@ void ActuatorFarm::init(Parameter *para, GridProvider *gridProvider, CudaMemoryM this->streamIndex = para->getStreamManager()->registerAndLaunchStream(CudaStreamIndex::ActuatorFarm); } -void ActuatorFarm::interact(Parameter *para, CudaMemoryManager *cudaManager, int level, unsigned int t) +void ActuatorFarm::interact(Parameter* para, CudaMemoryManager* cudaManager, int level, unsigned int t) { - if (level != this->level) return; + if (level != this->level) + return; cudaStream_t stream = para->getStreamManager()->getStream(CudaStreamIndex::ActuatorFarm, this->streamIndex); - if (useHostArrays) cudaManager->cudaCopyBladeCoordsHtoD(this); + if (useHostArrays) + cudaManager->cudaCopyBladeCoordsHtoD(this); - if(this->writeOutput && ((t-this->tStartOut) % this->tOut == 0)) - { - if(!useHostArrays){ + if (this->writeOutput && ((t - this->tStartOut) % this->tOut == 0)) { + if (!useHostArrays) { cudaManager->cudaCopyBladeCoordsDtoH(this); cudaManager->cudaCopyBladeVelocitiesDtoH(this); cudaManager->cudaCopyBladeForcesDtoH(this); @@ -208,15 +219,17 @@ void ActuatorFarm::interact(Parameter *para, CudaMemoryManager *cudaManager, int this->bladeIndicesD}; vf::cuda::CudaGrid bladeGrid = vf::cuda::CudaGrid(para->getParH(level)->numberofthreads, this->numberOfGridNodes); - interpolateVelocities<<< bladeGrid.grid, bladeGrid.threads, 0, stream >>>(gridData, turbineData, bladeData); + interpolateVelocities<<<bladeGrid.grid, bladeGrid.threads, 0, stream>>>(gridData, turbineData, bladeData); cudaStreamSynchronize(stream); - if (useHostArrays) cudaManager->cudaCopyBladeVelocitiesDtoH(this); - + if (useHostArrays) + cudaManager->cudaCopyBladeVelocitiesDtoH(this); + this->updateForcesAndCoordinates(); this->swapDeviceArrays(); - if (useHostArrays) cudaManager->cudaCopyBladeForcesHtoD(this); + if (useHostArrays) + cudaManager->cudaCopyBladeForcesHtoD(this); vf::cuda::CudaGrid sphereGrid = vf::cuda::CudaGrid(para->getParH(level)->numberofthreads, this->numberOfIndices); @@ -224,7 +237,7 @@ void ActuatorFarm::interact(Parameter *para, CudaMemoryManager *cudaManager, int cudaStreamSynchronize(stream); } -void ActuatorFarm::free(Parameter *para, CudaMemoryManager *cudaManager) +void ActuatorFarm::free(Parameter* para, CudaMemoryManager* cudaManager) { cudaManager->cudaFreeBladeGeometries(this); cudaManager->cudaFreeBladeCoords(this); @@ -234,13 +247,13 @@ void ActuatorFarm::free(Parameter *para, CudaMemoryManager *cudaManager) cudaManager->cudaFreeSphereIndices(this); } -void ActuatorFarm::getTaggedFluidNodes(Parameter *para, GridProvider *gridProvider) +void ActuatorFarm::getTaggedFluidNodes(Parameter* para, GridProvider* gridProvider) { std::vector<uint> indicesInSphere(this->boundingSphereIndicesH, this->boundingSphereIndicesH + this->numberOfIndices); gridProvider->tagFluidNodeIndices(indicesInSphere, CollisionTemplate::AllFeatures, this->level); } -void ActuatorFarm::initTurbineGeometries(CudaMemoryManager *cudaManager) +void ActuatorFarm::initTurbineGeometries(CudaMemoryManager* cudaManager) { this->numberOfGridNodes = this->numberOfTurbines * this->numberOfNodesPerTurbine; @@ -255,19 +268,18 @@ void ActuatorFarm::initTurbineGeometries(CudaMemoryManager *cudaManager) this->factorGaussian = pow(this->smearingWidth * sqrt(cPi), -c3o1) / this->forceRatio; } -void ActuatorFarm::initBladeCoords(CudaMemoryManager *cudaManager) +void ActuatorFarm::initBladeCoords(CudaMemoryManager* cudaManager) { cudaManager->cudaAllocBladeCoords(this); - for (uint turbine = 0; turbine < this->numberOfTurbines; turbine++) - { - for (uint blade = 0; blade < ActuatorFarm::numberOfBlades; blade++) - { + for (uint turbine = 0; turbine < this->numberOfTurbines; turbine++) { + for (uint blade = 0; blade < ActuatorFarm::numberOfBlades; blade++) { const real local_azimuth = this->azimuths[turbine] + blade * c2Pi / ActuatorFarm::numberOfBlades; for (uint bladeNode = 0; bladeNode < this->numberOfNodesPerBlade; bladeNode++) { - const uint node = calcNodeIndexInBladeArrays({ turbine, blade, bladeNode }, this->numberOfNodesPerBlade, ActuatorFarm::numberOfBlades); - + const uint node = calcNodeIndexInBladeArrays({ turbine, blade, bladeNode }, this->numberOfNodesPerBlade, + ActuatorFarm::numberOfBlades); + real x, y, z; rotateFromBladeToGlobal(c0o1, c0o1, this->bladeRadii[bladeNode], x, y, z, local_azimuth); bladeCoordsXH[node] = x + this->turbinePosXH[turbine]; @@ -283,7 +295,7 @@ void ActuatorFarm::initBladeCoords(CudaMemoryManager *cudaManager) cudaManager->cudaCopyBladeCoordsHtoD(this); } -void ActuatorFarm::initBladeVelocities(CudaMemoryManager *cudaManager) +void ActuatorFarm::initBladeVelocities(CudaMemoryManager* cudaManager) { cudaManager->cudaAllocBladeVelocities(this); @@ -298,7 +310,7 @@ void ActuatorFarm::initBladeVelocities(CudaMemoryManager *cudaManager) cudaManager->cudaCopyBladeVelocitiesHtoD(this); } -void ActuatorFarm::initBladeForces(CudaMemoryManager *cudaManager) +void ActuatorFarm::initBladeForces(CudaMemoryManager* cudaManager) { cudaManager->cudaAllocBladeForces(this); @@ -313,7 +325,7 @@ void ActuatorFarm::initBladeForces(CudaMemoryManager *cudaManager) cudaManager->cudaCopyBladeForcesHtoD(this); } -void ActuatorFarm::initBladeIndices(Parameter *para, CudaMemoryManager *cudaManager) +void ActuatorFarm::initBladeIndices(Parameter* para, CudaMemoryManager* cudaManager) { cudaManager->cudaAllocBladeIndices(this); @@ -322,12 +334,13 @@ void ActuatorFarm::initBladeIndices(Parameter *para, CudaMemoryManager *cudaMana cudaManager->cudaCopyBladeIndicesHtoD(this); } -void ActuatorFarm::initBoundingSpheres(Parameter *para, CudaMemoryManager *cudaManager) +void ActuatorFarm::initBoundingSpheres(Parameter* para, CudaMemoryManager* cudaManager) { std::vector<int> nodesInSpheres; const real sphereRadius = getBoundingSphereRadius(this->diameter, this->smearingWidth); const real sphereRadiusSqrd = sphereRadius * sphereRadius; - const uint minimumNumberOfNodesPerSphere = (uint)(c4o3 * cPi * pow(sphereRadius - this->deltaX, c3o1) / pow(this->deltaX, c3o1)); + const uint minimumNumberOfNodesPerSphere = + (uint)(c4o3 * cPi * pow(sphereRadius - this->deltaX, c3o1) / pow(this->deltaX, c3o1)); for (uint turbine = 0; turbine < this->numberOfTurbines; turbine++) { @@ -348,7 +361,8 @@ void ActuatorFarm::initBoundingSpheres(Parameter *para, CudaMemoryManager *cudaM } if (nodesInThisSphere < minimumNumberOfNodesPerSphere) { - VF_LOG_CRITICAL("Found only {} nodes in bounding sphere of turbine no. {}, expected at least {}!", nodesInThisSphere, turbine, minimumNumberOfNodesPerSphere); + VF_LOG_CRITICAL("Found only {} nodes in bounding sphere of turbine no. {}, expected at least {}!", + nodesInThisSphere, turbine, minimumNumberOfNodesPerSphere); throw std::runtime_error("ActuatorFarm::initBoundingSpheres: Turbine bounding sphere partially out of domain."); } } @@ -360,41 +374,46 @@ void ActuatorFarm::initBoundingSpheres(Parameter *para, CudaMemoryManager *cudaM cudaManager->cudaCopySphereIndicesHtoD(this); } -void ActuatorFarm::setAllBladeCoords(const real * _bladeCoordsX, const real * _bladeCoordsY, const real * _bladeCoordsZ) +void ActuatorFarm::setAllBladeCoords(const real* _bladeCoordsX, const real* _bladeCoordsY, const real* _bladeCoordsZ) { std::copy_n(_bladeCoordsX, this->numberOfGridNodes, this->bladeCoordsXH); std::copy_n(_bladeCoordsY, this->numberOfGridNodes, this->bladeCoordsYH); std::copy_n(_bladeCoordsZ, this->numberOfGridNodes, this->bladeCoordsZH); } -void ActuatorFarm::setAllBladeVelocities(const real *_bladeVelocitiesX, const real *_bladeVelocitiesY, const real *_bladeVelocitiesZ) +void ActuatorFarm::setAllBladeVelocities(const real* _bladeVelocitiesX, const real* _bladeVelocitiesY, + const real* _bladeVelocitiesZ) { std::copy_n(_bladeVelocitiesX, this->numberOfGridNodes, this->bladeVelocitiesXH); std::copy_n(_bladeVelocitiesY, this->numberOfGridNodes, this->bladeVelocitiesYH); std::copy_n(_bladeVelocitiesZ, this->numberOfGridNodes, this->bladeVelocitiesZH); } -void ActuatorFarm::setAllBladeForces(const real *_bladeForcesX, const real *_bladeForcesY, const real *_bladeForcesZ) +void ActuatorFarm::setAllBladeForces(const real* _bladeForcesX, const real* _bladeForcesY, const real* _bladeForcesZ) { std::copy_n(_bladeForcesX, this->numberOfGridNodes, this->bladeForcesXH); std::copy_n(_bladeForcesY, this->numberOfGridNodes, this->bladeForcesYH); std::copy_n(_bladeForcesZ, this->numberOfGridNodes, this->bladeForcesZH); } -void ActuatorFarm::setTurbineBladeCoords(uint turbine, const real *_bladeCoordsX, const real *_bladeCoordsY, const real *_bladeCoordsZ) + +void ActuatorFarm::setTurbineBladeCoords(uint turbine, const real* _bladeCoordsX, const real* _bladeCoordsY, + const real* _bladeCoordsZ) { std::copy_n(_bladeCoordsX, this->numberOfNodesPerTurbine, &this->bladeCoordsXH[turbine * this->numberOfNodesPerTurbine]); std::copy_n(_bladeCoordsY, this->numberOfNodesPerTurbine, &this->bladeCoordsYH[turbine * this->numberOfNodesPerTurbine]); std::copy_n(_bladeCoordsZ, this->numberOfNodesPerTurbine, &this->bladeCoordsZH[turbine * this->numberOfNodesPerTurbine]); } -void ActuatorFarm::setTurbineBladeVelocities(uint turbine, const real *_bladeVelocitiesX, const real *_bladeVelocitiesY, const real *_bladeVelocitiesZ) +void ActuatorFarm::setTurbineBladeVelocities(uint turbine, const real* _bladeVelocitiesX, const real* _bladeVelocitiesY, + const real* _bladeVelocitiesZ) { std::copy_n(_bladeVelocitiesX, this->numberOfNodesPerTurbine, &this->bladeVelocitiesXH[turbine * this->numberOfNodesPerTurbine]); std::copy_n(_bladeVelocitiesY, this->numberOfNodesPerTurbine, &this->bladeVelocitiesYH[turbine * this->numberOfNodesPerTurbine]); std::copy_n(_bladeVelocitiesZ, this->numberOfNodesPerTurbine, &this->bladeVelocitiesZH[turbine * this->numberOfNodesPerTurbine]); } -void ActuatorFarm::setTurbineBladeForces(uint turbine, const real *_bladeForcesX, const real *_bladeForcesY, const real *_bladeForcesZ) +void ActuatorFarm::setTurbineBladeForces(uint turbine, const real* _bladeForcesX, const real* _bladeForcesY, + const real* _bladeForcesZ) { std::copy_n(_bladeForcesX, this->numberOfNodesPerTurbine, &this->bladeForcesXH[turbine * this->numberOfNodesPerTurbine]); std::copy_n(_bladeForcesY, this->numberOfNodesPerTurbine, &this->bladeForcesYH[turbine * this->numberOfNodesPerTurbine]); @@ -433,8 +452,9 @@ void ActuatorFarm::write(const std::string &filename) }; std::vector<UbTupleFloat3> nodes(numberOfGridNodes); std::vector<std::vector<double>> nodeData(6); - for(auto& data : nodeData) data.resize(numberOfGridNodes); - for(uint i=0; i<numberOfGridNodes; i++){ + for (auto& data : nodeData) + data.resize(numberOfGridNodes); + for (uint i = 0; i < numberOfGridNodes; i++) { nodes[i] = UbTupleFloat3(this->bladeCoordsXH[i], this->bladeCoordsYH[i], this->bladeCoordsZH[i]); nodeData[0][i] = this->bladeVelocitiesXH[i]; nodeData[1][i] = this->bladeVelocitiesYH[i]; diff --git a/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarmStandalone.cu b/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarmStandalone.cu index ef2c2eca0a7d8e772a64c9e8930c115355a85146..6e577d1f06967ccf0e80db9474b9b7c3b8861c50 100644 --- a/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarmStandalone.cu +++ b/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarmStandalone.cu @@ -51,12 +51,12 @@ using namespace vf::basics::constant; - -std::vector<real> ActuatorFarmStandalone::computeBladeRadii(const real diameter, const uint numberOfNodesPerBlade){ +std::vector<real> ActuatorFarmStandalone::computeBladeRadii(const real diameter, const uint numberOfNodesPerBlade) +{ const real dr = c1o2 * diameter / numberOfNodesPerBlade; - std::vector<real>blade_radii(numberOfNodesPerBlade); - for(uint node=0; node<numberOfNodesPerBlade; node++) - blade_radii[node] = dr * (node+0.5); + std::vector<real> blade_radii(numberOfNodesPerBlade); + for (uint node = 0; node < numberOfNodesPerBlade; node++) + blade_radii[node] = dr * (node + 0.5); return blade_radii; } @@ -67,24 +67,22 @@ void ActuatorFarmStandalone::updateForcesAndCoordinates() const real c0 = 20 * c1o10; const real delta_azimuth = c2Pi / this->numberOfBlades; - for (uint turbine = 0; turbine < this->numberOfTurbines; turbine++) - { + for (uint turbine = 0; turbine < this->numberOfTurbines; turbine++) { const real rotor_speed = this->rotorSpeeds[turbine]; const real azimuth_old = this->azimuths[turbine]; const real azimuth_new = azimuth_old + deltaT * rotor_speed; this->azimuths[turbine] = azimuth_new > c2Pi ? azimuth_new - c2Pi : azimuth_new; - for (uint blade = 0; blade < this->numberOfBlades; blade++) - { - const real local_azimuth_new = azimuth_new + blade*delta_azimuth; + for (uint blade = 0; blade < this->numberOfBlades; blade++) { + const real local_azimuth_new = azimuth_new + blade * delta_azimuth; real last_node_radius = c0o1; real current_node_radius = c0o1; real next_node_radius = this->bladeRadii[0]; - for (uint bladeNode = 0; bladeNode < this->numberOfNodesPerBlade; bladeNode++) - { - const uint node = calcNodeIndexInBladeArrays({ turbine, blade, bladeNode }, this->numberOfNodesPerBlade, this->numberOfBlades); + for (uint bladeNode = 0; bladeNode < this->numberOfNodesPerBlade; bladeNode++) { + const uint node = calcNodeIndexInBladeArrays({ turbine, blade, bladeNode }, this->numberOfNodesPerBlade, + this->numberOfBlades); real u_rel, v_rel, w_rel; rotateFromGlobalToBlade(u_rel, v_rel, w_rel, @@ -92,10 +90,11 @@ void ActuatorFarmStandalone::updateForcesAndCoordinates() this->bladeVelocitiesYH[node], this->bladeVelocitiesZH[node], azimuth_old + delta_azimuth); - + last_node_radius = current_node_radius; current_node_radius = next_node_radius; - next_node_radius = bladeNode < this->numberOfNodesPerBlade-1 ? this->bladeRadii[bladeNode+1] : this->diameter * c1o2; + next_node_radius = + bladeNode < this->numberOfNodesPerBlade - 1 ? this->bladeRadii[bladeNode + 1] : this->diameter * c1o2; const real dr = c1o2 * (next_node_radius - last_node_radius); @@ -123,4 +122,3 @@ void ActuatorFarmStandalone::updateForcesAndCoordinates() } } } -