diff --git a/Python/actuator_line/actuator_line.py b/Python/actuator_line/actuator_line.py
index 7a29b1c3ad066193e1743953ae50b80b2e60dafe..ecd0fe0602bba83275798928fabce9339f20763e 100644
--- a/Python/actuator_line/actuator_line.py
+++ b/Python/actuator_line/actuator_line.py
@@ -1,5 +1,4 @@
#%%
-#%%
import numpy as np
from pathlib import Path
from mpi4py import MPI
@@ -73,8 +72,9 @@ t_out_probe = config.get_float_value("tOutProbe")
length = np.array([6,4,1])*boundary_layer_height
dx = boundary_layer_height/nodes_per_height
dt = dx * mach / (np.sqrt(3) * velocity)
-velocity_LB = velocity * dt / dx # LB units
-viscosity_LB = viscosity * dt / (dx * dx) # LB units
+velocity_ratio = dx/dt
+velocity_LB = velocity / velocity_ratio # LB units
+viscosity_LB = viscosity / (velocity_ratio * dx) # LB units
pressure_gradient = u_star * u_star / boundary_layer_height
pressure_gradient_LB = pressure_gradient * (dt*dt)/dx
@@ -85,8 +85,7 @@ logger.vf_log_info(f"viscosity [10^8 dx^2/dt] = {viscosity_LB*1e8}")
logger.vf_log_info(f"u* /(dx/dt) = {u_star*dt/dx}")
logger.vf_log_info(f"dpdx = {pressure_gradient}")
logger.vf_log_info(f"dpdx /(dx/dt^2) = {pressure_gradient_LB}")
-
-#%%
+
#%%
para.set_output_prefix(sim_name)
@@ -127,16 +126,20 @@ if read_precursor:
bc_factory.set_stress_boundary_condition(gpu.StressBC.StressPressureBounceBack)
bc_factory.set_slip_boundary_condition(gpu.SlipBC.SlipBounceBack)
bc_factory.set_pressure_boundary_condition(gpu.PressureBC.OutflowNonReflective)
-bc_factory.set_precursor_boundary_condition(gpu.PrecursorBC.DistributionsPrecursor if use_distributions else gpu.PrecursorBC.VelocityPrecursor)
+if read_precursor:
+ bc_factory.set_precursor_boundary_condition(gpu.PrecursorBC.DistributionsPrecursor if use_distributions else gpu.PrecursorBC.VelocityPrecursor)
para.set_outflow_pressure_correction_factor(0.0);
#%%
-def init_func(coord_x, coord_y, coord_z):
- return [
- 0.0,
- (u_star/0.4 * np.log(np.maximum(coord_z,z0)/z0) + 2.0*np.sin(np.pi*16*coord_x/length[0])*np.sin(np.pi*8*coord_z/boundary_layer_height)/(np.square(coord_z/boundary_layer_height)+1)) * dt / dx,
- 2.0*np.sin(np.pi*16.*coord_x/length[0])*np.sin(np.pi*8.*coord_z/boundary_layer_height)/(np.square(coord_z/boundary_layer_height)+1.) * dt / dx,
- 8.0*u_star/0.4*(np.sin(np.pi*8.0*coord_y/boundary_layer_height)*np.sin(np.pi*8.0*coord_z/boundary_layer_height)+np.sin(np.pi*8.0*coord_x/length[0]))/(np.square(length[2]/2.0-coord_z)+1.) * dt / dx]
-para.set_initial_condition(init_func)
+# don't use python init functions, they are very slow! Just kept as an example.
+# Define lambda in bindings and set it here.
+# def init_func(coord_x, coord_y, coord_z):
+# return [
+# 0.0,
+# (u_star/0.4 * np.log(np.maximum(coord_z,z0)/z0) + 2.0*np.sin(np.pi*16*coord_x/length[0])*np.sin(np.pi*8*coord_z/boundary_layer_height)/(np.square(coord_z/boundary_layer_height)+1)) * dt / dx,
+# 2.0*np.sin(np.pi*16.*coord_x/length[0])*np.sin(np.pi*8.*coord_z/boundary_layer_height)/(np.square(coord_z/boundary_layer_height)+1.) * dt / dx,
+# 8.0*u_star/0.4*(np.sin(np.pi*8.0*coord_y/boundary_layer_height)*np.sin(np.pi*8.0*coord_z/boundary_layer_height)+np.sin(np.pi*8.0*coord_x/length[0]))/(np.square(length[2]/2.0-coord_z)+1.) * dt / dx]
+# para.set_initial_condition(init_func)
+para.set_initial_condition_perturbed_log_law(u_star, z0, length[0], length[2], boundary_layer_height, velocity_ratio)
#%%
turb_pos = np.array([3,3,3])*turbine_diameter
@@ -145,7 +148,7 @@ density = 1.225
level = 0
n_blades = 3
n_blade_nodes = 32
-alm = gpu.ActuatorLine(n_blades, density, n_blade_nodes, epsilon, *turb_pos, turbine_diameter, level, dt, dx)
+alm = gpu.ActuatorLine(n_blades, density, n_blade_nodes, epsilon, *turb_pos, turbine_diameter, level, dt, dx, True)
para.add_actuator(alm)
#%%
planar_average_probe = gpu.probes.PlanarAverageProbe("horizontalPlanes", para.get_output_path(), 0, int(t_start_tmp_averaging/dt), int(t_averaging/dt) , int(t_start_out_probe/dt), int(t_out_probe/dt), 'z')
diff --git a/Python/actuator_line/configActuatorLinetxt b/Python/actuator_line/configActuatorLine.txt
similarity index 96%
rename from Python/actuator_line/configActuatorLinetxt
rename to Python/actuator_line/configActuatorLine.txt
index 45f5d4b8adfaa5ed0f8cb4155782478984bfae7c..4b38bc41c1d5b510e7b262423fff861dc1a9c030 100644
--- a/Python/actuator_line/configActuatorLinetxt
+++ b/Python/actuator_line/configActuatorLine.txt
@@ -23,8 +23,9 @@ Ma = 0.1
nz = 96
bodyForce = true
-UseAMD = true
-SGSconstant = 0.2
+SGSconstant = 0.333
+TurbulenceModel = QR
+
QuadricLimiterP = 100000.0
QuadricLimiterM = 100000.0
QuadricLimiterD = 100000.0
diff --git a/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp b/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
index 302bc69a429d342f5ad4b82757c76df8a7df2029..038f444f4489afcffecd91e29028d8f6154a456e 100644
--- a/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
+++ b/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
@@ -188,10 +188,10 @@ void multipleLevel(const std::string& configPath)
uint nBlades = 3;
uint nBladeNodes = 32;
- SPtr<ActuatorLine> actuator_line =SPtr<ActuatorLine>( new ActuatorLine(nBlades, density, nBladeNodes, epsilon, turbPos[0], turbPos[1], turbPos[2], reference_diameter, level, dt, dx) );
+ SPtr<ActuatorLine> actuator_line = std::make_shared<ActuatorLine>(nBlades, density, nBladeNodes, epsilon, turbPos[0], turbPos[1], turbPos[2], reference_diameter, level, dt, dx, true);
para->addActuator( actuator_line );
- SPtr<PointProbe> pointProbe = SPtr<PointProbe>( new PointProbe("pointProbe", para->getOutputPath(), 100, 1, 500, 100) );
+ SPtr<PointProbe> pointProbe = std::make_shared<PointProbe>("pointProbe", para->getOutputPath(), 100, 1, 500, 100);
std::vector<real> probeCoordsX = {reference_diameter,2*reference_diameter,5*reference_diameter};
std::vector<real> probeCoordsY = {3*reference_diameter,3*reference_diameter,3*reference_diameter};
std::vector<real> probeCoordsZ = {3*reference_diameter,3*reference_diameter,3*reference_diameter};
@@ -202,7 +202,7 @@ void multipleLevel(const std::string& configPath)
pointProbe->addStatistic(Statistic::Variances);
para->addProbe( pointProbe );
- SPtr<PlaneProbe> planeProbe = SPtr<PlaneProbe>( new PlaneProbe("planeProbe", para->getOutputPath(), 100, 500, 100, 100) );
+ SPtr<PlaneProbe> planeProbe = std::make_shared<PlaneProbe>("planeProbe", para->getOutputPath(), 100, 500, 100, 100);
planeProbe->setProbePlane(5*reference_diameter, 0, 0, dx, L_y, L_z);
planeProbe->addStatistic(Statistic::Means);
para->addProbe( planeProbe );
diff --git a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.cu b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.cu
index f2257edb0b017b0e253e970cd48c12ff5120066e..60dd7d3b581a102ad7b9c77f9eb6fb9a56f64bd7 100644
--- a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.cu
+++ b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.cu
@@ -188,9 +188,9 @@ __global__ void applyBodyForces(real* gridCoordsX, real* gridCoordsY, real* grid
}
}
- atomicAdd(&gridForcesX[gridIndex], gridForceX_RF);
- atomicAdd(&gridForcesY[gridIndex], gridForceY_RF);
- atomicAdd(&gridForcesZ[gridIndex], gridForceZ_RF);
+ gridForcesX[gridIndex] = gridForceX_RF;
+ gridForcesY[gridIndex] = gridForceY_RF;
+ gridForcesZ[gridIndex] = gridForceZ_RF;
}
@@ -210,7 +210,7 @@ void ActuatorLine::interact(Parameter* para, CudaMemoryManager* cudaMemoryManage
{
if (level != this->level) return;
- cudaMemoryManager->cudaCopyBladeCoordsHtoD(this);
+ if(useHostArrays) cudaMemoryManager->cudaCopyBladeCoordsHtoD(this);
vf::cuda::CudaGrid bladeGrid = vf::cuda::CudaGrid(para->getParH(level)->numberofthreads, this->nNodes);
@@ -225,11 +225,11 @@ void ActuatorLine::interact(Parameter* para, CudaMemoryManager* cudaMemoryManage
this->turbinePosX, this->turbinePosY, this->turbinePosZ,
this->bladeIndicesD, para->getVelocityRatio(), this->invDeltaX);
- cudaMemoryManager->cudaCopyBladeVelocitiesDtoH(this);
+ if(useHostArrays) cudaMemoryManager->cudaCopyBladeVelocitiesDtoH(this);
this->calcBladeForces();
- cudaMemoryManager->cudaCopyBladeForcesHtoD(this);
+ if(useHostArrays) cudaMemoryManager->cudaCopyBladeForcesHtoD(this);
vf::cuda::CudaGrid sphereGrid = vf::cuda::CudaGrid(para->getParH(level)->numberofthreads, this->nIndices);
diff --git a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.h b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.h
index 73bbc5e4eef79cc2f399d3a3d61426be8798486f..a441387512cc86e83453d9a4689d541b17dfde0f 100644
--- a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.h
+++ b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorLine.h
@@ -22,7 +22,8 @@ public:
const real _diameter,
int _level,
const real _deltaT,
- const real _deltaX
+ const real _deltaX,
+ const bool _useHostArrays
) : nBlades(_nBlades),
density(_density),
nBladeNodes(_nBladeNodes),
@@ -30,6 +31,7 @@ public:
turbinePosX(_turbinePosX), turbinePosY(_turbinePosY), turbinePosZ(_turbinePosZ),
diameter(_diameter),
level(_level),
+ useHostArrays(_useHostArrays),
PreCollisionInteractor()
{
this->deltaT = _deltaT*exp2(-this->level);
@@ -124,6 +126,7 @@ public:
uint* boundingSphereIndicesD;
private:
+ const bool useHostArrays;
const real density;
real turbinePosX, turbinePosY, turbinePosZ;
real omega, azimuth, yaw, deltaT, deltaX, invDeltaX, forceRatio, factorGaussian, invEpsilonSqrd;