Various bug fixes for AeroDyn when modeling MHK turbines

modules/aerodyn/src/AeroDyn.f90

@@ -1976,18 +1976,27 @@ subroutine AD_CalcWind(t, u, FLowField, p, m, o, Inflow, ErrStat, ErrMsg)
    ! OLAF points
    if (allocated(o%WakeLocationPoints) .and. allocated(Inflow%InflowWakeVel)) then
       ! If rotor is MHK, add water depth to z coordinate
-      if (p%FVW%MHK > 0) then
+      if (p%FVW%MHK /= MHK_None) then
          PosOffset = [0.0_ReKi, 0.0_ReKi, p%FVW%WtrDpth]
       else
          PosOffset = 0.0_ReKi
       end if
 
-      call IfW_FlowField_GetVelAcc(FlowField, StartNode, t, &
-                                   o%WakeLocationPoints, &
-                                   Inflow%InflowWakeVel, &
-                                   NoAcc, ErrStat2, ErrMsg2, &
-                                   BoxExceedAllow=.true., PosOffset=PosOffset)
-      if(Failed()) return
+      if (p%FVW%MHK /= MHK_None .and. p%CompSeaSt) then ! MHK turbines with waves
+         call WaveField_GetWaveVelAcc_AD(p%WaveField, m%WaveField_m, &
+                                         StartNode, t, &
+                                         o%WakeLocationPoints, &
+                                         Inflow%InflowWakeVel, &
+                                         NoAcc, ErrStat2, ErrMsg2)
+         if(Failed()) return
+      else
+         call IfW_FlowField_GetVelAcc(FlowField, StartNode, t, &
+                                      o%WakeLocationPoints, &
+                                      Inflow%InflowWakeVel, &
+                                      NoAcc, ErrStat2, ErrMsg2, &
+                                      BoxExceedAllow=.true., PosOffset=PosOffset)
+         if(Failed()) return
+      end if
       StartNode = StartNode + size(o%WakeLocationPoints)
    end if
 
