Skip to content

Commit

Permalink
Merge branch 'f/AD-shear' of https://github.com/ebranlard/openfast in…
Browse files Browse the repository at this point in the history 
…to f/AD-shear
  • Loading branch information
Emmanuel Branlard committed May 6, 2024
2 parents be84c59 + c228bd3 commit 791956d
Show file tree
Hide file tree
Showing 9 changed files with 579 additions and 953 deletions.
4 changes: 2 additions & 2 deletions docs/source/user/aerodyn/input.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ When ``Wake_Mod`` is set to 3, the free vortex wake model is used, also referred
:numref:`OLAF`). ``Wake_Mod`` cannot be set to 3 during linearization analyses.

.. note::
Link to old inputs: The previous input `WakeMod` is removed, `WakeMod=2` used to mean DBEMT, but this now controled using `DBEMT_Mod`.
Link to old inputs: The previous input `WakeMod` is removed, `WakeMod=2` used to mean DBEMT, but this now controlled using `DBEMT_Mod`.
`WakeMod=2` is a placeholder for future induction calculation method.


Expand Down Expand Up @@ -328,7 +328,7 @@ Most ``UA_Mod`` will require `AoA34` to be set to true. But when using quasi-ste
This feature is currently not implemented due to a lag between the `dev` and `dev-unstable` branch.

.. note::
Link to previous inputs: `AFAeroMod=1` implies `AoA34=False`. But to have a fair comparison between steady and unsteady aerodynamics model, it would be best to set `AoA34=True` when when doing quasi-steady aero.
Link to previous inputs: `AFAeroMod=1` implies `AoA34=False`. But to have a fair comparison between steady and unsteady aerodynamics model, it would be best to set `AoA34=True` when doing quasi-steady aero.



Expand Down
36 changes: 20 additions & 16 deletions modules/aerodyn/python-lib/aerodyn_inflow_library.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@
c_byte,
c_int,
c_double,
c_float,
c_float,
c_char,
c_char_p,
c_wchar,
Expand Down Expand Up @@ -98,7 +98,7 @@ def __init__(self, library_path):
self.WtrDpth = 0.0 # Water depth (m)
self.MSL2SWL = 0.0 # Offset between still-water level and mean sea level (m) [positive upward]

# flags
# flags
self.storeHHVel = 1 # 0=false, 1=true
self.transposeDCM= 1 # 0=false, 1=true
self.debuglevel = 0 # 0-4 levels
Expand Down Expand Up @@ -149,6 +149,7 @@ def __init__(self, library_path):
self.numMeshPts = 1
self.initMeshPos = np.zeros(shape=(self.numMeshPts,3),dtype=c_float ) # Nx3 array [x,y,z]
self.initMeshOrient = np.zeros(shape=(self.numMeshPts,9),dtype=c_double) # Nx9 array [r11,r12,r13,r21,r22,r23,r31,r32,r33]
self.meshPtToBladeNum = np.zeros(shape=(self.numMeshPts),dtype=c_int) # Nx1 array [blade number]

# OutRootName
# If HD writes a file (echo, summary, or other), use this for the
Expand All @@ -161,11 +162,11 @@ def _initialize_routines(self):
self.ADI_C_PreInit.argtypes = [
POINTER(c_int), # numTurbines
POINTER(c_int), # transposeDCM
POINTER(c_int), # debuglevel
POINTER(c_int), # debuglevel
POINTER(c_int), # ErrStat_C
POINTER(c_char) # ErrMsg_C
]
self.ADI_C_PreInit.restype = c_int
self.ADI_C_PreInit.restype = c_int

# setup one rotor
self.ADI_C_SetupRotor.argtypes = [
Expand All @@ -182,6 +183,7 @@ def _initialize_routines(self):
POINTER(c_int), # numMeshPts
POINTER(c_float), # initMeshPos_flat
POINTER(c_double), # initMeshOrient_flat
POINTER(c_int), # meshPtToBladeNum
POINTER(c_int), # ErrStat_C
POINTER(c_char) # ErrMsg_C
]
Expand All @@ -194,7 +196,7 @@ def _initialize_routines(self):
POINTER(c_int), # IfW input file passed as string
POINTER(c_char_p), # IfW input file as string
POINTER(c_int), # IfW input file string length
POINTER(c_char), # OutRootName
POINTER(c_char), # OutRootName
POINTER(c_float), # gravity
POINTER(c_float), # defFldDens
POINTER(c_float), # defKinVisc
Expand All @@ -203,9 +205,9 @@ def _initialize_routines(self):
POINTER(c_float), # defPvap
POINTER(c_float), # WtrDpth
POINTER(c_float), # MSL2SWL
POINTER(c_int), # InterpOrder
POINTER(c_int), # InterpOrder
POINTER(c_double), # dt
POINTER(c_double), # tmax
POINTER(c_double), # tmax
POINTER(c_int), # storeHHVel
POINTER(c_int), # WrVTK
POINTER(c_int), # WrVTK_Type
Expand All @@ -219,16 +221,16 @@ def _initialize_routines(self):
POINTER(c_int), # ErrStat_C
POINTER(c_char) # ErrMsg_C
]
self.ADI_C_Init.restype = c_int
self.ADI_C_Init.restype = c_int

