Skip to content

Commit

Permalink
Merge pull request #16 from jjonkman/f/ff-curl-review
Browse files Browse the repository at this point in the history 
Review of FAST.Farm Changes Regarding Polar to Rectangular Grid
  • Loading branch information
@ebranlard
ebranlard authored Jun 30, 2021
2 parents 1df2245 + 132dc3d commit 492cb4f
Show file tree
Hide file tree
Showing 4 changed files with 15 additions and 46 deletions.
15 changes: 4 additions & 11 deletions glue-codes/fast-farm/src/FASTWrapper.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -535,7 +535,6 @@ SUBROUTINE FWrap_CalcOutput(p, u, y, m, ErrStat, ErrMsg)
REAL(R8Ki) :: theta(3)
REAL(R8Ki) :: orientation(3,3)
REAL(ReKi) :: tmp_sz_z, tmp_sz_y
REAL(ReKi) :: R_rotor ! Rotor radius

INTEGER(IntKi) :: j, k ! loop counters

Expand Down Expand Up @@ -625,12 +624,6 @@ SUBROUTINE FWrap_CalcOutput(p, u, y, m, ErrStat, ErrMsg)
call setErrStat(ErrStat2,ErrMsg2,ErrStat2,ErrMsg,RoutineName)
if (ErrStat >= AbortErrLev) return
end do

! Rotor radius
R_rotor = 0.0_ReKi
do k=1,size(m%ADRotorDisk) ! loop on blades
R_rotor = max(R_rotor, TwoNorm(m%ADRotorDisk(k)%Position(:,p%nr) - p0) )
enddo

! --- Ct and Cq on polar grid (goes beyond rotor radius)
if (EqualRealNos(y%DiskAvg_Vx_Rel,0.0_ReKi)) then
Expand Down Expand Up @@ -658,14 +651,14 @@ SUBROUTINE FWrap_CalcOutput(p, u, y, m, ErrStat, ErrMsg)
do k=1,size(m%ADRotorDisk) ! loop on blades force contribution
num = num + dot_product(y%xHat_Disk, m%ADRotorDisk(k)%Moment(:,j) )
end do
y%AzimAvg_Cq(j) = num / (denom * R_rotor)
y%AzimAvg_Cq(j) = 2.0_ReKi*num / (denom * y%D_rotor)
end do

end if

! --- Variables needed to orient wake planes in "skew" coordinate system
! chi_skew and psi_skew
y%chi_skew = Calc_Chi0(m%Turbine%AD%m%rotors(1)%V_diskAvg, m%turbine%AD%m%rotors(1)%V_dot_x) * R2D ! AeroDyn_IO
y%chi_skew = Calc_Chi0(m%Turbine%AD%m%rotors(1)%V_diskAvg, m%turbine%AD%m%rotors(1)%V_dot_x) ! AeroDyn_IO

! TODO place me in an AeroDyn Function like Calc_Chi0
! Construct y_hat, orthogonal to x_hat when its z component is neglected (in a projected horizontal plane)
Expand All @@ -678,8 +671,8 @@ SUBROUTINE FWrap_CalcOutput(p, u, y, m, ErrStat, ErrMsg)
zHat_plane(1:3) = zHat_plane/TwoNorm(zHat_plane)

zHat_Disk = m%Turbine%AD%Input(1)%rotors(1)%HubMotion%Orientation(3,:,1) ! TODO TODO, shoudn't rotate
tmp_sz_y = dot_product(zHat_Disk,yHat_plane)
tmp_sz_z = -1.0* dot_product(zHat_Disk,zHat_plane)
tmp_sz_y = dot_product(zHat_Disk,yHat_plane)
tmp_sz_z = -1.0_ReKi* dot_product(zHat_Disk,zHat_plane)
if ( EqualRealNos(tmp_sz_y,0.0_ReKi) .and. EqualRealNos(tmp_sz_z,0.0_ReKi) ) then
y%psi_skew = 0.0_ReKi
else
Expand Down
16 changes: 5 additions & 11 deletions modules/awae/src/AWAE.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -210,8 +210,6 @@ subroutine LowResGridCalcOutput(n, u, p, y, m, errStat, errMsg)
real(ReKi) :: xplane_sq, yplane_sq, xysq_Z(3), xzplane_X(3)
real(ReKi), ALLOCATABLE :: tmp_xhat_plane(:,:), tmp_yhat_plane(:,:), tmp_zhat_plane(:,:)
real(ReKi), ALLOCATABLE :: tmp_Vx_wake(:), tmp_Vz_wake(:), tmp_Vy_wake(:)
real(ReKi) :: Vx_debug
integer(IntKi) :: ILo
integer(IntKi) :: maxPln, tmpPln, maxN_wake
integer(IntKi) :: i,np1,errStat2
character(*), parameter :: RoutineName = 'LowResGridCalcOutput'
Expand Down Expand Up @@ -242,7 +240,7 @@ subroutine LowResGridCalcOutput(n, u, p, y, m, errStat, errMsg)

