Skip to content

Commit

Permalink
Merge pull request #1715 AD: Initial AeroMap changes for AeroDyn and …
Browse files Browse the repository at this point in the history 
…misc UA/DBEMT changes

* If statements in AeroDyn to perform the linearization differently when the AeroMap flag is set
* Uses UA%lin_xIndx in UA and AeroDyn to store the indices of the states in the operating point vector x_op
* Uses utimes in DBEMT in cases the inputs are at multiple times (and not just dt and t+dt).
* Uses double precision variables for some of the critical variables of UA
* Prepares the ground for a followup pull request that will affect the test results


Co-authored-by: bjonkman <[email protected]>
  • Loading branch information
@ebranlard @bjonkman
ebranlard and bjonkman authored Aug 3, 2023
2 parents 10c9e33 + ee8f54a commit 2199ad3
Show file tree
Hide file tree
Showing 9 changed files with 985 additions and 384 deletions.
810 changes: 471 additions & 339 deletions modules/aerodyn/src/AeroDyn.f90
View file

Large diffs are not rendered by default.

3 changes: 3 additions & 0 deletions modules/aerodyn/src/AeroDyn_Registry.txt
View file
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
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
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
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
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

0 comments on commit 2199ad3

Please sign in to comment.