#self.ADI_C_ReInit.argtypes = [
# POINTER(c_double), # t_initial
# POINTER(c_double), # t_initial
# POINTER(c_double), # dt
# POINTER(c_double), # tmax
# POINTER(c_double), # tmax
# POINTER(c_int), # ErrStat_C
# POINTER(c_char) # ErrMsg_C
#]
#self.ADI_C_ReInit.restype = c_int
#self.ADI_C_ReInit.restype = c_int

self.ADI_C_SetRotorMotion.argtypes = [
POINTER(c_int), # iturb
Expand Down Expand Up @@ -322,6 +324,7 @@ def adi_setuprotor(self,iturb,isHAWT,turbRefPos):
initRootOrient_flat_c = self.flatOrientArr(self._initNumBlades, self.numBlades,self.initRootOrient, 'RootOrient')
initMeshPos_flat_c = self.flatPosArr( self._initNumMeshPts,self.numMeshPts,self.initMeshPos, 'MeshPos')
initMeshOrient_flat_c = self.flatOrientArr(self._initNumMeshPts,self.numMeshPts,self.initMeshOrient,'MeshOrient')
initMeshPtToBladeNum_flat_c = (c_int * len(self.meshPtToBladeNum))(*self.meshPtToBladeNum)

self.ADI_C_SetupRotor(
c_int(iturb), # IN: iturb -- current turbine number
Expand All @@ -337,6 +340,7 @@ def adi_setuprotor(self,iturb,isHAWT,turbRefPos):
byref(c_int(self.numMeshPts)), # IN: number of mesh points expected
initMeshPos_flat_c, # IN: initMeshPos -- initial node positions in flat array of 3*numMeshPts
initMeshOrient_flat_c, # IN: initMeshOrient -- initial node orientation DCMs in flat array of 9*numMeshPts
initMeshPtToBladeNum_flat_c, # IN: initMeshPtToBladeNum -- initial mesh point to blade number mapping
byref(self.error_status_c), # OUT: ErrStat_C
self.error_message_c # OUT: ErrMsg_C
)
Expand Down Expand Up @@ -402,7 +406,7 @@ def adi_init(self, AD_input_string_array, IfW_input_string_array):
)

self.check_error()

# Initialize output channels
self.numChannels = self._numChannels_c.value

Expand Down Expand Up @@ -497,7 +501,7 @@ def adi_getrotorloads(self, iturb, meshFrcMom):
)

self.check_error()

## Reshape Force/Moment into [N,6]
count = 0
for j in range(0,self.numMeshPts):
Expand All @@ -518,7 +522,7 @@ def adi_calcOutput(self, time, outputChannelValues):