@@ -2034,7 +2043,7 @@ subroutine AD_CalcWind_Rotor(t, u, FlowField, p, p_AD, m, RotInflow, StartNode,
          call WaveField_GetWaveVelAcc_AD(p_AD%WaveField, m%WaveField_m, StartNode, t, &
             real(u%HubMotion%TranslationDisp + u%HubMotion%Position, ReKi), &
             RotInflow%InflowOnHub, NoAcc, ErrStat2, ErrMsg2)
-         if(Failed()) return 
+         if(Failed()) return
       else
          RotInflow%InflowOnHub = 0.0_ReKi
       end if
@@ -2085,7 +2094,7 @@ subroutine AD_CalcWind_Rotor(t, u, FlowField, p, p_AD, m, RotInflow, StartNode,
          call IfW_FlowField_GetVelAcc(FlowField, StartNode, t, &
             real(u%HubMotion%TranslationDisp + u%HubMotion%Position, ReKi), &
             RotInflow%InflowOnHub, NoAcc, ErrStat2, ErrMsg2, PosOffset=PosOffset)
-         if(Failed()) return 
+         if(Failed()) return
       else
          RotInflow%InflowOnHub = 0.0_ReKi
       end if
@@ -3208,7 +3217,7 @@ end subroutine SetDisturbedInflow
 !! Loop on blade nodes and computed a weighted sector average inflow at each node
 subroutine SetSectAvgInflow(t, p, p_AD, u, RotInflow, m, errStat, errMsg)
    real(DbKi),                   intent(in   )  :: t                      !< Current simulation time in seconds
-   type(RotParameterType),       intent(in   )  :: p                      !< AD parameters
+   type(RotParameterType),       intent(in   )  :: p                      !< Rotor parameters
    type(AD_ParameterType),       intent(in   )  :: p_AD                   !< AD parameters
    type(RotInputType),           intent(in   )  :: u                      !< AD Inputs at Time
    type(RotInflowType),          intent(in   )  :: RotInflow              !< Rotor inflow at Time
@@ -3227,6 +3236,7 @@ subroutine SetSectAvgInflow(t, p, p_AD, u, RotInflow, m, errStat, errMsg)
    real(ReKi)              :: e_t(3)      !< Polar unit vector perpendicular to rHA_perp ("e_theta")
    real(ReKi)              :: psi         !< Azimuthal offset in the current sector, runs from -psi_bwd to psi_fwd
    real(ReKi)              :: dpsi        !< Azimuthal increment
+   real(ReKi)              :: PosOffset(3)!< IfW position offset for MHK turbines
    real(ReKi), allocatable :: SectPos(:,:)!< Points used to define a given sector (for a given blade node A)
    real(ReKi), allocatable :: SectVel(:,:)!< Inflow velocity at a given sector (Undisturbed and then disturbed)
    real(ReKi), allocatable :: SectAcc(:,:)!< Inflow velocity at a given sector (Undisturbed and then disturbed)
@@ -3293,7 +3303,21 @@ subroutine SetSectAvgInflow(t, p, p_AD, u, RotInflow, m, errStat, errMsg)
 
          ! --- Inflow on sector points
          ! Undisturbed
-         call IfW_FlowField_GetVelAcc(p_AD%FlowField, 1, t, SectPos, SectVel, SectAcc, errStat=errStat2, errMsg=errMsg2); if(Failed()) return
+         if (p%MHK /= MHK_None .and. p_AD%CompSeaSt) then ! MHK turbines with waves
+            call WaveField_GetWaveVelAcc_AD(p_AD%WaveField, m%WaveField_m, 1_IntKi, t, &
+                                            SectPos, SectVel, SectAcc, ErrStat2, ErrMsg2)
+            if(Failed()) return
+         else
+            if (p%MHK /= MHK_None) then
+               PosOffset = [0.0_ReKi, 0.0_ReKi, p%WtrDpth]
+            else
+               PosOffset = 0.0_ReKi
+            end if
+            call IfW_FlowField_GetVelAcc(p_AD%FlowField, 1_IntKi, t, &
+                                         SectPos, SectVel, SectAcc, errStat=errStat2, errMsg=errMsg2, &
+                                         PosOffset=PosOffset); if(Failed()) return
+         end if
+
          ! --- Option 1 Disturbed inflow Before averaging - SectVel is modified in place
          !if (p%TwrPotent /= TwrPotent_none .or. p%TwrShadow /= TwrShadow_none) then
          !   call TwrInflArray(p, u, RotInflow, m, SectPos, SectVel, errStat2, errMsg2); if(Failed()) return

modules/aerodyn/src/AeroDyn_Registry.txt

@@ -517,6 +517,7 @@ typedef	^	RotMiscVarType	ReKi	TFinVrel_i 	3	-	-	"Relative velocity at the refere
 typedef	^	RotMiscVarType	ReKi	TFinSTV_i 	3	-	-	"Structural velocity at the reference point of the fin in the inertial system"
 typedef	^	RotMiscVarType	ReKi	TFinF_i 	3	-	-	"Forces at the reference point of the fin in the inertial system"
 typedef	^	RotMiscVarType	ReKi	TFinM_i 	3	-	-	"Moments at the reference point of the fin in the inertial system"
+typedef	^	RotMiscVarType	GridInterp_MiscVarType	WaveField_m	-	-	-	"misc var information from the GridInterp module"
 
 typedef	^	MiscVarType	RotMiscVarType	rotors	{:}	-	-	"MiscVars for each rotor"	-
 typedef	^	MiscVarType	FVW_InputType	FVW_u	:	-	-	"Inputs to the FVW module"	-
@@ -525,7 +526,7 @@ typedef	^	MiscVarType	FVW_MiscVarType	FVW	-	-	-	"MiscVars from the FVW module"	-
 typedef	^	MiscVarType	ReKi            WindPos	{:}{:} - -	"XYZ coordinates to query for wind velocity/acceleration"	-
 typedef	^	MiscVarType	ReKi            WindVel	{:}{:} - -	"XYZ components of wind velocity"	-
 typedef	^	MiscVarType	ReKi            WindAcc	{:}{:} - -	"XYZ components of wind acceleration"	-
-typedef ^   MiscVarType GridInterp_MiscVarType WaveField_m - - - "misc var information from the SeaState WaveField module"   -
+typedef	^	MiscVarType	GridInterp_MiscVarType	WaveField_m	-	-	-	"misc var information from the GridInterp module"	-
 typedef	^	MiscVarType	AD_InflowType	Inflow	{:}	-	-	"Inflow storage (size of u for history of inputs)"	-
 typedef	^	MiscVarType	AD_InputType	u_perturb	-	-	-	"input perturbation for linearization"	-
 typedef	^	MiscVarType	AD_OutputType	y_lin	-	-	-	"output perturbation for linearization"	-

modules/aerodyn/src/AeroDyn_Types.f90

@@ -557,6 +557,7 @@ MODULE AeroDyn_Types
     REAL(ReKi) , DIMENSION(1:3)  :: TFinSTV_i = 0.0_ReKi      !< Structural velocity at the reference point of the fin in the inertial system [-]
     REAL(ReKi) , DIMENSION(1:3)  :: TFinF_i = 0.0_ReKi      !< Forces at the reference point of the fin in the inertial system [-]
     REAL(ReKi) , DIMENSION(1:3)  :: TFinM_i = 0.0_ReKi      !< Moments at the reference point of the fin in the inertial system [-]
+    TYPE(GridInterp_MiscVarType)  :: WaveField_m      !< misc var information from the GridInterp module [-]
   END TYPE RotMiscVarType
 ! =======================
 ! =========  AD_MiscVarType  =======
@@ -568,7 +569,7 @@ MODULE AeroDyn_Types
     REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE  :: WindPos      !< XYZ coordinates to query for wind velocity/acceleration [-]
     REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE  :: WindVel      !< XYZ components of wind velocity [-]
     REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE  :: WindAcc      !< XYZ components of wind acceleration [-]
-    TYPE(GridInterp_MiscVarType)  :: WaveField_m      !< misc var information from the SeaState WaveField module [-]
+    TYPE(GridInterp_MiscVarType)  :: WaveField_m      !< misc var information from the GridInterp module [-]
     TYPE(AD_InflowType) , DIMENSION(:), ALLOCATABLE  :: Inflow      !< Inflow storage (size of u for history of inputs) [-]
     TYPE(AD_InputType)  :: u_perturb      !< input perturbation for linearization [-]
     TYPE(AD_OutputType)  :: y_lin      !< output perturbation for linearization [-]
@@ -5553,6 +5554,9 @@ subroutine AD_CopyRotMiscVarType(SrcRotMiscVarTypeData, DstRotMiscVarTypeData, C
    DstRotMiscVarTypeData%TFinSTV_i = SrcRotMiscVarTypeData%TFinSTV_i
    DstRotMiscVarTypeData%TFinF_i = SrcRotMiscVarTypeData%TFinF_i
    DstRotMiscVarTypeData%TFinM_i = SrcRotMiscVarTypeData%TFinM_i
+   call GridInterp_CopyMisc(SrcRotMiscVarTypeData%WaveField_m, DstRotMiscVarTypeData%WaveField_m, CtrlCode, ErrStat2, ErrMsg2)
+   call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
+   if (ErrStat >= AbortErrLev) return
 end subroutine
 
 subroutine AD_DestroyRotMiscVarType(RotMiscVarTypeData, ErrStat, ErrMsg)
@@ -5776,6 +5780,8 @@ subroutine AD_DestroyRotMiscVarType(RotMiscVarTypeData, ErrStat, ErrMsg)
    call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
    call NWTC_Library_DestroyMeshMapType(RotMiscVarTypeData%T_P_2_T_L, ErrStat2, ErrMsg2)
    call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
+   call GridInterp_DestroyMisc(RotMiscVarTypeData%WaveField_m, ErrStat2, ErrMsg2)
+   call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
 end subroutine
 
 subroutine AD_PackRotMiscVarType(RF, Indata)
@@ -5916,6 +5922,7 @@ subroutine AD_PackRotMiscVarType(RF, Indata)
    call RegPack(RF, InData%TFinSTV_i)
    call RegPack(RF, InData%TFinF_i)
    call RegPack(RF, InData%TFinM_i)
+   call GridInterp_PackMisc(RF, InData%WaveField_m) 
    if (RegCheckErr(RF, RoutineName)) return
 end subroutine
 
@@ -6083,6 +6090,7 @@ subroutine AD_UnPackRotMiscVarType(RF, OutData)
    call RegUnpack(RF, OutData%TFinSTV_i); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%TFinF_i); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%TFinM_i); if (RegCheckErr(RF, RoutineName)) return
+   call GridInterp_UnpackMisc(RF, OutData%WaveField_m) ! WaveField_m 
 end subroutine
 
 subroutine AD_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg)

