Skip to content

AD: Initial AeroMap changes for AeroDyn and misc UA/DBEMT changes #1715

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ebranlard merged 7 commits into from
Aug 3, 2023
Merged

AD: Initial AeroMap changes for AeroDyn and misc UA/DBEMT changes #1715

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
0777e23
AD: using utimes for DBEMT
ebranlard Aug 2, 2023
69c8b90
AD: preparing UA lin indexing
ebranlard Aug 2, 2023
7256277
AD: preparing for aeromap
ebranlard Aug 2, 2023
9e77cdb
AD: aeromap and using ua lin idx
ebranlard Aug 2, 2023
1e81918
AD: some UA vars in R8
ebranlard Aug 3, 2023
1652796
AD: reverting to size 31 for du (no DiskAvgVel for now)
ebranlard Aug 3, 2023
ee8f54a
AD: setting UA_NumLinStates to 0 when no cont-states
ebranlard Aug 3, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
810 changes: 471 additions & 339 deletions modules/aerodyn/src/AeroDyn.f90
View file
Open in desktop

Large diffs are not rendered by default.

3 changes: 3 additions & 0 deletions modules/aerodyn/src/AeroDyn_Registry.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,13 +84,15 @@ typedef ^ RotInitInputType R8Ki NacellePosition {3} - - "X-Y-Z reference positio
typedef ^ RotInitInputType R8Ki NacelleOrientation {3}{3} - - "DCM reference orientation of nacelle" -
typedef ^ RotInitInputType IntKi AeroProjMod - 1 - "Flag to switch between different projection models" -
typedef ^ RotInitInputType IntKi AeroBEM_Mod - -1 - "Flag to switch between different BEM Model" -
typedef ^ RotInitInputType ReKi RotSpeed - - 0 "Rotor speed used when AeroDyn is computing aero maps" "rad/s"

typedef ^ InitInputType RotInitInputType rotors {:} - - "Init Input Types for rotors" -
typedef ^ InitInputType CHARACTER(1024) InputFile - - - "Name of the input file" -
typedef ^ InitInputType CHARACTER(1024) RootName - - - "RootName for writing output files" -
typedef ^ InitInputType LOGICAL UsePrimaryInputFile - .TRUE. - "Read input file instead of passed data" -
typedef ^ InitInputType FileInfoType PassedPrimaryInputData - - - "Primary input file as FileInfoType (set by driver/glue code)" -
typedef ^ InitInputType Logical Linearize - .FALSE. - "Flag that tells this module if the glue code wants to linearize." -
typedef ^ InitInputType LOGICAL CompAeroMaps - .FALSE. - "flag to determine if AeroDyn is computing aero maps (true) or running a normal simulation (false)" -
typedef ^ InitInputType ReKi Gravity - - - "Gravity force" Nm/s^2
typedef ^ InitInputType IntKi MHK - - - "MHK turbine type switch" -
typedef ^ InitInputType ReKi defFldDens - - - "Default fluid density from the driver; may be overwritten" kg/m^3
Expand Down Expand Up @@ -421,6 +423,7 @@ typedef ^ RotInputType ReKi InflowOnTower {:}{:} - - "U,V,W at nodes on the towe
typedef ^ RotInputType ReKi InflowOnHub {3} - - "U,V,W at hub" m/s
typedef ^ RotInputType ReKi InflowOnNacelle {3} - - "U,V,W at nacelle" m/s
typedef ^ RotInputType ReKi InflowOnTailFin {3} - - "U,V,W at tailfin" m/s
typedef ^ RotInputType ReKi AvgDiskVel {3} - 0.0 "disk-averaged U,V,W" m/s
typedef ^ RotInputType ReKi UserProp {:}{:} - - "Optional user property for interpolating airfoils (per element per blade)" -