# Run ADI_C_CalcOutput
self.ADI_C_CalcOutput(
byref(c_double(time)), # IN: time at which to calculate output forces
byref(c_double(time)), # IN: time at which to calculate output forces
outputChannelValues_c, # OUT: output channel values as described in input file
byref(self.error_status_c), # OUT: ErrStat_C
self.error_message_c # OUT: ErrMsg_C
Expand All @@ -535,7 +539,7 @@ def adi_updateStates(self, time, timeNext):

# Run AeroDyn_Inflow_UpdateStates_c
self.ADI_C_UpdateStates(
byref(c_double(time)), # IN: time at which to calculate output forces
byref(c_double(time)), # IN: time at which to calculate output forces
byref(c_double(timeNext)), # IN: time T+dt we are stepping to
byref(self.error_status_c), # OUT: ErrStat_C
self.error_message_c # OUT: ErrMsg_C
Expand Down Expand Up @@ -662,7 +666,7 @@ def check_init_hubroot(self):
print("Expecting a 9 element array for initNacelleOrient DCM [r11,r12,r13,r21,r22,r23,r31,r32,r33]")
self.adi_end()
raise Exception("\nAeroDyn terminated prematurely.")


def check_init_mesh(self):
#print("shape of initMeshPos ", self.initMeshPos.shape)
Expand Down
6 changes: 4 additions & 2 deletions modules/aerodyn/src/AeroDyn.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -360,7 +360,6 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut
endif
enddo

print*,'AD_Init: AeroProjMod B:',AeroProjMod

! -----------------------------------------------------------------
! Read the AeroDyn blade files, or copy from passed input
Expand Down Expand Up @@ -4495,6 +4494,9 @@ SUBROUTINE Init_BEMTmodule( InputFileData, RotInputFileData, u_AD, u, p, p_AD, x
elseif (p%AeroProjMod==APM_LiftingLine) then
Label='Projection: Lifting Line'
else
! Normally we wouldn't want to do a print or STOP, but we
! should never get here unless a programmer made a mistake.
! I'll leave this as is for now. - ADP
print*,'Invalid projection method'
STOP
endif
Expand Down Expand Up @@ -4745,7 +4747,7 @@ SUBROUTINE TFin_CalcOutput(p, p_AD, u, m, y, ErrStat, ErrMsg )
V_ind = 0.0_ReKi
elseif(p%TFin%TFinIndMod==TFinIndMod_rotavg) then
! TODO TODO
print*,'TODO TailFin: compute rotor average induced velocity'
call WrScr('TODO TailFin: compute rotor average induced velocity')
V_ind = 0.0_ReKi
else
STOP ! Will never happen
Expand Down
47 changes: 25 additions & 22 deletions modules/aerodyn/src/AeroDyn_IO.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -659,6 +659,7 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade
integer(IntKi) :: CurLine !< current entry in FileInfo_In%Lines array
real(ReKi) :: TmpRe5(5) !< temporary 8 number array for reading values in
character(1024) :: sDummy !< temporary string
character(1024) :: tmpOutStr !< temporary string for writing to screen
logical :: wakeModProvided, frozenWakeProvided, skewModProvided, AFAeroModProvided, UAModProvided, isLegalComment, firstWarn !< Temporary for legacy purposes
logical :: AoA34_Missing
integer :: UAMod_Old
Expand Down Expand Up @@ -1177,19 +1178,20 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade
! pass comment lines
elseif (index(sDummy, 'SECTAVG')>1) then
read(sDummy, '(L1)') InputFileData%SectAvg
print*,' >>> SectAvg ',InputFileData%SectAvg
write(tmpOutStr,*) ' >>> SectAvg ',InputFileData%SectAvg
elseif (index(sDummy, 'SA_PSIBWD')>1) then
read(sDummy, *) InputFileData%SA_PsiBwd
print*,' >>> SA_PsiBwd ',InputFileData%SA_PsiBwd
write(tmpOutStr,*) ' >>> SA_PsiBwd ',InputFileData%SA_PsiBwd
elseif (index(sDummy, 'SA_PSIFWD')>1) then
read(sDummy, *) InputFileData%SA_PsiFwd
print*,' >>> SA_PsiFwd ',InputFileData%SA_PsiFwd
write(tmpOutStr,*) ' >>> SA_PsiFwd ',InputFileData%SA_PsiFwd
elseif (index(sDummy, 'SA_NPERSEC')>1) then
read(sDummy, *) InputFileData%SA_nPerSec
print*,' >>> SA_nPerSec ',InputFileData%SA_nPerSec
write(tmpOutStr,*) ' >>> SA_nPerSec ',InputFileData%SA_nPerSec
else
print*,'[WARN] AeroDyn Line ignored: '//trim(sDummy)
write(tmpOutStr,*) '[WARN] AeroDyn Line ignored: '//trim(sDummy)
endif
call WrScr(trim(tmpOutStr))
enddo


Expand Down Expand Up @@ -1271,6 +1273,7 @@ end function newInputMissing

!-------------------------------------------------------------------------------------------------
subroutine printNewOldInputs()
character(1024) :: tmpStr
! Temporary HACK, for WakeMod=10, 11 or 12 use AeroProjMod 2 (will trigger PolarBEM)
if (InputFileData%Wake_Mod==10) then
call WrScr('[WARN] Wake_Mod=10 is a temporary hack. Setting BEM_Mod to 0')
Expand All @@ -1285,23 +1288,23 @@ subroutine printNewOldInputs()
!====== Summary of new AeroDyn options ===============================================================
! NOTE: remove me in future release
call WrScr('-------------- New AeroDyn inputs (with new meaning):')
write (*,'(A20,I0)') 'Wake_Mod: ' , InputFileData%Wake_Mod
write (*,'(A20,I0)') 'BEM_Mod: ' , InputFileData%BEM_Mod
write (*,'(A20,L0)') 'SectAvg: ' , InputFileData%SectAvg
write (*,'(A20,I0)') 'SectAvgWeighting: ', InputFileData%SA_Weighting
write (*,'(A20,I0)') 'SectAvgNPoints: ', InputFileData%SA_nPerSec
write (*,'(A20,I0)') 'DBEMT_Mod:' , InputFileData%DBEMT_Mod
write (*,'(A20,I0)') 'Skew_Mod: ' , InputFileData%Skew_Mod
write (*,'(A20,L0)') 'SkewMomCorr:' , InputFileData%SkewMomCorr
write (*,'(A20,I0)') 'SkewRedistr_Mod:' , InputFileData%SkewRedistr_Mod
write (*,'(A20,L0)') 'AoA34: ' , InputFileData%AoA34
write (*,'(A20,I0)') 'UA_Mod: ' , InputFileData%UAMod
write (tmpStr,'(A20,I0)') 'Wake_Mod: ' , InputFileData%Wake_Mod; call WrScr(trim(tmpStr))
write (tmpStr,'(A20,I0)') 'BEM_Mod: ' , InputFileData%BEM_Mod; call WrScr(trim(tmpStr))
write (tmpStr,'(A20,L1)') 'SectAvg: ' , InputFileData%SectAvg; call WrScr(trim(tmpStr))
write (tmpStr,'(A20,I0)') 'SectAvgWeighting: ', InputFileData%SA_Weighting; call WrScr(trim(tmpStr))
write (tmpStr,'(A20,I0)') 'SectAvgNPoints: ', InputFileData%SA_nPerSec; call WrScr(trim(tmpStr))
write (tmpStr,'(A20,I0)') 'DBEMT_Mod:' , InputFileData%DBEMT_Mod; call WrScr(trim(tmpStr))
write (tmpStr,'(A20,I0)') 'Skew_Mod: ' , InputFileData%Skew_Mod; call WrScr(trim(tmpStr))
write (tmpStr,'(A20,L1)') 'SkewMomCorr:' , InputFileData%SkewMomCorr; call WrScr(trim(tmpStr))
write (tmpStr,'(A20,I0)') 'SkewRedistr_Mod:' , InputFileData%SkewRedistr_Mod; call WrScr(trim(tmpStr))
write (tmpStr,'(A20,L1)') 'AoA34: ' , InputFileData%AoA34; call WrScr(trim(tmpStr))
write (tmpStr,'(A20,I0)') 'UA_Mod: ' , InputFileData%UAMod; call WrScr(trim(tmpStr))
call WrScr('-------------- Old AeroDyn inputs:')
write (*,'(A20,I0)') 'WakeMod: ', InputFileData%WakeMod
write (*,'(A20,I0)') 'SkewMod: ', InputFileData%SkewMod
write (*,'(A20,I0)') 'AFAeroMod:', InputFileData%AFAeroMod
write (*,'(A20,L0)') 'FrozenWake:', InputFileData%FrozenWake
write (*,'(A20,I0)') 'UAMod: ', UAMod_Old
write (tmpStr,'(A20,I0)') 'WakeMod: ', InputFileData%WakeMod; call WrScr(trim(tmpStr))
write (tmpStr,'(A20,I0)') 'SkewMod: ', InputFileData%SkewMod; call WrScr(trim(tmpStr))
write (tmpStr,'(A20,I0)') 'AFAeroMod:', InputFileData%AFAeroMod; call WrScr(trim(tmpStr))
write (tmpStr,'(A20,L1)') 'FrozenWake:', InputFileData%FrozenWake; call WrScr(trim(tmpStr))
write (tmpStr,'(A20,I0)') 'UAMod: ', UAMod_Old; call WrScr(trim(tmpStr))
call WrScr('------------------------------------------------------')
end subroutine printNewOldInputs

Expand Down Expand Up @@ -2408,7 +2411,7 @@ subroutine calcCantAngle(f, xi,stencilSize,n,cantAngle)
call differ_stencil ( xi(i), 1, 2, xiIn, cx, info )
!call differ_stencil ( xi(i), 1, 2, fIn, cf, info )
if (info /= 0) then
print*,'Cant Calc failed at i=',i
call WrScr('Cant Calc failed at i='//trim(Num2LStr(i)))
else
cPrime(i) = 0.0
fPrime(i) = 0.0
Expand Down
Loading

0 comments on commit 791956d

Please sign in to comment.