Merge pull request OpenFAST#2444 from luwang00/b/ExctnDisp2_Init

HD: Reimplement the initialization of low-pass-filtered potential-flow body positions for ExctnDisp=2
andrew-platt · Oct 3, 2024 · 853ba8c · 853ba8c
2 parents 4da8b4e + be8a671
commit 853ba8c
Show file tree

Hide file tree

Showing 8 changed files with 34 additions and 60 deletions.
diff --git a/modules/hydrodyn/src/HydroDyn.f90 b/modules/hydrodyn/src/HydroDyn.f90
@@ -293,6 +293,7 @@ SUBROUTINE HydroDyn_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, I
             p%vecMultiplier          = InputFileData%vecMultiplier  ! Multiply all vectors and matrices row/column lengths by NBody
             InputFileData%WAMIT%NBodyMod = InputFileData%NBodyMod
             InputFileData%WAMIT%Gravity  = InitInp%Gravity
+            InputFileData%WAMIT%PlatformPos = InitInp%PlatformPos   ! Initial platform/HD origin position
             p%NBody                  = InputFileData%NBody
             p%NBodyMod               = InputFileData%NBodyMod
             call AllocAry( m%F_PtfmAdd, 6*InputFileData%NBody, "m%F_PtfmAdd", ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

diff --git a/modules/hydrodyn/src/HydroDyn.txt b/modules/hydrodyn/src/HydroDyn.txt
@@ -71,13 +71,13 @@ typedef   HydroDyn/HydroDyn            InitInputType                 CHARACTER(1
 typedef   ^                            ^                             LOGICAL                  UseInputFile                    -       .TRUE.        -         "Supplied by Driver:  .TRUE. if using a input file, .FALSE. if all inputs are being passed in by the caller"    -
 typedef   ^                            ^                             FileInfoType             PassedFileData                  -          -          -         "If we don't use the input file, pass everything through this"   -
 typedef   ^                            ^                             CHARACTER(1024)          OutRootName                     -          -          -         "Supplied by Driver:  The name of the root file (without extension) including the full path"    -
-typedef	  ^	                           ^                             Logical	              Linearize                       -       .FALSE.       -         "Flag that tells this module if the glue code wants to linearize."	-
+typedef	  ^	                       ^                             Logical	              Linearize                       -       .FALSE.       -         "Flag that tells this module if the glue code wants to linearize."	-
 typedef   ^                            ^                             ReKi                     Gravity                         -          -          -         "Supplied by Driver:  Gravitational acceleration"  "(m/s^2)"
 typedef   ^                            ^                             DbKi                     TMax                            -          -          -         "Supplied by Driver:  The total simulation time"    "(sec)"
 typedef   ^                            ^                             logical                  VisMeshes                       -       .false.       -         "Output visualization meshes" -
-#
 typedef   ^                            ^                             LOGICAL                  InvalidWithSSExctn              -          -          -         "Whether SeaState configuration is invalid with HydroDyn's state-space excitation (ExctnMod=2)" (-)
 typedef   ^                            ^                             SeaSt_WaveFieldType     *WaveField                       -          -          -         "Pointer to SeaState wave field" -
+typedef   ^                            ^                             ReKi                     PlatformPos                    {6}         -          -         "Initial platform position (6 DOFs)"
 #
 #
 # Define outputs from the initialization routine here:

diff --git a/modules/hydrodyn/src/HydroDyn_DriverCode.f90 b/modules/hydrodyn/src/HydroDyn_DriverCode.f90
@@ -240,23 +240,6 @@ PROGRAM HydroDynDriver
       CALL SetHDInputs(0.0_R8Ki, n, u(1), mappingData, drvrData, ErrStat, ErrMsg);   CALL CheckError()
    END IF
 
-   ! Set the initial low-pass-filtered displacements of potential-flow bodies if ExctnDisp = 2
-   IF ( p%PotMod == 1_IntKi ) THEN
-      IF ( p%WAMIT(1)%ExctnDisp == 2_IntKi ) THEN
-         IF (p%NBodyMod .EQ. 1_IntKi) THEN ! One instance of WAMIT with NBody
-            DO i = 1,p%NBody
-               xd%WAMIT(1)%BdyPosFilt(1,i,:) = u(1)%WAMITMesh%TranslationDisp(1,i)
-               xd%WAMIT(1)%BdyPosFilt(2,i,:) = u(1)%WAMITMesh%TranslationDisp(2,i)
-            END DO
-         ELSE IF (p%NBodyMod > 1_IntKi) THEN ! NBody instances of WAMIT with one body each
-            DO i = 1,p%NBody
-               xd%WAMIT(i)%BdyPosFilt(1,1,:) = u(1)%WAMITMesh%TranslationDisp(1,i)
-               xd%WAMIT(i)%BdyPosFilt(2,1,:) = u(1)%WAMITMesh%TranslationDisp(2,i)
-            END DO
-         END IF
-      END IF
-   END IF
-
    !...............................................................................................................................
    ! --- Linearization
    !...............................................................................................................................
@@ -348,6 +331,12 @@ subroutine SetHD_InitInputs()
 
    InitInData_HD%WaveField => InitOutData_SeaSt%WaveField
 
+   IF (( drvrData%PRPInputsMod /= 2 ) .AND. ( drvrData%PRPInputsMod >= 0 )) THEN
+      InitInData_HD%PlatformPos  = drvrData%uPRPInSteady
+   ELSE
+      InitInData_HD%PlatformPos  = drvrData%PRPin(1,1:6)
+   END IF
+
 end subroutine SetHD_InitInputs
 !----------------------------------------------------------------------------------------------------------------------------------
 subroutine CheckError()

diff --git a/modules/hydrodyn/src/HydroDyn_Types.f90 b/modules/hydrodyn/src/HydroDyn_Types.f90
@@ -94,6 +94,7 @@ MODULE HydroDyn_Types
     LOGICAL  :: VisMeshes = .false.      !< Output visualization meshes [-]
     LOGICAL  :: InvalidWithSSExctn = .false.      !< Whether SeaState configuration is invalid with HydroDyn's state-space excitation (ExctnMod=2) [(-)]
     TYPE(SeaSt_WaveFieldType) , POINTER :: WaveField => NULL()      !< Pointer to SeaState wave field [-]
+    REAL(ReKi) , DIMENSION(1:6)  :: PlatformPos = 0.0_ReKi      !< Initial platform position (6 DOFs) [-]
   END TYPE HydroDyn_InitInputType
 ! =======================
 ! =========  HydroDyn_InitOutputType  =======
@@ -594,7 +595,7 @@ subroutine HydroDyn_CopyInitInput(SrcInitInputData, DstInitInputData, CtrlCode,
    integer(IntKi),  intent(in   ) :: CtrlCode
    integer(IntKi),  intent(  out) :: ErrStat
    character(*),    intent(  out) :: ErrMsg
-   integer(B8Ki)                  :: LB(0), UB(0)
+   integer(B8Ki)                  :: LB(1), UB(1)
    integer(IntKi)                 :: ErrStat2
    character(ErrMsgLen)           :: ErrMsg2
    character(*), parameter        :: RoutineName = 'HydroDyn_CopyInitInput'
@@ -612,6 +613,7 @@ subroutine HydroDyn_CopyInitInput(SrcInitInputData, DstInitInputData, CtrlCode,
    DstInitInputData%VisMeshes = SrcInitInputData%VisMeshes
    DstInitInputData%InvalidWithSSExctn = SrcInitInputData%InvalidWithSSExctn
    DstInitInputData%WaveField => SrcInitInputData%WaveField
+   DstInitInputData%PlatformPos = SrcInitInputData%PlatformPos
 end subroutine
 
 subroutine HydroDyn_DestroyInitInput(InitInputData, ErrStat, ErrMsg)
@@ -650,14 +652,15 @@ subroutine HydroDyn_PackInitInput(RF, Indata)
          call SeaSt_WaveField_PackSeaSt_WaveFieldType(RF, InData%WaveField) 
       end if
    end if
+   call RegPack(RF, InData%PlatformPos)
    if (RegCheckErr(RF, RoutineName)) return
 end subroutine
 
 subroutine HydroDyn_UnPackInitInput(RF, OutData)
    type(RegFile), intent(inout)    :: RF
    type(HydroDyn_InitInputType), intent(inout) :: OutData
    character(*), parameter            :: RoutineName = 'HydroDyn_UnPackInitInput'
-   integer(B8Ki)   :: LB(0), UB(0)
+   integer(B8Ki)   :: LB(1), UB(1)
    integer(IntKi)  :: stat
    logical         :: IsAllocAssoc
    integer(B8Ki)   :: PtrIdx
@@ -690,6 +693,7 @@ subroutine HydroDyn_UnPackInitInput(RF, OutData)
    else
       OutData%WaveField => null()
    end if
+   call RegUnpack(RF, OutData%PlatformPos); if (RegCheckErr(RF, RoutineName)) return
 end subroutine
 
 subroutine HydroDyn_CopyInitOutput(SrcInitOutputData, DstInitOutputData, CtrlCode, ErrStat, ErrMsg)

diff --git a/modules/hydrodyn/src/WAMIT.f90 b/modules/hydrodyn/src/WAMIT.f90
@@ -164,7 +164,7 @@ SUBROUTINE WAMIT_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, ErrS
       REAL(ReKi), ALLOCATABLE                :: WAMITPer  (:)                        ! Period         components as ordered in the WAMIT output files (sec    )
       REAL(ReKi), ALLOCATABLE                :: WAMITWvDir(:)                        ! Wave direction components as ordered in the WAMIT output files (degrees)
 
-      INTEGER                                :: I,iGrid,iX,iY,iHdg,iBdy              ! Generic index
+      INTEGER                                :: I,iGrid,iX,iY,iHdg,iBdy,iStp         ! Generic index
       INTEGER                                :: InsertInd                            ! The lowest sorted index whose associated frequency component is higher than the current frequency component -- this is to sort the frequency components from lowest to highest
       INTEGER                                :: J                                    ! Generic index
       INTEGER                                :: K                                    ! Generic index
@@ -190,6 +190,7 @@ SUBROUTINE WAMIT_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, ErrS
       TYPE(FFT_DataType)                     :: FFT_Data                             ! the instance of the FFT module we're using
       integer(IntKi)                         :: iSub, jSub                           ! indices into the 6x6 sub-matrices used to redimensionalize the WAMIT data (Needed because NBodyMod=1 could have WAMIT matrices which are 6N x 6N)
       integer(IntKi)                         :: iBody                                ! WAMIT body index
+      real(ReKi)                             :: BdyPos0(3)                           ! Initial translational displacement of the WAMIT body
       real(R8Ki)                             :: orientation(3,3)                     ! Initial orientation of the WAMIT body 
       real(R8Ki)                             :: theta(3)                             ! Euler angle rotations of the WAMIT body
       real(ReKi)                             :: WaveNmbr                             ! Frequency-dependent wave number
@@ -1371,15 +1372,24 @@ SUBROUTINE WAMIT_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, ErrS
                   end if   
             end if
 
-            IF ( (p%ExctnMod>0) .AND. (p%ExctnDisp==2) ) THEN ! Allocate array for filtered potential-flow body positions
+            IF ( (p%ExctnMod>0) .AND. (p%ExctnDisp==2) ) THEN ! Allocate and initialize array for filtered potential-flow body positions
                p%ExctnFiltConst = exp(-2.0*Pi*p%ExctnCutOff * Interval)
                ALLOCATE ( xd%BdyPosFilt(1:2, 1:p%NBody, 1:3) , STAT=ErrStat2 )
                IF ( ErrStat2 /= 0 )  THEN
                   CALL SetErrStat( ErrID_Fatal, 'Error allocating memory for the BdyPosFilt array.', ErrStat, ErrMsg, RoutineName)
                   CALL Cleanup()
                   RETURN            
                END IF
-               xd%BdyPosFilt = 0.0_ReKi
+               orientation  = EulerConstructZYX(InitInp%PlatformPos(4:6));
+               DO iBdy = 1,p%NBody
+                  ! Initial WAMIT body position
+                  BdyPos0 = InitInp%PlatformPos(1:3) &
+                            + matmul((/InitInp%PtfmRefxt(iBdy),InitInp%PtfmRefyt(iBdy),InitInp%PtfmRefzt(iBdy)/),orientation) &
+                            - (/InitInp%PtfmRefxt(iBdy),InitInp%PtfmRefyt(iBdy),InitInp%PtfmRefzt(iBdy)/)
+                  DO iStp = 1,3
+                     xd%BdyPosFilt(1:2,iBdy,iStp) = BdyPos0(1:2)
+                  END DO
+               END DO
             END IF
 
          ENDSELECT   

diff --git a/modules/hydrodyn/src/WAMIT.txt b/modules/hydrodyn/src/WAMIT.txt
@@ -41,6 +41,7 @@ typedef   ^                            ^                             Conv_Rdtn_I
 typedef   ^                            ^                             SeaSt_WaveFieldType     *WaveField                       -          -        -         "Pointer to wave field"
 typedef   ^                            ^                             INTEGER                  PtfmYMod                        -          -        -         "Large yaw model" -
 typedef   ^                            ^                             ReKi                     PtfmRefY                        -          -        -         "Initial reference yaw offset" (rad)
+typedef   ^                            ^                             ReKi                     PlatformPos                     {6}        -        -         "Initial platform position (6 DOFs)"
 #
 #
 # Define outputs from the initialization routine here:

diff --git a/modules/hydrodyn/src/WAMIT_Types.f90 b/modules/hydrodyn/src/WAMIT_Types.f90
@@ -61,6 +61,7 @@ MODULE WAMIT_Types
     TYPE(SeaSt_WaveFieldType) , POINTER :: WaveField => NULL()      !< Pointer to wave field [-]
     INTEGER(IntKi)  :: PtfmYMod = 0_IntKi      !< Large yaw model [-]
     REAL(ReKi)  :: PtfmRefY = 0.0_ReKi      !< Initial reference yaw offset [(rad)]
+    REAL(ReKi) , DIMENSION(1:6)  :: PlatformPos = 0.0_ReKi      !< Initial platform position (6 DOFs) [-]
   END TYPE WAMIT_InitInputType
 ! =======================
 ! =========  WAMIT_ContinuousStateType  =======
@@ -262,6 +263,7 @@ subroutine WAMIT_CopyInitInput(SrcInitInputData, DstInitInputData, CtrlCode, Err
    DstInitInputData%WaveField => SrcInitInputData%WaveField
    DstInitInputData%PtfmYMod = SrcInitInputData%PtfmYMod
    DstInitInputData%PtfmRefY = SrcInitInputData%PtfmRefY
+   DstInitInputData%PlatformPos = SrcInitInputData%PlatformPos
 end subroutine
 
 subroutine WAMIT_DestroyInitInput(InitInputData, ErrStat, ErrMsg)
@@ -334,6 +336,7 @@ subroutine WAMIT_PackInitInput(RF, Indata)
    end if
    call RegPack(RF, InData%PtfmYMod)
    call RegPack(RF, InData%PtfmRefY)
+   call RegPack(RF, InData%PlatformPos)
    if (RegCheckErr(RF, RoutineName)) return
 end subroutine
 
@@ -387,6 +390,7 @@ subroutine WAMIT_UnPackInitInput(RF, OutData)
    end if
    call RegUnpack(RF, OutData%PtfmYMod); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%PtfmRefY); if (RegCheckErr(RF, RoutineName)) return
+   call RegUnpack(RF, OutData%PlatformPos); if (RegCheckErr(RF, RoutineName)) return
 end subroutine
 
 subroutine WAMIT_CopyContState(SrcContStateData, DstContStateData, CtrlCode, ErrStat, ErrMsg)

diff --git a/modules/openfast-library/src/FAST_Subs.f90 b/modules/openfast-library/src/FAST_Subs.f90
@@ -931,6 +931,7 @@ SUBROUTINE FAST_InitializeAll( t_initial, p_FAST, y_FAST, m_FAST, ED, SED, BD, S
       Init%InData_HD%OutRootName   = TRIM(p_FAST%OutFileRoot)//'.'//TRIM(y_FAST%Module_Abrev(Module_HD))
       Init%InData_HD%TMax          = p_FAST%TMax
       Init%InData_HD%Linearize     = p_FAST%Linearize
+      Init%InData_HD%PlatformPos   = Init%OutData_ED%PlatformPos ! Initial platform position; PlatformPos(1:3) is effectively the initial position of the HD origin
       if (p_FAST%WrVTK /= VTK_None) Init%InData_HD%VisMeshes=.true.
 
       ! if ( p_FAST%CompSeaSt == Module_SeaSt ) then  ! this is always true
@@ -1664,42 +1665,6 @@ SUBROUTINE FAST_InitializeAll( t_initial, p_FAST, y_FAST, m_FAST, ED, SED, BD, S
          ErrMsg = ""
       END IF
 
-   ! ----------------------------------------------------------------------------
-   ! Initialize low-pass-filtered displacements of HydroDyn potential-flow bodies
-   ! ----------------------------------------------------------------------------
-   IF ( (p_FAST%CompHydro == Module_HD) .AND. (HD%p%PotMod == 1_IntKi) ) THEN
-      IF ( HD%p%WAMIT(1)%ExctnDisp == 2_IntKi ) THEN
-         ! Set the initial displacement of ED%PlatformPtMesh here to use MeshMapping
-         ED%y%PlatformPtMesh%TranslationDisp(:,1) = Init%OutData_ED%PlatformPos(1:3)
-         CALL SmllRotTrans( 'initial platform rotation ', &
-                             REAL(Init%OutData_ED%PlatformPos(4),R8Ki), &
-                             REAL(Init%OutData_ED%PlatformPos(5),R8Ki), &
-                             REAL(Init%OutData_ED%PlatformPos(6),R8Ki), &
-                             ED%y%PlatformPtMesh%Orientation(:,:,1), '', ErrStat2, ErrMsg2 )
-         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
-         ED%y%PlatformPtMesh%TranslationDisp(1,1) = ED%y%PlatformPtMesh%TranslationDisp(1,1) + ED%y%PlatformPtMesh%Orientation(3,1,1) * ED%p%PtfmRefzt
-         ED%y%PlatformPtMesh%TranslationDisp(2,1) = ED%y%PlatformPtMesh%TranslationDisp(2,1) + ED%y%PlatformPtMesh%Orientation(3,2,1) * ED%p%PtfmRefzt
-         ED%y%PlatformPtMesh%TranslationDisp(3,1) = ED%y%PlatformPtMesh%TranslationDisp(3,1) + ED%y%PlatformPtMesh%Orientation(3,3,1) * ED%p%PtfmRefzt - ED%p%PtfmRefzt
-         CALL Transfer_PlatformMotion_to_HD( ED%y%PlatformPtMesh, HD%Input(1), MeshMapData, ErrStat2, ErrMsg2 )
-         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
-         IF (ErrStat >= AbortErrLev) THEN
-            CALL Cleanup()
-            RETURN
-         END IF
-         IF (HD%p%NBodyMod .EQ. 1_IntKi) THEN ! One instance of WAMIT with NBody
-            DO i = 1,HD%p%NBody
-               HD%xd(STATE_CURR)%WAMIT(1)%BdyPosFilt(1,i,:) = HD%Input(1)%WAMITMesh%TranslationDisp(1,i)
-               HD%xd(STATE_CURR)%WAMIT(1)%BdyPosFilt(2,i,:) = HD%Input(1)%WAMITMesh%TranslationDisp(2,i)
-            END DO
-         ELSE IF (HD%p%NBodyMod > 1_IntKi) THEN ! NBody instances of WAMIT with one body each
-            DO i = 1,HD%p%NBody
-               HD%xd(STATE_CURR)%WAMIT(i)%BdyPosFilt(1,1,:) = HD%Input(1)%WAMITMesh%TranslationDisp(1,i)
-               HD%xd(STATE_CURR)%WAMIT(i)%BdyPosFilt(2,1,:) = HD%Input(1)%WAMITMesh%TranslationDisp(2,i)
-            END DO
-         END IF
-      END IF
-   END IF
-
    ! -------------------------------------------------------------------------
    ! Initialize for linearization or computing aero maps:
    ! -------------------------------------------------------------------------