diff --git a/modules/beamdyn/src/BeamDyn.f90 b/modules/beamdyn/src/BeamDyn.f90 index ec6ab10191..4eaddf4b38 100644 --- a/modules/beamdyn/src/BeamDyn.f90 +++ b/modules/beamdyn/src/BeamDyn.f90 @@ -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) @@ -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 @@ -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 @@ -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) diff --git a/modules/elastodyn/src/ElastoDyn.f90 b/modules/elastodyn/src/ElastoDyn.f90 index 07e84f7fea..ffb058a9fb 100644 --- a/modules/elastodyn/src/ElastoDyn.f90 +++ b/modules/elastodyn/src/ElastoDyn.f90 @@ -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) @@ -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 @@ -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 @@ -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. @@ -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 diff --git a/modules/openfast-library/src/FAST_Lin.f90 b/modules/openfast-library/src/FAST_Lin.f90 index e50931cda9..63b06f2344 100644 --- a/modules/openfast-library/src/FAST_Lin.f90 +++ b/modules/openfast-library/src/FAST_Lin.f90 @@ -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 @@ -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 @@ -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 diff --git a/modules/subdyn/src/SubDyn.f90 b/modules/subdyn/src/SubDyn.f90 index f5ddeb4e12..248e3e7478 100644 --- a/modules/subdyn/src/SubDyn.f90 +++ b/modules/subdyn/src/SubDyn.f90 @@ -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 @@ -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 @@ -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. @@ -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