Homec4science
Diffusion wireles-gpu (master) Browse

wireles-gpu/src/turbine_gpu.cuf
61fde96d0eeemaster

Last Change
Raw File

turbine_gpu.cuf
View Options

module turbine_m
!===============================================================================
! compute turbine
!===============================================================================
! warning: chord and twist are hardcoded for wire01
! ---------------------------------------------------------------------------
! Coordinate system:
! Top view: BACK view:
! Z
! Y .-|-.
! \ | / | \
! \|_____ X | X-- |-->Y,
! \ \ | /
! |\ "-|-"
! |
! counter-clockwise for positive gamma and phi angle from y axis
! ---------------------------------------------------------------------------
! Staggering grid:
! --w-----w-- --w--
! | u | u | u |
! |--w-o|--w--|--w--|
! | u | u | u |
! --w-----w-- --w--
! | u | u | u |
! --w-----w-- --w--
use precision
use dimen
real(fp), dimension(n_phi), parameter :: cos_Tang = [(cos(real(idx3-1)*dang), idx3=1,n_phi)]
real(fp), dimension(n_phi), parameter :: sin_Tang = [(sin(real(idx3-1)*dang), idx3=1,n_phi)]
real(fp), dimension(n_r,n_turb), parameter :: dA = &
reshape([((r(idx2,idx1)*dang*dr(idx1), idx2=1,n_r), idx1=1,n_turb)],(/n_r,n_turb/))
integer, parameter, dimension(n_turb) :: turb_i = int(turb_x/dx)
integer, parameter, dimension(n_turb) :: turb_j = int(turb_y/dy)
integer, parameter, dimension(n_turb) :: turb_k = int(turb_z/dz)
! Hard coded wire-01 data, to be replaced
real(fp), dimension(16), parameter :: radius_dat = &
(/0.0,0.0075,0.0125,0.0175,0.0226,0.0275,0.0325,&
0.0375,0.0425,0.0476,0.0525,0.0575,0.0625,0.0675,0.0726,0.075/)
real(fp), dimension(16), parameter :: twist_dat = &
(/0.7016,0.7016,0.5525,0.4389,0.3559,0.2949,0.2458,&
0.2085,0.1797,0.1542,0.1339,0.1187,0.1051,0.0932,0.0831,0.0831/)*rad2deg
real(fp), dimension(16), parameter :: chord_dat = &
(/0.0138,0.0138,0.0163,0.0176,0.0186,0.0195,&
0.0193,0.0174,0.0155,0.0141,0.0127,0.0115,0.0104,0.0096,0.0088,0.0088/)
contains
Subroutine turbine_gpu(Fx,Fy,Fz,u,v,w,alpha_cl_cd,t,me,nall)
!---------------------------------------------------------------------------
! declaration
!---------------------------------------------------------------------------
implicit none
integer*4 :: t,me,nall
real(fp),dimension(nx,ny,nz2),intent(out),device :: fx,fy,fz
real(fp),dimension(nx,ny,nz2),intent(in),device :: u,v,w
real(fp),dimension(3,181),intent(in),device :: alpha_cl_cd
integer :: i,j,k
!---------------------------------------------------------------------------
! main code
!---------------------------------------------------------------------------
fx = 0._fp
fy = 0._fp
fz = 0._fp
!---------------------------------------------------------------------------
! compute models
if(nacelle_model == 1)then
call nacelle_gpu(fx,fy,fz,u,v,w,t,me,nall)
end if
if(turbine_model == 1)then
call admryaw_gpu(fx,fy,fz,u,v,w,alpha_cl_cd,t,me,nall)
end if
end Subroutine turbine_gpu
!===============================================================================
! compute ADMR for yaw turbine
! ===============================================================================
Subroutine admryaw_gpu(fx,fy,fz,u,v,w,alpha_cl_cd,t,me,nall)
!---------------------------------------------------------------------------
! declaration
!---------------------------------------------------------------------------
use cudafor
implicit none
real(fp),dimension(nx,ny,nz2),intent(inout),device :: fx,fy,fz
real(fp),dimension(nx,ny,nz2),intent(in),device :: u,v,w
real(fp),dimension(3,181),intent(in),device :: alpha_cl_cd
integer,intent(in) :: t,me,nall
! Grid relative coordinates to the hub-center
! CTx,CTy,CTzu,CTzw,CTru,CTrw
real(fp),dimension(:,:,:,:,:),allocatable,device :: coord
! interpoltaion coefficient
! uv node, w node
real(fp),dimension(:,:,:,:,:),allocatable,device :: coeff
real(fp),dimension(:,:,:,:),allocatable,device :: sin_CTangx,cos_CTangz
! blade element center relative coordinates to the hub-center
! mu_x,mu_y,mu_z
real(fp),dimension(:,:,:,:),allocatable,device :: mu
! Normalization for the smearing kernel
! sum_ker_uvv,sum_ker_w
real(fp),dimension(:,:,:,:),allocatable,device :: sum_ker
! Closest grid point for a blade element
! blade_i,blade_j,blade_k_uv,blade_k_w
integer,dimension(:,:,:,:),allocatable,device :: blade_ijk
! Force on the blade element
! CTFx,CTFt
real(fp),dimension(:,:,:,:),allocatable,device :: CTF
real(fp),dimension(:,:),allocatable,device :: twist,CTsoli
real(fp),dimension(n_turb),device :: inflow,omega
integer :: i,j,k,me_k,index,x_idx_global,hub_node
type(dim3) :: grid_be,tBlock_be,&
grid_full,tBlock_full
! kernel config for force computing on blade element
tBlock_be = dim3(8,8,1)
grid_be = dim3(n_phi/8,n_r/8,n_turb)
! kernel config for force prejection on cartesian grid
tBlock_full = dim3(ker_d_x,ker_d_y,ker_d_z)
grid_full = dim3(n_phi,n_r,n_turb)
save sin_CTangx,cos_CTangz,coeff,sum_ker,blade_ijk,CTF,twist,CTsoli
!---------------------------------------------------------------------------
! pre-process for interp and project
!---------------------------------------------------------------------------
if (t == 1) then
allocate(sin_CTangx(mask_x,mask_y,mask_z,n_turb))
allocate(cos_CTangz(mask_x,mask_y,mask_z,n_turb))
allocate(coeff(n_phi,n_r,n_turb,2,8))
allocate(sum_ker(n_phi,n_r,n_turb,2))
allocate(blade_ijk(n_phi,n_r,n_turb,4))
allocate(CTF(n_phi,n_r,n_turb,2))
allocate(twist(n_r,n_turb))
allocate(CTsoli(n_r,n_turb))
! On the axi-symmetric blade elements:
! Pre-compute the interpolation coefficient on the blade element
call force_compute_pre<<<grid_be,tBlock_be>>>(coeff,blade_ijk,twist,CTsoli)
! Pre-compute some smearing parameter
sum_ker = 0
call force_projection_pre<<<grid_full,tBlock_full>>>(sum_ker)
end if
! ! Compute the force on the blade elements
call force_compute<<<grid_be,tBlock_be>>>(CTF,&
u,v,w,inflow,omega,&
coeff,blade_ijk,twist,CTsoli,&
alpha_cl_cd)
! Distribute the force from the blade element to the cartesian grid
call force_projection<<<grid_full,tBlock_full>>>(fx,fy,fz,CTF,sum_ker)
end Subroutine admryaw_gpu
! Subroutine nacelle_gpu(Fx,Fy,Fz,u,v,w,t,me,nall)
! !---------------------------------------------------------------------------
! ! declaration
! !---------------------------------------------------------------------------
! implicit none
! include 'dimen.h'
! real*8, dimension(nx,ny,nz2), intent(inout) :: fx,fy,fz
! real*8, dimension(nx,ny,nz2), intent(in) :: u,v,w
! integer*4, intent(in) :: t
! integer :: i,j,k,i0
! real*8 :: hub_kk,sigma,tmp1,tmp2,tmp3
! real*8 :: CTU,CTFx
! real*8, dimension(n_turb) :: sum_ker
! ! csts
! real*8, dimension(n_turb), parameter :: dA = &
! [(pi*nac_r(i_turb)**2_fp, i_turb=1,n_turb)]
! real*8, dimension(n_turb) :: hub_k,lz1,lz2
! integer, dimension(n_turb) :: me_flag,me_k
! real*8, dimension(ny,nz2,n_turb) :: ker,flag_j
! save hub_k,lz1,lz2,me_k,me_flag,ker,flag_j
! !---------------------------------------------------------------------------
! ! pre-process for interp and project
! !---------------------------------------------------------------------------
! if (t == 1) then
! !-----------------------------------------------------------------------
! ! interp
! me_flag = 0
! do i_turb=1,n_turb
! hub_kk= real(turb_z(i_turb)/dz)-0.5_fp ! global uv node index of hub
! me_k(i_turb)=hub_kk/Nzb ! processor of hub
! ! if (me == me_k(i_turb)) then
! me_flag(i_turb) = 1
! hub_kk = hub_kk-Nzb*me+2._fp ! local uv node index of hub
! hub_k(i_turb) = int(hub_kk)
! lz1(i_turb) = real(hub_kk-hub_k(i_turb))
! lz2(i_turb) = (1._fp_fp-lz1(i_turb))
! ! endif
! enddo
! !-----------------------------------------------------------------------
! ! project
! ker = 0._fp
! flag_j = 0
! !sigma = 0.75_fp*max(dy,dz)
! sigma = 0.5_fp*max(dy,dz)
! tmp1 = 1._fp/(sigma*dsqrt(pi))**2_fp
! tmp2 = -1._fp/sigma**2_fp
! do i_turb=1,n_turb
! do k=2,nzb+1
! do j=1,ny
! tmp3 = (real(turb_j(i_turb)-j)*dy)**2_fp &
! + (real(me*nzb+k-1)*dz-turb_z(i_turb))**2_fp
! if( sqrt(tmp3) <= 4._fp*sigma ) then
! ker(j,k,i_turb) = tmp1*dexp(tmp2*tmp3)
! endif
! end do
! end do
! sum_ker(i_turb) = sum(ker(:,:,i_turb))
! enddo
! ! call mpi_allreduce(mpi_in_place, sum_ker(1),&
! ! n_turb,mpi_double_precision,&
! ! mpi_sum,nall,ierr)
! do i_turb=1,n_turb
! do k=2,nzb+1
! do j=1,ny
! ker(j,k,i_turb) = ker(j,k,i_turb)/sum_ker(i_turb)
! if (ker(j,k,i_turb) >= 1E-15) then
! flag_j(j,k,i_turb) =1
! endif
! end do
! end do
! end do
! endif
! !---------------------------------------------------------------------------
! ! compute
! !---------------------------------------------------------------------------
! !---------------------------------------------------------------------------
! ! compute force (and vel)
! do i_turb=1,n_turb
! i = turb_i(i_turb)
! i0 = 4._fp*turb_r(i_turb)/dx
! if (me_flag(i_turb) == 1) then
! !-------------------------------------------------------------------
! ! compute vel
! CTU = u(i-i0,turb_j(i_turb),hub_k(i_turb))*lz1(i_turb) + &
! u(i-i0,turb_j(i_turb),hub_k(i_turb)+1)*lz2(i_turb)
! !-------------------------------------------------------------------
! ! compute force from drag law
! CTFx = 0.5_fp*(CTU**2_fp)*nac_c(i_turb)*dA(i_turb)/(dx*dy*dz)
! ! thrust = CTFx*(dx*dy*dz)*z_i**2
! endif
! !-----------------------------------------------------------------------
! ! project
! do k=2,nzb+1
! do j=1,ny
! if (flag_j(j,k,i_turb) == 1) then
! Fx(i,j,k) = Fx(i,j,k) + CTFx*ker(j,k,i_turb)
! ! print *, 'nac: me,i,j,k,Fx(i,j,k):',me,i,j,k,Fx(i,j,k)
! endif
! end do
! end do
! enddo
! end Subroutine nacelle_gpu
!===============================================================================
! funct lib
!===============================================================================
! compute_ang
attributes(device) subroutine compute_ang(angle,coord_y,coord_z)
implicit none
real(fp),intent(in) :: coord_y,coord_z
real(fp),intent(out) :: angle
real(fp) :: gau,RR
RR = sqrt(coord_y**2+coord_z**2)
if (RR <= 1E-16) then
angle = 0._fp
else
gau = abs(coord_z/sqrt(coord_y**2+coord_z**2))
if (coord_z >= 0._fp)then
if(coord_y >= 0._fp) angle = asin(gau)
if(coord_y < 0._fp) angle = PI-abs(asin(gau))
else
if(coord_y < 0._fp) angle = PI+abs(asin(gau))
if(coord_y >= 0._fp) angle = 2._fp*PI-asin(gau)
end if
end if
end subroutine compute_ang
attributes(global) subroutine force_compute_pre(coeff,blade_ijk,twist,CTsoli)
implicit none
integer,dimension(n_phi,n_r,n_turb,4),intent(out) :: blade_ijk
real(fp),dimension(n_phi,n_r,n_turb,2,8) ,intent(out) :: coeff
real(fp),dimension(n_r,n_turb),intent(out) :: twist,CTsoli
real(fp),dimension(8) :: c_w_1d,c_uv_1d
real(fp) :: blade_ii,blade_jj,blade_kk
real(fp) :: tmp3,chordl
integer :: i,j,k,i_turb,i_phi,i_r,me_k
i_phi = (blockIdx%x - 1) * blockDim%x + threadIdx%x
i_r = (blockIdx%y - 1) * blockDim%y + threadIdx%y
i_turb = (blockIdx%z - 1) * blockDim%z + threadIdx%z
!-----------------------------------------------------------
! pre-process for interp
! WiRE-01
call itp1D( radius_dat,chord_dat,size(radius_dat),r(i_r,i_turb)*z_i,chordl )
! Chord length
chordl = chordl/z_i
! Twist angle
call itp1D( radius_dat,twist_dat,size(radius_dat),r(i_r,i_turb)*z_i,twist(i_r,i_turb) )
! Solidity
CTsoli(i_r,i_turb) = 3._fp*chordl/(2._fp*pi*r(i_r,i_turb))
! blade element centre location
blade_ii = turb_i(i_turb)-0.5_fp-r(i_r,i_turb)*sin(yaw_angle(i_turb))*cos_Tang(i_phi)/dx
blade_jj = turb_j(i_turb)-0.5_fp+r(i_r,i_turb)*cos(yaw_angle(i_turb))*cos_Tang(i_phi)/dy
blade_kk = turb_k(i_turb)-0.5_fp+r(i_r,i_turb)*sin_Tang(i_phi)/dz
! print *,int(floor(blade_ii)),int(floor(blade_jj)),int(floor(blade_kk))
! Global index of the left_lower grid to a blade element centre
blade_ijk(i_phi,i_r,i_turb,1) = int(floor(blade_ii))
blade_ijk(i_phi,i_r,i_turb,2) = int(floor(blade_jj))
blade_ijk(i_phi,i_r,i_turb,3) = int(floor(blade_kk))
! Trilinear interpolation coefficients
! Interpolation at o for w node needs to be treated with care
call get_interp_coeff(c_uv_1d,blade_ii,blade_jj,blade_kk)
if ((blade_kk - int(floor(blade_kk))) < 0.5_fp) then
blade_ijk(i_phi,i_r,i_turb,4) = int(floor(blade_kk))
call get_interp_coeff(c_w_1d,blade_ii,blade_jj,blade_kk)
else ! For the case when interpolation is in the upper layer
blade_ijk(i_phi,i_r,i_turb,4) = int(floor(blade_kk))+1
call get_interp_coeff(c_w_1d,blade_ii,blade_jj,blade_kk+1)
end if
coeff(i_phi,i_r,i_turb,1,:) = c_uv_1d(:)
coeff(i_phi,i_r,i_turb,2,:) = c_w_1d(:)
end subroutine force_compute_pre
attributes(global) subroutine force_projection_pre(sum_ker)
implicit none
real(fp),dimension(n_phi,n_r,n_turb,2),intent(inout) :: sum_ker
real(fp) :: sigma,tmp1,tmp2,tmp3,ker_uv,ker_w
real(fp) :: mu_x,mu_y,mu_z
real(fp) :: coord_x,coord_y,coord_z_uv,coord_z_w,coord_r_uv,coord_r_w
integer :: i,j,k,i_turb,i_phi,i_r,istat
i = threadIdx%x - int(0.5_fp*ker_d_x)
j = threadIdx%y - int(0.5_fp*ker_d_y)
k = threadIdx%z - int(0.5_fp*ker_d_z)
i_phi = blockIdx%x
i_r = blockIdx%y
i_turb = blockIdx%z
!-----------------------------------------------------------
! pre-process for force projection
! Blade element relative coordinates to the hub
if (i_r >= n_r_s) then
mu_x = (turb_i(i_turb)+0.5_fp)*dx+r(i_r,i_turb)*cos_Tang(i_phi)*sin(yaw_angle(i_turb))
mu_y = (turb_j(i_turb)+0.5_fp)*dy+r(i_r,i_turb)*cos_Tang(i_phi)*cos(yaw_angle(i_turb))
mu_z = (turb_k(i_turb)+0.5_fp)*dz+r(i_r,i_turb)*sin_Tang(i_phi)
sigma = sqrt((r(i_r,i_turb)*dang)**2+dr(i_turb)**2)!2*dx!1.5*sqrt((r(i_r,i_turb)*dang)**2+dr(i_turb)**2)
tmp1 = 1._fp/(sigma*sqrt(pi))**3
tmp2 = -1._fp/(2._fp*sigma**2)
tmp3 = (i*dx)**2+(j*dy)**2+(k*dz)**2
ker_uv = tmp1*exp(tmp2*tmp3)
tmp3 = (i*dx)**2+(j*dy)**2+(k*dz-0.5_fp*dz)**2
ker_w = tmp1*exp(tmp2*tmp3)
!-----------------------------------------------------------
! Loop through all kernel grid points to sum the smearing coeff.
! Also, we don't pre-computed the smearing coeff. due to
! memeory limit (avoiding saving a 6d array)
istat = atomicAdd(sum_ker(i_phi,i_r,i_turb,1),ker_uv)
istat = atomicAdd(sum_ker(i_phi,i_r,i_turb,2),ker_w )
end if
call syncthreads()
end subroutine force_projection_pre
attributes(global) subroutine force_compute(CTF,&
u,v,w,inflow,omega,&
coeff,blade_ijk,twist,CTsoli,&
alpha_cl_cd)
implicit none
real(fp),dimension(n_phi,n_r,n_turb,2),intent(out) :: CTF
real(fp),dimension(n_turb),intent(out) :: omega
real(fp),dimension(nx,ny,nz2),intent(in) :: u,v,w
real(fp),dimension(n_turb),intent(in) :: inflow
real(fp),dimension(n_phi,n_r,n_turb,2,8),intent(in) :: coeff
integer,dimension(n_phi,n_r,n_turb,4),intent(in) :: blade_ijk
real(fp),dimension(n_r,n_turb),intent(in) :: twist,CTsoli
real(fp),dimension(3,181),intent(in) :: alpha_cl_cd
real(fp),dimension(8) :: c_w_1d,c_uv_1d
real(fp) :: CTphi,sin_CTphi,cos_CTphi
real(fp) :: CTU,CTV,CTW,CTUrel,AoA
real(fp) :: CTCL,CTCD,CTU_origin,CTV_origin,CTW_origin
real(fp) :: f1,f2
real(fp) :: temp
integer :: i_turb,i_phi,i_r
i_phi = (blockIdx%x - 1) * blockDim%x + threadIdx%x
i_r = (blockIdx%y - 1) * blockDim%y + threadIdx%y
i_turb = (blockIdx%z - 1) * blockDim%z + threadIdx%z
if (i_r >= n_r_s) then
c_uv_1d = coeff(i_phi,i_r,i_turb,1,:)
c_w_1d = coeff(i_phi,i_r,i_turb,2,:)
! TODO: read inflow
omega(i_turb) = 3.8*5.0/(turb_r(i_turb)*z_i)
call compute_vel_3D(CTU_origin,u,&
blade_ijk(i_phi,i_r,i_turb,1),&
blade_ijk(i_phi,i_r,i_turb,2),&
blade_ijk(i_phi,i_r,i_turb,3),c_uv_1d)
call compute_vel_3D(CTV_origin,v,&
blade_ijk(i_phi,i_r,i_turb,1),&
blade_ijk(i_phi,i_r,i_turb,2),&
blade_ijk(i_phi,i_r,i_turb,3),c_uv_1d)
call compute_vel_3D(CTW_origin,w,&
blade_ijk(i_phi,i_r,i_turb,1),&
blade_ijk(i_phi,i_r,i_turb,2),&
blade_ijk(i_phi,i_r,i_turb,4),c_w_1d)
!-----------------------------------------------------------
! compute flow
! Project from the original grid-align reference to the disk normal reference
CTU = CTU_origin*cos(yaw_angle(i_turb))+CTV_origin*sin(yaw_angle(i_turb))
CTV = CTV_origin*cos(yaw_angle(i_turb))-CTU_origin*sin(yaw_angle(i_turb))
CTW = CTW_origin
! TODO: Check it
CTV = CTW*cos_Tang(i_phi)-CTV*sin_Tang(i_phi) + omega(i_turb)*(r(i_r,i_turb)*z_i)
if (CTV > 0) then ! exclude the negative incoming flow
CTphi = atan(CTU/CTV)
sin_CTphi = sin(CTphi)
cos_CTphi = cos(CTphi)
! AoA in degree
AoA = CTphi*rad2deg-twist(i_r,i_turb)
CTUrel = CTU/sin_CTphi
!-----------------------------------------------------------
! compute cl/cd from AoA (assume one Re)
if (AoA < -90._fp) then
AoA = -90._fp
elseif (AoA > 90._fp) then
AoA = 90._fp
endif
call compute_cl_cd(CTCL,CTCD,AOA,alpha_cl_cd)
!-----------------------------------------------------------
! compute BEM
! Tip loss factor
f1 = 0.5_fp*3._fp*(turb_r(i_turb)-r(i_r,i_turb))/(r(i_r,i_turb)*sin_CTphi)
f2 = 2._fp/pi*acos(exp(-f1))
if (AoA > -10 .and. AoA < 50) then ! Impose a stronger restriction
CTF(i_phi,i_r,i_turb,1) = 0.5_fp*(CTUrel**2)*CTsoli(i_r,i_turb)&
*dA(i_r,i_turb)/(dx*dy*dz)&
*(CTCL*cos_CTphi+CTCD*sin_CTphi)*f2
CTF(i_phi,i_r,i_turb,2) = 0.5_fp*(CTUrel**2)*CTsoli(i_r,i_turb)&
*dA(i_r,i_turb)/(dx*dy*dz)&
*(CTCL*sin_CTphi-CTCD*cos_CTphi)*f2
end if
end if
end if
end subroutine force_compute
attributes(global) subroutine force_projection(fx,fy,fz,CTF,sum_ker)
implicit none
real(fp),dimension(nx,ny,nz2),intent(inout) ::Fx,Fy,Fz
real(fp),dimension(n_phi,n_r,n_turb,2),intent(in) :: CTF
real(fp),dimension(n_phi,n_r,n_turb,2),intent(in) :: sum_ker
real(fp) :: sin_CTangx,cos_CTangz
real(fp) :: tmp3_a,sigma,tmp1,tmp2,tmp3
real(fp) :: ker_uv,ker_w
real(fp) :: coord_x,coord_y,coord_z_uv,coord_z_w,r_uv,r_w
real(fp) :: angle_uv,angle_z
real(fp) :: mu_x,mu_y,mu_z
real(fp) :: CTFy,CTFx_origin,CTFy_origin,CTFz_origin
integer :: i,j,k,i_turb,i_phi,i_r,x_idx_global,y_idx_global,z_idx_global,istat
! Center the 0
i = threadIdx%x - int(0.5*ker_d_x)
j = threadIdx%y - int(0.5*ker_d_y)
k = threadIdx%z - int(0.5*ker_d_z)
i_phi = blockIdx%x
i_r = blockIdx%y
i_turb = blockIdx%z
if (i_r >= n_r_s) then
mu_x = (turb_i(i_turb)+0.5)*dx + r(i_r,i_turb)*cos_Tang(i_phi)*sin(yaw_angle(i_turb))
mu_y = (turb_j(i_turb)+0.5)*dy + r(i_r,i_turb)*cos_Tang(i_phi)*cos(yaw_angle(i_turb))
mu_z = (turb_k(i_turb)+0.5 + 1)*dz + r(i_r,i_turb)*sin_Tang(i_phi) ! why +1: z = 2 to nzb+1
! global index for kernel
x_idx_global = nint(mu_x/dx) + i
y_idx_global = nint(mu_y/dy) + j
z_idx_global = nint(mu_z/dz) + k
! distance to hub cell
coord_x = (x_idx_global - turb_i(i_turb))*dx
coord_y = (y_idx_global - turb_j(i_turb))*dy
coord_z_uv = (z_idx_global - turb_k(i_turb) - 1._fp)*dz
coord_z_w = (z_idx_global - turb_k(i_turb) - 1.5_fp)*dz
! Distance to the hub
r_uv = sqrt(coord_x**2 + coord_y**2 + coord_z_uv**2)
r_w = sqrt(coord_x**2 + coord_y**2 + coord_z_w**2)
sin_CTangx = coord_z_uv/(r_uv+1e-16)
cos_CTangz = coord_y/(r_w+1e-16)
sigma = 1._fp*sqrt((r(i_r,i_turb)*dang)**2+dr(i_turb)**2)!2*dx!1.5*sqrt((r(i_r,i_turb)*dang)**2+dr(i_turb)**2)
tmp1 = 1._fp/(sigma*sqrt(pi))**3
tmp2 = -1._fp/(2._fp*sigma**2)
tmp3 = (i*dx)**2+(j*dy)**2+(k*dz)**2
ker_uv = tmp1*exp(tmp2*tmp3)/sum_ker(i_phi,i_r,i_turb,1)
tmp3 = (i*dx)**2+(j*dy)**2+(k*dz - 0.5_fp*dz)**2
ker_w = tmp1*exp(tmp2*tmp3)/sum_ker(i_phi,i_r,i_turb,2)
! ! ! Rotate the force from the disk-normal reference to the original reference
CTFy = - CTF(i_phi,i_r,i_turb,2)*sin_CTangx
CTFx_origin = CTF(i_phi,i_r,i_turb,1)*cos(yaw_angle(i_turb))-CTFy*sin(yaw_angle(i_turb))
CTFy_origin = CTF(i_phi,i_r,i_turb,1)*sin(yaw_angle(i_turb))+CTFy*cos(yaw_angle(i_turb))
CTFz_origin = CTF(i_phi,i_r,i_turb,2)*cos_CTangz
! ! Project the force
istat = atomicAdd(Fx(x_idx_global,y_idx_global,z_idx_global),CTFx_origin*ker_uv)
istat = atomicAdd(Fy(x_idx_global,y_idx_global,z_idx_global),CTFy_origin*ker_uv)
istat = atomicAdd(Fz(x_idx_global,y_idx_global,z_idx_global),CTFz_origin*ker_w)
end if
end subroutine force_projection
attributes(device) subroutine get_interp_coeff(c,blade_ii,blade_jj,blade_kk)
implicit none
real(fp),dimension(8),intent(out) :: c
real(fp),intent(in) :: blade_ii,blade_jj,blade_kk
real(fp),dimension(3) :: blade_point
integer :: i,j,k,counter
real(fp),dimension(3,2) :: l
blade_point = (/blade_ii,blade_jj,blade_kk/)
! Get the location in the cube
do i = 1,3
l(i,1) = 1._fp-(dble(blade_point(i)-int(floor(blade_point(i)))))
l(i,2) = dble(blade_point(i)-int(floor(blade_point(i))))
! if (i == 1) print *,i,l(i,1),l(i,2)
end do
! Get the coefficient
counter = 1
do i = 1,2
do j = 1,2
do k = 1,2
c(counter) = l(1,i)*l(2,j)*l(3,k)
counter = counter + 1
end do
end do
end do
end subroutine get_interp_coeff
attributes(device) subroutine compute_vel_3D(u_interp,u,blade_i,blade_j,blade_k,c)
implicit none
real(fp),intent(out) :: u_interp
real(fp),dimension(nx,ny,nz2),intent(in) :: u
integer,intent(in) :: blade_i,blade_j,blade_k
integer :: i,j,k,counter
real(fp),dimension(8),intent(in) :: c
! To be rewritten using dot_product
counter = 1
u_interp = 0
do i = 0,1
do j = 0,1
do k = 0,1
u_interp = u_interp + c(counter) * u(blade_i+i,blade_j+j,blade_k+k)
counter = counter + 1
end do
end do
end do
end subroutine compute_vel_3D
attributes(device) subroutine compute_cl_cd(CTCL,CTCD,x,alpha_cl_cd)
implicit none
real(fp) x,compute_cl,CTCL,CTCD
real(fp),dimension(3,181) :: alpha_cl_cd
integer i,j
i=floor(x)+91
j=ceiling(x)+91
CTCL = dble(x-floor(x))*alpha_cl_cd(2,j)+dble(ceiling(x)-x)*alpha_cl_cd(2,i)
CTCD = dble(x-floor(x))*alpha_cl_cd(3,j)+dble(ceiling(x)-x)*alpha_cl_cd(3,i)
end subroutine compute_cl_cd
attributes(device) subroutine itp1D(xData,yData,array_size,xVal,yVal)
implicit none
integer,intent(in) :: array_size
real(fp),intent(in) :: xData(array_size),yData(array_size)
real(fp),intent(in) :: xVal
real(fp),intent(out) :: yVal
integer :: dataIndex
real(fp) :: minXdata,maxXdata,xRange,weight
minXData = xData(1)
maxXData = xData(size(xData))
if (xVal<minXData .or. xVal>maxXData) then
print *,"Error in interpolation"
stop
else
do dataIndex = 1,array_size-1
if (xVal>xData(dataIndex) .and. xVal<xData(dataIndex+1)) then
weight = (xVal-xData(dataIndex))/(xData(dataIndex+1)-xData(dataIndex))
yVal = (1._fp-weight)*yData(dataIndex) + weight*yData(dataIndex+1)
exit
end if
end do
end if
end subroutine itp1D
end module turbine_m
c4science · Help