Merge pull request #3306 from hkross/bug/MHK_Inflow_Calls

MHK inflow calls

Author: andrew-platt

docs/source/user/inflowwind/input.rst

@@ -117,26 +117,26 @@ be performed in a different order than if both angles are specified in the same
 
 .. _inflow_superposition:
 
-Superposition of Wave and Current Inflow
-========================================
-For MHK turbines, wave and current velocities and accelerations are superimposed (i.e., summed) such that all submerged components are exposed
+Superposition/Coupling of Wave and Current Inflow
+=================================================
+For MHK turbines, wave and current velocities and accelerations are superimposed (i.e., summed) or coupled (see :ref:`sea-waves`) such that all submerged components are exposed
 to the same inflow field. Both AeroDyn and HydroDyn can query SeaState for wave field information. SeaState then queries InflowWind for the current
-field, sums the velocities and accelerations, and returns the superimposed flow field information. This has several implications for modeling
+field, sums or couples the velocities and accelerations, and returns the flow field information. This has several implications for modeling
 MHK turbines, which are listed below. Note that dynamic pressure contributions from InflowWind are neglected.
 
 When modeling a rotor or rotor/tower only (i.e., hydrodynamics modeled in AeroDyn only):
 
 - SeaState must be used when defining a flow field with waves
-- Current definition in SeaState must always be set to 0
-- If SeaState is activated, InflowWind must also be activated, though the current can be set to 0
 - InflowWind must be used when defining a flow field with currents
+- Current definition in SeaState must always be set to 0
 - For combined wave and current flow fields, SeaState will query InflowWind
 
 When modeling a rotor or rotor/tower and support structure (i.e., hydrodynamics modeled in AeroDyn and HydroDyn):
 
 - SeaState must always be used, even when defining a flow field with no waves
-- Current definition in SeaState must always be set to 0
-- If SeaState is activated, InflowWind must also be activated, though the current can be set to 0
 - InflowWind must be used when defining a flow field with currents
+- Current definition in SeaState must always be set to 0
 - For current only cases, set the SeaState wave field to 0; current information will be passed through SeaState from InflowWind
-- For combined wave and current flow fields, SeaState will query InflowWind
+- For combined wave and current flow fields, SeaState will query InflowWind
+
+Wave and current coupling is only possible when running the OpenFAST glue code or the AeroDyn driver. This feature is not supported by the HydroDyn or MoorDyn drivers.

docs/source/user/seastate/input_files.rst

@@ -107,6 +107,8 @@ When setting up the wave grid, it is necessary to make sure the wave grid is lar
 
 OpenFAST precomputes and saves the wave-field velocity, acceleration, dynamic pressure, and wave elevation at the start of the simulation. Generating and maintaining the wave grid can be memory intensive for long simulations. Users should set the wave grid to be no larger or finer than necessary to reduce memory use. Reducing **WaveTMax** or increasing **WaveDT** (see WAVES section below) also reduces memory use. For long crested waves (no directional spreading) aligned with the *X*-direction (or *Y*-direction), **NY** (or **NX**) can be reduced to the minimum allowed value of 2 to save memory.
 
+.. _sea-waves:
+
 Waves
 -----
 
@@ -211,7 +213,8 @@ time-averaged current velocity at the still water level. For applicable **WindTy
 InflowWind, users should ensure that the flow-field grid from InflowWind reaches the still 
 water level. **WvCrntMod** has no effect when **WaveMod** = 0 or 6, or when there is no 
 current from either SeaState (**CurrMod** = 0) or InflowWind if simulating marine 
-hydrokinetic turbines.
+hydrokinetic turbines. See :ref:`inflow_superposition`
+for additional context around wave-current coupling when simulating MHK turbines.
 
 **WaveTMax** sets the length of the incident wave kinematics time
 series, but it also determines the frequency step used in the inverse

modules/aerodyn/src/AeroDyn.f90

@@ -1412,7 +1412,7 @@ subroutine Init_RotInflow( p, RotInflow, errStat, ErrMsg )
       if (Failed()) return
       RotInflow%Blade(k)%InflowVel = 0.0_ReKi
 
