Skip to content
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

Fix for linear trim solution #1158

Merged
merged 2 commits into from
Jul 12, 2022
Merged

Fix for linear trim solution #1158

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ea9c6f3
Update orientation differences in linear trim solution
bjonkman Jun 14, 2022
8fcbdaf
Merge remote-tracking branch 'NREL/dev' into b/TrimLinear
bjonkman Jun 20, 2022
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
10 changes: 6 additions & 4 deletions modules/beamdyn/src/BeamDyn.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6519,7 +6519,7 @@ SUBROUTINE BD_JacobianPConstrState( t, u, p, x, xd, z, OtherState, y, m, ErrStat
END SUBROUTINE BD_JacobianPConstrState
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Routine to pack the data structures representing the operating points into arrays for linearization.
SUBROUTINE BD_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op, y_op, x_op, dx_op, xd_op, z_op, NeedLogMap )
SUBROUTINE BD_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op, y_op, x_op, dx_op, xd_op, z_op, NeedPackedOrient )

REAL(DbKi), INTENT(IN ) :: t !< Time in seconds at operating point
TYPE(BD_InputType), INTENT(INOUT) :: u !< Inputs at operating point (may change to inout if a mesh copy is required)
Expand All @@ -6538,7 +6538,7 @@ SUBROUTINE BD_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op,
REAL(ReKi), ALLOCATABLE, OPTIONAL, INTENT(INOUT) :: dx_op(:) !< values of first time derivatives of linearized continuous states
REAL(ReKi), ALLOCATABLE, OPTIONAL, INTENT(INOUT) :: xd_op(:) !< values of linearized discrete states
REAL(ReKi), ALLOCATABLE, OPTIONAL, INTENT(INOUT) :: z_op(:) !< values of linearized constraint states
LOGICAL, OPTIONAL, INTENT(IN ) :: NeedLogMap !< whether a y_op values should contain log maps instead of full orientation matrices
LOGICAL, OPTIONAL, INTENT(IN ) :: NeedPackedOrient !< whether a y_op values should contain 3-value representation instead of full orientation matrices

INTEGER(IntKi) :: index, i, dof
INTEGER(IntKi) :: nu
Expand Down Expand Up @@ -6584,8 +6584,9 @@ SUBROUTINE BD_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op,


IF ( PRESENT( y_op ) ) THEN
if (present(NeedLogMap)) then
ReturnSmallAngle = NeedLogMap
! Only the y operating points need to potentially return a smaller array than the "normal" call to this return. In the trim solution, we use a smaller array for y.
if (present(NeedPackedOrient)) then
ReturnSmallAngle = NeedPackedOrient
else
ReturnSmallAngle = .false.
end if
Expand All @@ -6598,6 +6599,7 @@ SUBROUTINE BD_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op,
if (ErrStat >= AbortErrLev) return
end if

if (ReturnSmallAngle) y_op = 0.0_ReKi ! initialize in case we are returning packed orientations and don't fill the entire array

index = 1
call PackLoadMesh(y%ReactionForce, y_op, index)
Expand Down
19 changes: 10 additions & 9 deletions modules/elastodyn/src/ElastoDyn.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11668,7 +11668,7 @@ SUBROUTINE Compute_dY(p, y_p, y_m, delta, dY)
END SUBROUTINE Compute_dY
!----------------------------------------------------------------------------------------------------------------------------------
!> Routine to pack the data structures representing the operating points into arrays for linearization.
SUBROUTINE ED_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op, y_op, x_op, dx_op, xd_op, z_op, NeedLogMap )
SUBROUTINE ED_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op, y_op, x_op, dx_op, xd_op, z_op, NeedPackedOrient )

REAL(DbKi), INTENT(IN ) :: t !< Time in seconds at operating point
TYPE(ED_InputType), INTENT(IN ) :: u !< Inputs at operating point (may change to inout if a mesh copy is required)
Expand All @@ -11687,7 +11687,7 @@ SUBROUTINE ED_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op,
REAL(ReKi), ALLOCATABLE, OPTIONAL, INTENT(INOUT) :: dx_op(:) !< values of first time derivatives of linearized continuous states
REAL(ReKi), ALLOCATABLE, OPTIONAL, INTENT(INOUT) :: xd_op(:) !< values of linearized discrete states
REAL(ReKi), ALLOCATABLE, OPTIONAL, INTENT(INOUT) :: z_op(:) !< values of linearized constraint states
LOGICAL, OPTIONAL, INTENT(IN ) :: NeedLogMap !< whether a y_op values should contain log maps instead of full orientation matrices
LOGICAL, OPTIONAL, INTENT(IN ) :: NeedPackedOrient !< whether a y_op values should contain 3-value representation instead of full orientation matrices