modules/seastate/src/SeaSt_WaveField.f90

@@ -684,7 +684,7 @@ end subroutine WaveField_GetWaveKin
 
 
 ! This subroutine is intended for AeroDyn when modeling MHK turbines
-SUBROUTINE WaveField_GetWaveVelAcc_AD( WaveField, WaveField_m, StartNode, Time, pos, FV, FA, ErrStat, ErrMsg )
+SUBROUTINE WaveField_GetWaveVelAcc_AD( WaveField, WaveField_m, StartNode, Time, pos, FV, FA, ErrStat, ErrMsg, BoxExceedAllow )
    type(SeaSt_WaveFieldType),          intent(in   ) :: WaveField
    type(GridInterp_MiscVarType),       intent(inout) :: WaveField_m
    integer(IntKi),                     intent(in   ) :: StartNode
@@ -694,6 +694,7 @@ SUBROUTINE WaveField_GetWaveVelAcc_AD( WaveField, WaveField_m, StartNode, Time,
    real(ReKi),       allocatable,      intent(inout) :: FA(:,:)
    integer(IntKi),                     intent(  out) :: ErrStat  ! Error status of the operation
    character(*),                       intent(  out) :: ErrMsg   ! Error message if errStat /= ErrID_None
+   logical,          optional,         intent(in   ) :: BoxExceedAllow
    integer(IntKi),   allocatable                     :: nodeInWater(:)
    integer(IntKi)                                    :: NumPoints, i
    real(SiKi)                                        :: FV_node(3), FA_node(3)
@@ -734,7 +735,12 @@ SUBROUTINE WaveField_GetWaveVelAcc_AD( WaveField, WaveField_m, StartNode, Time,
       IF (getAcc) THEN
          ALLOCATE(FA_DC( 3, NumPoints ), STAT=ErrStat2); if (FailedMsg('Error allocating FA_DC')) return;
       END IF
-      CALL IfW_FlowField_GetVelAcc(WaveField%CurrField, StartNode, Time, pos, FV_DC, FA_DC, ErrStat2, ErrMsg2, PosOffset=PosOffset); if (Failed()) return;
+
+      IF (PRESENT(BoxExceedAllow)) THEN
+         CALL IfW_FlowField_GetVelAcc(WaveField%CurrField, StartNode, Time, pos, FV_DC, FA_DC, ErrStat2, ErrMsg2, BoxExceedAllow=BoxExceedAllow, PosOffset=PosOffset); if (Failed()) return
+      ELSE
+         CALL IfW_FlowField_GetVelAcc(WaveField%CurrField, StartNode, Time, pos, FV_DC, FA_DC, ErrStat2, ErrMsg2, PosOffset=PosOffset); if (Failed()) return
+      END IF
 
       ! Add contributions from IfW current field if node is in water
       DO i = 1, NumPoints

reg_tests/CTestList.cmake

@@ -455,7 +455,6 @@ ad_regression("ad_MultipleHAWT"             "aerodyn;bem")
 ad_regression("ad_QuadRotor_OLAF"           "aerodyn;bem")
 ad_regression("ad_VerticalAxis_OLAF"        "aerodyn;bem")
 ad_regression("ad_MHK_RM1_Fixed"            "aerodyn;bem;mhk")
-ad_regression("ad_MHK_RM1_Fixed_IfW"        "aerodyn;bem;mhk")
 ad_regression("ad_MHK_RM1_Floating"         "aerodyn;bem;mhk")
 ad_regression("ad_BAR_CombinedCases"        "aerodyn;bem") # NOTE: doing BAR at the end to avoid copy errors
 ad_regression("ad_BAR_OLAF"                 "aerodyn;bem")