-      if (p%MHK > 0) then
+      if (p%MHK /= MHK_None) then
          call AllocAry( RotInflow%Blade(k)%InflowAcc, 3_IntKi, p%NumBlNds, 'RotInflow%Blade(k)%InflowAcc', ErrStat2, ErrMsg2 )
          if (Failed()) return
          RotInflow%Blade(k)%InflowAcc = 0.0_ReKi
@@ -1422,7 +1422,7 @@ subroutine Init_RotInflow( p, RotInflow, errStat, ErrMsg )
    call AllocAry( RotInflow%Tower%InflowVel, 3_IntKi, p%NumTwrNds, 'RotInflow%Tower%InflowVel', ErrStat2, ErrMsg2 ) ! could be size zero
    if (Failed()) return
 
-   if (p%MHK > 0) then
+   if (p%MHK /= MHK_None) then
       call AllocAry( RotInflow%Tower%InflowAcc, 3_IntKi, p%NumTwrNds, 'RotInflow%Tower%InflowAcc', ErrStat2, ErrMsg2 ) ! could be size zero
       if (Failed()) return
    end if
@@ -2031,7 +2031,7 @@ subroutine AD_CalcWind_Rotor(t, u, FlowField, p, p_AD, m, RotInflow, StartNode,
    if (.not. associated(FlowField)) return  ! use the initial (or input) values for these inputs
 
    ! If rotor is MHK, add water depth to z coordinate
-   if (p%MHK > 0) then
+   if (p%MHK /= MHK_None) then
       PosOffset = [0.0_ReKi, 0.0_ReKi, p%WtrDpth]
    else
       PosOffset = 0.0_ReKi
@@ -4541,7 +4541,7 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, calcCrvAngle, ErrSt
       ! .............................
       ! check tower mesh data:
       ! .............................
-   if (InputFileData%TwrPotent /= TwrPotent_none .or. InputFileData%TwrShadow /= TwrShadow_none .or. InputFileData%TwrAero /= TwrAero_none .or. InitInp%MHK > 0) then
+   if (InputFileData%TwrPotent /= TwrPotent_none .or. InputFileData%TwrShadow /= TwrShadow_none .or. InputFileData%TwrAero /= TwrAero_none .or. InitInp%MHK /= MHK_None) then
       do iR = 1,size(NumBl)
          if (InputFileData%rotors(iR)%NumTwrNds <= 0) cycle !bjj: this could be removed since the loops here already take into account the number of tower nodes
 

modules/aerodyn/src/AeroDyn_Driver_Subs.f90

@@ -1411,7 +1411,7 @@ subroutine ValidateInputs(dvr, errStat, errMsg)
 
    if (dvr%MHK /= MHK_None .and. dvr%MHK /= MHK_FixedBottom .and. dvr%MHK /= MHK_Floating) call SetErrStat(ErrID_Fatal, 'MHK switch must be 0, 1, or 2.', ErrStat, ErrMsg, RoutineName)
 
-   if (dvr%MHK /= MHK_None .and. dvr%SS_InitInp%CompSeaSt == 1 .and. dvr%IW_InitInp%CompInflow /= 1) call SetErrStat( ErrID_Fatal, 'InflowWind must be activated for MHK turbines when SeaState is used.', ErrStat, ErrMsg, RoutineName )
+   if (dvr%MHK /= MHK_None .and. dvr%SS_InitInp%CompSeaSt == 1 .and. dvr%IW_InitInp%CompInflow == 0 .and. dvr%IW_InitInp%HWindSpeed > 0) call SetErrStat( ErrID_Fatal, 'Steady Wind option in AeroDyn driver cannot be used for MHK turbines with SeaState.', ErrStat, ErrMsg, RoutineName )
 
    if (dvr%MHK == MHK_None .and. dvr%SS_InitInp%CompSeaSt /= 0) call SetErrStat( ErrID_Fatal, 'SeaState cannot be used with wind turbines.', ErrStat, ErrMsg, RoutineName )
 

modules/seastate/src/SeaSt_WaveField.f90

@@ -168,6 +168,16 @@ SUBROUTINE WaveField_GetNodeWaveKin( WaveField, WaveField_m, Time, pos, forceNod
    WaveElev = WaveField_GetNodeTotalWaveElev(WaveField, WaveField_m, Time, pos, ErrStat2, ErrMsg2, Elev1=WaveElev1, Elev2=WaveElev2)
    if (Failed()) return
 
+   ! Check if point is below the seabed
+   if (pos(3)<-WaveField%EffWtrDpth) then
+      nodeInWater = 1_IntKi  ! Prevent problems with HydroDyn logic
+      FV(:)       = 0.0_SiKi
+      FA(:)       = 0.0_SiKi
+      FDynP       = 0.0_SiKi
+      FAMCF(:)    = 0.0_SiKi
+      return
+   end if
+
    IF (WaveField%WaveStMod == 0) THEN ! No wave stretching
 
       IF ( pos(3) <= 0.0_ReKi) THEN ! Node is at or below the SWL
@@ -182,10 +192,10 @@ SUBROUTINE WaveField_GetNodeWaveKin( WaveField, WaveField_m, Time, pos, forceNod
          END IF
       ELSE ! Node is above the SWL
          nodeInWater = 0_IntKi
-         FV(:)       = 0.0
-         FA(:)       = 0.0
-         FDynP       = 0.0
-         FAMCF(:)    = 0.0
+         FV(:)       = 0.0_SiKi
+         FA(:)       = 0.0_SiKi
+         FDynP       = 0.0_SiKi
+         FAMCF(:)    = 0.0_SiKi
       END IF
 
    ELSE ! Wave stretching enabled
@@ -251,10 +261,10 @@ SUBROUTINE WaveField_GetNodeWaveKin( WaveField, WaveField_m, Time, pos, forceNod
       ELSE ! Node is out of water - zero-out all wave dynamics
 
          nodeInWater = 0_IntKi
-         FV(:)       = 0.0
-         FA(:)       = 0.0
-         FDynP       = 0.0
-         FAMCF(:)    = 0.0
+         FV(:)       = 0.0_SiKi
+         FA(:)       = 0.0_SiKi
+         FDynP       = 0.0_SiKi
+         FAMCF(:)    = 0.0_SiKi
 
       END IF ! If node is in or out of water
 
@@ -315,6 +325,13 @@ SUBROUTINE WaveField_GetDynP( WaveField, WaveField_m, Time, pos, forceNodeInWate
    ! Wave elevation (Calls WaveField_Interp_Setup3D internally so WaveField_Interp_3D_vec can be used below)
    WaveElev  = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, pos, ErrStat2, ErrMsg2 ); if (Failed()) return;
 
+   ! Check if point is below the seabed
+   if (pos(3)<-WaveField%EffWtrDpth) then
+      nodeInWater = 1_IntKi  ! Prevent problems with HydroDyn logic
+      FDynP       = 0.0_SiKi
+      return
+   end if
+
    IF (WaveField%WaveStMod == 0) THEN ! No wave stretching
 
       IF ( pos(3) <= 0.0_ReKi) THEN ! Node is at or below the SWL
@@ -324,7 +341,7 @@ SUBROUTINE WaveField_GetDynP( WaveField, WaveField_m, Time, pos, forceNodeInWate
          FDynP = GridInterp4D   ( WaveField%WaveDynP, WaveField_m )
       ELSE ! Node is above the SWL
          nodeInWater = 0_IntKi
-         FDynP       = 0.0
+         FDynP       = 0.0_SiKi
       END IF
 
    ELSE ! Wave stretching enabled
@@ -365,7 +382,7 @@ SUBROUTINE WaveField_GetDynP( WaveField, WaveField_m, Time, pos, forceNodeInWate
       ELSE ! Node is out of water - zero-out all wave dynamics
 
          nodeInWater = 0_IntKi
-         FDynP       = 0.0
+         FDynP       = 0.0_SiKi
 
       END IF ! If node is in or out of water
    END IF ! If wave stretching is on or off
@@ -407,6 +424,13 @@ SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNod
    ! Wave elevation (Calls WaveField_Interp_Setup3D internally so WaveField_Interp_3D_vec can be used below)
    WaveElev  = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, pos, ErrStat2, ErrMsg2 ); if (Failed()) return;
 
+   ! Check if point is below the seabed
+   if (pos(3)<-WaveField%EffWtrDpth) then
+      nodeInWater = 1_IntKi  ! Prevent problems with HydroDyn logic
+      FV(:)       = 0.0_SiKi
+      return
+   end if
+
    IF (WaveField%WaveStMod == 0) THEN ! No wave stretching
 
       IF ( pos(3) <= 0.0_ReKi) THEN ! Node is at or below the SWL
@@ -416,7 +440,7 @@ SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNod
          FV(:) = GridInterp4DVec( WaveField%WaveVel,  WaveField_m )
       ELSE ! Node is above the SWL
          nodeInWater = 0_IntKi
-         FV(:)       = 0.0
+         FV(:)       = 0.0_SiKi
       END IF
 
    ELSE ! Wave stretching enabled
@@ -463,7 +487,7 @@ SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNod
       ELSE ! Node is out of water - zero-out all wave dynamics
 
          nodeInWater = 0_IntKi
-         FV(:)       = 0.0
+         FV(:)       = 0.0_SiKi
 
       END IF ! If node is in or out of water
 
@@ -520,7 +544,15 @@ SUBROUTINE WaveField_GetNodeWaveVelAcc( WaveField, WaveField_m, Time, pos, force
 
    ! Wave elevation
    WaveElev  = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, pos, ErrStat2, ErrMsg2 ); if (Failed()) return;
-    
+
+   ! Check if point is below the seabed
+   if (pos(3)<-WaveField%EffWtrDpth) then
+      nodeInWater = 1_IntKi  ! Prevent problems with HydroDyn logic
+      FV(:)       = 0.0_SiKi
+      FA(:)       = 0.0_SiKi
+      return
+   end if
+
    IF (WaveField%WaveStMod == 0) THEN ! No wave stretching
 
       IF ( pos(3) <= 0.0_ReKi) THEN ! Node is at or below the SWL
@@ -531,8 +563,8 @@ SUBROUTINE WaveField_GetNodeWaveVelAcc( WaveField, WaveField_m, Time, pos, force
          FA(:) = GridInterp4DVec( WaveField%WaveAcc,  WaveField_m )
       ELSE ! Node is above the SWL
          nodeInWater = 0_IntKi
-         FV(:)       = 0.0
-         FA(:)       = 0.0
+         FV(:)       = 0.0_SiKi
+         FA(:)       = 0.0_SiKi
       END IF
 
    ELSE ! Wave stretching enabled
@@ -582,8 +614,8 @@ SUBROUTINE WaveField_GetNodeWaveVelAcc( WaveField, WaveField_m, Time, pos, force
       ELSE ! Node is out of water - zero-out all wave dynamics
 
          nodeInWater = 0_IntKi
-         FV(:)       = 0.0
-         FA(:)       = 0.0
+         FV(:)       = 0.0_SiKi
+         FA(:)       = 0.0_SiKi
       END IF ! If node is in or out of water
 
    END IF ! If wave stretching is on or off

modules/seastate/src/SeaState_Input.f90

@@ -925,6 +925,11 @@ subroutine SeaStateInput_ProcessInitData( InitInp, p, InputFileData, ErrStat, Er
 
 
       ! CurrMod - Current profile model switch
+   if ( ( InputFileData%Current%CurrMod /= 0 ) .AND. ( InitInp%MHK /= MHK_None ) ) then
+      call SetErrStat( ErrID_Fatal,'CurrMod must be set to 0 for an MHK turbine.',ErrStat,ErrMsg,RoutineName)
+      return
+   end if
+   
    if ( InitInp%hasCurrField ) then
       call SetErrStat( ErrID_Warn,'Expecting current field from InflowWind. Setting CurrMod to 0.',ErrStat,ErrMsg,RoutineName)
       InputFileData%Current%CurrMod = 0
@@ -940,10 +945,6 @@ subroutine SeaStateInput_ProcessInitData( InitInp, p, InputFileData, ErrStat, Er
       return
    end if
 
-   ! if ( ( InputFileData%Current%CurrMod /= 0 ) .AND. ( InitInp%MHK /= MHK_None ) ) then
-   !    call SetErrStat( ErrID_Fatal,'CurrMod must be set to 0 for an MHK turbine.',ErrStat,ErrMsg,RoutineName)
-   !    return
-   ! end if
 
 
    if ( InputFileData%Current%CurrMod == 1 )  then  ! .TRUE if we have standard current.