Skip to content

Commit

Permalink
Linear Trim Solution Improvements
Browse files Browse the repository at this point in the history 
- The operating point used for trim solutions cannot use GetSmllRotAngs since we can't assume all angles in the simulation are less than 25 degrees. It reverted back to log maps for this calculation.
- I removed the rotational accelerations from the packed arrays for trim solution. These can be noisy outputs, and we are already making sure that the position and velocity values are converging.
- I updated the argument names in the `GetOP` routines for the structural codes that output motions.

This is a follow-up to OpenFAST#1158.
  • Loading branch information
@bjonkman
bjonkman committed Oct 6, 2022
1 parent a992760 commit b023a36
Show file tree
Hide file tree
Showing 5 changed files with 72 additions and 58 deletions.
16 changes: 8 additions & 8 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, NeedPackedOrient )
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, NeedTrimOP )

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 ) :: NeedPackedOrient !< whether a y_op values should contain 3-value representation instead of full orientation matrices
LOGICAL, OPTIONAL, INTENT(IN ) :: NeedTrimOP !< whether a y_op values should contain values for trim solution (3-value representation instead of full orientation matrices, no rotation acc)

INTEGER(IntKi) :: index, i, dof
INTEGER(IntKi) :: nu
Expand All @@ -6547,7 +6547,7 @@ SUBROUTINE BD_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op,
CHARACTER(ErrMsgLen) :: ErrMsg2
CHARACTER(*), PARAMETER :: RoutineName = 'BD_GetOP'
LOGICAL :: FieldMask(FIELDMASK_SIZE)
LOGICAL :: ReturnSmallAngle
LOGICAL :: ReturnTrimOP
TYPE(BD_ContinuousStateType) :: dx ! derivative of continuous states at operating point


Expand Down Expand Up @@ -6585,10 +6585,10 @@ SUBROUTINE BD_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op,

IF ( PRESENT( y_op ) ) THEN
! 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
if (present(NeedTrimOP)) then
ReturnTrimOP = NeedTrimOP
else
ReturnSmallAngle = .false.
ReturnTrimOP = .false.
end if

if (.not. allocated(y_op)) then
Expand All @@ -6599,7 +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
if (ReturnTrimOP) 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 All @@ -6611,7 +6611,7 @@ SUBROUTINE BD_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%BldMotion, y_op, index, FieldMask=FieldMask, UseSmlAngle=ReturnSmallAngle)
call PackMotionMesh(y%BldMotion, y_op, index, FieldMask=FieldMask, TrimOP=ReturnTrimOP)

index = index - 1
do i=1,p%NumOuts + p%BldNd_TotNumOuts
Expand Down
32 changes: 18 additions & 14 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, NeedPackedOrient )
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, NeedTrimOP )

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 ) :: NeedPackedOrient !< whether a y_op values should contain 3-value representation instead of full orientation matrices
LOGICAL, OPTIONAL, INTENT(IN ) :: NeedTrimOP !< whether a y_op values should contain values for trim solution (3-value representation instead of full orientation matrices, no rotation acc)



Expand All @@ -11696,7 +11696,7 @@ SUBROUTINE ED_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op,
INTEGER(IntKi) :: ErrStat2
CHARACTER(ErrMsgLen) :: ErrMsg2
CHARACTER(*), PARAMETER :: RoutineName = 'ED_GetOP'
LOGICAL :: ReturnLogMap
LOGICAL :: ReturnTrimOP
TYPE(ED_ContinuousStateType) :: dx !< derivative of continuous states at operating point
LOGICAL :: Mask(FIELDMASK_SIZE) !< flags to determine if this field is part of the packing

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

!..................................
IF ( PRESENT( y_op ) ) THEN
if (present(NeedPackedOrient)) then
ReturnLogMap = NeedPackedOrient
if (present(NeedTrimOP)) then
ReturnTrimOP = NeedTrimOP
else
ReturnLogMap = .false.
ReturnTrimOP = .false.
end if

