Skip to content

Commit

Permalink
Update orientation differences in linear trim solution (#1158)
Browse files Browse the repository at this point in the history 
- Changed name of optional parameter in`ModName_GetOP` routine from `NeedLogMap` to `NeedPackedOrient` in an attempt to clarify what that does (returns the orientation matrix as a 3-valued vector instead of 9-valued one)
- Made the use of this value consistent in ElastoDyn (all meshes with orientation fields that are packed in the GetOP routine need to have the same value for their  call to `PackMotionMesh`)
- Added the optional `NeedPackedOrient` input to `SD_GetOP`
- Made ED, BD, and SD each pack their arrays the same way in `FAST_DiffInterpOutputs`; otherwise, the linear trim solution does not difference all of the outputs correctly to determine if it has converged
  • Loading branch information
@bjonkman
bjonkman authored Jul 12, 2022
1 parent 9bac504 commit 5206df0
Show file tree
Hide file tree
Showing 4 changed files with 36 additions and 19 deletions.
10 changes: 6 additions & 4 deletions modules/beamdyn/src/BeamDyn.f90
View file
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
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
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
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

0 comments on commit 5206df0

Please sign in to comment.