diff --git a/modules/openfast-library/src/FAST_Registry.txt b/modules/openfast-library/src/FAST_Registry.txt
index e11ad319c4..08822a8140 100644
--- a/modules/openfast-library/src/FAST_Registry.txt
+++ b/modules/openfast-library/src/FAST_Registry.txt
@@ -76,8 +76,9 @@ typedef ^ FAST_VTK_SurfaceType SiKi GroundRad - - - "radius for plotting circle
typedef ^ FAST_VTK_SurfaceType SiKi NacelleBox {3}{8} - - "X-Y-Z locations of 8 points that define the nacelle box, relative to the nacelle position" m
typedef ^ FAST_VTK_SurfaceType SiKi TowerRad {:} - - "radius of each ED tower node" m
typedef ^ FAST_VTK_SurfaceType IntKi NWaveElevPts {2} - - "number of points for wave elevation visualization" -
-typedef ^ FAST_VTK_SurfaceType SiKi WaveElevXY {:}{:} - - "X-Y locations for WaveElev output (for visualization). First dimension is the X (1) and Y (2) coordinate. Second dimension is the point number." "m,-"
-typedef ^ FAST_VTK_SurfaceType SiKi WaveElev {:}{:} - - "wave elevation at WaveElevXY; first dimension is time step; second dimension is point number" "m,-"
+typedef ^ FAST_VTK_SurfaceType SiKi WaveElevVisX {:} - - "X locations for WaveElev output (for visualization)." "m,-"
+typedef ^ FAST_VTK_SurfaceType SiKi WaveElevVisY {:} - - "Y locations for WaveElev output (for visualization)." "m,-"
+typedef ^ FAST_VTK_SurfaceType SiKi WaveElevVisGrid {:}{:}{:} - - "wave elevation at WaveElevVis{XY}; first dimension is time step; second/third dimensions are grid of elevations" "m,-"
typedef ^ FAST_VTK_SurfaceType FAST_VTK_BLSurfaceType BladeShape {:} - - "AirfoilCoords for each blade" m
typedef ^ FAST_VTK_SurfaceType SiKi MorisonVisRad {:} - - "radius of each Morison node" m
diff --git a/modules/openfast-library/src/FAST_Subs.f90 b/modules/openfast-library/src/FAST_Subs.f90
index 1732086c5f..0954c7b2e6 100644
--- a/modules/openfast-library/src/FAST_Subs.f90
+++ b/modules/openfast-library/src/FAST_Subs.f90
@@ -776,20 +776,6 @@ SUBROUTINE FAST_InitializeAll( t_initial, p_FAST, y_FAST, m_FAST, ED, BD, SrvD,
END IF
- ! ........................
- ! set some VTK parameters required before SeaState init (so we can get wave elevations for visualization)
- ! ........................
-
- ! get wave elevation data for visualization
- if ( p_FAST%WrVTK > VTK_None ) then
- call SetVTKParameters_B4SeaSt(p_FAST, Init%OutData_ED, Init%InData_SeaSt, BD, ErrStat2, ErrMsg2)
- CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
- IF (ErrStat >= AbortErrLev) THEN
- CALL Cleanup()
- RETURN
- END IF
- end if
-
! ........................
! initialize SeaStates
! ........................
@@ -816,9 +802,12 @@ SUBROUTINE FAST_InitializeAll( t_initial, p_FAST, y_FAST, m_FAST, ED, BD, SrvD,
Init%InData_SeaSt%WaveFieldMod = p_FAST%WaveFieldMod
Init%InData_SeaSt%PtfmLocationX = p_FAST%TurbinePos(1)
Init%InData_SeaSt%PtfmLocationY = p_FAST%TurbinePos(2)
-
+
Init%InData_SeaSt%TMax = p_FAST%TMax
+ ! wave field visualization
+ if (p_FAST%WrVTK == VTK_Animate .and. p_FAST%VTK_Type == VTK_Surf) Init%InData_SeaSt%SurfaceVis = .true.
+
CALL SeaSt_Init( Init%InData_SeaSt, SeaSt%Input(1), SeaSt%p, SeaSt%x(STATE_CURR), SeaSt%xd(STATE_CURR), SeaSt%z(STATE_CURR), &
SeaSt%OtherSt(STATE_CURR), SeaSt%y, SeaSt%m, p_FAST%dt_module( MODULE_SeaSt ), Init%OutData_SeaSt, ErrStat2, ErrMsg2 )
CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
@@ -826,7 +815,15 @@ SUBROUTINE FAST_InitializeAll( t_initial, p_FAST, y_FAST, m_FAST, ED, BD, SrvD,
p_FAST%ModuleInitialized(Module_SeaSt) = .TRUE.
CALL SetModuleSubstepTime(Module_SeaSt, p_FAST, y_FAST, ErrStat2, ErrMsg2)
CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
-
+
+ if (allocated(Init%OutData_SeaSt%WaveElevVisGrid)) then
+ p_FAST%VTK_surface%NWaveElevPts(1) = size(Init%OutData_SeaSt%WaveElevVisX)
+ p_FAST%VTK_surface%NWaveElevPts(2) = size(Init%OutData_SeaSt%WaveElevVisX)
+ else
+ p_FAST%VTK_surface%NWaveElevPts(1) = 0
+ p_FAST%VTK_surface%NWaveElevPts(2) = 0
+ endif
+
IF (ErrStat >= AbortErrLev) THEN
CALL Cleanup()
RETURN
@@ -3795,75 +3792,6 @@ end subroutine cleanup
!...............................................................................................................................
END SUBROUTINE FAST_ReadSteadyStateFile
!----------------------------------------------------------------------------------------------------------------------------------
-!> This subroutine sets up some of the information needed for plotting VTK surfaces. It initializes only the data needed before
-!! SeaState initialization. (SeaSt needs some of this data so it can return the wave elevation data we want.)
-SUBROUTINE SetVTKParameters_B4SeaSt(p_FAST, InitOutData_ED, InitInData_SeaSt, BD, ErrStat, ErrMsg)
-
- TYPE(FAST_ParameterType), INTENT(INOUT) :: p_FAST !< The parameters of the glue code
- TYPE(ED_InitOutputType), INTENT(IN ) :: InitOutData_ED !< The initialization output from structural dynamics module
- TYPE(SeaSt_InitInputType), INTENT(INOUT) :: InitInData_SeaSt !< The initialization input to SeaState
- TYPE(BeamDyn_Data), INTENT(IN ) :: BD !< BeamDyn data
- INTEGER(IntKi), INTENT( OUT) :: ErrStat !< Error status of the operation
- CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None
-
-
- REAL(SiKi) :: BladeLength, HubRad, Width, WidthBy2
- REAL(SiKi) :: dx, dy
- INTEGER(IntKi) :: i, j, n
- INTEGER(IntKi) :: ErrStat2
- CHARACTER(ErrMsgLen) :: ErrMsg2
- CHARACTER(*), PARAMETER :: RoutineName = 'SetVTKParameters_B4SeaSt'
-
-
- ErrStat = ErrID_None
- ErrMsg = ""
-
- ! Get radius for ground (blade length + hub radius):
- if ( p_FAST%CompElast == Module_BD ) then
- BladeLength = TwoNorm(BD%y(1)%BldMotion%Position(:,1) - BD%y(1)%BldMotion%Position(:,BD%y(1)%BldMotion%Nnodes))
- HubRad = InitOutData_ED%HubRad
- else
- BladeLength = InitOutData_ED%BladeLength
- HubRad = InitOutData_ED%HubRad
- end if
- p_FAST%VTK_Surface%HubRad = HubRad
- p_FAST%VTK_Surface%GroundRad = BladeLength + HubRad
-
- !........................................................................................................
- ! We don't use the rest of this routine for stick-figure output
- if (p_FAST%VTK_Type /= VTK_Surf) return
- !........................................................................................................
-
- ! initialize wave elevation data:
- if ( p_FAST%CompSeaSt == Module_SeaSt ) then
-
- p_FAST%VTK_surface%NWaveElevPts(1) = 25
- p_FAST%VTK_surface%NWaveElevPts(2) = 25
-
- call allocAry( InitInData_SeaSt%WaveElevXY, 2, p_FAST%VTK_surface%NWaveElevPts(1)*p_FAST%VTK_surface%NWaveElevPts(2), 'WaveElevXY', ErrStat2, ErrMsg2)
- call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
- if (ErrStat >= AbortErrLev) return
-
- Width = p_FAST%VTK_Surface%GroundRad * VTK_GroundFactor
- if (p_FAST%MHK /= MHK_None) Width = Width * 5.0_SiKi
- dx = Width / (p_FAST%VTK_surface%NWaveElevPts(1) - 1)
- dy = Width / (p_FAST%VTK_surface%NWaveElevPts(2) - 1)
-
- WidthBy2 = Width / 2.0_SiKi
- n = 1
- do i=1,p_FAST%VTK_surface%NWaveElevPts(1)
- do j=1,p_FAST%VTK_surface%NWaveElevPts(2)
- InitInData_SeaSt%WaveElevXY(1,n) = dx*(i-1) - WidthBy2 ! SeaSt takes p_FAST%TurbinePos into account already
- InitInData_SeaSt%WaveElevXY(2,n) = dy*(j-1) - WidthBy2
- n = n+1
- end do
- end do
-
- end if
-
-
-END SUBROUTINE SetVTKParameters_B4SeaSt
-!----------------------------------------------------------------------------------------------------------------------------------
!> This subroutine sets up the information needed for plotting VTK surfaces.
SUBROUTINE SetVTKParameters(p_FAST, InitOutData_ED, InitOutData_AD, InitInData_SeaSt, InitOutData_SeaSt, InitOutData_HD, ED, BD, AD, HD, ErrStat, ErrMsg)
@@ -3883,6 +3811,7 @@ SUBROUTINE SetVTKParameters(p_FAST, InitOutData_ED, InitOutData_AD, InitInData_S
REAL(SiKi) :: RefPoint(3), RefLengths(2)
REAL(SiKi) :: x, y
REAL(SiKi) :: TwrDiam_top, TwrDiam_base, TwrRatio, TwrLength
+ REAL(SiKi) :: BladeLength, HubRad
INTEGER(IntKi) :: topNode, baseNode
INTEGER(IntKi) :: NumBl, k, Indx
LOGICAL :: UseADtwr
@@ -3921,10 +3850,18 @@ SUBROUTINE SetVTKParameters(p_FAST, InitOutData_ED, InitOutData_AD, InitInData_S
NumBl = InitOutData_ED%NumBl
! initialize the vtk data
-
p_FAST%VTK_Surface%NumSectors = 25
- ! NOTE: we set p_FAST%VTK_Surface%GroundRad and p_FAST%VTK_Surface%HubRad in SetVTKParameters_B4SeaSt
+ ! Get radius for ground (blade length + hub radius):
+ if ( p_FAST%CompElast == Module_BD ) then
+ BladeLength = TwoNorm(BD%y(1)%BldMotion%Position(:,1) - BD%y(1)%BldMotion%Position(:,BD%y(1)%BldMotion%Nnodes))
+ HubRad = InitOutData_ED%HubRad
+ else
+ BladeLength = InitOutData_ED%BladeLength
+ HubRad = InitOutData_ED%HubRad
+ end if
+ p_FAST%VTK_Surface%HubRad = HubRad
+ p_FAST%VTK_Surface%GroundRad = BladeLength + HubRad
! write the ground or seabed reference polygon:
RefPoint = p_FAST%TurbinePos
@@ -4070,13 +4007,18 @@ SUBROUTINE SetVTKParameters(p_FAST, InitOutData_ED, InitOutData_AD, InitInData_S
!.......................
!bjj: interpolate here instead of each time step?
- if ( allocated(InitOutData_SeaSt%WaveElevSeries) ) then
- call move_alloc( InitInData_SeaSt%WaveElevXY, p_FAST%VTK_Surface%WaveElevXY )
- call move_alloc( InitOutData_SeaSt%WaveElevSeries, p_FAST%VTK_Surface%WaveElev )
+ if ( allocated(InitOutData_SeaSt%WaveElevVisGrid) ) then
+print*,'Storing Wave surface visualization'
+ call move_alloc( InitOutData_SeaSt%WaveElevVisX, p_FAST%VTK_Surface%WaveElevVisX )
+ call move_alloc( InitOutData_SeaSt%WaveElevVisY, p_FAST%VTK_Surface%WaveElevVisY )
+ call move_alloc( InitOutData_SeaSt%WaveElevVisGrid,p_FAST%VTK_Surface%WaveElevVisGrid )
! put the following lines in loops to avoid stack-size issues:
- do k=1,size(p_FAST%VTK_Surface%WaveElevXY,2)
- p_FAST%VTK_Surface%WaveElevXY(:,k) = p_FAST%VTK_Surface%WaveElevXY(:,k) + p_FAST%TurbinePos(1:2)
+ do k=1,size(p_FAST%VTK_Surface%WaveElevVisX)
+ p_FAST%VTK_Surface%WaveElevVisX(k) = p_FAST%VTK_Surface%WaveElevVisX(k) + p_FAST%TurbinePos(1)
+ end do
+ do k=1,size(p_FAST%VTK_Surface%WaveElevVisY)
+ p_FAST%VTK_Surface%WaveElevVisY(k) = p_FAST%VTK_Surface%WaveElevVisY(k) + p_FAST%TurbinePos(2)
end do
end if
@@ -6375,7 +6317,7 @@ SUBROUTINE WrVTK_Surfaces(t_global, p_FAST, y_FAST, MeshMapData, ED, BD, AD, IfW
! Ground (written at initialization)
! Wave elevation
- if ( allocated( p_FAST%VTK_Surface%WaveElev ) ) call WrVTK_WaveElev( t_global, p_FAST, y_FAST, SeaSt)
+ if ( allocated( p_FAST%VTK_Surface%WaveElevVisGrid ) ) call WrVTK_WaveElevVisGrid( t_global, p_FAST, y_FAST, SeaSt)
if (allocated(ED%Input)) then
! Nacelle
@@ -6480,7 +6422,7 @@ SUBROUTINE WrVTK_Surfaces(t_global, p_FAST, y_FAST, MeshMapData, ED, BD, AD, IfW
END SUBROUTINE WrVTK_Surfaces
!----------------------------------------------------------------------------------------------------------------------------------
!> This subroutine writes the wave elevation data for a given time step
-SUBROUTINE WrVTK_WaveElev(t_global, p_FAST, y_FAST, SeaSt)
+SUBROUTINE WrVTK_WaveElevVisGrid(t_global, p_FAST, y_FAST, SeaSt)
REAL(DbKi), INTENT(IN ) :: t_global !< Current global time
TYPE(FAST_ParameterType), INTENT(IN ) :: p_FAST !< Parameters for the glue code
@@ -6499,10 +6441,10 @@ SUBROUTINE WrVTK_WaveElev(t_global, p_FAST, y_FAST, SeaSt)
CHARACTER(1024) :: Tstr
INTEGER(IntKi) :: ErrStat2
CHARACTER(ErrMsgLen) :: ErrMsg2
- CHARACTER(*),PARAMETER :: RoutineName = 'WrVTK_WaveElev'
+ CHARACTER(*),PARAMETER :: RoutineName = 'WrVTK_WaveElevVisGrid'
- NumberOfPoints = size(p_FAST%VTK_surface%WaveElevXY,2)
+ NumberOfPoints = p_FAST%VTK_surface%NWaveElevPts(1) * p_FAST%VTK_surface%NWaveElevPts(2)
! I'm going to make triangles for now. we should probably just make this a structured file at some point
NumberOfPolys = ( p_FAST%VTK_surface%NWaveElevPts(1) - 1 ) * &
( p_FAST%VTK_surface%NWaveElevPts(2) - 1 ) * 2
@@ -6520,49 +6462,47 @@ SUBROUTINE WrVTK_WaveElev(t_global, p_FAST, y_FAST, SeaSt)
if (ErrStat2 >= AbortErrLev) return
! points (nodes, augmented with NumSegments):
- WRITE(Un,'(A)') ' '
- WRITE(Un,'(A)') ' '
+ WRITE(Un,'(A)') ' '
+ WRITE(Un,'(A)') ' '
- ! I'm not going to interpolate in time; I'm just going to get the index of the closest wave time value
- t = REAL(t_global,SiKi)
- call GetWaveElevIndx( t, SeaSt%p%WaveField%WaveTime, y_FAST%VTK_LastWaveIndx )
+ ! I'm not going to interpolate in time; I'm just going to get the index of the closest wave time value
+ t = REAL(t_global,SiKi)
+ call GetWaveElevIndx( t, SeaSt%p%WaveField%WaveTime, y_FAST%VTK_LastWaveIndx )
- n = 1
- do ix=1,p_FAST%VTK_surface%NWaveElevPts(1)
- do iy=1,p_FAST%VTK_surface%NWaveElevPts(2)
- WRITE(Un,VTK_AryFmt) p_FAST%VTK_surface%WaveElevXY(:,n), p_FAST%VTK_surface%WaveElev(y_FAST%VTK_LastWaveIndx,n)
- n = n+1
- end do
+ do ix=1,p_FAST%VTK_surface%NWaveElevPts(1)
+ do iy=1,p_FAST%VTK_surface%NWaveElevPts(2)
+ WRITE(Un,VTK_AryFmt) p_FAST%VTK_surface%WaveElevVisX(ix), p_FAST%VTK_surface%WaveElevVisX(iy), p_FAST%VTK_surface%WaveElevVisGrid(y_FAST%VTK_LastWaveIndx,ix,iy)
end do
+ end do
- WRITE(Un,'(A)') ' '
- WRITE(Un,'(A)') ' '
+ WRITE(Un,'(A)') ' '
+ WRITE(Un,'(A)') ' '
- WRITE(Un,'(A)') ' '
- WRITE(Un,'(A)') ' '
+ WRITE(Un,'(A)') ' '
+ WRITE(Un,'(A)') ' '
- do ix=1,p_FAST%VTK_surface%NWaveElevPts(1)-1
- do iy=1,p_FAST%VTK_surface%NWaveElevPts(2)-1
- n = p_FAST%VTK_surface%NWaveElevPts(1)*(ix-1)+iy - 1 ! points start at 0
+ do ix=1,p_FAST%VTK_surface%NWaveElevPts(1)-1
+ do iy=1,p_FAST%VTK_surface%NWaveElevPts(2)-1
+ n = p_FAST%VTK_surface%NWaveElevPts(1)*(ix-1)+iy - 1 ! points start at 0
- WRITE(Un,'(3(i7))') n, n+1, n+p_FAST%VTK_surface%NWaveElevPts(2)
- WRITE(Un,'(3(i7))') n+1, n+1+p_FAST%VTK_surface%NWaveElevPts(2), n+p_FAST%VTK_surface%NWaveElevPts(2)
+ WRITE(Un,'(3(i7))') n, n+1, n+p_FAST%VTK_surface%NWaveElevPts(2)
+ WRITE(Un,'(3(i7))') n+1, n+1+p_FAST%VTK_surface%NWaveElevPts(2), n+p_FAST%VTK_surface%NWaveElevPts(2)
- end do
end do
- WRITE(Un,'(A)') ' '
+ end do
+ WRITE(Un,'(A)') ' '
- WRITE(Un,'(A)') ' '
- do n=1,NumberOfPolys
- WRITE(Un,'(i7)') 3*n
- end do
- WRITE(Un,'(A)') ' '
- WRITE(Un,'(A)') ' '
+ WRITE(Un,'(A)') ' '
+ do n=1,NumberOfPolys
+ WRITE(Un,'(i7)') 3*n
+ end do
+ WRITE(Un,'(A)') ' '
+ WRITE(Un,'(A)') ' '
- call WrVTK_footer( Un )
+ call WrVTK_footer( Un )
-END SUBROUTINE WrVTK_WaveElev
+END SUBROUTINE WrVTK_WaveElevVisGrid
!----------------------------------------------------------------------------------------------------------------------------------
!> This function returns the index, Ind, of the XAry closest to XValIn, where XAry is assumed to be periodic. It starts
!! searching at the value of Ind from a previous step.
diff --git a/modules/openfast-library/src/FAST_Types.f90 b/modules/openfast-library/src/FAST_Types.f90
index 8c5b43ef30..69cf370a60 100644
--- a/modules/openfast-library/src/FAST_Types.f90
+++ b/modules/openfast-library/src/FAST_Types.f90
@@ -92,8 +92,9 @@ MODULE FAST_Types
REAL(SiKi) , DIMENSION(1:3,1:8) :: NacelleBox = 0.0_R4Ki !< X-Y-Z locations of 8 points that define the nacelle box, relative to the nacelle position [m]
REAL(SiKi) , DIMENSION(:), ALLOCATABLE :: TowerRad !< radius of each ED tower node [m]
INTEGER(IntKi) , DIMENSION(1:2) :: NWaveElevPts = 0_IntKi !< number of points for wave elevation visualization [-]
- REAL(SiKi) , DIMENSION(:,:), ALLOCATABLE :: WaveElevXY !< X-Y locations for WaveElev output (for visualization). First dimension is the X (1) and Y (2) coordinate. Second dimension is the point number. [m,-]
- REAL(SiKi) , DIMENSION(:,:), ALLOCATABLE :: WaveElev !< wave elevation at WaveElevXY; first dimension is time step; second dimension is point number [m,-]
+ REAL(SiKi) , DIMENSION(:), ALLOCATABLE :: WaveElevVisX !< X locations for WaveElev output (for visualization). [m,-]
+ REAL(SiKi) , DIMENSION(:), ALLOCATABLE :: WaveElevVisY !< Y locations for WaveElev output (for visualization). [m,-]
+ REAL(SiKi) , DIMENSION(:,:,:), ALLOCATABLE :: WaveElevVisGrid !< wave elevation at WaveElevVis{XY}; first dimension is time step; second/third dimensions are grid of elevations [m,-]
TYPE(FAST_VTK_BLSurfaceType) , DIMENSION(:), ALLOCATABLE :: BladeShape !< AirfoilCoords for each blade [m]
REAL(SiKi) , DIMENSION(:), ALLOCATABLE :: MorisonVisRad !< radius of each Morison node [m]
END TYPE FAST_VTK_SurfaceType
@@ -896,8 +897,8 @@ subroutine FAST_CopyVTK_SurfaceType(SrcVTK_SurfaceTypeData, DstVTK_SurfaceTypeDa
integer(IntKi), intent(in ) :: CtrlCode
integer(IntKi), intent( out) :: ErrStat
character(*), intent( out) :: ErrMsg
- integer(B8Ki) :: i1, i2
- integer(B8Ki) :: LB(2), UB(2)
+ integer(B8Ki) :: i1, i2, i3
+ integer(B8Ki) :: LB(3), UB(3)
integer(IntKi) :: ErrStat2
character(ErrMsgLen) :: ErrMsg2
character(*), parameter :: RoutineName = 'FAST_CopyVTK_SurfaceType'
@@ -920,29 +921,41 @@ subroutine FAST_CopyVTK_SurfaceType(SrcVTK_SurfaceTypeData, DstVTK_SurfaceTypeDa
DstVTK_SurfaceTypeData%TowerRad = SrcVTK_SurfaceTypeData%TowerRad
end if
DstVTK_SurfaceTypeData%NWaveElevPts = SrcVTK_SurfaceTypeData%NWaveElevPts
- if (allocated(SrcVTK_SurfaceTypeData%WaveElevXY)) then
- LB(1:2) = lbound(SrcVTK_SurfaceTypeData%WaveElevXY, kind=B8Ki)
- UB(1:2) = ubound(SrcVTK_SurfaceTypeData%WaveElevXY, kind=B8Ki)
- if (.not. allocated(DstVTK_SurfaceTypeData%WaveElevXY)) then
- allocate(DstVTK_SurfaceTypeData%WaveElevXY(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2)
+ if (allocated(SrcVTK_SurfaceTypeData%WaveElevVisX)) then
+ LB(1:1) = lbound(SrcVTK_SurfaceTypeData%WaveElevVisX, kind=B8Ki)
+ UB(1:1) = ubound(SrcVTK_SurfaceTypeData%WaveElevVisX, kind=B8Ki)
+ if (.not. allocated(DstVTK_SurfaceTypeData%WaveElevVisX)) then
+ allocate(DstVTK_SurfaceTypeData%WaveElevVisX(LB(1):UB(1)), stat=ErrStat2)
if (ErrStat2 /= 0) then
- call SetErrStat(ErrID_Fatal, 'Error allocating DstVTK_SurfaceTypeData%WaveElevXY.', ErrStat, ErrMsg, RoutineName)
+ call SetErrStat(ErrID_Fatal, 'Error allocating DstVTK_SurfaceTypeData%WaveElevVisX.', ErrStat, ErrMsg, RoutineName)
return
end if
end if
- DstVTK_SurfaceTypeData%WaveElevXY = SrcVTK_SurfaceTypeData%WaveElevXY
+ DstVTK_SurfaceTypeData%WaveElevVisX = SrcVTK_SurfaceTypeData%WaveElevVisX
end if
- if (allocated(SrcVTK_SurfaceTypeData%WaveElev)) then
- LB(1:2) = lbound(SrcVTK_SurfaceTypeData%WaveElev, kind=B8Ki)
- UB(1:2) = ubound(SrcVTK_SurfaceTypeData%WaveElev, kind=B8Ki)
- if (.not. allocated(DstVTK_SurfaceTypeData%WaveElev)) then
- allocate(DstVTK_SurfaceTypeData%WaveElev(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2)
+ if (allocated(SrcVTK_SurfaceTypeData%WaveElevVisY)) then
+ LB(1:1) = lbound(SrcVTK_SurfaceTypeData%WaveElevVisY, kind=B8Ki)
+ UB(1:1) = ubound(SrcVTK_SurfaceTypeData%WaveElevVisY, kind=B8Ki)
+ if (.not. allocated(DstVTK_SurfaceTypeData%WaveElevVisY)) then
+ allocate(DstVTK_SurfaceTypeData%WaveElevVisY(LB(1):UB(1)), stat=ErrStat2)
if (ErrStat2 /= 0) then
- call SetErrStat(ErrID_Fatal, 'Error allocating DstVTK_SurfaceTypeData%WaveElev.', ErrStat, ErrMsg, RoutineName)
+ call SetErrStat(ErrID_Fatal, 'Error allocating DstVTK_SurfaceTypeData%WaveElevVisY.', ErrStat, ErrMsg, RoutineName)
return
end if
end if
- DstVTK_SurfaceTypeData%WaveElev = SrcVTK_SurfaceTypeData%WaveElev
+ DstVTK_SurfaceTypeData%WaveElevVisY = SrcVTK_SurfaceTypeData%WaveElevVisY
+ end if
+ if (allocated(SrcVTK_SurfaceTypeData%WaveElevVisGrid)) then
+ LB(1:3) = lbound(SrcVTK_SurfaceTypeData%WaveElevVisGrid, kind=B8Ki)
+ UB(1:3) = ubound(SrcVTK_SurfaceTypeData%WaveElevVisGrid, kind=B8Ki)
+ if (.not. allocated(DstVTK_SurfaceTypeData%WaveElevVisGrid)) then
+ allocate(DstVTK_SurfaceTypeData%WaveElevVisGrid(LB(1):UB(1),LB(2):UB(2),LB(3):UB(3)), stat=ErrStat2)
+ if (ErrStat2 /= 0) then
+ call SetErrStat(ErrID_Fatal, 'Error allocating DstVTK_SurfaceTypeData%WaveElevVisGrid.', ErrStat, ErrMsg, RoutineName)
+ return
+ end if
+ end if
+ DstVTK_SurfaceTypeData%WaveElevVisGrid = SrcVTK_SurfaceTypeData%WaveElevVisGrid
end if
if (allocated(SrcVTK_SurfaceTypeData%BladeShape)) then
LB(1:1) = lbound(SrcVTK_SurfaceTypeData%BladeShape, kind=B8Ki)
@@ -978,8 +991,8 @@ subroutine FAST_DestroyVTK_SurfaceType(VTK_SurfaceTypeData, ErrStat, ErrMsg)
type(FAST_VTK_SurfaceType), intent(inout) :: VTK_SurfaceTypeData
integer(IntKi), intent( out) :: ErrStat
character(*), intent( out) :: ErrMsg
- integer(B8Ki) :: i1, i2
- integer(B8Ki) :: LB(2), UB(2)
+ integer(B8Ki) :: i1, i2, i3
+ integer(B8Ki) :: LB(3), UB(3)
integer(IntKi) :: ErrStat2
character(ErrMsgLen) :: ErrMsg2
character(*), parameter :: RoutineName = 'FAST_DestroyVTK_SurfaceType'
@@ -988,11 +1001,14 @@ subroutine FAST_DestroyVTK_SurfaceType(VTK_SurfaceTypeData, ErrStat, ErrMsg)
if (allocated(VTK_SurfaceTypeData%TowerRad)) then
deallocate(VTK_SurfaceTypeData%TowerRad)
end if
- if (allocated(VTK_SurfaceTypeData%WaveElevXY)) then
- deallocate(VTK_SurfaceTypeData%WaveElevXY)
+ if (allocated(VTK_SurfaceTypeData%WaveElevVisX)) then
+ deallocate(VTK_SurfaceTypeData%WaveElevVisX)
+ end if
+ if (allocated(VTK_SurfaceTypeData%WaveElevVisY)) then
+ deallocate(VTK_SurfaceTypeData%WaveElevVisY)
end if
- if (allocated(VTK_SurfaceTypeData%WaveElev)) then
- deallocate(VTK_SurfaceTypeData%WaveElev)
+ if (allocated(VTK_SurfaceTypeData%WaveElevVisGrid)) then
+ deallocate(VTK_SurfaceTypeData%WaveElevVisGrid)
end if
if (allocated(VTK_SurfaceTypeData%BladeShape)) then
LB(1:1) = lbound(VTK_SurfaceTypeData%BladeShape, kind=B8Ki)
@@ -1012,8 +1028,8 @@ subroutine FAST_PackVTK_SurfaceType(RF, Indata)
type(RegFile), intent(inout) :: RF
type(FAST_VTK_SurfaceType), intent(in) :: InData
character(*), parameter :: RoutineName = 'FAST_PackVTK_SurfaceType'
- integer(B8Ki) :: i1, i2
- integer(B8Ki) :: LB(2), UB(2)
+ integer(B8Ki) :: i1, i2, i3
+ integer(B8Ki) :: LB(3), UB(3)
if (RF%ErrStat >= AbortErrLev) return
call RegPack(RF, InData%NumSectors)
call RegPack(RF, InData%HubRad)
@@ -1021,8 +1037,9 @@ subroutine FAST_PackVTK_SurfaceType(RF, Indata)
call RegPack(RF, InData%NacelleBox)
call RegPackAlloc(RF, InData%TowerRad)
call RegPack(RF, InData%NWaveElevPts)
- call RegPackAlloc(RF, InData%WaveElevXY)
- call RegPackAlloc(RF, InData%WaveElev)
+ call RegPackAlloc(RF, InData%WaveElevVisX)
+ call RegPackAlloc(RF, InData%WaveElevVisY)
+ call RegPackAlloc(RF, InData%WaveElevVisGrid)
call RegPack(RF, allocated(InData%BladeShape))
if (allocated(InData%BladeShape)) then
call RegPackBounds(RF, 1, lbound(InData%BladeShape, kind=B8Ki), ubound(InData%BladeShape, kind=B8Ki))
@@ -1040,8 +1057,8 @@ subroutine FAST_UnPackVTK_SurfaceType(RF, OutData)
type(RegFile), intent(inout) :: RF
type(FAST_VTK_SurfaceType), intent(inout) :: OutData
character(*), parameter :: RoutineName = 'FAST_UnPackVTK_SurfaceType'
- integer(B8Ki) :: i1, i2
- integer(B8Ki) :: LB(2), UB(2)
+ integer(B8Ki) :: i1, i2, i3
+ integer(B8Ki) :: LB(3), UB(3)
integer(IntKi) :: stat
logical :: IsAllocAssoc
if (RF%ErrStat /= ErrID_None) return
@@ -1051,8 +1068,9 @@ subroutine FAST_UnPackVTK_SurfaceType(RF, OutData)
call RegUnpack(RF, OutData%NacelleBox); if (RegCheckErr(RF, RoutineName)) return
call RegUnpackAlloc(RF, OutData%TowerRad); if (RegCheckErr(RF, RoutineName)) return
call RegUnpack(RF, OutData%NWaveElevPts); if (RegCheckErr(RF, RoutineName)) return
- call RegUnpackAlloc(RF, OutData%WaveElevXY); if (RegCheckErr(RF, RoutineName)) return
- call RegUnpackAlloc(RF, OutData%WaveElev); if (RegCheckErr(RF, RoutineName)) return
+ call RegUnpackAlloc(RF, OutData%WaveElevVisX); if (RegCheckErr(RF, RoutineName)) return
+ call RegUnpackAlloc(RF, OutData%WaveElevVisY); if (RegCheckErr(RF, RoutineName)) return
+ call RegUnpackAlloc(RF, OutData%WaveElevVisGrid); if (RegCheckErr(RF, RoutineName)) return
if (allocated(OutData%BladeShape)) deallocate(OutData%BladeShape)
call RegUnpack(RF, IsAllocAssoc); if (RegCheckErr(RF, RoutineName)) return
if (IsAllocAssoc) then
diff --git a/modules/seastate/src/SeaState.f90 b/modules/seastate/src/SeaState.f90
index 9d56069af5..0e3e27033d 100644
--- a/modules/seastate/src/SeaState.f90
+++ b/modules/seastate/src/SeaState.f90
@@ -381,35 +381,11 @@ SUBROUTINE SeaSt_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, Init
! If requested, output wave elevation data for VTK visualization
-
- IF (ALLOCATED(InitInp%WaveElevXY)) THEN
- ! maybe instead of getting these requested points, we just output the grid that SeaState is generated on?
- ALLOCATE(InitOut%WaveElevSeries( 0:p%WaveField%NStepWave, 1:SIZE(InitInp%WaveElevXY, DIM=2)),STAT=ErrStat2)
- if (ErrStat2 /= 0) then
- CALL SetErrStat(ErrID_Fatal,"Error allocating InitOut%WaveElevSeries.",ErrStat,ErrMsg,RoutineName)
- CALL CleanUp()
- RETURN
- end if
-
- do it = 1,size(p%WaveField%WaveTime)
- do i = 1, size(InitOut%WaveElevSeries,DIM=2)
- InitOut%WaveElevSeries(it,i) = SeaSt_Interp_3D( real(p%WaveField%WaveTime(it),DbKi), real(InitInp%WaveElevXY(:,i),ReKi), p%WaveField%WaveElev1, p%WaveField%seast_interp_p, m%seast_interp_m%FirstWarn_Clamp, ErrStat2, ErrMsg2 )
- call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
- end do
- end do
-
- if (allocated(p%WaveField%WaveElev2)) then
- do it = 1,size(p%WaveField%WaveTime)
- do i = 1, size(InitOut%WaveElevSeries,DIM=2)
- TmpElev = SeaSt_Interp_3D( real(p%WaveField%WaveTime(it),DbKi), real(InitInp%WaveElevXY(:,i),ReKi), p%WaveField%WaveElev2, p%WaveField%seast_interp_p, m%seast_interp_m%FirstWarn_Clamp, ErrStat2, ErrMsg2 )
- call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
- InitOut%WaveElevSeries(it,i) = InitOut%WaveElevSeries(it,i) + TmpElev
- end do
- end do
- end if
-
-
- ENDIF
+ if (InitInp%SurfaceVis) then
+ call SurfaceVisGenerate(ErrStat2, ErrMsg2)
+ call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
+ if (ErrStat >= AbortErrLev) return
+ endif
IF ( InitInp%hasIce ) THEN
@@ -468,6 +444,88 @@ SUBROUTINE CleanUp()
END SUBROUTINE CleanUp
!................................
+ subroutine SurfaceVisGenerate(ErrStat3, ErrMsg3)
+ integer(IntKi), intent( out) :: ErrStat3
+ character(ErrMsgLen),intent( out) :: ErrMsg3
+ integer(IntKi) :: Nx,Ny,i1,i2
+ real(SiKi) :: HWidX, HWidY, dx, dy, TmpElev
+ real(ReKi) :: loc(2) ! location (x,y)
+ integer(IntKi) :: ErrStat4
+ character(ErrMsgLen) :: ErrMsg4
+ character(*), parameter :: RtnName="SurfaceVisGenerate"
+
+ ErrStat3 = ErrID_None
+ ErrMsg3 = ""
+
+ ! Grid half width from the WaveField
+ HWidX = (real(p%WaveField%SeaSt_Interp_p%n(2)-1,SiKi)) / 2.0_SiKi * p%WaveField%SeaSt_Interp_p%delta(2)
+ HWidY = (real(p%WaveField%SeaSt_Interp_p%n(3)-1,SiKi)) / 2.0_SiKi * p%WaveField%SeaSt_Interp_p%delta(3)
+
+ if ((InitInp%SurfaceVisNx <= 0) .or. (InitInp%SurfaceVisNy <= 0))then ! use the SeaState points exactly
+ ! Set number of points to the number of seastate grid points in each direction
+ Nx = p%WaveField%SeaSt_Interp_p%n(2)
+ Ny = p%WaveField%SeaSt_Interp_p%n(3)
+ dx = p%WaveField%SeaSt_Interp_p%delta(2)
+ dy = p%WaveField%SeaSt_Interp_p%delta(3)
+ call SetErrStat(ErrID_Info,"Setting wavefield visualization grid to "//trim(Num2LStr(Nx))//" x "//trim(Num2LStr(Ny))//"points",ErrStat3,ErrMsg3,RoutineName)
+ elseif ((InitInp%SurfaceVisNx < 3) .or. (InitInp%SurfaceVisNx < 3)) then ! Set to 3 for minimum
+ Nx = 3
+ Ny = 3
+ dx = HWidX
+ dy = HWidY
+ call SetErrStat(ErrID_Warn,"Setting wavefield visualization grid to 3 points in each direction",ErrStat3,ErrMsg3,RoutineName)
+ else ! Specified number of points
+ Nx = InitInp%SurfaceVisNx
+ Ny = InitInp%SurfaceVisNy
+ dx = 2.0_SiKi * HWidX / (real(Nx,SiKi)-1)
+ dy = 2.0_SiKi * HWidY / (real(Ny,SiKi)-1)
+ endif
+
+ ! allocate arrays
+ call AllocAry(InitOut%WaveElevVisX,Nx,"InitOut%NWaveElevVisX",ErrStat4,ErrMsg4)
+ call SetErrStat(ErrStat4,ErrMsg4,ErrStat3,ErrMsg3,RtnName)
+ call AllocAry(InitOut%WaveElevVisY,Ny,"InitOut%NWaveElevVisY",ErrStat4,ErrMsg4)
+ call SetErrStat(ErrStat4,ErrMsg4,ErrStat3,ErrMsg3,RtnName)
+ allocate(InitOut%WaveElevVisGrid( 0:size(p%WaveField%WaveTime),Nx,Ny ),STAT=ErrStat4)
+ if (ErrStat4 /= 0) then
+ CALL SetErrStat(ErrID_Fatal,"Error allocating InitOut%WaveElevVisGrid.",ErrStat3,ErrMsg3,RoutineName)
+ return
+ end if
+
+ ! Populate the arrays
+ do i1=1,Nx
+ InitOut%WaveElevVisX(i1) = -HWidX + real(i1-1,SiKi)*dx
+ enddo
+ do i2=1,Ny
+ InitOut%WaveElevVisY(i2) = -HWidY + real(i2-1,SiKi)*dy
+ enddo
+
+!FIXME: calculate from the FFT of the data.
+ do it = 0,size(p%WaveField%WaveTime)-1
+ do i1 = 1, nx
+ loc(1) = InitOut%WaveElevVisX(i1)
+ do i2 = 1, ny
+ loc(2) = InitOut%WaveElevVisX(i2)
+ InitOut%WaveElevVisGrid(it,i1,i2) = SeaSt_Interp_3D( real(p%WaveField%WaveTime(it),DbKi), real(loc,ReKi), p%WaveField%WaveElev1, p%WaveField%seast_interp_p, m%seast_interp_m%FirstWarn_Clamp, ErrStat4, ErrMsg4 )
+ call SetErrStat( ErrStat4, ErrMsg4, ErrStat3, ErrMsg3, RoutineName )
+ enddo
+ end do
+ end do
+
+ if (allocated(p%WaveField%WaveElev2)) then
+ do it = 0,size(p%WaveField%WaveTime)-1
+ do i1 = 1, nx
+ loc(1) = InitOut%WaveElevVisX(i1)
+ do i2 = 1, ny
+ loc(2) = InitOut%WaveElevVisX(i2)
+ TmpElev = SeaSt_Interp_3D( real(p%WaveField%WaveTime(it),DbKi), real(loc,ReKi), p%WaveField%WaveElev2, p%WaveField%seast_interp_p, m%seast_interp_m%FirstWarn_Clamp, ErrStat4, ErrMsg4 )
+ call SetErrStat( ErrStat4, ErrMsg4, ErrStat3, ErrMsg3, RoutineName )
+ InitOut%WaveElevVisGrid(it,i1,i2) = InitOut%WaveElevVisGrid(it,i1,i2) + TmpElev
+ end do
+ end do
+ end do
+ end if
+ end subroutine SurfaceVisGenerate
END SUBROUTINE SeaSt_Init
!----------------------------------------------------------------------------------------------------------------------------------
diff --git a/modules/seastate/src/SeaState.txt b/modules/seastate/src/SeaState.txt
index 03e659fac1..1ef0d93440 100644
--- a/modules/seastate/src/SeaState.txt
+++ b/modules/seastate/src/SeaState.txt
@@ -73,13 +73,15 @@ typedef ^ ^ ReKi def
typedef ^ ^ ReKi defWtrDpth - - - "Default water depth from the driver; may be overwritten " "m"
typedef ^ ^ ReKi defMSL2SWL - - - "Default mean sea level to still water level from the driver; may be overwritten" "m"
typedef ^ ^ DbKi TMax - - - "Supplied by Driver: The total simulation time" "(sec)"
-typedef ^ ^ SiKi WaveElevXY {:}{:} - - "Supplied by Driver: X-Y locations for WaveElevation output (for visualization). First dimension is the X (1) and Y (2) coordinate. Second dimension is the point number." "m,-"
typedef ^ ^ INTEGER WaveFieldMod - - - "Wave field handling (-) (switch) 0: use individual SeaState inputs without adjustment, 1: adjust wave phases based on turbine offsets from farm origin" -
typedef ^ ^ ReKi PtfmLocationX - - - "Supplied by Driver: X coordinate of platform location in the wave field" "m"
typedef ^ ^ ReKi PtfmLocationY - - - "Supplied by Driver: Y coordinate of platform location in the wave field" "m"
typedef ^ ^ IntKi WrWvKinMod - 0 - "0,1, or 2 indicating whether we are going to write out kinematics files. [ignored if WaveMod = 6, if 1 or 2 then files are written using the outrootname]" -
typedef ^ ^ LOGICAL HasIce - - - "Supplied by Driver: Whether this simulation has ice loading (flag)" -
typedef ^ ^ Logical Linearize - .FALSE. - "Flag that tells this module if the glue code wants to linearize." -
+typedef ^ ^ Logical SurfaceVis - .FALSE. - "Turn on grid surface visualization outputs" -
+typedef ^ ^ IntKi SurfaceVisNx - 0 - "Number of points in X direction to output for visualization grid. Use 0 or negative to set to SeaState resolution." -
+typedef ^ ^ IntKi SurfaceVisNy - 0 - "Number of points in Y direction to output for visualization grid. Use 0 or negative to set to SeaState resolution." -
#
#
@@ -89,8 +91,11 @@ typedef ^ InitOutputType CHARACTER(ChanLen) Wri
typedef ^ ^ CHARACTER(ChanLen) WriteOutputUnt {:} - - "The is the list of all HD-related output channel unit strings (includes all sub-module channels)" -
typedef ^ ^ ProgDesc Ver - - - "Version of SeaState"
typedef ^ ^ LOGICAL InvalidWithSSExctn - - - "Whether SeaState configuration is invalid with HydroDyn's state-space excitation (ExctnMod=2)" (-)
-typedef ^ ^ SiKi WaveElevSeries {:}{:} - - "Wave elevation time-series at each of the points given by WaveElevXY. First dimension is the timestep. Second dimension is XY point number corresponding to second dimension of WaveElevXY." (m)
+typedef ^ ^ SiKi WaveElevVisX {:} - - "X locations of grid output" "m,-"
+typedef ^ ^ SiKi WaveElevVisY {:} - - "X locations of grid output" "m,-"
+typedef ^ ^ SiKi WaveElevVisGrid {:}{:}{:} - - "Wave elevation time-series at each of the points given by WaveElevXY. First dimension is the timestep. Second/third dimensions are the grid of points." (m)
typedef ^ ^ SeaSt_WaveFieldType *WaveField - - - "Pointer to wave field"
+
#
#
# ..... States ....................................................................................................................
diff --git a/modules/seastate/src/SeaState_DriverCode.f90 b/modules/seastate/src/SeaState_DriverCode.f90
index 276352fbd3..885a1f50e7 100644
--- a/modules/seastate/src/SeaState_DriverCode.f90
+++ b/modules/seastate/src/SeaState_DriverCode.f90
@@ -42,11 +42,9 @@ program SeaStateDriver
integer :: WrWvKinMod
integer :: NSteps
real(DbKi) :: TimeInterval
- logical :: WaveElevSeriesFlag !< Should we put together a wave elevation series and save it to file?
- real(ReKi) :: WaveElevdX !< Spacing in the X direction for wave elevation series (m)
- real(ReKi) :: WaveElevdY !< Spacing in the Y direction for the wave elevation series (m)
- integer(IntKi) :: WaveElevNX !< Number of points in the X direction for the wave elevation series (-)
- integer(IntKi) :: WaveElevNY !< Number of points in the X direction for the wave elevation series (-)
+ logical :: WaveElevVis !< Should we put together a wave elevation series and save it to file?
+ integer(IntKi) :: WaveElevVisNx !< Number of points in the X direction for the wave elevation series (-)
+ integer(IntKi) :: WaveElevVisNy !< Number of points in the X direction for the wave elevation series (-)
end type SeaSt_Drvr_InitInput
! -----------------------------------------------------------------------------------
@@ -178,28 +176,14 @@ program SeaStateDriver
!-------------------------------------------------------------------------------------
! Setup the arrays for the wave elevation timeseries if requested by the driver input file
- !if ( drvrInitInp%WaveElevSeriesFlag ) then
- ! ALLOCATE ( InitInData%WaveElevXY(2,drvrInitInp%WaveElevNX*drvrInitInp%WaveElevNY), STAT=ErrStat )
- ! if ( ErrStat >= ErrID_Fatal ) then
- ! call SeaSt_End( u(1), p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg )
- ! if ( ErrStat /= ErrID_None ) then
- ! call WrScr( ErrMsg )
- ! end if
- ! stop
- ! end if
- !
- ! ! Set the values
- ! n = 0 ! Dummy counter we are using to get the current point number
- ! do I = 0,drvrInitInp%WaveElevNX-1
- ! do J = 0, drvrInitInp%WaveElevNY-1
- ! n = n+1
- ! ! X dimension
- ! InitInData%WaveElevXY(1,n) = drvrInitInp%WaveElevDX*(I - 0.5*(drvrInitInp%WaveElevNX-1))
- ! ! Y dimension
- ! InitInData%WaveElevXY(2,n) = drvrInitInp%WaveElevDY*(J - 0.5*(drvrInitInp%WaveElevNY-1))
- ! ENDDO
- ! ENDDO
- !endif
+ if ( drvrInitInp%WaveElevVis ) then
+ InitInData%SurfaceVis = .true.
+!FIXME: enable this when we can use an arbitrary number of points from the FFT of the data.
+ !InitInData%SurfaceVisNx = drvrInitInp%WaveElevVisNx ! Number of points in X
+ !InitInData%SurfaceVisNy = drvrInitInp%WaveElevVisNy ! Number of points in Y
+ InitInData%SurfaceVisNx = 0 ! use the WaveField grid resolution
+ InitInData%SurfaceVisNy = 0 ! use the WaveField grid resolution
+ endif
! Initialize the module
Interval = drvrInitInp%TimeInterval
@@ -217,7 +201,7 @@ program SeaStateDriver
! Write the gridded wave elevation data to a file
- if ( drvrInitInp%WaveElevSeriesFlag ) call WaveElevGrid_Output (drvrInitInp, InitInData, InitOutData, p, ErrStat, ErrMsg)
+ if ( drvrInitInp%WaveElevVis ) call WaveElevGrid_Output (drvrInitInp, InitInData, InitOutData, p, ErrStat, ErrMsg)
if (errStat >= AbortErrLev) then
! Clean up and exit
call SeaSt_DvrCleanup()
@@ -434,13 +418,13 @@ SUBROUTINE ReadDriverInputFile( inputFile, InitInp, ErrStat, ErrMsg )
if (InitInp%Echo .and. UnEchoLocal>0) close(UnEchoLocal)
close( UnIn )
return
- end if
+ end if
end if
!-------------------------------------------------------------------------------------------------
! Environmental conditions section
!-------------------------------------------------------------------------------------------------
-
+
! Header
call ReadCom( UnIn, FileName, 'Environmental conditions header', ErrStat, ErrMsg, UnEchoLocal )
@@ -499,7 +483,7 @@ SUBROUTINE ReadDriverInputFile( inputFile, InitInp, ErrStat, ErrMsg )
! Header
call ReadCom( UnIn, FileName, 'SeaState header', ErrStat, ErrMsg, UnEchoLocal )
-
+
if ( ErrStat >= AbortErrLev ) then
if (InitInp%Echo .and. UnEchoLocal>0) close(UnEchoLocal)
close( UnIn )
@@ -580,21 +564,39 @@ SUBROUTINE ReadDriverInputFile( inputFile, InitInp, ErrStat, ErrMsg )
!> Header
call ReadCom( UnIn, FileName, 'Waves multipoint elevation output header', ErrStat, ErrMsg, UnEchoLocal )
-
+
if ( ErrStat >= AbortErrLev ) then
- if (InitInp%Echo .and. UnEchoLocal>0) close(UnEchoLocal)
- close( UnIn )
+ if (InitInp%Echo .and. UnEchoLocal>0) close(UnEchoLocal)
+ close( UnIn )
return
end if
!> WaveElevSeriesFlag -- are we doing multipoint wave elevation output?
- call ReadVar ( UnIn, FileName, InitInp%WaveElevSeriesFlag, 'WaveElevSeriesFlag', 'WaveElevSeriesFlag', ErrStat, ErrMsg )
+ call ReadVar ( UnIn, FileName, InitInp%WaveElevVis, 'WaveElevVis', 'WaveElevVis', ErrStat, ErrMsg )
if ( ErrStat >= AbortErrLev ) then
if (InitInp%Echo .and. UnEchoLocal>0) close(UnEchoLocal)
close( UnIn )
return
end if
+!FIXME: enable this when we can use an arbitrary number of points from the FFT of the data.
+! !> WaveElevVisNx -- number of points in X if visualizing
+! call ReadVar ( UnIn, FileName, InitInp%WaveElevVisNx, 'WaveElevVisNX', 'WaveElevVisNx', ErrStat, ErrMsg )
+! if ( ErrStat >= AbortErrLev ) then
+! if (InitInp%Echo .and. UnEchoLocal>0) close(UnEchoLocal)
+! close( UnIn )
+! return
+! end if
+!
+! !> WaveElevVisNy -- number of points in Y if visualizing
+! call ReadVar ( UnIn, FileName, InitInp%WaveElevVisNy, 'WaveElevVisNy', 'WaveElevVisNy', ErrStat, ErrMsg )
+! if ( ErrStat >= AbortErrLev ) then
+! if (InitInp%Echo .and. UnEchoLocal>0) close(UnEchoLocal)
+! close( UnIn )
+! return
+! end if
+
+
if (InitInp%Echo .and. UnEchoLocal>0) close(UnEchoLocal)
close( UnIn )
@@ -667,15 +669,11 @@ SUBROUTINE WaveElevGrid_Output (drvrInitInp, SeaStateInitInp, SeaStateInitOut, S
write (WaveElevFileUn,'(A)', IOSTAT=ErrStatTmp ) NewLine
write (WaveElevFileUn,'(A8,F10.3)', IOSTAT=ErrStatTmp ) '# Time: ',SeaState_p%WaveField%WaveTime(I)
! Now output the X,Y, Elev info for this timestep
- do j=1,SeaState_p%NGrid(1)
- xpos = -SeaState_p%deltaGrid(1)*(SeaState_p%NGrid(1)-1)/2.0 + (J-1)*SeaState_p%deltaGrid(1)
+ do j=1,size(SeaStateInitOut%WaveElevVisX)
+ xpos = SeaStateInitOut%WaveElevVisX(j)
do k=1, SeaState_p%NGrid(2)
- ypos = -SeaState_p%deltaGrid(2)*(SeaState_p%NGrid(2)-1)/2.0 + (K-1)*SeaState_p%deltaGrid(2)
- if (allocated(SeaState_p%WaveField%WaveElev2)) then
- WaveElev = SeaState_p%WaveField%WaveElev1(I,J,K) + SeaState_p%WaveField%WaveElev2(I,J,K)
- else
- WaveElev = SeaState_p%WaveField%WaveElev1(I,J,K)
- end if
+ ypos = SeaStateInitOut%WaveElevVisY(k)
+ WaveElev = SeaStateInitOut%WaveElevVisGrid(i,j,k)
write (WaveElevFileUn,WaveElevFmt, IOSTAT=ErrStatTmp ) xpos, ypos, WaveElev
end do
end do
diff --git a/modules/seastate/src/SeaState_Types.f90 b/modules/seastate/src/SeaState_Types.f90
index ff4581d486..f2807fdb04 100644
--- a/modules/seastate/src/SeaState_Types.f90
+++ b/modules/seastate/src/SeaState_Types.f90
@@ -94,13 +94,15 @@ MODULE SeaState_Types
REAL(ReKi) :: defWtrDpth = 0.0_ReKi !< Default water depth from the driver; may be overwritten [m]
REAL(ReKi) :: defMSL2SWL = 0.0_ReKi !< Default mean sea level to still water level from the driver; may be overwritten [m]
REAL(DbKi) :: TMax = 0.0_R8Ki !< Supplied by Driver: The total simulation time [(sec)]
- REAL(SiKi) , DIMENSION(:,:), ALLOCATABLE :: WaveElevXY !< Supplied by Driver: X-Y locations for WaveElevation output (for visualization). First dimension is the X (1) and Y (2) coordinate. Second dimension is the point number. [m,-]
INTEGER(IntKi) :: WaveFieldMod = 0_IntKi !< Wave field handling (-) (switch) 0: use individual SeaState inputs without adjustment, 1: adjust wave phases based on turbine offsets from farm origin [-]
REAL(ReKi) :: PtfmLocationX = 0.0_ReKi !< Supplied by Driver: X coordinate of platform location in the wave field [m]
REAL(ReKi) :: PtfmLocationY = 0.0_ReKi !< Supplied by Driver: Y coordinate of platform location in the wave field [m]
INTEGER(IntKi) :: WrWvKinMod = 0 !< 0,1, or 2 indicating whether we are going to write out kinematics files. [ignored if WaveMod = 6, if 1 or 2 then files are written using the outrootname] [-]
LOGICAL :: HasIce = .false. !< Supplied by Driver: Whether this simulation has ice loading (flag) [-]
LOGICAL :: Linearize = .FALSE. !< Flag that tells this module if the glue code wants to linearize. [-]
+ LOGICAL :: SurfaceVis = .FALSE. !< Turn on grid surface visualization outputs [-]
+ INTEGER(IntKi) :: SurfaceVisNx = 0 !< Number of points in X direction to output for visualization grid. Use 0 or negative to set to SeaState resolution. [-]
+ INTEGER(IntKi) :: SurfaceVisNy = 0 !< Number of points in Y direction to output for visualization grid. Use 0 or negative to set to SeaState resolution. [-]
END TYPE SeaSt_InitInputType
! =======================
! ========= SeaSt_InitOutputType =======
@@ -109,7 +111,9 @@ MODULE SeaState_Types
CHARACTER(ChanLen) , DIMENSION(:), ALLOCATABLE :: WriteOutputUnt !< The is the list of all HD-related output channel unit strings (includes all sub-module channels) [-]
TYPE(ProgDesc) :: Ver !< Version of SeaState [-]
LOGICAL :: InvalidWithSSExctn = .false. !< Whether SeaState configuration is invalid with HydroDyn's state-space excitation (ExctnMod=2) [(-)]
- REAL(SiKi) , DIMENSION(:,:), ALLOCATABLE :: WaveElevSeries !< Wave elevation time-series at each of the points given by WaveElevXY. First dimension is the timestep. Second dimension is XY point number corresponding to second dimension of WaveElevXY. [(m)]
+ REAL(SiKi) , DIMENSION(:), ALLOCATABLE :: WaveElevVisX !< X locations of grid output [m,-]
+ REAL(SiKi) , DIMENSION(:), ALLOCATABLE :: WaveElevVisY !< X locations of grid output [m,-]
+ REAL(SiKi) , DIMENSION(:,:,:), ALLOCATABLE :: WaveElevVisGrid !< Wave elevation time-series at each of the points given by WaveElevXY. First dimension is the timestep. Second/third dimensions are the grid of points. [(m)]
TYPE(SeaSt_WaveFieldType) , POINTER :: WaveField => NULL() !< Pointer to wave field [-]
END TYPE SeaSt_InitOutputType
! =======================
@@ -445,7 +449,6 @@ subroutine SeaSt_CopyInitInput(SrcInitInputData, DstInitInputData, CtrlCode, Err
integer(IntKi), intent(in ) :: CtrlCode
integer(IntKi), intent( out) :: ErrStat
character(*), intent( out) :: ErrMsg
- integer(B8Ki) :: LB(2), UB(2)
integer(IntKi) :: ErrStat2
character(ErrMsgLen) :: ErrMsg2
character(*), parameter :: RoutineName = 'SeaSt_CopyInitInput'
@@ -462,24 +465,15 @@ subroutine SeaSt_CopyInitInput(SrcInitInputData, DstInitInputData, CtrlCode, Err
DstInitInputData%defWtrDpth = SrcInitInputData%defWtrDpth
DstInitInputData%defMSL2SWL = SrcInitInputData%defMSL2SWL
DstInitInputData%TMax = SrcInitInputData%TMax
- if (allocated(SrcInitInputData%WaveElevXY)) then
- LB(1:2) = lbound(SrcInitInputData%WaveElevXY, kind=B8Ki)
- UB(1:2) = ubound(SrcInitInputData%WaveElevXY, kind=B8Ki)
- if (.not. allocated(DstInitInputData%WaveElevXY)) then
- allocate(DstInitInputData%WaveElevXY(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2)
- if (ErrStat2 /= 0) then
- call SetErrStat(ErrID_Fatal, 'Error allocating DstInitInputData%WaveElevXY.', ErrStat, ErrMsg, RoutineName)
- return
- end if
- end if
- DstInitInputData%WaveElevXY = SrcInitInputData%WaveElevXY
- end if
DstInitInputData%WaveFieldMod = SrcInitInputData%WaveFieldMod
DstInitInputData%PtfmLocationX = SrcInitInputData%PtfmLocationX
DstInitInputData%PtfmLocationY = SrcInitInputData%PtfmLocationY
DstInitInputData%WrWvKinMod = SrcInitInputData%WrWvKinMod
DstInitInputData%HasIce = SrcInitInputData%HasIce
DstInitInputData%Linearize = SrcInitInputData%Linearize
+ DstInitInputData%SurfaceVis = SrcInitInputData%SurfaceVis
+ DstInitInputData%SurfaceVisNx = SrcInitInputData%SurfaceVisNx
+ DstInitInputData%SurfaceVisNy = SrcInitInputData%SurfaceVisNy
end subroutine
subroutine SeaSt_DestroyInitInput(InitInputData, ErrStat, ErrMsg)
@@ -493,9 +487,6 @@ subroutine SeaSt_DestroyInitInput(InitInputData, ErrStat, ErrMsg)
ErrMsg = ''
call NWTC_Library_DestroyFileInfoType(InitInputData%PassedFileData, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
- if (allocated(InitInputData%WaveElevXY)) then
- deallocate(InitInputData%WaveElevXY)
- end if
end subroutine
subroutine SeaSt_PackInitInput(RF, Indata)
@@ -512,13 +503,15 @@ subroutine SeaSt_PackInitInput(RF, Indata)
call RegPack(RF, InData%defWtrDpth)
call RegPack(RF, InData%defMSL2SWL)
call RegPack(RF, InData%TMax)
- call RegPackAlloc(RF, InData%WaveElevXY)
call RegPack(RF, InData%WaveFieldMod)
call RegPack(RF, InData%PtfmLocationX)
call RegPack(RF, InData%PtfmLocationY)
call RegPack(RF, InData%WrWvKinMod)
call RegPack(RF, InData%HasIce)
call RegPack(RF, InData%Linearize)
+ call RegPack(RF, InData%SurfaceVis)
+ call RegPack(RF, InData%SurfaceVisNx)
+ call RegPack(RF, InData%SurfaceVisNy)
if (RegCheckErr(RF, RoutineName)) return
end subroutine
@@ -526,9 +519,6 @@ subroutine SeaSt_UnPackInitInput(RF, OutData)
type(RegFile), intent(inout) :: RF
type(SeaSt_InitInputType), intent(inout) :: OutData
character(*), parameter :: RoutineName = 'SeaSt_UnPackInitInput'
- integer(B8Ki) :: LB(2), UB(2)
- integer(IntKi) :: stat
- logical :: IsAllocAssoc
if (RF%ErrStat /= ErrID_None) return
call RegUnpack(RF, OutData%InputFile); if (RegCheckErr(RF, RoutineName)) return
call RegUnpack(RF, OutData%UseInputFile); if (RegCheckErr(RF, RoutineName)) return
@@ -539,13 +529,15 @@ subroutine SeaSt_UnPackInitInput(RF, OutData)
call RegUnpack(RF, OutData%defWtrDpth); if (RegCheckErr(RF, RoutineName)) return
call RegUnpack(RF, OutData%defMSL2SWL); if (RegCheckErr(RF, RoutineName)) return
call RegUnpack(RF, OutData%TMax); if (RegCheckErr(RF, RoutineName)) return
- call RegUnpackAlloc(RF, OutData%WaveElevXY); if (RegCheckErr(RF, RoutineName)) return
call RegUnpack(RF, OutData%WaveFieldMod); if (RegCheckErr(RF, RoutineName)) return
call RegUnpack(RF, OutData%PtfmLocationX); if (RegCheckErr(RF, RoutineName)) return
call RegUnpack(RF, OutData%PtfmLocationY); if (RegCheckErr(RF, RoutineName)) return
call RegUnpack(RF, OutData%WrWvKinMod); if (RegCheckErr(RF, RoutineName)) return
call RegUnpack(RF, OutData%HasIce); if (RegCheckErr(RF, RoutineName)) return
call RegUnpack(RF, OutData%Linearize); if (RegCheckErr(RF, RoutineName)) return
+ call RegUnpack(RF, OutData%SurfaceVis); if (RegCheckErr(RF, RoutineName)) return
+ call RegUnpack(RF, OutData%SurfaceVisNx); if (RegCheckErr(RF, RoutineName)) return
+ call RegUnpack(RF, OutData%SurfaceVisNy); if (RegCheckErr(RF, RoutineName)) return
end subroutine
subroutine SeaSt_CopyInitOutput(SrcInitOutputData, DstInitOutputData, CtrlCode, ErrStat, ErrMsg)
@@ -554,7 +546,7 @@ subroutine SeaSt_CopyInitOutput(SrcInitOutputData, DstInitOutputData, CtrlCode,
integer(IntKi), intent(in ) :: CtrlCode
integer(IntKi), intent( out) :: ErrStat
character(*), intent( out) :: ErrMsg
- integer(B8Ki) :: LB(2), UB(2)
+ integer(B8Ki) :: LB(3), UB(3)
integer(IntKi) :: ErrStat2
character(ErrMsgLen) :: ErrMsg2
character(*), parameter :: RoutineName = 'SeaSt_CopyInitOutput'
@@ -588,17 +580,41 @@ subroutine SeaSt_CopyInitOutput(SrcInitOutputData, DstInitOutputData, CtrlCode,
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
if (ErrStat >= AbortErrLev) return
DstInitOutputData%InvalidWithSSExctn = SrcInitOutputData%InvalidWithSSExctn
- if (allocated(SrcInitOutputData%WaveElevSeries)) then
- LB(1:2) = lbound(SrcInitOutputData%WaveElevSeries, kind=B8Ki)
- UB(1:2) = ubound(SrcInitOutputData%WaveElevSeries, kind=B8Ki)
- if (.not. allocated(DstInitOutputData%WaveElevSeries)) then
- allocate(DstInitOutputData%WaveElevSeries(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2)
+ if (allocated(SrcInitOutputData%WaveElevVisX)) then
+ LB(1:1) = lbound(SrcInitOutputData%WaveElevVisX, kind=B8Ki)
+ UB(1:1) = ubound(SrcInitOutputData%WaveElevVisX, kind=B8Ki)
+ if (.not. allocated(DstInitOutputData%WaveElevVisX)) then
+ allocate(DstInitOutputData%WaveElevVisX(LB(1):UB(1)), stat=ErrStat2)
+ if (ErrStat2 /= 0) then
+ call SetErrStat(ErrID_Fatal, 'Error allocating DstInitOutputData%WaveElevVisX.', ErrStat, ErrMsg, RoutineName)
+ return
+ end if
+ end if
+ DstInitOutputData%WaveElevVisX = SrcInitOutputData%WaveElevVisX
+ end if
+ if (allocated(SrcInitOutputData%WaveElevVisY)) then
+ LB(1:1) = lbound(SrcInitOutputData%WaveElevVisY, kind=B8Ki)
+ UB(1:1) = ubound(SrcInitOutputData%WaveElevVisY, kind=B8Ki)
+ if (.not. allocated(DstInitOutputData%WaveElevVisY)) then
+ allocate(DstInitOutputData%WaveElevVisY(LB(1):UB(1)), stat=ErrStat2)
+ if (ErrStat2 /= 0) then
+ call SetErrStat(ErrID_Fatal, 'Error allocating DstInitOutputData%WaveElevVisY.', ErrStat, ErrMsg, RoutineName)
+ return
+ end if
+ end if
+ DstInitOutputData%WaveElevVisY = SrcInitOutputData%WaveElevVisY
+ end if
+ if (allocated(SrcInitOutputData%WaveElevVisGrid)) then
+ LB(1:3) = lbound(SrcInitOutputData%WaveElevVisGrid, kind=B8Ki)
+ UB(1:3) = ubound(SrcInitOutputData%WaveElevVisGrid, kind=B8Ki)
+ if (.not. allocated(DstInitOutputData%WaveElevVisGrid)) then
+ allocate(DstInitOutputData%WaveElevVisGrid(LB(1):UB(1),LB(2):UB(2),LB(3):UB(3)), stat=ErrStat2)
if (ErrStat2 /= 0) then
- call SetErrStat(ErrID_Fatal, 'Error allocating DstInitOutputData%WaveElevSeries.', ErrStat, ErrMsg, RoutineName)
+ call SetErrStat(ErrID_Fatal, 'Error allocating DstInitOutputData%WaveElevVisGrid.', ErrStat, ErrMsg, RoutineName)
return
end if
end if
- DstInitOutputData%WaveElevSeries = SrcInitOutputData%WaveElevSeries
+ DstInitOutputData%WaveElevVisGrid = SrcInitOutputData%WaveElevVisGrid
end if
DstInitOutputData%WaveField => SrcInitOutputData%WaveField
end subroutine
@@ -620,8 +636,14 @@ subroutine SeaSt_DestroyInitOutput(InitOutputData, ErrStat, ErrMsg)
end if
call NWTC_Library_DestroyProgDesc(InitOutputData%Ver, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
- if (allocated(InitOutputData%WaveElevSeries)) then
- deallocate(InitOutputData%WaveElevSeries)
+ if (allocated(InitOutputData%WaveElevVisX)) then
+ deallocate(InitOutputData%WaveElevVisX)
+ end if
+ if (allocated(InitOutputData%WaveElevVisY)) then
+ deallocate(InitOutputData%WaveElevVisY)
+ end if
+ if (allocated(InitOutputData%WaveElevVisGrid)) then
+ deallocate(InitOutputData%WaveElevVisGrid)
end if
nullify(InitOutputData%WaveField)
end subroutine
@@ -636,7 +658,9 @@ subroutine SeaSt_PackInitOutput(RF, Indata)
call RegPackAlloc(RF, InData%WriteOutputUnt)
call NWTC_Library_PackProgDesc(RF, InData%Ver)
call RegPack(RF, InData%InvalidWithSSExctn)
- call RegPackAlloc(RF, InData%WaveElevSeries)
+ call RegPackAlloc(RF, InData%WaveElevVisX)
+ call RegPackAlloc(RF, InData%WaveElevVisY)
+ call RegPackAlloc(RF, InData%WaveElevVisGrid)
call RegPack(RF, associated(InData%WaveField))
if (associated(InData%WaveField)) then
call RegPackPointer(RF, c_loc(InData%WaveField), PtrInIndex)
@@ -651,7 +675,7 @@ subroutine SeaSt_UnPackInitOutput(RF, OutData)
type(RegFile), intent(inout) :: RF
type(SeaSt_InitOutputType), intent(inout) :: OutData
character(*), parameter :: RoutineName = 'SeaSt_UnPackInitOutput'
- integer(B8Ki) :: LB(2), UB(2)
+ integer(B8Ki) :: LB(3), UB(3)
integer(IntKi) :: stat
logical :: IsAllocAssoc
integer(B8Ki) :: PtrIdx
@@ -661,7 +685,9 @@ subroutine SeaSt_UnPackInitOutput(RF, OutData)
call RegUnpackAlloc(RF, OutData%WriteOutputUnt); if (RegCheckErr(RF, RoutineName)) return
call NWTC_Library_UnpackProgDesc(RF, OutData%Ver) ! Ver
call RegUnpack(RF, OutData%InvalidWithSSExctn); if (RegCheckErr(RF, RoutineName)) return
- call RegUnpackAlloc(RF, OutData%WaveElevSeries); if (RegCheckErr(RF, RoutineName)) return
+ call RegUnpackAlloc(RF, OutData%WaveElevVisX); if (RegCheckErr(RF, RoutineName)) return
+ call RegUnpackAlloc(RF, OutData%WaveElevVisY); if (RegCheckErr(RF, RoutineName)) return
+ call RegUnpackAlloc(RF, OutData%WaveElevVisGrid); if (RegCheckErr(RF, RoutineName)) return
if (associated(OutData%WaveField)) deallocate(OutData%WaveField)
call RegUnpack(RF, IsAllocAssoc); if (RegCheckErr(RF, RoutineName)) return
if (IsAllocAssoc) then