if (.not. allocated(y_op)) then
Expand All @@ -11774,7 +11774,7 @@ 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
if (ReturnTrimOP) y_op = 0.0_ReKi ! initialize in case we are returning packed orientations and don't fill the entire array


Mask = .false.
Expand All @@ -11785,24 +11785,28 @@ 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=ReturnLogMap)
call PackMotionMesh(y%BladeLn2Mesh(k), y_op, index, TrimOP=ReturnTrimOP)
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=ReturnLogMap)
call PackMotionMesh(y%PlatformPtMesh, y_op, index, TrimOP=ReturnTrimOP)
call PackMotionMesh(y%TowerLn2Mesh, y_op, index, TrimOP=ReturnTrimOP)
call PackMotionMesh(y%HubPtMotion, y_op, index, TrimOP=ReturnTrimOP)

do k=1,p%NumBl
call PackMotionMesh(y%BladeRootMotion(k), y_op, index, UseSmlAngle=ReturnLogMap)
call PackMotionMesh(y%BladeRootMotion(k), y_op, index, TrimOP=ReturnTrimOP)
end do
call PackMotionMesh(y%NacelleMotion, y_op, index, UseSmlAngle=ReturnLogMap)
call PackMotionMesh(y%PlatformPtMesh, y_op, index, TrimOP=ReturnTrimOP)
call PackMotionMesh(y%TowerLn2Mesh, y_op, index, TrimOP=ReturnTrimOP)
call PackMotionMesh(y%HubPtMotion, y_op, index, FieldMask=Mask, TrimOP=ReturnTrimOP)
call PackMotionMesh(y%NacelleMotion, y_op, index, TrimOP=ReturnTrimOP)

y_op(index) = y%Yaw ; index = index + 1
y_op(index) = y%YawRate ; index = index + 1
y_op(index) = y%HSS_Spd

do i=1,p%NumOuts + p%BldNd_TotNumOuts
y_op(i+index) = y%WriteOutput(i)
end do
end do

END IF

Expand Down
54 changes: 32 additions & 22 deletions modules/nwtc-library/src/ModMesh.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3072,24 +3072,24 @@ SUBROUTINE PackMotionMesh_Names(M, MeshName, Names, indx_first, FieldMask)


END SUBROUTINE PackMotionMesh_Names
!...............................................................................................................................
!> This subroutine returns the operating point values of the mesh fields. It assumes all fields marked
!! by FieldMask are allocated; Some fields may be allocated by the ModMesh module and not used in
!! the linearization procedure, thus I am not using the check if they are allocated to determine if they should be included.
SUBROUTINE PackMotionMesh(M, Ary, indx_first, FieldMask, UseSmlAngle)
SUBROUTINE PackMotionMesh(M, Ary, indx_first, FieldMask, TrimOP)

TYPE(MeshType) , INTENT(IN ) :: M !< Motion mesh
REAL(ReKi) , INTENT(INOUT) :: Ary(:) !< array to pack this mesh into
INTEGER(IntKi) , INTENT(INOUT) :: indx_first !< index into Ary; gives location of next array position to fill
LOGICAL, OPTIONAL , INTENT(IN ) :: FieldMask(FIELDMASK_SIZE) !< flags to determine if this field is part of the packing
LOGICAL, OPTIONAL , INTENT(IN ) :: UseSmlAngle !< flag to determine if the orientation should be packed as a DCM or a log map
LOGICAL, OPTIONAL , INTENT(IN ) :: TrimOP !< flag to determine if the orientation should be packed as a DCM or a log map


! local variables:
INTEGER(IntKi) :: i, j, k
LOGICAL :: Mask(FIELDMASK_SIZE) !< flags to determine if this field is part of the packing
LOGICAL :: OutputSmlAngle
!REAL(R8Ki) :: logmap(3) !< array to pack logmaps into
REAL(R8Ki) :: angles(3) !< array to pack logmaps into
LOGICAL :: PackForTrimSolution
REAL(R8Ki) :: logmap(3) !< array to pack dcm vector representation (logmaps) into
INTEGER(IntKi) :: ErrStat2
CHARACTER(ErrMsgLen) :: ErrMsg2