Expand Down Expand Up @@ -11747,8 +11747,8 @@ SUBROUTINE ED_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op,

!..................................
IF ( PRESENT( y_op ) ) THEN
if (present(NeedLogMap)) then
ReturnLogMap = NeedLogMap
if (present(NeedPackedOrient)) then
ReturnLogMap = NeedPackedOrient
else
ReturnLogMap = .false.
end if
Expand All @@ -11774,7 +11774,8 @@ SUBROUTINE ED_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op,
if (ErrStat>=AbortErrLev) return
end if


if (ReturnLogMap) y_op = 0.0_ReKi ! initialize in case we are returning packed orientations and don't fill the entire array


Mask = .false.
Mask(MASKID_TRANSLATIONDISP) = .true.
Expand All @@ -11784,16 +11785,16 @@ SUBROUTINE ED_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op,
index = 1
if (allocated(y%BladeLn2Mesh)) then
do k=1,p%NumBl
call PackMotionMesh(y%BladeLn2Mesh(k), y_op, index, UseSmlAngle=.false.)
call PackMotionMesh(y%BladeLn2Mesh(k), y_op, index, UseSmlAngle=ReturnLogMap)
end do
end if
call PackMotionMesh(y%PlatformPtMesh, y_op, index, UseSmlAngle=ReturnLogMap)
call PackMotionMesh(y%TowerLn2Mesh, y_op, index, UseSmlAngle=ReturnLogMap)
call PackMotionMesh(y%HubPtMotion, y_op, index, FieldMask=Mask, UseSmlAngle=.false.)
call PackMotionMesh(y%HubPtMotion, y_op, index, FieldMask=Mask, UseSmlAngle=ReturnLogMap)
do k=1,p%NumBl
call PackMotionMesh(y%BladeRootMotion(k), y_op, index, UseSmlAngle=.false.)
call PackMotionMesh(y%BladeRootMotion(k), y_op, index, UseSmlAngle=ReturnLogMap)
end do
call PackMotionMesh(y%NacelleMotion, y_op, index, UseSmlAngle=.false.)
call PackMotionMesh(y%NacelleMotion, y_op, index, UseSmlAngle=ReturnLogMap)

y_op(index) = y%Yaw ; index = index + 1
y_op(index) = y%YawRate ; index = index + 1
Expand Down
6 changes: 3 additions & 3 deletions modules/openfast-library/src/FAST_Lin.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6405,7 +6405,7 @@ SUBROUTINE FAST_DiffInterpOutputs( psi_target, p_FAST, y_FAST, m_FAST, ED, BD, S
CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName )

call ED_GetOP( t_global, ED%Input(1), ED%p, ED%x(STATE_CURR), ED%xd(STATE_CURR), ED%z(STATE_CURR), ED%OtherSt(STATE_CURR), &
ED%y_interp, ED%m, ErrStat2, ErrMsg2, y_op=y_FAST%Lin%Modules(Module_ED)%Instance(1)%op_y, NeedLogMap=.true.)
ED%y_interp, ED%m, ErrStat2, ErrMsg2, y_op=y_FAST%Lin%Modules(Module_ED)%Instance(1)%op_y, NeedPackedOrient=.true.)
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)

! BeamDyn
Expand All @@ -6417,7 +6417,7 @@ SUBROUTINE FAST_DiffInterpOutputs( psi_target, p_FAST, y_FAST, m_FAST, ED, BD, S
CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName )

call BD_GetOP( t_global, BD%Input(1,k), BD%p(k), BD%x(k,STATE_CURR), BD%xd(k,STATE_CURR), BD%z(k,STATE_CURR), BD%OtherSt(k,STATE_CURR), &
BD%y_interp(k), BD%m(k), ErrStat2, ErrMsg2, y_op=y_FAST%Lin%Modules(Module_BD)%Instance(k)%op_y, NeedLogMap=.false.)
BD%y_interp(k), BD%m(k), ErrStat2, ErrMsg2, y_op=y_FAST%Lin%Modules(Module_BD)%Instance(k)%op_y, NeedPackedOrient=.true.)
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
END DO ! k=p_FAST%nBeams

Expand Down Expand Up @@ -6478,7 +6478,7 @@ SUBROUTINE FAST_DiffInterpOutputs( psi_target, p_FAST, y_FAST, m_FAST, ED, BD, S
CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName )

