Skip to content

Commit

Permalink
Lin: fix indexing in AD15 lin, and add access to IfW OP
Browse files Browse the repository at this point in the history
  • Loading branch information
@andrew-platt
andrew-platt committed Jan 29, 2024
1 parent f0e1d89 commit 617c39e
Show file tree
Hide file tree
Showing 5 changed files with 32 additions and 94 deletions.
37 changes: 21 additions & 16 deletions modules/aerodyn/src/AeroDyn.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ module AeroDyn
use UnsteadyAero
use FVW
use FVW_Subs, only: FVW_AeroOuts
use IfW_FlowField, only: IfW_FlowField_GetVelAcc
use IfW_FlowField, only: IfW_FlowField_GetVelAcc, IfW_UniformWind_GetOP

implicit none
private
Expand Down Expand Up @@ -1193,6 +1193,7 @@ subroutine Init_u( u, p, p_AD, InputFileData, MHK, WtrDpth, InitInp, errStat, er
,TranslationVel = .true. &
,RotationVel = .true. &
,TranslationAcc = .true. &
,RotationAcc = .true. &
)
call SetErrStat( errStat2, errMsg2, errStat, errMsg, RoutineName )

Expand Down Expand Up @@ -6071,6 +6072,7 @@ SUBROUTINE RotGetOP( t, u, RotInflow, p, p_AD, x, xd, z, OtherState, y, m, ErrSt
CHARACTER(*), PARAMETER :: RoutineName = 'AD_GetOP'
LOGICAL :: FieldMask(FIELDMASK_SIZE)
TYPE(RotContinuousStateType) :: dxdt
real(ReKi) :: OP_out(3) !< operating point of wind (HWindSpeed, PLexp, and AngleH)


! Initialize ErrStat
Expand Down Expand Up @@ -6196,11 +6198,18 @@ SUBROUTINE RotGetOP( t, u, RotInflow, p, p_AD, x, xd, z, OtherState, y, m, ErrSt
end do

!------------------------------
! Extended inputs
! Extended inputs -- Linearization is only possible with Steady or Uniform Wind, so take advantage of that here
! Module/Mesh/Field: HWindSpeed = 37
! Module/Mesh/Field: PLexp = 38
! Module/Mesh/Field: PropagationDir = 39
!FIXME: Extended inputs
call IfW_UniformWind_GetOP(p_AD%FlowField%Uniform, t, .false. , OP_out)
! HWindSpeed
u_op(index) = OP_out(1); index = index + 1
! PLexp
u_op(index) = OP_out(2); index = index + 1
! PropagationDir (include AngleH in calculation if any)
u_op(index) = OP_out(3) + p_AD%FlowField%PropagationDir; index = index + 1

end if
END IF

Expand Down Expand Up @@ -6658,7 +6667,7 @@ SUBROUTINE Init_Jacobian_u( InputFileData, p, p_AD, u, InitOut, ErrStat, ErrMsg)
p%Jac_u_idxStartList%Tower = index
call SetJac_u_idx(6,9,u%TowerMotion%NNodes,index)
! Perturbations
p%du(5) = perturb_t
p%du(6) = perturb_t
p%du(7) = perturb
p%du(8) = perturb_t
p%du(9) = perturb_t
Expand Down Expand Up @@ -6707,7 +6716,7 @@ SUBROUTINE Init_Jacobian_u( InputFileData, p, p_AD, u, InitOut, ErrStat, ErrMsg)
p%Jac_u_idxStartList%Blade = index
call SetJac_u_idx(13,18,u%BladeMotion(1)%NNodes,index)
if (p%NumBl_Lin > 1) call SetJac_u_idx(19,24,u%BladeMotion(2)%NNodes,index)
if (p%NumBl_Lin > 2) call SetJac_u_idx(15,30,u%BladeMotion(3)%NNodes,index)
if (p%NumBl_Lin > 2) call SetJac_u_idx(25,30,u%BladeMotion(3)%NNodes,index)
! Perturbations
do k=1,p%NumBl_Lin
p%du(13 + (k-1)*6) = perturb_b(k)
Expand Down Expand Up @@ -6798,18 +6807,14 @@ SUBROUTINE Init_Jacobian_u( InputFileData, p, p_AD, u, InitOut, ErrStat, ErrMsg)
! Module/Mesh/Field: PLexp = 38
! Module/Mesh/Field: PropagationDir = 39
p%Jac_u_idxStartList%Extended = index
p%Jac_u_indx(index,1)=37; p%Jac_u_indx(index,2)=1; p%Jac_u_indx(index,3)=1; index = index + 1
p%Jac_u_indx(index,1)=38; p%Jac_u_indx(index,2)=1; p%Jac_u_indx(index,3)=1; index = index + 1
p%Jac_u_indx(index,1)=39; p%Jac_u_indx(index,2)=1; p%Jac_u_indx(index,3)=1; index = index + 1
p%Jac_u_indx(index,1)=37; p%Jac_u_indx(index,2)=1; p%Jac_u_indx(index,3)=1; InitOut%LinNames_u(index) = 'Extended input: horizontal wind speed (steady/uniform wind), m/s'; index=index+1
p%Jac_u_indx(index,1)=38; p%Jac_u_indx(index,2)=1; p%Jac_u_indx(index,3)=1; InitOut%LinNames_u(index) = 'Extended input: vertical power-law shear exponent, -'; index=index+1
p%Jac_u_indx(index,1)=39; p%Jac_u_indx(index,2)=1; p%Jac_u_indx(index,3)=1; InitOut%LinNames_u(index) = 'Extended input: propagation direction, rad'; index=index+1
! Perturbations
p%du(37) = perturb
p%du(38) = perturb
p%du(39) = perturb
! Names
InitOut%LinNames_u(index) = 'Extended input: horizontal wind speed (steady/uniform wind), m/s'; index=index+1
InitOut%LinNames_u(index) = 'Extended input: vertical power-law shear exponent, -'; index=index+1
InitOut%LinNames_u(index) = 'Extended input: propagation direction, rad'; index=index+1


end if ! .not. compAeroMaps

contains
Expand Down Expand Up @@ -7063,7 +7068,7 @@ SUBROUTINE Perturb_u( p, n, perturb_sign, u, du )
case(15); u%BladeMotion(1)%TranslationVel( fieldIndx,node) = u%BladeMotion(1)%TranslationVel(fieldIndx,node) + du * perturb_sign
case(16); u%BladeMotion(1)%RotationVel( fieldIndx,node) = u%BladeMotion(1)%RotationVel(fieldIndx,node) + du * perturb_sign
case(17); u%BladeMotion(1)%TranslationAcc( fieldIndx,node) = u%BladeMotion(1)%TranslationAcc(fieldIndx,node) + du * perturb_sign
case(18); u%BladeMotion(1)%RotationAcc( fieldIndx,node) = u%BladeMotion(1)%TranslationAcc(fieldIndx,node) + du * perturb_sign
case(18); u%BladeMotion(1)%RotationAcc( fieldIndx,node) = u%BladeMotion(1)%RotationAcc( fieldIndx,node) + du * perturb_sign

! Blade 2
! Module/Mesh/Field: u%BladeMotion(2)%TranslationDisp = 19;
Expand All @@ -7077,7 +7082,7 @@ SUBROUTINE Perturb_u( p, n, perturb_sign, u, du )
case(21); u%BladeMotion(2)%TranslationVel( fieldIndx,node) = u%BladeMotion(2)%TranslationVel(fieldIndx,node) + du * perturb_sign
case(22); u%BladeMotion(2)%RotationVel( fieldIndx,node) = u%BladeMotion(2)%RotationVel(fieldIndx,node) + du * perturb_sign
case(23); u%BladeMotion(2)%TranslationAcc( fieldIndx,node) = u%BladeMotion(2)%TranslationAcc(fieldIndx,node) + du * perturb_sign
case(24); u%BladeMotion(2)%RotationAcc( fieldIndx,node) = u%BladeMotion(2)%TranslationAcc(fieldIndx,node) + du * perturb_sign
case(24); u%BladeMotion(2)%RotationAcc( fieldIndx,node) = u%BladeMotion(2)%RotationAcc( fieldIndx,node) + du * perturb_sign

! Blade 3
! Module/Mesh/Field: u%BladeMotion(3)%TranslationDisp = 25;
Expand All @@ -7091,7 +7096,7 @@ SUBROUTINE Perturb_u( p, n, perturb_sign, u, du )
case(27); u%BladeMotion(3)%TranslationVel( fieldIndx,node) = u%BladeMotion(3)%TranslationVel(fieldIndx,node) + du * perturb_sign
case(28); u%BladeMotion(3)%RotationVel( fieldIndx,node) = u%BladeMotion(3)%RotationVel(fieldIndx,node) + du * perturb_sign
case(29); u%BladeMotion(3)%TranslationAcc( fieldIndx,node) = u%BladeMotion(3)%TranslationAcc(fieldIndx,node) + du * perturb_sign
case(30); u%BladeMotion(3)%RotationAcc( fieldIndx,node) = u%BladeMotion(3)%TranslationAcc(fieldIndx,node) + du * perturb_sign
case(30); u%BladeMotion(3)%RotationAcc( fieldIndx,node) = u%BladeMotion(3)%RotationAcc( fieldIndx,node) + du * perturb_sign

! TailFin
! Module/Mesh/Field: u%TFinMotion%TranslationDisp = 31;
Expand Down
3 changes: 2 additions & 1 deletion modules/inflowwind/src/IfW_FlowField.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -716,7 +716,7 @@ subroutine IfW_UniformWind_GetOP(UF, t, InterpCubic, OP_out)
type(UniformFieldType), intent(IN) :: UF !< Parameters
real(DbKi), intent(IN) :: t !< Current simulation time in seconds
logical, intent(in) :: InterpCubic !< flag for using cubic interpolation
real(ReKi), intent(OUT) :: OP_out(2) !< operating point (HWindSpeed and PLexp
real(ReKi), intent(OUT) :: OP_out(3) !< operating point (HWindSpeed, PLexp, and AngleH)

type(UniformField_Interp) :: op ! interpolated values of InterpParams

Expand All @@ -729,6 +729,7 @@ subroutine IfW_UniformWind_GetOP(UF, t, InterpCubic, OP_out)

OP_out(1) = op%VelH
OP_out(2) = op%ShrV
OP_out(3) = op%AngleH

end subroutine

Expand Down
10 changes: 6 additions & 4 deletions modules/inflowwind/src/InflowWind.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -814,8 +814,10 @@ SUBROUTINE InflowWind_JacobianPInput( t, u, p, x, xd, z, OtherState, y, m, ErrSt
IF ( PRESENT( dYdu ) ) THEN

! If dYdu is allocated, make sure it is the correct size
if (size(dYdu,1) /= NumExtendedIO + p%NumOuts) deallocate (dYdu)
if (size(dYdu,2) /= NumExtendedIO) deallocate (dYdu)
if (allocated(dYdu)) then
if (size(dYdu,1) /= NumExtendedIO + p%NumOuts) deallocate (dYdu)
if (size(dYdu,2) /= NumExtendedIO) deallocate (dYdu)
endif

! Calculate the partial derivative of the output functions (Y) with respect to the inputs (u) here:
! - inputs are extended inputs only
Expand Down Expand Up @@ -1094,8 +1096,8 @@ SUBROUTINE InflowWind_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMs
if (ErrStat >= AbortErrLev) return
end if

call IfW_UniformWind_GetOP( p%FlowField%Uniform, t, p%FlowField%VelInterpCubic, u_op(1:2) )
u_op(3) = p%FlowField%PropagationDir
call IfW_UniformWind_GetOP( p%FlowField%Uniform, t, p%FlowField%VelInterpCubic, u_op )
u_op(3) = p%FlowField%PropagationDir + u_op(3) ! include the AngleH from Uniform Wind input files
end if

if ( PRESENT( y_op ) ) then
Expand Down
72 changes: 1 addition & 71 deletions modules/openfast-library/src/FAST_Lin.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1510,13 +1510,7 @@ SUBROUTINE Glue_Jacobians( p_FAST, y_FAST, m_FAST, ED, BD, SrvD, AD, IfW, ExtInf
end do


!............
! \f$ \frac{\partial U_\Lambda^{IfW}}{\partial u^{AD}} \end{bmatrix} = \f$ (dUdu block row 1=IfW)
!............
IF (p_FAST%CompInflow == MODULE_IfW .and. p_FAST%CompAero == MODULE_AD) THEN
call Linear_IfW_InputSolve_du_AD( p_FAST, y_FAST, AD%Input(1), dUdu )
end if ! we're using the InflowWind module


!............
! \f$ \frac{\partial U_\Lambda^{SrvD}}{\partial u^{SrvD}} \end{bmatrix} = \f$ (dUdu block row 2=SrvD)
!............
Expand Down Expand Up @@ -1702,70 +1696,6 @@ SUBROUTINE Glue_Jacobians( p_FAST, y_FAST, m_FAST, ED, BD, SrvD, AD, IfW, ExtInf

END SUBROUTINE Glue_Jacobians


!----------------------------------------------------------------------------------------------------------------------------------
!> This routine forms the dU^{IfW}/du^{AD} block of dUdu. (i.e., how do changes in the AD inputs affect IfW inputs?)
SUBROUTINE Linear_IfW_InputSolve_du_AD( p_FAST, y_FAST, u_AD, dUdu )
!FIXME: this routine needs revision (info passed is no longer correct)
TYPE(FAST_ParameterType), INTENT(IN ) :: p_FAST !< FAST parameter data
TYPE(FAST_OutputFileType), INTENT(IN ) :: y_FAST !< FAST output data (for linearization)
TYPE(AD_InputType), INTENT(IN) :: u_AD !< The input meshes (already calculated) from AeroDyn
REAL(R8Ki), INTENT(INOUT) :: dUdu(:,:) !< Jacobian matrix of which we are computing the dU^(IfW)/du^(AD) block


INTEGER(IntKi) :: i, j, k ! loop counters
INTEGER(IntKi) :: i2, j2 ! loop counters
INTEGER(IntKi) :: AD_Start_Bl ! starting index of dUdu (column) where AD blade motion inputs are located
INTEGER(IntKi) :: Node ! InflowWind node number


! compare with IfW_InputSolve():

Node = 0 !InflowWind node
if (p_FAST%CompServo == MODULE_SrvD) Node = Node + 1

IF (p_FAST%CompAero == MODULE_AD) THEN

! blades:
AD_Start_Bl = y_FAST%Lin%Modules(MODULE_AD)%Instance(1)%LinStartIndx(LIN_INPUT_COL) &
+ u_AD%rotors(1)%TowerMotion%NNodes * 9 & ! 3 fields (MASKID_TRANSLATIONDISP,MASKID_Orientation,MASKID_TRANSLATIONVel) with 3 components
+ u_AD%rotors(1)%HubMotion%NNodes * 9 ! 3 fields (MASKID_TRANSLATIONDISP,MASKID_Orientation,MASKID_RotationVel) with 3 components

do k = 1,size(u_AD%rotors(1)%BladeRootMotion)
AD_Start_Bl = AD_Start_Bl + u_AD%rotors(1)%BladeRootMotion(k)%NNodes * 3 ! 1 field (MASKID_Orientation) with 3 components
end do
! next is u_AD%rotors(1)%BladeMotion(k):

DO K = 1,SIZE(u_AD%rotors(1)%BladeMotion)
DO J = 1,u_AD%rotors(1)%BladeMotion(k)%Nnodes
Node = Node + 1 ! InflowWind node
do i=1,3 !XYZ components of this node
i2 = y_FAST%Lin%Modules(MODULE_IfW)%Instance(1)%LinStartIndx(LIN_INPUT_COL) + (Node-1)*3 + i - 1
j2 = AD_Start_Bl + (j-1)*3 + i - 1
dUdu( i2, j2 ) = -1.0_R8Ki
end do
END DO !J = 1,p%BldNodes ! Loop through the blade nodes / elements

! get starting AD index of BladeMotion for next blade
AD_Start_Bl = AD_Start_Bl + u_AD%rotors(1)%BladeMotion(k)%Nnodes * 9 ! 3 fields (MASKID_TRANSLATIONDISP,MASKID_Orientation,MASKID_TRANSLATIONVel) with 3 components
END DO !K = 1,p%NumBl

! tower:
DO J=1,u_AD%rotors(1)%TowerMotion%nnodes
Node = Node + 1
do i=1,3 !XYZ components of this node
i2 = y_FAST%Lin%Modules(MODULE_IfW)%Instance(1)%LinStartIndx(LIN_INPUT_COL) + (Node-1)*3 + i - 1
j2 = y_FAST%Lin%Modules(MODULE_AD )%Instance(1)%LinStartIndx(LIN_INPUT_COL) + (j-1)*3 + i - 1
dUdu( i2, j2 ) = -1.0_R8Ki
end do
END DO

! HubPosition and HubOrientation from ElastoDyn are missing from this
!FIXME: add in the extended inputs here
END IF
END SUBROUTINE Linear_IfW_InputSolve_du_AD


!----------------------------------------------------------------------------------------------------------------------------------
!> This routine forms the dU^{ED}/du^{BD} and dU^{ED}/du^{AD} blocks (ED row) of dUdu. (i.e., how do changes in the AD and BD inputs affect the ED inputs?)
SUBROUTINE Linear_ED_InputSolve_du( p_FAST, y_FAST, SrvD, u_ED, y_ED, y_AD, u_AD, p_AD, BD, HD, SD, MAPp, MD, MeshMapData, dUdu, ErrStat, ErrMsg )
Expand Down
4 changes: 2 additions & 2 deletions reg_tests/CTestList.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -336,8 +336,8 @@ of_regression_py("EllipticalWing_OLAF_py" "openfast;fastlib;p
of_regression_aeroacoustic("IEA_LB_RWT-AeroAcoustics" "openfast;aerodyn15;aeroacoustics")

# Linearized OpenFAST regression tests
# of_regression_linear("Fake5MW_AeroLin_B1_UA4_DBEMT3" "openfast;linear;elastodyn") #Also: aerodyn
# of_regression_linear("Fake5MW_AeroLin_B3_UA6" "openfast;linear;elastodyn") #Also: aerodyn
of_regression_linear("Fake5MW_AeroLin_B1_UA4_DBEMT3" "openfast;linear;elastodyn") #Also: aerodyn
of_regression_linear("Fake5MW_AeroLin_B3_UA6" "openfast;linear;elastodyn") #Also: aerodyn
of_regression_linear("WP_Stationary_Linear" "openfast;linear;elastodyn")
of_regression_linear("Ideal_Beam_Fixed_Free_Linear" "openfast;linear;beamdyn")
of_regression_linear("Ideal_Beam_Free_Free_Linear" "openfast;linear;beamdyn")
Expand Down

0 comments on commit 617c39e

Please sign in to comment.