! Loop over the entire grid of low resolution ambient wind data to compute:
! 1) the disturbed flow at each point and 2) the averaged disturbed velocity of each wake plane
!$OMP PARALLEL DO PRIVATE(nx_low,ny_low,nz_low, nXYZ_low, n_wake, xhatBar_plane, x_end_plane,nt,np,ILo,x_start_plane,delta,deltad,Vx_debug,&
!$OMP PARALLEL DO PRIVATE(nx_low,ny_low,nz_low, nXYZ_low, n_wake, xhatBar_plane, x_end_plane,nt,np,x_start_plane,delta,deltad,&
!$OMP& p_tmp_plane,tmp_vec,r_vec_plane,y_tmp_plane,z_tmp_plane,xhatBar_plane_norm, Vx_wake_tmp, Vr_wake_tmp,&
!$OMP& nw,Vr_term,Vx_term,tmp_x,tmp_y,tmp_z,&
!$OMP& xHat_plane, yHat_plane, zHat_plane,&
Expand Down Expand Up @@ -277,8 +275,6 @@ subroutine LowResGridCalcOutput(n, u, p, y, m, errStat, errMsg)
x_end_plane = tmp_x + tmp_y + tmp_z

do np = 0, maxPln
! Reset interpolation counter
ILo = 0 ! TODO, curently not used
np1 = np + 1
! Construct the endcaps of the current wake plane volume
x_start_plane = x_end_plane
Expand Down Expand Up @@ -333,7 +329,7 @@ subroutine LowResGridCalcOutput(n, u, p, y, m, errStat, errMsg)
! Increment number of wakes contributing to current grid point
n_wake = n_wake + 1

! Store unit vectors for projection, TODO Compute Vr_term directly
! Store unit vectors for projection
tmp_xhat_plane(:,n_wake) = xHat_plane
tmp_yhat_plane(:,n_wake) = yHat_plane
tmp_zhat_plane(:,n_wake) = zHat_plane
Expand Down Expand Up @@ -427,7 +423,7 @@ subroutine LowResGridCalcOutput(n, u, p, y, m, errStat, errMsg)

Vsum_low = 0.0_ReKi
iwsum = 0
n_r_polar = FLOOR((p%C_Meander*u%D_wake(np,nt))/(2.0_ReKi*p%dpol)) ! TODO change me Dwake*sqrt(2)?
n_r_polar = FLOOR((p%C_Meander*u%D_wake(np,nt))/(2.0_ReKi*p%dpol))

do nr = 0,n_r_polar

Expand Down Expand Up @@ -567,7 +563,7 @@ subroutine HighResGridCalcOutput(n, u, p, y, m, errStat, errMsg)
real(ReKi) :: p_tmp_plane(3)
real(ReKi) :: tmp_vec(3)
real(ReKi) :: delta, deltad
integer(IntKi) :: np1, ILo
integer(IntKi) :: np1
integer(IntKi) :: maxPln
integer(IntKi) :: n_high_low
character(*), parameter :: RoutineName = 'HighResGridCalcOutput'
Expand Down Expand Up @@ -612,8 +608,6 @@ subroutine HighResGridCalcOutput(n, u, p, y, m, errStat, errMsg)
x_end_plane = dot_product(u%xhat_plane(:,0,nt2), (p%Grid_high(:,nXYZ_high,nt) - u%p_plane(:,0,nt2)) )

do np = 0, maxPln !p%NumPlanes-2
! Reset interpolation counter
ILo = 0 ! TODO, curently not used
np1 = np + 1
! Construct the endcaps of the current wake plane volume
x_start_plane = x_end_plane
Expand Down Expand Up @@ -663,7 +657,7 @@ subroutine HighResGridCalcOutput(n, u, p, y, m, errStat, errMsg)
! Increment number of wakes contributing to current grid point
n_wake = n_wake + 1

! Store unit vectors for projection, TODO Compute Vr_term directly
! Store unit vectors for projection
m%xhat_plane(:,n_wake) = xHat_plane
m%yhat_plane(:,n_wake) = yHat_plane
m%zhat_plane(:,n_wake) = zHat_plane
Expand Down
6 changes: 3 additions & 3 deletions modules/awae/src/AWAE_Registry.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -212,8 +212,8 @@ typedef ^ OutputType ReKi Vx_wind_disk {:} -
# Define inputs that are contained on the mesh here:
typedef ^ InputType ReKi xhat_plane {:}{:}{:} - - "Orientations of wake planes, normal to wake planes, for each turbine" -
typedef ^ InputType ReKi p_plane {:}{:}{:} - - "Center positions of wake planes for each turbine" m
typedef ^ InputType ReKi Vx_wake {:}{:}{:}{:} - - "Axial wake velocity deficit at wake planes, distributed acrros the plane, for each turbine (ny,nz,np,nWT)" m/s
typedef ^ InputType ReKi Vy_wake {:}{:}{:}{:} - - "Transverse horizonal wake velocity deficit at wake planes, distributed acrros the plane, for each turbine (ny,nz,np,nWT)" m/s
typedef ^ InputType ReKi Vz_wake {:}{:}{:}{:} - - "Transverse nominally vertical wake velocity deficit at wake planes, distributed acrros the plane, for each turbine (ny,nz,np,nWT)" m/s
typedef ^ InputType ReKi Vx_wake {:}{:}{:}{:} - - "Axial wake velocity deficit at wake planes, distributed across the plane, for each turbine (ny,nz,np,nWT)" m/s
typedef ^ InputType ReKi Vy_wake {:}{:}{:}{:} - - "Transverse horizonal wake velocity deficit at wake planes, distributed across the plane, for each turbine (ny,nz,np,nWT)" m/s
typedef ^ InputType ReKi Vz_wake {:}{:}{:}{:} - - "Transverse nominally vertical wake velocity deficit at wake planes, distributed across the plane, for each turbine (ny,nz,np,nWT)" m/s
typedef ^ InputType ReKi D_wake {:}{:} - - "Wake diameters at wake planes for each turbine" m

24 changes: 3 additions & 21 deletions modules/wakedynamics/src/WakeDynamics.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1103,25 +1103,6 @@ subroutine WD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg )
call Axisymmetric2Cartesian(y%Vx_wake(:,i), y%Vr_wake(:,i), p%r, p%y, p%z, y%Vx_wake2(:,:,i), y%Vy_wake2(:,:,i), y%Vz_wake2(:,:,i))
enddo

! DEBUG check
do i = 0, min(n+1,p%NumPlanes-1)
if( any(abs(y%Vx_wake2(0:,0,i)-y%Vx_wake(:,i))>1e-20)) then
print*,'Problem Vx',i
print*,'y%Vx_wake2(0:,0,i)',y%Vx_wake2(0:,0,i)
print*,'y%Vx_wake(:,i)',y%Vx_wake(:,i)
STOP
endif
!y%Vy_wake2(0:,0,:)
if (any(abs(y%Vy_wake2(0:,0,i)-y%Vr_wake(:,i))>1e-9)) then
print*,'Problem Vr',i
print*,'y%Vx_wake2(0:,0,i)',y%Vy_wake2(0:,0,i)
print*,'y%Vx_wake(:,i)',y%Vr_wake(:,i)
STOP
endif
enddo




end subroutine WD_CalcOutput

Expand Down Expand Up @@ -1211,8 +1192,9 @@ subroutine InitStatesWithInputs(numPlanes, numRadii, u, p, xd, m, errStat, errMs

xd%Vx_rel_disk_filt = u%Vx_rel_disk

! Initialze Ct_azavg_filt and Vx_wake; Vr_wake is already initialized to zero, so, we don't need to do that here.
xd%Ct_azavg_filt (:) = u%Ct_azavg(:)
! Initialze Ct_azavg_filt, Cq_azavg_filt, and Vx_wake; Vr_wake is already initialized to zero, so, we don't need to do that here.
xd%Ct_azavg_filt (:) = u%Ct_azavg(:)
xd%Cq_azavg_filt (:) = u%Cq_azavg(:)

call NearWakeCorrection( xd%Ct_azavg_filt, xd%Cq_azavg_filt, xd%Vx_rel_disk_filt, p, m, xd%Vx_wake(:,0), xd%D_rotor_filt(0), errStat, errMsg )
xd%Vx_wake(:,1) = xd%Vx_wake(:,0)
Expand Down

0 comments on commit 492cb4f

Please sign in to comment.