call SD_GetOP( t_global, SD%Input(1), SD%p, SD%x(STATE_CURR), SD%xd(STATE_CURR), SD%z(STATE_CURR), SD%OtherSt(STATE_CURR), &
SD%y_interp, SD%m, ErrStat2, ErrMsg2, y_op=y_FAST%Lin%Modules(Module_SD)%Instance(1)%op_y)
SD%y_interp, SD%m, ErrStat2, ErrMsg2, y_op=y_FAST%Lin%Modules(Module_SD)%Instance(1)%op_y, NeedPackedOrient=.true.)
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
ELSE IF ( p_FAST%CompSub == Module_ExtPtfm ) THEN
END IF ! SubDyn/ExtPtfm_MCKF
Expand Down
20 changes: 17 additions & 3 deletions modules/subdyn/src/SubDyn.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2187,7 +2187,7 @@ SUBROUTINE SD_JacobianPConstrState( t, u, p, x, xd, z, OtherState, y, m, ErrStat
END SUBROUTINE SD_JacobianPConstrState
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Routine to pack the data structures representing the operating points into arrays for linearization.
SUBROUTINE SD_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op, y_op, x_op, dx_op, xd_op, z_op )
SUBROUTINE SD_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op, y_op, x_op, dx_op, xd_op, z_op, NeedPackedOrient )
REAL(DbKi), INTENT(IN ) :: t !< Time in seconds at operating point
TYPE(SD_InputType), INTENT(INOUT) :: u !< Inputs at operating point (may change to inout if a mesh copy is required)
TYPE(SD_ParameterType), INTENT(IN ) :: p !< Parameters
Expand All @@ -2205,8 +2205,11 @@ SUBROUTINE SD_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op,
REAL(ReKi), ALLOCATABLE, OPTIONAL, INTENT(INOUT) :: dx_op(:) !< values of first time derivatives of linearized continuous states
REAL(ReKi), ALLOCATABLE, OPTIONAL, INTENT(INOUT) :: xd_op(:) !< values of linearized discrete states
REAL(ReKi), ALLOCATABLE, OPTIONAL, INTENT(INOUT) :: z_op(:) !< values of linearized constraint states
LOGICAL, OPTIONAL, INTENT(IN ) :: NeedPackedOrient !< whether a y_op values should contain 3-value representation instead of full orientation matrices

! Local
INTEGER(IntKi) :: idx, i
LOGICAL :: ReturnPackedOrientation
INTEGER(IntKi) :: nu
INTEGER(IntKi) :: ny
INTEGER(IntKi) :: ErrStat2
Expand All @@ -2232,11 +2235,21 @@ SUBROUTINE SD_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op,
call PackMotionMesh(u%TPMesh, u_op, idx, FieldMask=FieldMask, UseSmlAngle=.true.)
call PackLoadMesh(u%LMesh, u_op, idx)
END IF

IF ( PRESENT( y_op ) ) THEN
ny = p%Jac_ny + y%Y2Mesh%NNodes * 6 + y%Y3Mesh%NNodes * 6 ! Jac_ny has 3 orientation angles, but the OP needs the full 9 elements of the DCM (thus 6 more per node)
if (.not. allocated(y_op)) then
call AllocAry(y_op, ny, 'y_op', ErrStat2, ErrMsg2); if(Failed()) return
end if

if (present(NeedPackedOrient)) then
ReturnPackedOrientation = NeedPackedOrient
else
ReturnPackedOrientation = .false.
end if

if (ReturnPackedOrientation) y_op = 0.0_ReKi ! initialize in case we are returning packed orientations and don't fill the entire array

idx = 1
call PackLoadMesh(y%Y1Mesh, y_op, idx)
FieldMask = .false.
Expand All @@ -2246,13 +2259,14 @@ SUBROUTINE SD_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op,
FieldMask(MASKID_RotationVel) = .true.
FieldMask(MASKID_TranslationAcc) = .true.
FieldMask(MASKID_RotationAcc) = .true.
call PackMotionMesh(y%Y2Mesh, y_op, idx, FieldMask=FieldMask, UseSmlAngle=.true.)
call PackMotionMesh(y%Y3Mesh, y_op, idx, FieldMask=FieldMask, UseSmlAngle=.true.)
call PackMotionMesh(y%Y2Mesh, y_op, idx, FieldMask=FieldMask, UseSmlAngle=ReturnPackedOrientation)
call PackMotionMesh(y%Y3Mesh, y_op, idx, FieldMask=FieldMask, UseSmlAngle=ReturnPackedOrientation)
idx = idx - 1
do i=1,p%NumOuts
y_op(i+idx) = y%WriteOutput(i)
end do
END IF

IF ( PRESENT( x_op ) ) THEN
if (.not. allocated(x_op)) then
call AllocAry(x_op, p%Jac_nx*2,'x_op',ErrStat2,ErrMsg2); if (Failed()) return
Expand Down