From 617c39e5ce991a0a9db24c984645a94f9787a8be Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Mon, 29 Jan 2024 13:39:30 -0700 Subject: [PATCH] Lin: fix indexing in AD15 lin, and add access to IfW OP --- modules/aerodyn/src/AeroDyn.f90 | 37 +++++++----- modules/inflowwind/src/IfW_FlowField.f90 | 3 +- modules/inflowwind/src/InflowWind.f90 | 10 ++-- modules/openfast-library/src/FAST_Lin.f90 | 72 +---------------------- reg_tests/CTestList.cmake | 4 +- 5 files changed, 32 insertions(+), 94 deletions(-) diff --git a/modules/aerodyn/src/AeroDyn.f90 b/modules/aerodyn/src/AeroDyn.f90 index 3fe3b3cbb2..19af647680 100644 --- a/modules/aerodyn/src/AeroDyn.f90 +++ b/modules/aerodyn/src/AeroDyn.f90 @@ -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 @@ -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 ) @@ -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 @@ -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 @@ -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 @@ -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) @@ -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 @@ -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; @@ -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; @@ -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; diff --git a/modules/inflowwind/src/IfW_FlowField.f90 b/modules/inflowwind/src/IfW_FlowField.f90 index 4232dca0c2..dd3b0c31c6 100644 --- a/modules/inflowwind/src/IfW_FlowField.f90 +++ b/modules/inflowwind/src/IfW_FlowField.f90 @@ -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 @@ -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 diff --git a/modules/inflowwind/src/InflowWind.f90 b/modules/inflowwind/src/InflowWind.f90 index 6d2320d46f..8899965c79 100644 --- a/modules/inflowwind/src/InflowWind.f90 +++ b/modules/inflowwind/src/InflowWind.f90 @@ -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 @@ -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 diff --git a/modules/openfast-library/src/FAST_Lin.f90 b/modules/openfast-library/src/FAST_Lin.f90 index 6dabb3d5ac..7035d89d07 100644 --- a/modules/openfast-library/src/FAST_Lin.f90 +++ b/modules/openfast-library/src/FAST_Lin.f90 @@ -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) !............ @@ -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 ) diff --git a/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index bbe017d88f..c75d50398b 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -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")