Expand Down
40 changes: 40 additions & 0 deletions modules/aerodyn/src/AeroDyn_Types.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,7 @@ MODULE AeroDyn_Types
REAL(R8Ki) , DIMENSION(1:3,1:3) :: NacelleOrientation !< DCM reference orientation of nacelle [-]
INTEGER(IntKi) :: AeroProjMod = 1 !< Flag to switch between different projection models [-]
INTEGER(IntKi) :: AeroBEM_Mod = -1 !< Flag to switch between different BEM Model [-]
REAL(ReKi) :: RotSpeed !< Rotor speed used when AeroDyn is computing aero maps [rad/s]
END TYPE RotInitInputType
! =======================
! ========= AD_InitInputType =======
Expand All @@ -111,6 +112,7 @@ MODULE AeroDyn_Types
LOGICAL :: UsePrimaryInputFile = .TRUE. !< Read input file instead of passed data [-]
TYPE(FileInfoType) :: PassedPrimaryInputData !< Primary input file as FileInfoType (set by driver/glue code) [-]
LOGICAL :: Linearize = .FALSE. !< Flag that tells this module if the glue code wants to linearize. [-]
LOGICAL :: CompAeroMaps = .FALSE. !< flag to determine if AeroDyn is computing aero maps (true) or running a normal simulation (false) [-]
REAL(ReKi) :: Gravity !< Gravity force [Nm/s^2]
INTEGER(IntKi) :: MHK !< MHK turbine type switch [-]
REAL(ReKi) :: defFldDens !< Default fluid density from the driver; may be overwritten [kg/m^3]
Expand Down Expand Up @@ -457,6 +459,7 @@ MODULE AeroDyn_Types
REAL(ReKi) , DIMENSION(1:3) :: InflowOnHub !< U,V,W at hub [m/s]
REAL(ReKi) , DIMENSION(1:3) :: InflowOnNacelle !< U,V,W at nacelle [m/s]
REAL(ReKi) , DIMENSION(1:3) :: InflowOnTailFin !< U,V,W at tailfin [m/s]
REAL(ReKi) , DIMENSION(1:3) :: AvgDiskVel !< disk-averaged U,V,W [m/s]
REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: UserProp !< Optional user property for interpolating airfoils (per element per blade) [-]
END TYPE RotInputType
! =======================
Expand Down Expand Up @@ -1436,6 +1439,7 @@ SUBROUTINE AD_CopyRotInitInputType( SrcRotInitInputTypeData, DstRotInitInputType
DstRotInitInputTypeData%NacelleOrientation = SrcRotInitInputTypeData%NacelleOrientation
DstRotInitInputTypeData%AeroProjMod = SrcRotInitInputTypeData%AeroProjMod
DstRotInitInputTypeData%AeroBEM_Mod = SrcRotInitInputTypeData%AeroBEM_Mod
DstRotInitInputTypeData%RotSpeed = SrcRotInitInputTypeData%RotSpeed
END SUBROUTINE AD_CopyRotInitInputType

SUBROUTINE AD_DestroyRotInitInputType( RotInitInputTypeData, ErrStat, ErrMsg, DEALLOCATEpointers )
Expand Down Expand Up @@ -1519,6 +1523,7 @@ SUBROUTINE AD_PackRotInitInputType( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat,
Db_BufSz = Db_BufSz + SIZE(InData%NacelleOrientation) ! NacelleOrientation
Int_BufSz = Int_BufSz + 1 ! AeroProjMod
Int_BufSz = Int_BufSz + 1 ! AeroBEM_Mod
Re_BufSz = Re_BufSz + 1 ! RotSpeed
IF ( Re_BufSz .GT. 0 ) THEN
ALLOCATE( ReKiBuf( Re_BufSz ), STAT=ErrStat2 )
IF (ErrStat2 /= 0) THEN
Expand Down Expand Up @@ -1617,6 +1622,8 @@ SUBROUTINE AD_PackRotInitInputType( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat,
Int_Xferred = Int_Xferred + 1
IntKiBuf(Int_Xferred) = InData%AeroBEM_Mod
Int_Xferred = Int_Xferred + 1
ReKiBuf(Re_Xferred) = InData%RotSpeed
Re_Xferred = Re_Xferred + 1
END SUBROUTINE AD_PackRotInitInputType

SUBROUTINE AD_UnPackRotInitInputType( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMsg )
Expand Down Expand Up @@ -1737,6 +1744,8 @@ SUBROUTINE AD_UnPackRotInitInputType( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrSt
Int_Xferred = Int_Xferred + 1
OutData%AeroBEM_Mod = IntKiBuf(Int_Xferred)
Int_Xferred = Int_Xferred + 1
OutData%RotSpeed = ReKiBuf(Re_Xferred)
Re_Xferred = Re_Xferred + 1
END SUBROUTINE AD_UnPackRotInitInputType

SUBROUTINE AD_CopyInitInput( SrcInitInputData, DstInitInputData, CtrlCode, ErrStat, ErrMsg )
Expand Down Expand Up @@ -1777,6 +1786,7 @@ SUBROUTINE AD_CopyInitInput( SrcInitInputData, DstInitInputData, CtrlCode, ErrSt
CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg,RoutineName)
IF (ErrStat>=AbortErrLev) RETURN
DstInitInputData%Linearize = SrcInitInputData%Linearize
DstInitInputData%CompAeroMaps = SrcInitInputData%CompAeroMaps
DstInitInputData%Gravity = SrcInitInputData%Gravity
DstInitInputData%MHK = SrcInitInputData%MHK
DstInitInputData%defFldDens = SrcInitInputData%defFldDens
Expand Down Expand Up @@ -1900,6 +1910,7 @@ SUBROUTINE AD_PackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg
DEALLOCATE(Int_Buf)
END IF
Int_BufSz = Int_BufSz + 1 ! Linearize
Int_BufSz = Int_BufSz + 1 ! CompAeroMaps
Re_BufSz = Re_BufSz + 1 ! Gravity
Int_BufSz = Int_BufSz + 1 ! MHK
Re_BufSz = Re_BufSz + 1 ! defFldDens
Expand Down Expand Up @@ -2017,6 +2028,8 @@ SUBROUTINE AD_PackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg
ENDIF
IntKiBuf(Int_Xferred) = TRANSFER(InData%Linearize, IntKiBuf(1))
Int_Xferred = Int_Xferred + 1
IntKiBuf(Int_Xferred) = TRANSFER(InData%CompAeroMaps, IntKiBuf(1))
Int_Xferred = Int_Xferred + 1
ReKiBuf(Re_Xferred) = InData%Gravity
Re_Xferred = Re_Xferred + 1
IntKiBuf(Int_Xferred) = InData%MHK
Expand Down Expand Up @@ -2172,6 +2185,8 @@ SUBROUTINE AD_UnPackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, Err
IF(ALLOCATED(Int_Buf)) DEALLOCATE(Int_Buf)
OutData%Linearize = TRANSFER(IntKiBuf(Int_Xferred), OutData%Linearize)
Int_Xferred = Int_Xferred + 1
OutData%CompAeroMaps = TRANSFER(IntKiBuf(Int_Xferred), OutData%CompAeroMaps)
Int_Xferred = Int_Xferred + 1
OutData%Gravity = ReKiBuf(Re_Xferred)
Re_Xferred = Re_Xferred + 1
OutData%MHK = IntKiBuf(Int_Xferred)
Expand Down Expand Up @@ -15894,6 +15909,7 @@ SUBROUTINE AD_CopyRotInputType( SrcRotInputTypeData, DstRotInputTypeData, CtrlCo
DstRotInputTypeData%InflowOnHub = SrcRotInputTypeData%InflowOnHub
DstRotInputTypeData%InflowOnNacelle = SrcRotInputTypeData%InflowOnNacelle
DstRotInputTypeData%InflowOnTailFin = SrcRotInputTypeData%InflowOnTailFin
DstRotInputTypeData%AvgDiskVel = SrcRotInputTypeData%AvgDiskVel
IF (ALLOCATED(SrcRotInputTypeData%UserProp)) THEN
i1_l = LBOUND(SrcRotInputTypeData%UserProp,1)
i1_u = UBOUND(SrcRotInputTypeData%UserProp,1)
Expand Down Expand Up @@ -16127,6 +16143,7 @@ SUBROUTINE AD_PackRotInputType( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, Err
Re_BufSz = Re_BufSz + SIZE(InData%InflowOnHub) ! InflowOnHub
Re_BufSz = Re_BufSz + SIZE(InData%InflowOnNacelle) ! InflowOnNacelle
Re_BufSz = Re_BufSz + SIZE(InData%InflowOnTailFin) ! InflowOnTailFin
Re_BufSz = Re_BufSz + SIZE(InData%AvgDiskVel) ! AvgDiskVel
Int_BufSz = Int_BufSz + 1 ! UserProp allocated yes/no
IF ( ALLOCATED(InData%UserProp) ) THEN
Int_BufSz = Int_BufSz + 2*2 ! UserProp upper/lower bounds for each dimension
Expand Down Expand Up @@ -16410,6 +16427,10 @@ SUBROUTINE AD_PackRotInputType( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, Err
ReKiBuf(Re_Xferred) = InData%InflowOnTailFin(i1)
Re_Xferred = Re_Xferred + 1
END DO
DO i1 = LBOUND(InData%AvgDiskVel,1), UBOUND(InData%AvgDiskVel,1)
ReKiBuf(Re_Xferred) = InData%AvgDiskVel(i1)
Re_Xferred = Re_Xferred + 1
END DO
IF ( .NOT. ALLOCATED(InData%UserProp) ) THEN
IntKiBuf( Int_Xferred ) = 0
Int_Xferred = Int_Xferred + 1
Expand Down Expand Up @@ -16802,6 +16823,12 @@ SUBROUTINE AD_UnPackRotInputType( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat,
OutData%InflowOnTailFin(i1) = ReKiBuf(Re_Xferred)
Re_Xferred = Re_Xferred + 1
END DO
i1_l = LBOUND(OutData%AvgDiskVel,1)
i1_u = UBOUND(OutData%AvgDiskVel,1)
DO i1 = LBOUND(OutData%AvgDiskVel,1), UBOUND(OutData%AvgDiskVel,1)
OutData%AvgDiskVel(i1) = ReKiBuf(Re_Xferred)
Re_Xferred = Re_Xferred + 1
END DO
IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN ! UserProp not allocated
Int_Xferred = Int_Xferred + 1
ELSE
Expand Down Expand Up @@ -18311,6 +18338,12 @@ SUBROUTINE AD_Input_ExtrapInterp1(u1, u2, tin, u_out, tin_out, ErrStat, ErrMsg )
END DO
ENDDO
DO i01 = LBOUND(u_out%rotors,1),UBOUND(u_out%rotors,1)
DO i1 = LBOUND(u_out%rotors(i01)%AvgDiskVel,1),UBOUND(u_out%rotors(i01)%AvgDiskVel,1)
b = -(u1%rotors(i01)%AvgDiskVel(i1) - u2%rotors(i01)%AvgDiskVel(i1))
u_out%rotors(i01)%AvgDiskVel(i1) = u1%rotors(i01)%AvgDiskVel(i1) + b * ScaleFactor
END DO
ENDDO
DO i01 = LBOUND(u_out%rotors,1),UBOUND(u_out%rotors,1)
IF (ALLOCATED(u_out%rotors(i01)%UserProp) .AND. ALLOCATED(u1%rotors(i01)%UserProp)) THEN
DO i2 = LBOUND(u_out%rotors(i01)%UserProp,2),UBOUND(u_out%rotors(i01)%UserProp,2)
DO i1 = LBOUND(u_out%rotors(i01)%UserProp,1),UBOUND(u_out%rotors(i01)%UserProp,1)
Expand Down Expand Up @@ -18469,6 +18502,13 @@ SUBROUTINE AD_Input_ExtrapInterp2(u1, u2, u3, tin, u_out, tin_out, ErrStat, ErrM
END DO
ENDDO
DO i01 = LBOUND(u_out%rotors,1),UBOUND(u_out%rotors,1)
DO i1 = LBOUND(u_out%rotors(i01)%AvgDiskVel,1),UBOUND(u_out%rotors(i01)%AvgDiskVel,1)
b = (t(3)**2*(u1%rotors(i01)%AvgDiskVel(i1) - u2%rotors(i01)%AvgDiskVel(i1)) + t(2)**2*(-u1%rotors(i01)%AvgDiskVel(i1) + u3%rotors(i01)%AvgDiskVel(i1)))* scaleFactor
c = ( (t(2)-t(3))*u1%rotors(i01)%AvgDiskVel(i1) + t(3)*u2%rotors(i01)%AvgDiskVel(i1) - t(2)*u3%rotors(i01)%AvgDiskVel(i1) ) * scaleFactor
u_out%rotors(i01)%AvgDiskVel(i1) = u1%rotors(i01)%AvgDiskVel(i1) + b + c * t_out
END DO
ENDDO
DO i01 = LBOUND(u_out%rotors,1),UBOUND(u_out%rotors,1)
IF (ALLOCATED(u_out%rotors(i01)%UserProp) .AND. ALLOCATED(u1%rotors(i01)%UserProp)) THEN
DO i2 = LBOUND(u_out%rotors(i01)%UserProp,2),UBOUND(u_out%rotors(i01)%UserProp,2)
DO i1 = LBOUND(u_out%rotors(i01)%UserProp,1),UBOUND(u_out%rotors(i01)%UserProp,1)
Expand Down
2 changes: 1 addition & 1 deletion modules/aerodyn/src/BEMT.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -882,7 +882,7 @@ subroutine BEMT_UpdateStates( t, n, u, utimes, p, x, xd, z, OtherState, AFInfo,
!........................
do j = 1,p%numBlades
do i = 1,p%numBladeNodes
call DBEMT_UpdateStates(i, j, t, n, m%u_DBEMT, p%DBEMT, x%DBEMT, OtherState%DBEMT, m%DBEMT, errStat2, errMsg2)
call DBEMT_UpdateStates(i, j, t, n, m%u_DBEMT, uTimes, p%DBEMT, x%DBEMT, OtherState%DBEMT, m%DBEMT, errStat2, errMsg2)
if (ErrStat2 /= ErrID_None) then
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName//trim(NodeText(i,j)))
if (errStat >= AbortErrLev) return
Expand Down
10 changes: 4 additions & 6 deletions modules/aerodyn/src/DBEMT.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
! [2] R. Damiani, J.Jonkman
! DBEMT Theory Rev. 3
! Unpublished
!
module DBEMT

use NWTC_Library
Expand Down Expand Up @@ -328,13 +329,14 @@ end subroutine DBEMT_InitStates
!!----------------------------------------------------------------------------------------------------------------------------------
!> Loose coupling routine for solving for constraint states, integrating continuous states, and updating discrete and other states.
!! Continuous, constraint, discrete, and other states are updated for t + Interval
subroutine DBEMT_UpdateStates( i, j, t, n, u, p, x, OtherState, m, errStat, errMsg )
subroutine DBEMT_UpdateStates( i, j, t, n, u, uTimes, p, x, OtherState, m, errStat, errMsg )
!..................................................................................................................................
integer(IntKi), intent(in ) :: i !< blade node counter
integer(IntKi), intent(in ) :: j !< blade counter
real(DbKi), intent(in ) :: t !< Current simulation time in seconds
integer(IntKi), intent(in ) :: n !< Current simulation time step n = 0,1,...
type(DBEMT_InputType), intent(in ) :: u(2) !< Inputs at t and t+dt
real(DbKi), intent(in ) :: uTimes(2) ! Times associated with u(:), in seconds
type(DBEMT_ParameterType), intent(in ) :: p !< Parameters
type(DBEMT_ContinuousStateType), intent(inout) :: x !< Input: Continuous states at t;
!! Output: Continuous states at t + Interval
Expand All @@ -346,7 +348,6 @@ subroutine DBEMT_UpdateStates( i, j, t, n, u, p, x, OtherState, m, errStat, errM
! local variables
real(ReKi) :: A, B, C0, k_tau, C0_2 ! tau1_plus1, C_tau1, C, K1
integer(IntKi) :: indx
real(DbKi) :: utimes(2)

TYPE(DBEMT_ElementInputType) :: u_elem(2) !< Inputs at utimes

Expand All @@ -364,9 +365,6 @@ subroutine DBEMT_UpdateStates( i, j, t, n, u, p, x, OtherState, m, errStat, errM
call DBEMT_InitStates( i, j, u(1), p, x, OtherState )

if (p%DBEMT_Mod == DBEMT_cont_tauConst) then ! continuous formulation:
utimes(1) = t
utimes(2) = t + p%dt

u_elem(1) = u(1)%element(i,j)
u_elem(2) = u(2)%element(i,j)
call DBEMT_ABM4( i, j, t, n, u_elem, utimes, p, x, OtherState, m, ErrStat, ErrMsg )
Expand Down Expand Up @@ -924,4 +922,4 @@ subroutine DBEMT_End( u, p, x, OtherState, m, ErrStat, ErrMsg )

END SUBROUTINE DBEMT_End

end module DBEMT
end module DBEMT
5 changes: 2 additions & 3 deletions modules/aerodyn/src/FVW_Subs.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -462,7 +462,6 @@ logical function have_nan(p, m, x, z, u, label)
character(len=*), intent(in) :: label !< label for print
integer :: iW
have_nan=.False.
!bjj: If we used Is_NaN (or a version of it for this data type) instead of isnan, I'd get fewer compiler warnings that "Fortran 2003 does not allow this intrinsic procedure."
do iW = 1,size(p%W)
if (any(isnan(x%W(iW)%r_NW))) then
print*,trim(label),'NaN in W(iW)%r_NW'//trim(num2lstr(iW))
Expand All @@ -481,11 +480,11 @@ logical function have_nan(p, m, x, z, u, label)
have_nan=.True.
endif
if (any(isnan(x%W(iW)%Eps_NW))) then
print*,trim(label),'NaN in E_NW'
print*,trim(label),'NaN in E_NW'//trim(num2lstr(iW))
have_nan=.True.
endif
if (any(isnan(x%W(iW)%Eps_FW))) then
print*,trim(label),'NaN in E_FW'
print*,trim(label),'NaN in E_FW'//trim(num2lstr(iW))
have_nan=.True.
endif
if (any(isnan(z%W(iW)%Gamma_LL))) then
Expand Down
Loading