Expand All @@ -3100,6 +3100,12 @@ SUBROUTINE PackMotionMesh(M, Ary, indx_first, FieldMask, UseSmlAngle)
Mask = .true.
end if

if (present(TrimOP)) then
PackForTrimSolution = TrimOP
else
PackForTrimSolution = .false.
end if


if (Mask(MASKID_TRANSLATIONDISP)) then
do i=1,M%NNodes
Expand All @@ -3111,19 +3117,12 @@ SUBROUTINE PackMotionMesh(M, Ary, indx_first, FieldMask, UseSmlAngle)
end if

if (Mask(MASKID_ORIENTATION)) then
if (present(UseSmlAngle)) then
OutputSmlAngle = UseSmlAngle
else
OutputSmlAngle = .false.
end if

if (OutputSmlAngle) then
if (PackForTrimSolution) then
do i=1,M%NNodes
!call DCM_logMap(M%Orientation(:,:,i), logmap, ErrStat2, ErrMsg2)
angles = GetSmllRotAngs ( M%Orientation(:,:,i), ErrStat2, ErrMsg2 )
call DCM_logMap(M%Orientation(:,:,i), logmap, ErrStat2, ErrMsg2) !NOTE: we cannot use GetSmllRotAngs because we CANNOT assume that all DCMs in the code are small.
do k=1,3
!Ary(indx_first) = logmap(k)
Ary(indx_first) = angles(k)
Ary(indx_first) = logmap(k)
indx_first = indx_first + 1
end do
end do
Expand All @@ -3149,6 +3148,7 @@ SUBROUTINE PackMotionMesh(M, Ary, indx_first, FieldMask, UseSmlAngle)
end if

if (Mask(MASKID_ROTATIONVEL)) then

do i=1,M%NNodes
do j=1,3
Ary(indx_first) = M%RotationVel(j,i)
Expand All @@ -3162,17 +3162,27 @@ SUBROUTINE PackMotionMesh(M, Ary, indx_first, FieldMask, UseSmlAngle)
do j=1,3
Ary(indx_first) = M%TranslationAcc(j,i)
indx_first = indx_first + 1
end do
end do
end do
end if

if (Mask(MASKID_ROTATIONACC)) then
do i=1,M%NNodes
do j=1,3
Ary(indx_first) = M%RotationAcc(j,i)
indx_first = indx_first + 1
end do
end do
if (PackForTrimSolution) then ! these are difficult to converge in a trim solution
do i=1,M%NNodes
do j=1,3
Ary(indx_first) = 0.0_ReKi
indx_first = indx_first + 1
end do
end do
else
do i=1,M%NNodes
do j=1,3
Ary(indx_first) = M%RotationAcc(j,i)
indx_first = indx_first + 1
end do
end do
end if

end if


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 @@ -6725,7 +6725,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, NeedPackedOrient=.true.)
ED%y_interp, ED%m, ErrStat2, ErrMsg2, y_op=y_FAST%Lin%Modules(Module_ED)%Instance(1)%op_y, NeedTrimOP=.true.)
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)

! BeamDyn
Expand All @@ -6737,7 +6737,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, NeedPackedOrient=.true.)
BD%y_interp(k), BD%m(k), ErrStat2, ErrMsg2, y_op=y_FAST%Lin%Modules(Module_BD)%Instance(k)%op_y, NeedTrimOP=.true.)
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
END DO ! k=p_FAST%nBeams

Expand Down Expand Up @@ -6798,7 +6798,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, NeedPackedOrient=.true.)
SD%y_interp, SD%m, ErrStat2, ErrMsg2, y_op=y_FAST%Lin%Modules(Module_SD)%Instance(1)%op_y, NeedTrimOP=.true.)
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
ELSE IF ( p_FAST%CompSub == Module_ExtPtfm ) THEN
END IF ! SubDyn/ExtPtfm_MCKF
Expand Down
Loading

0 comments on commit b023a36

Please sign in to comment.