diff --git a/pythonbindings/pyfluids-stubs/gpu/__init__.pyi b/pythonbindings/pyfluids-stubs/gpu/__init__.pyi
index 8f78fd3387dc59c94ddf2315635e674a7330dae9..15386f7b0df4d2d06e01f8621805af8790c32442 100644
--- a/pythonbindings/pyfluids-stubs/gpu/__init__.pyi
+++ b/pythonbindings/pyfluids-stubs/gpu/__init__.pyi
@@ -38,10 +38,11 @@ from typing import Callable, ClassVar, List, Optional
from typing import overload, Union
import numpy as np
import numpy.typing as npt
-import basics
+import basics, parallel
from . import grid_generator as grid_generator
from . import probes as probes
+from . import Kernel as Kernel
class PreCollisionInteractor:
def __init__(self, *args, **kwargs) -> None: ...
@@ -51,10 +52,8 @@ class FileCollection:
def __init__(self, *args, **kwargs) -> None: ...
class ActuatorFarm(PreCollisionInteractor):
- def __init__(self, number_of_blades_per_turbine: int, density: float, number_of_nodes_per_blade: int, epsilon: float, level: int, delta_t: float, delta_x: float, use_host_arrays: bool) -> None: ...
- def add_turbine(self, posX: float, posY: float, posZ: float, diameter: float, omega: float, azimuth: float, yaw: float, bladeRadii: List[float]) -> None: ...
- def calc_blade_forces(self) -> None: ...
- def get_all_azimuths(self) -> npt.NDArray[np.float32]: ...
+ def __init__(self, diameter: float, blade_radii: npt.NDArray[np.float32],turbine_positions_x: npt.NDArray[np.float32], turbine_positions_y: npt.NDArray[np.float32], turbine_positions_z: npt.NDArray[np.float32], density: float, smearing_width: float, level: int, delta_t: float, delta_x: float, use_host_arrays: bool) -> None: ...
+ def update_forces_and_coordinates(self) -> None: ...
def get_all_blade_coords_x(self) -> npt.NDArray[np.float32]: ...
def get_all_blade_coords_x_device(self) -> int: ...
def get_all_blade_coords_y(self) -> npt.NDArray[np.float32]: ...
@@ -67,20 +66,15 @@ class ActuatorFarm(PreCollisionInteractor):
def get_all_blade_forces_y_device(self) -> int: ...
def get_all_blade_forces_z(self) -> npt.NDArray[np.float32]: ...
def get_all_blade_forces_z_device(self) -> int: ...
- def get_all_blade_radii(self) -> npt.NDArray[np.float32]: ...
- def get_all_blade_radii_device(self) -> int: ...
def get_all_blade_velocities_x(self) -> npt.NDArray[np.float32]: ...
def get_all_blade_velocities_x_device(self) -> int: ...
def get_all_blade_velocities_y(self) -> npt.NDArray[np.float32]: ...
def get_all_blade_velocities_y_device(self) -> int: ...
def get_all_blade_velocities_z(self) -> npt.NDArray[np.float32]: ...
def get_all_blade_velocities_z_device(self) -> int: ...
- def get_all_omegas(self) -> npt.NDArray[np.float32]: ...
def get_all_turbine_pos_x(self) -> npt.NDArray[np.float32]: ...
def get_all_turbine_pos_y(self) -> npt.NDArray[np.float32]: ...
def get_all_turbine_pos_z(self) -> npt.NDArray[np.float32]: ...
- def get_all_yaws(self) -> npt.NDArray[np.float32]: ...
- def get_turbine_azimuth(self, turbine: int) -> float: ...
def get_turbine_blade_coords_x(self, turbine: int) -> npt.NDArray[np.float32]: ...
def get_turbine_blade_coords_x_device(self, turbine: int) -> int: ...
def get_turbine_blade_coords_y(self, turbine: int) -> npt.NDArray[np.float32]: ...
@@ -93,29 +87,20 @@ class ActuatorFarm(PreCollisionInteractor):
def get_turbine_blade_forces_y_device(self, turbine: int) -> int: ...
def get_turbine_blade_forces_z(self, turbine: int) -> npt.NDArray[np.float32]: ...
def get_turbine_blade_forces_z_device(self, turbine: int) -> int: ...
- def get_turbine_blade_radii(self, turbine: int) -> npt.NDArray[np.float32]: ...
- def get_turbine_blade_radii_device(self, turbine: int) -> int: ...
def get_turbine_blade_velocities_x(self, turbine: int) -> npt.NDArray[np.float32]: ...
def get_turbine_blade_velocities_x_device(self, turbine: int) -> int: ...
def get_turbine_blade_velocities_y(self, turbine: int) -> npt.NDArray[np.float32]: ...
def get_turbine_blade_velocities_y_device(self, turbine: int) -> int: ...
def get_turbine_blade_velocities_z(self, turbine: int) -> npt.NDArray[np.float32]: ...
def get_turbine_blade_velocities_z_device(self, turbine: int) -> int: ...
- def get_turbine_omega(self, turbine: int) -> float: ...
- def get_turbine_pos(self, turbine: int) -> npt.NDArray[np.float32]: ...
- def get_turbine_yaw(self, turbine: int) -> float: ...
- def set_all_azimuths(self, azimuths: npt.NDArray[np.float32]) -> None: ...
def set_all_blade_coords(self, blade_coords_x: npt.NDArray[np.float32], blade_coords_y: npt.NDArray[np.float32], blade_coords_z: npt.NDArray[np.float32]) -> None: ...
def set_all_blade_forces(self, blade_forces_x: npt.NDArray[np.float32], blade_forces_y: npt.NDArray[np.float32], blade_forces_z: npt.NDArray[np.float32]) -> None: ...
def set_all_blade_velocities(self, blade_velocities_x: npt.NDArray[np.float32], blade_velocities_y: npt.NDArray[np.float32], blade_velocities_z: npt.NDArray[np.float32]) -> None: ...
- def set_all_omegas(self, omegas: npt.NDArray[np.float32]) -> None: ...
- def set_all_yaws(self, yaws: npt.NDArray[np.float32]) -> None: ...
- def set_turbine_azimuth(self, turbine: int, azimuth: float) -> None: ...
def set_turbine_blade_coords(self, turbine: int, blade_coords_x: npt.NDArray[np.float32], blade_coords_y: npt.NDArray[np.float32], blade_coords_z: npt.NDArray[np.float32]) -> None: ...
def set_turbine_blade_forces(self, turbine: int, blade_forces_x: npt.NDArray[np.float32], blade_forces_y: npt.NDArray[np.float32], blade_forces_z: npt.NDArray[np.float32]) -> None: ...
def set_turbine_blade_velocities(self, turbine: int, blade_velocities_x: npt.NDArray[np.float32], blade_velocities_y: npt.NDArray[np.float32], blade_velocities_z: npt.NDArray[np.float32]) -> None: ...
- def set_turbine_omega(self, turbine: int, omega: float) -> None: ...
- def set_turbine_yaw(self, turbine: int, yaw: float) -> None: ...
+ def set_turbine_azimuth(self, turbine: int, azimuth: float) -> None: ...
+ def enable_output(self, output_name: str, t_start_out: int, t_out: int) -> None: ...
@property
def delta_t(self) -> float: ...
@property
@@ -133,6 +118,8 @@ class ActuatorFarm(PreCollisionInteractor):
@property
def number_of_turbines(self) -> int: ...
+class ActuatorFarmStandalone(ActuatorFarm):
+ def __init__(self, diameter: float, number_of_nodes_per_blade: int, turbine_positions_x: npt.NDArray[np.float32], turbine_positions_y: npt.NDArray[np.float32], turbine_positions_z: npt.NDArray[np.float32], rotor_speeds: npt.NDArray[np.float32], density: float, smearing_width: float, level: int, delta_t: float, delta_x: float) -> None: ...
class BoundaryConditionFactory:
def __init__(self) -> None: ...
@@ -145,14 +132,6 @@ class BoundaryConditionFactory:
def set_velocity_boundary_condition(self, boundary_condition_type: VelocityBC) -> None: ...
-class MpiCommunicator:
- def __init__(self, *args, **kwargs) -> None: ...
- @staticmethod
- def get_instance() -> MpiCommunicator: ...
- def get_number_of_process(self) -> int: ...
- def get_pid(self) -> int: ...
-
-
class CudaMemoryManager:
def __init__(self, parameter: Parameter) -> None: ...
@@ -176,7 +155,7 @@ class FileType:
class GridProvider:
def __init__(self, *args, **kwargs) -> None: ...
@staticmethod
- def make_grid_generator(builder: grid_generator.GridBuilder, para: Parameter, cuda_memory_manager: CudaMemoryManager, communicator: MpiCommunicator) -> GridProvider: ...
+ def make_grid_generator(builder: grid_generator.GridBuilder, para: Parameter, cuda_memory_manager: CudaMemoryManager, communicator: parallel.Communicator) -> GridProvider: ...
class MultipleGridBuilder:
def __init__(self) -> None: ...
@@ -273,7 +252,7 @@ class Parameter:
def set_initial_condition_perturbed_log_law(self, u_star: float, z0: float, length_x: float, length_z: float, height: float, velocity_ratio: float) -> None: ...
def set_initial_condition_uniform(self, velocity_x: float, velocity_y: float, velocity_z: float) -> None: ...
def set_is_body_force(self, is_body_force: bool) -> None: ...
- def set_main_kernel(self, kernel: str) -> None: ...
+ def configure_main_kernel(self, kernel: str) -> None: ...
def set_max_dev(self, max_dev: int) -> None: ...
def set_max_level(self, number_of_levels: int) -> None: ...
def set_outflow_pressure_correction_factor(self, correction_factor: float) -> None: ...
diff --git a/pythonbindings/src/gpu/submodules/actuator_farm.cpp b/pythonbindings/src/gpu/submodules/actuator_farm.cpp
index ee5275483f85559a15aefa769ba30a87af993355..c057d6bb15cc750adf8f6266fba7d27aede03b61 100644
--- a/pythonbindings/src/gpu/submodules/actuator_farm.cpp
+++ b/pythonbindings/src/gpu/submodules/actuator_farm.cpp
@@ -33,17 +33,18 @@
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
-#include <gpu/core/PreCollisionInteractor/ActuatorFarm.h>
-#include <gpu/core/PreCollisionInteractor/PreCollisionInteractor.h>
+#include <gpu/VirtualFluids_GPU/PreCollisionInteractor/Actuator/ActuatorFarm.h>
+#include <gpu/VirtualFluids_GPU/PreCollisionInteractor/Actuator/ActuatorFarmStandalone.h>
+#include <gpu/VirtualFluids_GPU/PreCollisionInteractor/PreCollisionInteractor.h>
class PyActuatorFarm : public ActuatorFarm
{
public:
using ActuatorFarm::ActuatorFarm; // Inherit constructors
- void calcBladeForces() override
+ void updateForcesAndCoordinates() override
{
- PYBIND11_OVERRIDE_NAME(void, ActuatorFarm, "calc_blade_forces", calcBladeForces);
+ PYBIND11_OVERRIDE_PURE_NAME(void, ActuatorFarm, "update_forces_and_coordinates", updateForcesAndCoordinates);
}
};
@@ -62,24 +63,31 @@ namespace actuator_farm
using arr = py::array_t<real, py::array::c_style>;
py::class_<ActuatorFarm, PreCollisionInteractor, PyActuatorFarm, std::shared_ptr<ActuatorFarm>>(parentModule, "ActuatorFarm", py::dynamic_attr())
- .def(py::init< const uint,
+ .def(py::init< const real,
+ const std::vector<real>,
+ const std::vector<real>,
+ const std::vector<real>,
+ const std::vector<real>,
const real,
- const uint,
const real,
- int,
+ const int,
const real,
const real,
const bool>(),
- py::arg("number_of_blades_per_turbine"),
+ py::arg("diameter"),
+ py::arg("blade_radii"),
+ py::arg("turbine_positions_x"),
+ py::arg("turbine_positions_y"),
+ py::arg("turbine_positions_z"),
py::arg("density"),
- py::arg("number_of_nodes_per_blade"),
- py::arg("epsilon"),
+ py::arg("smearing_width"),
py::arg("level"),
py::arg("delta_t"),
py::arg("delta_x"),
py::arg("use_host_arrays"))
- .def_property_readonly("number_of_turbines", &ActuatorFarm::getNumberOfTurbines)
.def_property_readonly("number_of_nodes_per_blade", &ActuatorFarm::getNumberOfNodesPerBlade)
+ .def_property_readonly("number_of_turbines", &ActuatorFarm::getNumberOfTurbines)
+ .def_property_readonly("number_of_nodes_per_turbine", &ActuatorFarm::getNumberOfNodesPerTurbine)
.def_property_readonly("number_of_blades_per_turbine", &ActuatorFarm::getNumberOfBladesPerTurbine)
.def_property_readonly("number_of_grid_nodes", &ActuatorFarm::getNumberOfGridNodes)
.def_property_readonly("number_of_indices", &ActuatorFarm::getNumberOfIndices)
@@ -87,22 +95,13 @@ namespace actuator_farm
.def_property_readonly("delta_t", &ActuatorFarm::getDeltaT)
.def_property_readonly("delta_x", &ActuatorFarm::getDeltaX)
- .def("add_turbine", &ActuatorFarm::addTurbine, py::arg("posX"), py::arg("posY"), py::arg("posZ"), py::arg("diameter"), py::arg("omega"), py::arg("azimuth"), py::arg("yaw"), py::arg("bladeRadii"))
-
.def("get_turbine_pos", [](ActuatorFarm& al, uint turbine){
real position[3] = {al.getTurbinePosX(turbine), al.getTurbinePosY(turbine), al.getTurbinePosZ(turbine)}; return arr(3, position);
}, py::arg("turbine"))
- .def("get_turbine_azimuth", &ActuatorFarm::getTurbineAzimuth, py::arg("turbine"))
- .def("get_turbine_yaw", &ActuatorFarm::getTurbineYaw, py::arg("turbine"))
- .def("get_turbine_omega", &ActuatorFarm::getTurbineOmega, py::arg("turbine"))
- .def("get_all_azimuths", [](ActuatorFarm& al){ return arr(al.getNumberOfTurbines(), al.getAllAzimuths()); } )
- .def("get_all_yaws", [](ActuatorFarm& al){ return arr(al.getNumberOfTurbines(), al.getAllYaws()); } )
- .def("get_all_omegas", [](ActuatorFarm& al){ return arr(al.getNumberOfTurbines(), al.getAllOmegas()); } )
.def("get_all_turbine_pos_x", [](ActuatorFarm& al){ return arr(al.getNumberOfTurbines(), al.getAllTurbinePosX()); } )
.def("get_all_turbine_pos_y", [](ActuatorFarm& al){ return arr(al.getNumberOfTurbines(), al.getAllTurbinePosY()); } )
.def("get_all_turbine_pos_z", [](ActuatorFarm& al){ return arr(al.getNumberOfTurbines(), al.getAllTurbinePosZ()); } )
- .def("get_all_blade_radii", [](ActuatorFarm& al){ return arr({al.getNumberOfTurbines(), al.getNumberOfNodesPerBlade()}, al.getAllBladeRadii()); } )
.def("get_all_blade_coords_x", [](ActuatorFarm& al){ return arr({al.getNumberOfTurbines(), al.getNumberOfBladesPerTurbine(), al.getNumberOfNodesPerBlade()}, al.getAllBladeCoordsX()); } )
.def("get_all_blade_coords_y", [](ActuatorFarm& al){ return arr({al.getNumberOfTurbines(), al.getNumberOfBladesPerTurbine(), al.getNumberOfNodesPerBlade()}, al.getAllBladeCoordsY()); } )
.def("get_all_blade_coords_z", [](ActuatorFarm& al){ return arr({al.getNumberOfTurbines(), al.getNumberOfBladesPerTurbine(), al.getNumberOfNodesPerBlade()}, al.getAllBladeCoordsZ()); } )
@@ -113,7 +112,6 @@ namespace actuator_farm
.def("get_all_blade_forces_y", [](ActuatorFarm& al){ return arr({al.getNumberOfTurbines(), al.getNumberOfBladesPerTurbine(), al.getNumberOfNodesPerBlade()}, al.getAllBladeForcesY()); } )
.def("get_all_blade_forces_z", [](ActuatorFarm& al){ return arr({al.getNumberOfTurbines(), al.getNumberOfBladesPerTurbine(), al.getNumberOfNodesPerBlade()}, al.getAllBladeForcesZ()); } )
- .def("get_turbine_blade_radii", [](ActuatorFarm& al, uint turbine){ return arr(al.getNumberOfNodesPerBlade(), al.getTurbineBladeRadiiDevice(turbine)); } , py::arg("turbine"))
.def("get_turbine_blade_coords_x", [](ActuatorFarm& al, uint turbine){ return arr({al.getNumberOfBladesPerTurbine(), al.getNumberOfNodesPerBlade()}, al.getTurbineBladeCoordsXDevice(turbine)); }, py::arg("turbine") )
.def("get_turbine_blade_coords_y", [](ActuatorFarm& al, uint turbine){ return arr({al.getNumberOfBladesPerTurbine(), al.getNumberOfNodesPerBlade()}, al.getTurbineBladeCoordsYDevice(turbine)); }, py::arg("turbine") )
.def("get_turbine_blade_coords_z", [](ActuatorFarm& al, uint turbine){ return arr({al.getNumberOfBladesPerTurbine(), al.getNumberOfNodesPerBlade()}, al.getTurbineBladeCoordsZDevice(turbine)); }, py::arg("turbine") )
@@ -124,7 +122,6 @@ namespace actuator_farm
.def("get_turbine_blade_forces_y", [](ActuatorFarm& al, uint turbine){ return arr({al.getNumberOfBladesPerTurbine(), al.getNumberOfNodesPerBlade()}, al.getTurbineBladeForcesYDevice(turbine)); }, py::arg("turbine") )
.def("get_turbine_blade_forces_z", [](ActuatorFarm& al, uint turbine){ return arr({al.getNumberOfBladesPerTurbine(), al.getNumberOfNodesPerBlade()}, al.getTurbineBladeForcesZDevice(turbine)); }, py::arg("turbine") )
- .def("get_all_blade_radii_device", [](ActuatorFarm& al) -> intptr_t { return arr_to_cp(al.getAllBladeRadiiDevice()); } )
.def("get_all_blade_coords_x_device", [](ActuatorFarm& al) -> intptr_t { return arr_to_cp(al.getAllBladeCoordsXDevice()); } )
.def("get_all_blade_coords_y_device", [](ActuatorFarm& al) -> intptr_t { return arr_to_cp(al.getAllBladeCoordsYDevice()); } )
.def("get_all_blade_coords_z_device", [](ActuatorFarm& al) -> intptr_t { return arr_to_cp(al.getAllBladeCoordsZDevice()); } )
@@ -135,7 +132,6 @@ namespace actuator_farm
.def("get_all_blade_forces_y_device", [](ActuatorFarm& al) -> intptr_t { return arr_to_cp(al.getAllBladeForcesYDevice()); } )
.def("get_all_blade_forces_z_device", [](ActuatorFarm& al) -> intptr_t { return arr_to_cp(al.getAllBladeForcesZDevice()); } )
- .def("get_turbine_blade_radii_device", [](ActuatorFarm& al, uint turbine) -> intptr_t { return arr_to_cp(al.getTurbineBladeRadiiDevice(turbine)); }, py::arg("turbine") )
.def("get_turbine_blade_coords_x_device", [](ActuatorFarm& al, uint turbine) -> intptr_t { return arr_to_cp(al.getTurbineBladeCoordsXDevice(turbine)); }, py::arg("turbine") )
.def("get_turbine_blade_coords_y_device", [](ActuatorFarm& al, uint turbine) -> intptr_t { return arr_to_cp(al.getTurbineBladeCoordsYDevice(turbine)); }, py::arg("turbine") )
.def("get_turbine_blade_coords_z_device", [](ActuatorFarm& al, uint turbine) -> intptr_t { return arr_to_cp(al.getTurbineBladeCoordsZDevice(turbine)); }, py::arg("turbine") )
@@ -146,14 +142,6 @@ namespace actuator_farm
.def("get_turbine_blade_forces_y_device", [](ActuatorFarm& al, uint turbine) -> intptr_t { return arr_to_cp(al.getTurbineBladeForcesYDevice(turbine)); }, py::arg("turbine") )
.def("get_turbine_blade_forces_z_device", [](ActuatorFarm& al, uint turbine) -> intptr_t { return arr_to_cp(al.getTurbineBladeForcesZDevice(turbine)); }, py::arg("turbine") )
- .def("set_all_azimuths", [](ActuatorFarm& al, arr azimuths){ al.setAllAzimuths(np_to_arr(azimuths)); }, py::arg("azimuths"))
- .def("set_all_yaws", [](ActuatorFarm& al, arr yaws){ al.setAllYaws(np_to_arr(yaws)); }, py::arg("yaws"))
- .def("set_all_omegas", [](ActuatorFarm& al, arr omegas){ al.setAllOmegas(np_to_arr(omegas)); }, py::arg("omegas"))
-
- .def("set_turbine_azimuth", &ActuatorFarm::setTurbineAzimuth, py::arg("turbine"), py::arg("azimuth"))
- .def("set_turbine_yaw", &ActuatorFarm::setTurbineYaw, py::arg("turbine"), py::arg("yaw"))
- .def("set_turbine_omega", &ActuatorFarm::setTurbineOmega, py::arg("turbine"), py::arg("omega"))
-
.def("set_all_blade_coords", [](ActuatorFarm& al, arr coordsX, arr coordsY, arr coordsZ){
al.setAllBladeCoords(np_to_arr(coordsX), np_to_arr(coordsY), np_to_arr(coordsZ));
}, py::arg("blade_coords_x"), py::arg("blade_coords_y"), py::arg("blade_coords_z") )
@@ -172,6 +160,34 @@ namespace actuator_farm
.def("set_turbine_blade_forces", [](ActuatorFarm& al, uint turbine, arr forcesX, arr forcesY, arr forcesZ){
al.setTurbineBladeForces(turbine, np_to_arr(forcesX), np_to_arr(forcesY), np_to_arr(forcesZ));
}, py::arg("turbine"), py::arg("blade_forces_x"), py::arg("blade_forces_y"), py::arg("blade_forces_z") )
- .def("calc_blade_forces", &ActuatorFarm::calcBladeForces);
+ .def("update_forces_and_coordinates", &ActuatorFarm::updateForcesAndCoordinates)
+ .def("enable_output", &ActuatorFarm::enableOutput, py::arg("output_name"), py::arg("t_start_out"), py::arg("t_out"))
+ .def("set_turbine_azimuth", &ActuatorFarm::setTurbineAzimuth, py::arg("turbine"), py::arg("azimuth"));
+
+ py::class_<ActuatorFarmStandalone, ActuatorFarm, std::shared_ptr<ActuatorFarmStandalone>>(parentModule, "ActuatorFarmStandalone")
+ .def(py::init< const real,
+ const uint,
+ const std::vector<real>,
+ const std::vector<real>,
+ const std::vector<real>,
+ const std::vector<real>,
+ const real,
+ const real,
+ const int,
+ const real,
+ const real>(),
+ py::arg("diameter"),
+ py::arg("number_of_nodes_per_blade"),
+ py::arg("turbine_positions_x"),
+ py::arg("turbine_positions_y"),
+ py::arg("turbine_positions_z"),
+ py::arg("rotor_speeds"),
+ py::arg("density"),
+ py::arg("smearing_width"),
+ py::arg("level"),
+ py::arg("delta_t"),
+ py::arg("delta_x"));
}
+
+
}
\ No newline at end of file
diff --git a/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarm.cu b/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarm.cu
index fac6b3cb70c8bb2694f89e886b9192dfeed86782..b0602d58a4aa3964d7ddc4edb12984b9fdcf01dc 100644
--- a/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarm.cu
+++ b/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarm.cu
@@ -42,7 +42,7 @@
#include <basics/constants/NumericConstants.h>
#include <basics/writer/WbWriterVtkXmlBinary.h>
-#include "gpu/core/GPU/GeometryUtils.h"
+#include "VirtualFluids_GPU/GPU/GeometryUtils.h"
#include "LBM/GPUHelperFunctions/KernelUtilities.h"
#include "Parameter/Parameter.h"
#include "Parameter/CudaStreamManager.h"
@@ -51,351 +51,187 @@
using namespace vf::basics::constant;
-__host__ __device__ __forceinline__ real distSqrd(real distX, real distY, real distZ)
+struct GridData
{
- return distX * distX + distY * distY + distZ * distZ;
-}
-
-void swapArrays(real* &arr1, real* &arr2)
-{
- real* tmp = arr1;
- arr1 = arr2;
- arr2 = tmp;
-}
-
-__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)
+ 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 inverseDeltaX, velocityRatio;
+};
+
+struct TurbineData
{
- 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);
-}
+ const real* posX, *posY, *posZ;
+ const uint numberOfTurbines;
+ const real smearingWidth, factorGaussian;
+};
-__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)
+struct ComponentData
{
+ const real referenceLength;
+ const uint numberOfNodesPerTurbine;
+ const real* coordsX, * coordsY, * coordsZ;
+ real* velocitiesX, * velocitiesY, * velocitiesZ;
+ const real* forcesX, * forcesY, * forcesZ;
+ uint* gridIndices;
+};
- 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,
- real* bladeCoordsX, real* bladeCoordsY, real* bladeCoordsZ,
- real* bladeVelocitiesX, real* bladeVelocitiesY, real* bladeVelocitiesZ,
- uint numberOfTurbines, uint numberOfBlades, uint numberOfBladeNodes,
- real* azimuths, real* yaws, real* omegas,
- real* turbPosX, real* turbPosY, real* turbPosZ,
- uint* bladeIndices, real velocityRatio, real invDeltaX)
+__global__ void interpolateVelocities(const GridData gridData, const TurbineData turbineData, ComponentData componentData)
{
-
- ////////////////////////////////////////////////////////////////////////////////
- //! - Get node index coordinates from threadIdx, blockIdx, blockDim and gridDim.
- //!
const unsigned nodeIndex = vf::gpu::getNodeIndex();
- if (nodeIndex >= numberOfBladeNodes * numberOfBlades * numberOfTurbines) return;
-
- auto [turbine, blade, bladeNodeIndexOnBlade] = calcTurbineBladeAndBladeNode(nodeIndex, numberOfBladeNodes, numberOfBlades);
-
- const real localAzimuth = azimuths[turbine] + blade * c2Pi / numberOfBlades;
- const real yaw = yaws[turbine];
+ if (nodeIndex >= componentData.numberOfNodesPerTurbine * turbineData.numberOfTurbines) return;
- real bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF;
- calcCoordinateOfBladeNodeInGridFrameOfReference(nodeIndex, turbine, localAzimuth, yaw,
- bladeCoordsX, bladeCoordsY, bladeCoordsZ,
- turbPosX, turbPosY, turbPosZ,
- bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF);
+ const real coordX = componentData.coordsX[nodeIndex];
+ const real coordY = componentData.coordsY[nodeIndex];
+ const real coordZ = componentData.coordsZ[nodeIndex];
uint k, ke, kn, kt;
uint kne, kte, ktn, ktne;
- k = findNearestCellBSW(bladeIndices[nodeIndex],
- gridCoordsX, gridCoordsY, gridCoordsZ,
- bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF,
- neighborsX, neighborsY, neighborsZ, neighborsWSB);
+ k = findNearestCellBSW(componentData.gridIndices[nodeIndex],
+ gridData.coordsX, gridData.coordsY, gridData.coordsZ,
+ coordX, coordY, coordZ,
+ gridData.neighborsX, gridData.neighborsY, gridData.neighborsZ, gridData.neighborsWSB);
- bladeIndices[nodeIndex] = k;
+ componentData.gridIndices[nodeIndex] = k;
- getNeighborIndicesOfBSW(k, ke, kn, kt, kne, kte, ktn, ktne, neighborsX, neighborsY, neighborsZ);
+ getNeighborIndicesOfBSW(k, ke, kn, kt, kne, kte, ktn, ktne, gridData.neighborsX, gridData.neighborsY, gridData.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]);
+ const real distX = gridData.inverseDeltaX * (coordX - gridData.coordsX[k]);
+ const real distY = gridData.inverseDeltaX * (coordY - gridData.coordsY[k]);
+ const real distZ = gridData.inverseDeltaX * (coordZ - gridData.coordsZ[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[nodeIndex] = bladeVelX_BF;
- bladeVelocitiesY[nodeIndex] = bladeVelY_BF + omegas[turbine] * bladeCoordsZ[nodeIndex];
- bladeVelocitiesZ[nodeIndex] = bladeVelZ_BF;
+ 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(real* gridCoordsX, real* gridCoordsY, real* gridCoordsZ,
- real* gridForcesX, real* gridForcesY, real* gridForcesZ,
- real* bladeCoordsX, real* bladeCoordsY, real* bladeCoordsZ,
- real* bladeForcesX, real* bladeForcesY, real* bladeForcesZ,
- const uint numberOfTurbines, const uint numberOfBlades, const uint numberOfBladeNodes,
- real* azimuths, real* yaws, real* diameters,
- real* turbPosX, real* turbPosY, real* turbPosZ,
- uint* gridIndices, uint nIndices,
- const real invEpsilonSqrd, const real factorGaussian)
+__global__ void applyBodyForces(GridData gridData, const TurbineData turbineData, const ComponentData componentData)
{
const uint index = vf::gpu::getNodeIndex();
- 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;
+ if (index >= gridData.nIndices) return;
- real dAzimuth = c2Pi / numberOfBlades;
+ const uint gridIndex = gridData.indices[index];
+ const real gridCoordX = gridData.coordsX[gridIndex];
+ const real gridCoordY = gridData.coordsY[gridIndex];
+ const real gridCoordZ = gridData.coordsZ[gridIndex];
- for (uint turbine = 0; turbine < numberOfTurbines; 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];
+ real gridForceX = c0o1;
+ real gridForceY = c0o1;
+ real gridForceZ = c0o1;
- if (distSqrd(gridCoordX_RF, gridCoordY_RF, gridCoordZ_RF) * invEpsilonSqrd > radius * radius * invEpsilonSqrd + c7o1) continue;
+ for (uint turbine = 0; turbine < turbineData.numberOfTurbines; turbine++) {
+ const real distToHubX = gridCoordX - turbineData.posX[turbine];
+ const real distToHubY = gridCoordY - turbineData.posY[turbine];
+ const real distToHubZ = gridCoordZ - turbineData.posZ[turbine];
- real azimuth = azimuths[turbine];
- real yaw = yaws[turbine];
+ if (!inBoundingSphere(distToHubX, distToHubY, distToHubZ, componentData.referenceLength, turbineData.smearingWidth)) continue;
- for (uint blade = 0; blade < numberOfBlades; blade++) {
- real localAzimuth = azimuth + blade * dAzimuth;
+ for (uint turbineNode = 0; turbineNode < componentData.numberOfNodesPerTurbine; turbineNode++) {
+ const uint node = turbine*componentData.numberOfNodesPerTurbine + turbineNode;
- real gridCoordX_BF, gridCoordY_BF, gridCoordZ_BF;
+ const real distX = componentData.coordsX[node] - gridCoordX;
+ const real distY = componentData.coordsY[node] - gridCoordY;
+ const real distZ = componentData.coordsZ[node] - gridCoordZ;
- rotateFromGlobalToBlade(gridCoordX_BF, gridCoordY_BF, gridCoordZ_BF,
- gridCoordX_RF, gridCoordY_RF, gridCoordZ_RF,
- localAzimuth, yaw);
-
- uint node;
- uint nextNode = calcNodeIndexInBladeArrays({ /*.turbine = */ turbine, /* .blade = */ blade, /*.bladeNode =*/0 }, numberOfBladeNodes, numberOfBlades);
-
- real last_z = c0o1;
- real current_z = c0o1;
- real next_z = bladeCoordsZ[nextNode];
-
- real x, y, dz, eta, forceX_RF, forceY_RF, forceZ_RF;
-
- for (uint bladeNode = 0; bladeNode < numberOfBladeNodes - 1; bladeNode++) {
- node = nextNode;
- nextNode = calcNodeIndexInBladeArrays({ /*.turbine = */ turbine, /* .blade = */ blade, /*.bladeNode =*/bladeNode + 1 }, numberOfBladeNodes, numberOfBlades);
-
- x = bladeCoordsX[node];
- y = bladeCoordsY[node];
- last_z = current_z;
- current_z = next_z;
- next_z = bladeCoordsZ[nextNode];
-
- 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;
- }
-
- // Handle last node separately
-
- node = nextNode;
-
- x = bladeCoordsX[node];
- y = bladeCoordsY[node];
- last_z = current_z;
- current_z = next_z;
-
- dz = c1o2 * (radius - 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;
+ const real eta = gaussianSmearing(distX, distY, distZ, turbineData.smearingWidth, turbineData.factorGaussian);
+ gridForceX += componentData.forcesX[node] * eta;
+ gridForceY += componentData.forcesY[node] * eta;
+ gridForceZ += componentData.forcesZ[node] * eta;
}
}
-
- gridForcesX[gridIndex] += gridForceX_RF;
- gridForcesY[gridIndex] += gridForceY_RF;
- gridForcesZ[gridIndex] += gridForceZ_RF;
+ gridData.fx[gridIndex] += gridForceX;
+ gridData.fy[gridIndex] += gridForceY;
+ gridData.fz[gridIndex] += gridForceZ;
}
-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)
+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.");
this->forceRatio = para->getForceRatio();
- this->initTurbineGeometries(cudaMemoryManager);
- this->initBladeCoords(cudaMemoryManager);
- this->initBladeIndices(para, cudaMemoryManager);
- this->initBladeVelocities(cudaMemoryManager);
- this->initBladeForces(cudaMemoryManager);
- this->initBoundingSpheres(para, cudaMemoryManager);
+ this->initTurbineGeometries(cudaManager);
+ this->initBladeCoords(cudaManager);
+ this->initBladeIndices(para, cudaManager);
+ this->initBladeVelocities(cudaManager);
+ this->initBladeForces(cudaManager);
+ this->initBoundingSpheres(para, cudaManager);
this->streamIndex = para->getStreamManager()->registerAndLaunchStream(CudaStreamIndex::ActuatorFarm);
}
-void ActuatorFarm::interact(Parameter *para, CudaMemoryManager *cudaMemoryManager, int level, unsigned int t)
+void ActuatorFarm::interact(Parameter *para, CudaMemoryManager *cudaManager, int level, unsigned int t)
{
if (level != this->level) return;
cudaStream_t stream = para->getStreamManager()->getStream(CudaStreamIndex::ActuatorFarm, this->streamIndex);
- if (useHostArrays) cudaMemoryManager->cudaCopyBladeCoordsHtoD(this);
+ if (useHostArrays) cudaManager->cudaCopyBladeCoordsHtoD(this);
- vf::cuda::CudaGrid bladeGrid = vf::cuda::CudaGrid(para->getParH(level)->numberofthreads, this->numberOfGridNodes);
-
- // if (t == 0) this->writeBladeCoordsToVtkForDebug("output/ALM/bladecoords_" + std::to_string(t));
+ if(this->writeOutput && ((t-this->tStartOut) % this->tOut == 0))
+ {
+ if(!useHostArrays){
+ cudaManager->cudaCopyBladeCoordsDtoH(this);
+ cudaManager->cudaCopyBladeVelocitiesDtoH(this);
+ cudaManager->cudaCopyBladeForcesDtoH(this);
+ }
+ this->write(this->getFilename(para, t));
+ }
- interpolateVelocities<<< bladeGrid.grid, bladeGrid.threads, 0, stream >>>(
+ const GridData gridData {
+ this->boundingSphereIndicesD, this->numberOfIndices,
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->bladeCoordsXDCurrentTimestep, this->bladeCoordsYDCurrentTimestep, this->bladeCoordsZDCurrentTimestep,
- this->bladeVelocitiesXDCurrentTimestep, this->bladeVelocitiesYDCurrentTimestep, this->bladeVelocitiesZDCurrentTimestep,
- this->numberOfTurbines, this->numberOfBlades, this->numberOfBladeNodes,
- this->azimuthsD, this->yawsD, this->omegasD,
+ para->getParD(this->level)->forceX_SP,para->getParD(this->level)->forceY_SP,para->getParD(this->level)->forceZ_SP,
+ this->invDeltaX, para->getVelocityRatio()};
+
+ const TurbineData turbineData {
this->turbinePosXD, this->turbinePosYD, this->turbinePosZD,
- this->bladeIndicesD, para->getVelocityRatio(), this->invDeltaX);
+ this->numberOfTurbines,
+ this->smearingWidth, this->factorGaussian};
+
+ const ComponentData bladeData {
+ this->diameter, this->numberOfNodesPerTurbine,
+ this->bladeCoordsXDCurrentTimestep, this->bladeCoordsYDCurrentTimestep, this->bladeCoordsZDCurrentTimestep,
+ this->bladeVelocitiesXDCurrentTimestep, this->bladeVelocitiesYDCurrentTimestep, this->bladeVelocitiesZDCurrentTimestep,
+ this->bladeForcesXDCurrentTimestep, this->bladeForcesYDCurrentTimestep, this->bladeForcesZDCurrentTimestep,
+ 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);
cudaStreamSynchronize(stream);
- if (useHostArrays) cudaMemoryManager->cudaCopyBladeVelocitiesDtoH(this);
- this->calcBladeForces();
+
+ if (useHostArrays) cudaManager->cudaCopyBladeVelocitiesDtoH(this);
+
+ this->updateForcesAndCoordinates();
this->swapDeviceArrays();
- if (useHostArrays) cudaMemoryManager->cudaCopyBladeForcesHtoD(this);
+ if (useHostArrays) cudaManager->cudaCopyBladeForcesHtoD(this);
vf::cuda::CudaGrid sphereGrid = vf::cuda::CudaGrid(para->getParH(level)->numberofthreads, this->numberOfIndices);
- applyBodyForces<<<sphereGrid.grid, sphereGrid.threads, 0, stream>>>(
- 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->bladeCoordsXDCurrentTimestep, this->bladeCoordsYDCurrentTimestep, this->bladeCoordsZDCurrentTimestep,
- this->bladeForcesXDCurrentTimestep, this->bladeForcesYDCurrentTimestep, this->bladeForcesZDCurrentTimestep,
- this->numberOfTurbines, this->numberOfBlades, this->numberOfBladeNodes,
- this->azimuthsD, this->yawsD, this->diametersD,
- this->turbinePosXD, this->turbinePosYD, this->turbinePosZD,
- this->boundingSphereIndicesD, this->numberOfIndices,
- this->invEpsilonSqrd, this->factorGaussian);
- cudaMemoryManager->cudaCopyBladeOrientationsHtoD(this);
+ applyBodyForces<<<sphereGrid.grid, sphereGrid.threads, 0, stream>>>(gridData, turbineData, bladeData);
cudaStreamSynchronize(stream);
}
-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 = 20 * c1o10;
- real c, Cn, Ct;
-
- for (uint turbine = 0; turbine < this->numberOfTurbines; turbine++)
- {
- real diameter = this->diametersH[turbine];
- for (uint blade = 0; blade < this->numberOfBlades; blade++)
- {
- for (uint bladeNode = 0; bladeNode < this->numberOfBladeNodes; bladeNode++)
- {
- uint node = calcNodeIndexInBladeArrays({ /*.turbine = */ turbine, /* .blade = */ blade, /*.bladeNode =*/bladeNode }, this->numberOfBladeNodes, this->numberOfBlades);
-
- 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);
- real fx = c1o2 * u_rel_sq * c * this->density * Cn;
- real fy = c1o2 * u_rel_sq * c * this->density * Ct;
- this->bladeForcesXH[node] = -fx;
- this->bladeForcesYH[node] = -fy;
- this->bladeForcesZH[node] = c0o1;
- // printf("u %f v %f fx %f fy %f \n", u_rel, v_rel, fx, fy);
- }
- }
- azimuthsH[turbine] = azimuthsH[turbine] + deltaT * omegasH[turbine];
- }
-}
-
-void ActuatorFarm::calcBladeForces()
+void ActuatorFarm::free(Parameter *para, CudaMemoryManager *cudaManager)
{
- this->calcForcesEllipticWing();
+ cudaManager->cudaFreeBladeGeometries(this);
+ cudaManager->cudaFreeBladeCoords(this);
+ cudaManager->cudaFreeBladeVelocities(this);
+ cudaManager->cudaFreeBladeForces(this);
+ cudaManager->cudaFreeBladeIndices(this);
+ cudaManager->cudaFreeSphereIndices(this);
}
void ActuatorFarm::getTaggedFluidNodes(Parameter *para, GridProvider *gridProvider)
@@ -404,111 +240,101 @@ void ActuatorFarm::getTaggedFluidNodes(Parameter *para, GridProvider *gridProvid
gridProvider->tagFluidNodeIndices(indicesInSphere, CollisionTemplate::AllFeatures, this->level);
}
-void ActuatorFarm::initTurbineGeometries(CudaMemoryManager *cudaMemoryManager)
+void ActuatorFarm::initTurbineGeometries(CudaMemoryManager *cudaManager)
{
- this->numberOfTurbines = uint(this->preInitDiameters.size());
- this->numberOfGridNodes = numberOfTurbines * numberOfBladeNodes * numberOfBlades;
+ this->numberOfGridNodes = this->numberOfTurbines * this->numberOfNodesPerTurbine;
- cudaMemoryManager->cudaAllocBladeGeometries(this);
- cudaMemoryManager->cudaAllocBladeOrientations(this);
+ cudaManager->cudaAllocBladeGeometries(this);
- for (uint turbine = 0; turbine < this->numberOfTurbines; turbine++) {
- for (uint node = 0; node < this->numberOfBladeNodes; node++) {
- this->bladeRadiiH[calcNodeIndexInBladeArrays({ /*.turbine = */ turbine, /* .blade = */ 0, /*.bladeNode =*/node }, numberOfBladeNodes, 1)] = 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);
- this->factorGaussian = pow(this->epsilon * sqrt(cPi), -c3o1) / this->forceRatio;
+ std::copy(initialTurbinePositionsX.begin(), initialTurbinePositionsX.end(), turbinePosXH);
+ std::copy(initialTurbinePositionsY.begin(), initialTurbinePositionsY.end(), turbinePosYH);
+ std::copy(initialTurbinePositionsZ.begin(), initialTurbinePositionsZ.end(), turbinePosZH);
+
+ cudaManager->cudaCopyBladeGeometriesHtoD(this);
+
+ this->factorGaussian = pow(this->smearingWidth * sqrt(cPi), -c3o1) / this->forceRatio;
}
-void ActuatorFarm::initBladeCoords(CudaMemoryManager *cudaMemoryManager)
+void ActuatorFarm::initBladeCoords(CudaMemoryManager *cudaManager)
{
- cudaMemoryManager->cudaAllocBladeCoords(this);
+ cudaManager->cudaAllocBladeCoords(this);
- for (uint turbine = 0; turbine < numberOfTurbines; turbine++)
+ for (uint turbine = 0; turbine < this->numberOfTurbines; turbine++)
{
- for (uint blade = 0; blade < this->numberOfBlades; blade++)
+ for (uint blade = 0; blade < ActuatorFarm::numberOfBlades; blade++)
{
- for (uint bladeNode = 0; bladeNode < this->numberOfBladeNodes; bladeNode++) {
- uint node = calcNodeIndexInBladeArrays({ /*.turbine = */ turbine, /* .blade = */ blade, /*.bladeNode =*/bladeNode }, this->numberOfBladeNodes, this->numberOfBlades);
-
- this->bladeCoordsXH[node] = c0o1;
- this->bladeCoordsYH[node] = c0o1;
- this->bladeCoordsZH[node] = this->bladeRadiiH[calcNodeIndexInBladeArrays({ /*.turbine = */ turbine, /* .blade = */ 0, /*.bladeNode =*/bladeNode }, numberOfBladeNodes, 1)];
+ 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);
+
+ real x, y, z;
+ rotateFromBladeToGlobal(c0o1, c0o1, this->bladeRadii[bladeNode], x, y, z, local_azimuth);
+ bladeCoordsXH[node] = x + this->turbinePosXH[turbine];
+ bladeCoordsYH[node] = y + this->turbinePosYH[turbine];
+ bladeCoordsZH[node] = z + this->turbinePosZH[turbine];
}
}
}
- cudaMemoryManager->cudaCopyBladeCoordsHtoD(this);
+ cudaManager->cudaCopyBladeCoordsHtoD(this);
swapArrays(this->bladeCoordsXDCurrentTimestep, this->bladeCoordsXDPreviousTimestep);
swapArrays(this->bladeCoordsYDCurrentTimestep, this->bladeCoordsYDPreviousTimestep);
swapArrays(this->bladeCoordsZDCurrentTimestep, this->bladeCoordsZDPreviousTimestep);
- cudaMemoryManager->cudaCopyBladeCoordsHtoD(this);
+ cudaManager->cudaCopyBladeCoordsHtoD(this);
}
-void ActuatorFarm::initBladeVelocities(CudaMemoryManager *cudaMemoryManager)
+void ActuatorFarm::initBladeVelocities(CudaMemoryManager *cudaManager)
{
- cudaMemoryManager->cudaAllocBladeVelocities(this);
+ cudaManager->cudaAllocBladeVelocities(this);
std::fill_n(this->bladeVelocitiesXH, this->numberOfGridNodes, c0o1);
std::fill_n(this->bladeVelocitiesYH, this->numberOfGridNodes, c0o1);
std::fill_n(this->bladeVelocitiesZH, this->numberOfGridNodes, c0o1);
- cudaMemoryManager->cudaCopyBladeVelocitiesHtoD(this);
+ cudaManager->cudaCopyBladeVelocitiesHtoD(this);
swapArrays(this->bladeVelocitiesXDCurrentTimestep, this->bladeVelocitiesXDPreviousTimestep);
swapArrays(this->bladeVelocitiesYDCurrentTimestep, this->bladeVelocitiesYDPreviousTimestep);
swapArrays(this->bladeVelocitiesZDCurrentTimestep, this->bladeVelocitiesZDPreviousTimestep);
- cudaMemoryManager->cudaCopyBladeVelocitiesHtoD(this);
+ cudaManager->cudaCopyBladeVelocitiesHtoD(this);
}
-void ActuatorFarm::initBladeForces(CudaMemoryManager *cudaMemoryManager)
+void ActuatorFarm::initBladeForces(CudaMemoryManager *cudaManager)
{
- cudaMemoryManager->cudaAllocBladeForces(this);
+ cudaManager->cudaAllocBladeForces(this);
std::fill_n(this->bladeForcesXH, this->numberOfGridNodes, c0o1);
std::fill_n(this->bladeForcesYH, this->numberOfGridNodes, c0o1);
std::fill_n(this->bladeForcesZH, this->numberOfGridNodes, c0o1);
- cudaMemoryManager->cudaCopyBladeForcesHtoD(this);
+ cudaManager->cudaCopyBladeForcesHtoD(this);
swapArrays(this->bladeForcesXDCurrentTimestep, this->bladeForcesXDPreviousTimestep);
swapArrays(this->bladeForcesYDCurrentTimestep, this->bladeForcesYDPreviousTimestep);
swapArrays(this->bladeForcesZDCurrentTimestep, this->bladeForcesZDPreviousTimestep);
- cudaMemoryManager->cudaCopyBladeForcesHtoD(this);
+ cudaManager->cudaCopyBladeForcesHtoD(this);
}
-void ActuatorFarm::initBladeIndices(Parameter *para, CudaMemoryManager *cudaMemoryManager)
+void ActuatorFarm::initBladeIndices(Parameter *para, CudaMemoryManager *cudaManager)
{
- cudaMemoryManager->cudaAllocBladeIndices(this);
+ cudaManager->cudaAllocBladeIndices(this);
std::fill_n(this->bladeIndicesH, this->numberOfGridNodes, 1);
- cudaMemoryManager->cudaCopyBladeIndicesHtoD(this);
+ cudaManager->cudaCopyBladeIndicesHtoD(this);
}
-void ActuatorFarm::initBoundingSpheres(Parameter *para, CudaMemoryManager *cudaMemoryManager)
+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));
for (uint turbine = 0; turbine < this->numberOfTurbines; turbine++) {
- real sphereRadius = c1o2 * this->diametersH[turbine] + c4o1 * this->epsilon;
-
- real posX = this->turbinePosXH[turbine];
- real posY = this->turbinePosYH[turbine];
- real posZ = this->turbinePosZH[turbine];
- real sphereRadiusSqrd = sphereRadius * sphereRadius;
+ const real posX = this->turbinePosXH[turbine];
+ const real posY = this->turbinePosYH[turbine];
+ const real posZ = this->turbinePosZH[turbine];
- uint minimumNumberOfNodesPerSphere =
- (uint)(c4o3 * cPi * pow(sphereRadius - this->deltaX, c3o1) / pow(this->deltaX, c3o1));
uint nodesInThisSphere = 0;
for (size_t pos = 1; pos <= para->getParH(this->level)->numberOfNodes; pos++) {
@@ -529,65 +355,50 @@ void ActuatorFarm::initBoundingSpheres(Parameter *para, CudaMemoryManager *cudaM
this->numberOfIndices = uint(nodesInSpheres.size());
- cudaMemoryManager->cudaAllocSphereIndices(this);
+ cudaManager->cudaAllocSphereIndices(this);
std::copy(nodesInSpheres.begin(), nodesInSpheres.end(), this->boundingSphereIndicesH);
- cudaMemoryManager->cudaCopySphereIndicesHtoD(this);
-}
-
-void ActuatorFarm::setAllAzimuths(real *_azimuths)
-{
- std::copy_n(_azimuths, this->numberOfTurbines, this->azimuthsH);
+ cudaManager->cudaCopySphereIndicesHtoD(this);
}
-void ActuatorFarm::setAllOmegas(real *_omegas)
-{
- std::copy_n(_omegas, this->numberOfTurbines, this->omegasH);
-}
-
-void ActuatorFarm::setAllYaws(real *_yaws)
-{
- std::copy_n(_yaws, this->numberOfTurbines, this->yawsH);
-}
-
-void ActuatorFarm::setAllBladeCoords(real *_bladeCoordsX, real *_bladeCoordsY, 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(real *_bladeVelocitiesX, real *_bladeVelocitiesY, 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(real *_bladeForcesX, real *_bladeForcesY, 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, real *_bladeCoordsX, real *_bladeCoordsY, real *_bladeCoordsZ)
+void ActuatorFarm::setTurbineBladeCoords(uint turbine, const real *_bladeCoordsX, const real *_bladeCoordsY, const real *_bladeCoordsZ)
{
- std::copy_n(_bladeCoordsX, numberOfBladeNodes * numberOfBlades, &this->bladeCoordsXH[turbine * numberOfBladeNodes * numberOfBlades]);
- std::copy_n(_bladeCoordsY, numberOfBladeNodes * numberOfBlades, &this->bladeCoordsYH[turbine * numberOfBladeNodes * numberOfBlades]);
- std::copy_n(_bladeCoordsZ, numberOfBladeNodes * numberOfBlades, &this->bladeCoordsZH[turbine * numberOfBladeNodes * numberOfBlades]);
+ 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, real *_bladeVelocitiesX, real *_bladeVelocitiesY, real *_bladeVelocitiesZ)
+void ActuatorFarm::setTurbineBladeVelocities(uint turbine, const real *_bladeVelocitiesX, const real *_bladeVelocitiesY, const real *_bladeVelocitiesZ)
{
- std::copy_n(_bladeVelocitiesX, numberOfBladeNodes * numberOfBlades, &this->bladeVelocitiesXH[turbine * numberOfBladeNodes * numberOfBlades]);
- std::copy_n(_bladeVelocitiesY, numberOfBladeNodes * numberOfBlades, &this->bladeVelocitiesYH[turbine * numberOfBladeNodes * numberOfBlades]);
- std::copy_n(_bladeVelocitiesZ, numberOfBladeNodes * numberOfBlades, &this->bladeVelocitiesZH[turbine * numberOfBladeNodes * numberOfBlades]);
+ 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, real *_bladeForcesX, real *_bladeForcesY, real *_bladeForcesZ)
+void ActuatorFarm::setTurbineBladeForces(uint turbine, const real *_bladeForcesX, const real *_bladeForcesY, const real *_bladeForcesZ)
{
- std::copy_n(_bladeForcesX, numberOfBladeNodes * numberOfBlades, &this->bladeForcesXH[turbine * numberOfBladeNodes * numberOfBlades]);
- std::copy_n(_bladeForcesY, numberOfBladeNodes * numberOfBlades, &this->bladeForcesYH[turbine * numberOfBladeNodes * numberOfBlades]);
- std::copy_n(_bladeForcesZ, numberOfBladeNodes * numberOfBlades, &this->bladeForcesZH[turbine * numberOfBladeNodes * numberOfBlades]);
+ std::copy_n(_bladeForcesX, this->numberOfNodesPerTurbine, &this->bladeForcesXH[turbine * this->numberOfNodesPerTurbine]);
+ std::copy_n(_bladeForcesY, this->numberOfNodesPerTurbine, &this->bladeForcesYH[turbine * this->numberOfNodesPerTurbine]);
+ std::copy_n(_bladeForcesZ, this->numberOfNodesPerTurbine, &this->bladeForcesZH[turbine * this->numberOfNodesPerTurbine]);
}
void ActuatorFarm::swapDeviceArrays()
@@ -605,25 +416,32 @@ void ActuatorFarm::swapDeviceArrays()
swapArrays(this->bladeForcesZDPreviousTimestep, this->bladeForcesZDCurrentTimestep);
}
-void ActuatorFarm::writeBladeCoordsToVtkForDebug(const std::string &filename)
+std::string ActuatorFarm::getFilename(Parameter* para, uint t) const
{
- std::vector<real> coordsX(numberOfGridNodes);
- std::vector<real> coordsY(numberOfGridNodes);
- std::vector<real> coordsZ(numberOfGridNodes);
-
- for (uint nodeIndex = 0; nodeIndex < numberOfGridNodes; nodeIndex++) {
- auto [turbine, blade, bladeNodeIndexOnBlade] = calcTurbineBladeAndBladeNode(nodeIndex, numberOfBladeNodes, numberOfBlades);
-
- const real localAzimuth = this->getTurbineAzimuth(turbine) + blade * c2Pi / numberOfBlades;
- const real yaw = this->getTurbineYaw(turbine);
- real bladeCoordX_GF;
- real bladeCoordY_GF;
- real bladeCoordZ_GF;
- calcCoordinateOfBladeNodeInGridFrameOfReference(nodeIndex, turbine, localAzimuth, yaw, this->bladeCoordsXH, this->bladeCoordsYH, this->bladeCoordsZH, this->turbinePosXH, this->turbinePosYH, this->turbinePosZH, bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF);
- coordsX[nodeIndex] = bladeCoordX_GF;
- coordsY[nodeIndex] = bladeCoordY_GF;
- coordsZ[nodeIndex] = bladeCoordZ_GF;
- }
+ return para->getOutputPath() + this->outputName + "_t_" + std::to_string(t);
+}
- WbWriterVtkXmlBinary::getInstance()->writePolyLines(filename, coordsX, coordsY, coordsZ);
-}
\ No newline at end of file
+void ActuatorFarm::write(const std::string &filename)
+{
+ std::vector<std::string> dataNames = {
+ "bladeVelocitiesX",
+ "bladeVelocitiesY",
+ "bladeVelocitiesZ",
+ "bladeForcesX",
+ "bladeForcesY",
+ "bladeForcesZ",
+ };
+ 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++){
+ nodes[i] = UbTupleFloat3(this->bladeCoordsXH[i], this->bladeCoordsYH[i], this->bladeCoordsZH[i]);
+ nodeData[0][i] = this->bladeVelocitiesXH[i];
+ nodeData[1][i] = this->bladeVelocitiesYH[i];
+ nodeData[2][i] = this->bladeVelocitiesZH[i];
+ nodeData[3][i] = this->bladeForcesXH[i];
+ nodeData[4][i] = this->bladeForcesYH[i];
+ nodeData[5][i] = this->bladeForcesZH[i];
+ }
+ WbWriterVtkXmlBinary::getInstance()->writeNodesWithNodeData(filename, nodes, dataNames, nodeData);
+}
diff --git a/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarm.h b/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarm.h
index 74f23f3e6a4dc9e3e798ca13aa3632171d0e3bd0..c94e9f885f2c6815880d14cc83f03d1811a6798d 100644
--- a/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarm.h
+++ b/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarm.h
@@ -1,146 +1,141 @@
#ifndef ActuatorFarm_H
#define ActuatorFarm_H
-#include "PreCollisionInteractor.h"
-#include <basics/PointerDefinitions.h>
+#include "PreCollisionInteractor/PreCollisionInteractor.h"
#include "basics/constants/NumericConstants.h"
+#include "basics/DataTypes.h"
#include <stdexcept>
-using namespace vf::basics::constant;
-
class Parameter;
class GridProvider;
using namespace vf::basics::constant;
+
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
+ const real diameter,
+ const std::vector<real> bladeRadii,
+ const std::vector<real> turbinePositionsX,
+ const std::vector<real> turbinePositionsY,
+ const std::vector<real> turbinePositionsZ,
+ const real density,
+ const real smearingWidth,
+ const int level,
+ const real deltaT,
+ const real deltaX,
+ const bool useHostArrays
) :
- numberOfBlades(_nBlades),
- density(_density),
- numberOfBladeNodes(_nBladeNodes),
- epsilon(_epsilon),
- level(_level),
- useHostArrays(_useHostArrays),
- numberOfTurbines(0),
- numberOfGridNodes(0),
- PreCollisionInteractor()
+ diameter(diameter),
+ bladeRadii(bladeRadii),
+ numberOfNodesPerBlade(static_cast<uint>(bladeRadii.size())),
+ numberOfNodesPerTurbine(numberOfNodesPerBlade*numberOfBlades),
+ numberOfTurbines(static_cast<uint>(turbinePositionsX.size())),
+ initialTurbinePositionsX(turbinePositionsX),
+ initialTurbinePositionsY(turbinePositionsY),
+ initialTurbinePositionsZ(turbinePositionsZ),
+ density(density),
+ smearingWidth(smearingWidth),
+ level(level),
+ useHostArrays(useHostArrays),
+ deltaT(deltaT*exp2(-level)),
+ deltaX(deltaX*exp2(-level)),
+ invDeltaX(c1o1/deltaX)
{
- this->deltaT = _deltaT*exp2(-this->level);
- this->deltaX = _deltaX*exp2(-this->level);
- this->invEpsilonSqrd = 1/(epsilon*epsilon);
- this->invDeltaX = c1o1/this->deltaX;
-
- if(this->epsilon<this->deltaX)
- throw std::runtime_error("ActuatorFarm::ActuatorFarm: epsilon needs to be larger than dx!");
+ if(this->smearingWidth < this->deltaX)
+ throw std::runtime_error("ActuatorFarm::ActuatorFarm: smearing width needs to be larger than dx!");
+ if(numberOfTurbines != turbinePositionsY.size() || numberOfTurbines != turbinePositionsZ.size())
+ throw std::runtime_error("ActuatorFarm::ActuatorFarm: turbine positions need to have the same length!");
+ azimuths = std::vector<real>(numberOfTurbines, 0.0);
}
~ActuatorFarm() override = default;
- 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) override;
void interact(Parameter* para, CudaMemoryManager* cudaManager, int level, uint t) override;
void free(Parameter* para, CudaMemoryManager* cudaManager) override;
void getTaggedFluidNodes(Parameter *para, GridProvider* gridProvider) override;
- void write(uint t);
-
- real getDensity(){ return this->density; };
- real getDeltaT(){ return this->deltaT; };
- real getDeltaX(){ return this->deltaX; };
-
- uint getNumberOfTurbines(){ return this->numberOfTurbines; };
- uint getNumberOfNodesPerBlade(){ return this->numberOfBladeNodes; };
- uint getNumberOfBladesPerTurbine(){ return this->numberOfBlades; };
-
- uint getNumberOfIndices(){ return this->numberOfIndices; };
- uint getNumberOfGridNodes(){ return this->numberOfGridNodes; };
-
- 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]; };
-
- real* getAllBladeRadii(){ return this->bladeRadiiH; };
- 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* getTurbineBladeRadii(uint turbine){ return &this->bladeRadiiH[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeCoordsX(uint turbine){ return &this->bladeCoordsXH[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeCoordsY(uint turbine){ return &this->bladeCoordsYH[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeCoordsZ(uint turbine){ return &this->bladeCoordsZH[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeVelocitiesX(uint turbine){ return &this->bladeVelocitiesXH[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeVelocitiesY(uint turbine){ return &this->bladeVelocitiesYH[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeVelocitiesZ(uint turbine){ return &this->bladeVelocitiesZH[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeForcesX(uint turbine){ return &this->bladeForcesXH[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeForcesY(uint turbine){ return &this->bladeForcesYH[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeForcesZ(uint turbine){ return &this->bladeForcesZH[turbine*numberOfBladeNodes*numberOfBlades]; };
-
- real* getAllBladeRadiiDevice(){ return this->bladeRadiiD; };
- real* getAllBladeCoordsXDevice(){ return this->bladeCoordsXDCurrentTimestep; };
- real* getAllBladeCoordsYDevice(){ return this->bladeCoordsYDCurrentTimestep; };
- real* getAllBladeCoordsZDevice(){ return this->bladeCoordsZDCurrentTimestep; };
- real* getAllBladeVelocitiesXDevice(){ return this->bladeVelocitiesXDCurrentTimestep; };
- real* getAllBladeVelocitiesYDevice(){ return this->bladeVelocitiesYDCurrentTimestep; };
- real* getAllBladeVelocitiesZDevice(){ return this->bladeVelocitiesZDCurrentTimestep; };
- real* getAllBladeForcesXDevice(){ return this->bladeForcesXDCurrentTimestep; };
- real* getAllBladeForcesYDevice(){ return this->bladeForcesYDCurrentTimestep; };
- real* getAllBladeForcesZDevice(){ return this->bladeForcesZDCurrentTimestep; };
-
- real* getTurbineBladeRadiiDevice(uint turbine){ return &this->bladeRadiiD[turbine*numberOfBladeNodes]; };
- real* getTurbineBladeCoordsXDevice(uint turbine){ return &this->bladeCoordsXDCurrentTimestep[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeCoordsYDevice(uint turbine){ return &this->bladeCoordsYDCurrentTimestep[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeCoordsZDevice(uint turbine){ return &this->bladeCoordsZDCurrentTimestep[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeVelocitiesXDevice(uint turbine){ return &this->bladeVelocitiesXDCurrentTimestep[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeVelocitiesYDevice(uint turbine){ return &this->bladeVelocitiesYDCurrentTimestep[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeVelocitiesZDevice(uint turbine){ return &this->bladeVelocitiesZDCurrentTimestep[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeForcesXDevice(uint turbine){ return &this->bladeForcesXDCurrentTimestep[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeForcesYDevice(uint turbine){ return &this->bladeForcesYDCurrentTimestep[turbine*numberOfBladeNodes*numberOfBlades]; };
- real* getTurbineBladeForcesZDevice(uint turbine){ return &this->bladeForcesZDCurrentTimestep[turbine*numberOfBladeNodes*numberOfBlades]; };
-
- void setAllAzimuths(real* _azimuth);
- void setAllOmegas(real* _omegas);
- void setAllYaws(real* yaws);
-
- void setTurbineAzimuth(uint turbine, real azimuth){ azimuthsH[turbine] = azimuth; };
- void setTurbineYaw(uint turbine, real yaw){ yawsH[turbine] = yaw; };
- void setTurbineOmega(uint turbine, real omega){ omegasH[turbine] = omega; };
-
- 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);
-
- virtual void calcBladeForces();
+ void enableOutput(const std::string outputName, uint tStart, uint tOut) {
+ this->outputName = outputName;
+ this->writeOutput = true;
+ this->tStartOut = tStart;
+ this->tOut = tOut;
+ }
+
+ void write(const std::string& filename);
+
+ real getDensity() const { return this->density; };
+ real getDeltaT() const { return this->deltaT; };
+ real getDeltaX() const { return this->deltaX; };
+
+ uint getNumberOfTurbines() const { return this->numberOfTurbines; };
+ uint getNumberOfNodesPerTurbine() const { return this->numberOfNodesPerTurbine; };
+ uint getNumberOfNodesPerBlade() const { return this->numberOfNodesPerBlade; };
+ uint getNumberOfBladesPerTurbine() const { return ActuatorFarm::numberOfBlades; };
+
+ uint getNumberOfIndices() const { return this->numberOfIndices; };
+ uint getNumberOfGridNodes() const { return this->numberOfGridNodes; };
+
+ real* getAllTurbinePosX() const { return turbinePosXH; };
+ real* getAllTurbinePosY() const { return turbinePosYH; };
+ real* getAllTurbinePosZ() const { return turbinePosZH; };
+
+ real getTurbinePosX(uint turbine) const { return turbinePosXH[turbine]; };
+ real getTurbinePosY(uint turbine) const { return turbinePosYH[turbine]; };
+ real getTurbinePosZ(uint turbine) const { return turbinePosZH[turbine]; };
+
+ real* getAllBladeCoordsX() const { return this->bladeCoordsXH; };
+ real* getAllBladeCoordsY() const { return this->bladeCoordsYH; };
+ real* getAllBladeCoordsZ() const { return this->bladeCoordsZH; };
+ real* getAllBladeVelocitiesX() const { return this->bladeVelocitiesXH; };
+ real* getAllBladeVelocitiesY() const { return this->bladeVelocitiesYH; };
+ real* getAllBladeVelocitiesZ() const { return this->bladeVelocitiesZH; };
+ real* getAllBladeForcesX() const { return this->bladeForcesXH; };
+ real* getAllBladeForcesY() const { return this->bladeForcesYH; };
+ real* getAllBladeForcesZ() const { return this->bladeForcesZH; };
+
+ real* getTurbineBladeCoordsX(uint turbine) const { return &this->bladeCoordsXH[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeCoordsY(uint turbine) const { return &this->bladeCoordsYH[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeCoordsZ(uint turbine) const { return &this->bladeCoordsZH[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeVelocitiesX(uint turbine) const { return &this->bladeVelocitiesXH[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeVelocitiesY(uint turbine) const { return &this->bladeVelocitiesYH[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeVelocitiesZ(uint turbine) const { return &this->bladeVelocitiesZH[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeForcesX(uint turbine) const { return &this->bladeForcesXH[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeForcesY(uint turbine) const { return &this->bladeForcesYH[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeForcesZ(uint turbine) const { return &this->bladeForcesZH[turbine*numberOfNodesPerTurbine]; };
+
+ real* getAllBladeCoordsXDevice() const { return this->bladeCoordsXDCurrentTimestep; };
+ real* getAllBladeCoordsYDevice() const { return this->bladeCoordsYDCurrentTimestep; };
+ real* getAllBladeCoordsZDevice() const { return this->bladeCoordsZDCurrentTimestep; };
+ real* getAllBladeVelocitiesXDevice() const { return this->bladeVelocitiesXDCurrentTimestep; };
+ real* getAllBladeVelocitiesYDevice() const { return this->bladeVelocitiesYDCurrentTimestep; };
+ real* getAllBladeVelocitiesZDevice() const { return this->bladeVelocitiesZDCurrentTimestep; };
+ real* getAllBladeForcesXDevice() const { return this->bladeForcesXDCurrentTimestep; };
+ real* getAllBladeForcesYDevice() const { return this->bladeForcesYDCurrentTimestep; };
+ real* getAllBladeForcesZDevice() const { return this->bladeForcesZDCurrentTimestep; };
+
+ real* getTurbineBladeCoordsXDevice(uint turbine) const { return &this->bladeCoordsXDCurrentTimestep[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeCoordsYDevice(uint turbine) const { return &this->bladeCoordsYDCurrentTimestep[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeCoordsZDevice(uint turbine) const { return &this->bladeCoordsZDCurrentTimestep[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeVelocitiesXDevice(uint turbine) const { return &this->bladeVelocitiesXDCurrentTimestep[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeVelocitiesYDevice(uint turbine) const { return &this->bladeVelocitiesYDCurrentTimestep[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeVelocitiesZDevice(uint turbine) const { return &this->bladeVelocitiesZDCurrentTimestep[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeForcesXDevice(uint turbine) const { return &this->bladeForcesXDCurrentTimestep[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeForcesYDevice(uint turbine) const { return &this->bladeForcesYDCurrentTimestep[turbine*numberOfNodesPerTurbine]; };
+ real* getTurbineBladeForcesZDevice(uint turbine) const { return &this->bladeForcesZDCurrentTimestep[turbine*numberOfNodesPerTurbine]; };
+
+ void setAllBladeCoords(const real* _bladeCoordsX, const real* _bladeCoordsY, const real* _bladeCoordsZ);
+ void setAllBladeVelocities(const real* _bladeVelocitiesX, const real* _bladeVelocitiesY, const real* _bladeVelocitiesZ);
+ void setAllBladeForces(const real* _bladeForcesX, const real* _bladeForcesY, const real* _bladeForcesZ);
+
+ void setTurbineBladeCoords(uint turbine, const real* _bladeCoordsX, const real* _bladeCoordsY, const real* _bladeCoordsZ);
+ void setTurbineBladeVelocities(uint turbine, const real* _bladeVelocitiesX, const real* _bladeVelocitiesY, const real* _bladeVelocitiesZ);
+ void setTurbineBladeForces(uint turbine, const real* _bladeForcesX, const real* _bladeForcesY, const real* _bladeForcesZ);
+
+ void setTurbineAzimuth(uint turbine, real azimuth){azimuths[turbine] = azimuth;}
+
+ virtual void updateForcesAndCoordinates()=0;
private:
void initTurbineGeometries(CudaMemoryManager* cudaManager);
@@ -149,20 +144,10 @@ private:
void initBladeVelocities(CudaMemoryManager* cudaManager);
void initBladeForces(CudaMemoryManager* cudaManager);
void initBladeIndices(Parameter* para, CudaMemoryManager* cudaManager);
-
- void calcForcesEllipticWing();
- void rotateBlades(real angle, uint turbineID);
-
- void writeBladeCoords(uint t);
- void writeBladeCoordsToVtkForDebug(const std::string& filename);
- void writeBladeForces(uint t);
- void writeBladeVelocities(uint t);
-
+ std::string getFilename(Parameter* para, uint t) const;
void swapDeviceArrays();
public:
- real* bladeRadiiH;
- real* bladeRadiiD;
real* bladeCoordsXH, * bladeCoordsYH, * bladeCoordsZH;
real* bladeCoordsXDPreviousTimestep, * bladeCoordsYDPreviousTimestep, * bladeCoordsZDPreviousTimestep;
real* bladeCoordsXDCurrentTimestep, * bladeCoordsYDCurrentTimestep, * bladeCoordsZDCurrentTimestep;
@@ -176,23 +161,30 @@ public:
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;
+ real* turbinePosXH, *turbinePosYH, *turbinePosZH;
+ real* turbinePosXD, *turbinePosYD, *turbinePosZD;
+
+protected:
+ static constexpr uint numberOfBlades{3};
+ std::vector<real> bladeRadii, initialTurbinePositionsX, initialTurbinePositionsY, initialTurbinePositionsZ;
+ std::vector<real> azimuths;
+ const real diameter;
const bool useHostArrays;
const real density;
- real deltaT, deltaX;
- const uint numberOfBladeNodes, numberOfBlades;
- uint numberOfTurbines;
- const real epsilon; // in m
+ const real deltaT, deltaX, invDeltaX;
+ const uint numberOfTurbines, numberOfNodesPerBlade, numberOfNodesPerTurbine;
+ const real smearingWidth; // in m
const int level;
- uint numberOfIndices;
- uint numberOfGridNodes;
- real forceRatio, factorGaussian, invEpsilonSqrd, invDeltaX;
+ uint numberOfIndices{0};
+ uint numberOfGridNodes{0};
+
+ real forceRatio, factorGaussian;
int streamIndex;
+
+ bool writeOutput{false};
+ std::string outputName;
+ uint tOut{0};
+ uint tStartOut{0};
};
#endif
diff --git a/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarmInlines.h b/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarmInlines.h
index cc151088db963f00d3d39af705c80d22337701c9..6a7113fa5f3a61a35598aa837a4cebe2975b2588 100644
--- a/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarmInlines.h
+++ b/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarmInlines.h
@@ -26,7 +26,7 @@
// You should have received a copy of the GNU General Public License along
// with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
-//! \file ActuatorFarm.cu
+//! \file ActuatorFarmInlines.cu
//! \ingroup PreCollisionInteractor
//! \author Henrik Asmuth, Henry Korb, Anna Wellmann
//======================================================================================
@@ -35,6 +35,10 @@
#define ACTUATOR_FARM_INLINES
#include "basics/DataTypes.h"
+#include "basics/constants/NumericConstants.h"
+#include "VirtualFluids_GPU/GPU/GeometryUtils.h"
+
+using namespace vf::basics::constant;
struct TurbineNodeIndex {
uint turbine;
@@ -57,11 +61,12 @@ __host__ __device__ __inline__ void calcTurbineBladeAndBladeNode(uint node, uint
{
// see https://git.rz.tu-bs.de/irmb/VirtualFluids_dev/-/merge_requests/248 for visualization
turbine = node / (numberOfNodesPerBlade * numberOfBlades);
- uint x_off = turbine * numberOfNodesPerBlade * numberOfBlades;
+ const uint x_off = turbine * numberOfNodesPerBlade * numberOfBlades;
blade = (node - x_off) / numberOfNodesPerBlade;
- uint y_off = numberOfNodesPerBlade * blade + x_off;
+ const uint y_off = numberOfNodesPerBlade * blade + x_off;
bladeNode = node - y_off;
}
+
__host__ __device__ __inline__ TurbineNodeIndex calcTurbineBladeAndBladeNode(uint node, uint numberOfNodesPerBlade, uint numberOfBlades)
{
uint turbine;
@@ -71,4 +76,47 @@ __host__ __device__ __inline__ TurbineNodeIndex calcTurbineBladeAndBladeNode(uin
return { /*.turbine = */ turbine, /*.blade = */ blade, /*.bladeNode = */ bladeNode }; // Designated initializers are a C++ 20 feature
}
+__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)
+{
+ rotateAboutX3D(azimuth, bladeCoordX_BF, bladeCoordY_BF, bladeCoordZ_BF, 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)
+{
+ invRotateAboutX3D(azimuth, bladeCoordX_GF, bladeCoordY_GF, bladeCoordZ_GF, bladeCoordX_BF, bladeCoordY_BF, bladeCoordZ_BF);
+}
+__host__ __device__ __inline__ real distSqrd(real distX, real distY, real distZ)
+{
+ return distX * distX + distY * distY + distZ * distZ;
+}
+
+__host__ __device__ __inline__ real getBoundingSphereRadius(real diameter, real smearingWidth)
+{
+ return c1o2 * diameter + c4o1 * smearingWidth;
+}
+
+__host__ __device__ __inline__ bool inBoundingSphere(real distX, real distY, real distZ, real diameter, real smearingWidth)
+{
+ const real boundingSphereRadius = getBoundingSphereRadius(diameter, smearingWidth);
+ return distSqrd(distX, distY, distZ) < boundingSphereRadius * boundingSphereRadius;
+}
+
+__host__ __device__ __inline__ real gaussianSmearing(real distX, real distY, real distZ, real epsilon, real factorGaussian)
+{
+ return factorGaussian * exp(-distSqrd(distX, distY, distZ) / (epsilon * epsilon));
+}
+
+__inline__ void swapArrays(real* &arr1, real* &arr2)
+{
+ real* tmp = arr1;
+ arr1 = arr2;
+ arr2 = tmp;
+}
+
#endif
diff --git a/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarmStandalone.cu b/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarmStandalone.cu
new file mode 100644
index 0000000000000000000000000000000000000000..228795b8d891bd09715f520be5ae981711c80fb7
--- /dev/null
+++ b/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarmStandalone.cu
@@ -0,0 +1,126 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// This file is part of VirtualFluids. VirtualFluids is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file ActuatorFarm.cu
+//! \ingroup PreCollisionInteractor
+//! \author Henrik Asmuth, Henry Korb
+//======================================================================================
+#include "ActuatorFarmInlines.h"
+#include "ActuatorFarmStandalone.h"
+
+#include <cuda.h>
+#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 "VirtualFluids_GPU/GPU/GeometryUtils.h"
+#include "LBM/GPUHelperFunctions/KernelUtilities.h"
+#include "Parameter/Parameter.h"
+#include "Parameter/CudaStreamManager.h"
+#include "DataStructureInitializer/GridProvider.h"
+#include "GPU/CudaMemoryManager.h"
+
+using namespace vf::basics::constant;
+
+
+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);
+ return blade_radii;
+}
+
+void ActuatorFarmStandalone::updateForcesAndCoordinates()
+{
+ const real lift_coefficient = c1o1;
+ const real drag_coefficient = c0o1;
+ const real c0 = 20 * c1o10;
+ const real delta_azimuth = c2Pi / this->numberOfBlades;
+
+ 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;
+
+ 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);
+
+ real u_rel, v_rel, w_rel;
+ rotateFromGlobalToBlade(u_rel, v_rel, w_rel,
+ this->bladeVelocitiesXH[node],
+ 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;
+
+ const real dr = c1o2 * (next_node_radius - last_node_radius);
+
+ v_rel += current_node_radius * rotor_speed;
+ const real u_rel_sq = u_rel * u_rel + v_rel * v_rel;
+ const real phi = atan2(u_rel, v_rel);
+
+ const real tmp = c4o1 * current_node_radius / this->diameter - c1o1;
+ const real chord = c0 * sqrt(c1o1 - tmp * tmp);
+ const real normal_coefficient = lift_coefficient * cos(phi) + drag_coefficient * sin(phi);
+ const real tangential_coefficient = lift_coefficient * sin(phi) - drag_coefficient * cos(phi);
+ const real fx = -c1o2 * u_rel_sq * chord * this->density * normal_coefficient * dr;
+ const real fy = -c1o2 * u_rel_sq * chord * this->density * tangential_coefficient * dr;
+
+ rotateFromBladeToGlobal(fx, fy, c0o1,
+ this->bladeForcesXH[node], this->bladeForcesYH[node], this->bladeForcesZH[node],
+ local_azimuth_new);
+ rotateFromBladeToGlobal(c0o1, c0o1, current_node_radius,
+ this->bladeCoordsXH[node], this->bladeCoordsYH[node], this->bladeCoordsZH[node],
+ local_azimuth_new);
+ bladeCoordsXH[node] += this->turbinePosXH[turbine];
+ bladeCoordsYH[node] += this->turbinePosYH[turbine];
+ bladeCoordsZH[node] += this->turbinePosZH[turbine];
+ }
+ }
+ }
+}
+
diff --git a/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarmStandalone.h b/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarmStandalone.h
new file mode 100644
index 0000000000000000000000000000000000000000..78f0ebf67e0cfd60683be7f36f25bf97eb839035
--- /dev/null
+++ b/src/gpu/core/PreCollisionInteractor/Actuator/ActuatorFarmStandalone.h
@@ -0,0 +1,34 @@
+#ifndef ActuatorFarmStandalone_H
+#define ActuatorFarmStandalone_H
+
+#include "ActuatorFarm.h"
+#include "basics/DataTypes.h"
+
+class ActuatorFarmStandalone : public ActuatorFarm
+{
+public:
+ ActuatorFarmStandalone(
+ const real diameter,
+ const uint numberOfNodesPerBlade,
+ const std::vector<real> turbinePositionsX,
+ const std::vector<real> turbinePositionsY,
+ const std::vector<real> turbinePositionsZ,
+ const std::vector<real> rotorSpeeds,
+ const real density,
+ const real smearingWidth,
+ const int level,
+ const real deltaT,
+ const real deltaX
+ ) : rotorSpeeds(rotorSpeeds),
+ ActuatorFarm(diameter, computeBladeRadii(diameter, numberOfNodesPerBlade), turbinePositionsX, turbinePositionsY, turbinePositionsZ, density, smearingWidth, level, deltaT, deltaX, true)
+ {}
+
+ void updateForcesAndCoordinates() override;
+ static std::vector<real> computeBladeRadii(const real diameter, const uint numberOfNodesPerBlade);
+
+private:
+ std::vector<real> rotorSpeeds;
+
+};
+
+#endif