Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

IfW: add VTK output of slice in XY to driver #1643

Closed
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 17 additions & 0 deletions modules/inflowwind/src/InflowWind_Driver.f90
Original file line number Diff line number Diff line change
Expand Up @@ -503,6 +503,23 @@ PROGRAM InflowWind_Driver
END IF
END IF


IF (SettingsFlags%XYslice) THEN
CALL IfW_WriteXYslice( InflowWind_p%FlowField, InflowWind_InitInp%RootName, Settings%XYslice_height, ErrStat, ErrMsg )
IF (ErrStat > ErrID_None) THEN
CALL WrScr( TRIM(ErrMsg) )
IF ( ErrStat >= AbortErrLev ) THEN
CALL DriverCleanup()
CALL ProgAbort( ErrMsg )
ELSEIF ( IfWDriver_Verbose >= 7_IntKi ) THEN
CALL WrScr(NewLine//' IfW_WriteXYslice returned: ErrStat: '//TRIM(Num2LStr(ErrStat)))
END IF
ELSE IF ( IfWDriver_Verbose >= 5_IntKi ) THEN
CALL WrScr(NewLine//'IfW_WriteXYslice CALL returned without errors.'//NewLine)
END IF
END IF


!--------------------------------------------------------------------------------------------------------------------------------
!-=-=- Other Setup -=-=-
!--------------------------------------------------------------------------------------------------------------------------------
Expand Down
81 changes: 81 additions & 0 deletions modules/inflowwind/src/InflowWind_Driver_Subs.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1308,6 +1308,40 @@ SUBROUTINE ReadDvrIptFile( DvrFileName, DvrFlags, DvrSettings, ProgInfo, ErrStat
ENDIF


!-------------------------------------------------------------------------------------------------
! XY slice output
!-------------------------------------------------------------------------------------------------

! Header
CALL ReadCom( UnIn, FileName,' XY slice output, comment line', ErrStatTmp, ErrMsgTmp, UnEchoLocal )
IF ( ErrStatTmp /= ErrID_None ) THEN
CALL SetErrStat(ErrID_Fatal,ErrMsgTmp,ErrStat,ErrMsg,RoutineName)
CALL CleanupEchoFile( EchoFileContents, UnEchoLocal )
CLOSE( UnIn )
RETURN
ENDIF


! XYslice -- Output a VTK slice in XY
CALL ReadVar( UnIn, FileName,DvrFlags%XYslice,'XYslice',' VTK slice in XY?', &
ErrStatTmp,ErrMsgTmp, UnEchoLocal )
IF ( ErrStatTmp /= ErrID_None ) THEN
CALL SetErrStat(ErrID_Fatal,ErrMsgTmp,ErrStat,ErrMsg,RoutineName)
CALL CleanupEchoFile( EchoFileContents, UnEchoLocal )
CLOSE( UnIn )
RETURN
ENDIF

! XYslice_height -- Height for XY slice
CALL ReadVar( UnIn, FileName,DvrSettings%XYslice_height,'XYslice_height',' VTK slice height', &
ErrStatTmp,ErrMsgTmp, UnEchoLocal )
IF ( ErrStatTmp /= ErrID_None ) THEN
CALL SetErrStat(ErrID_Fatal,ErrMsgTmp,ErrStat,ErrMsg,RoutineName)
CALL CleanupEchoFile( EchoFileContents, UnEchoLocal )
CLOSE( UnIn )
RETURN
ENDIF


! Close the echo and input file
CALL CleanupEchoFile( EchoFileContents, UnEchoLocal )
Expand Down Expand Up @@ -2623,6 +2657,49 @@ subroutine IfW_WriteVTK(FF, FileRootName, ErrStat, ErrMsg)
end subroutine IfW_WriteVTK


subroutine IfW_WriteXYslice(FF, FileRootName, XYslice_height, ErrStat, ErrMsg)
type(FlowFieldType), intent(in) :: FF !< Parameters
character(*), intent(in) :: FileRootName !< RootName for output files
real(ReKi), intent(in) :: XYslice_height
integer(IntKi), intent(out) :: ErrStat !< Error status of the operation
character(*), intent(out) :: ErrMsg !< Error message if ErrStat /= ErrID_None

character(*), parameter :: RoutineName = "IfW_WriteXYslice"
type(Grid3DFieldType) :: G3D
integer(IntKi) :: unit
integer(IntKi) :: ErrStat2
character(ErrMsgLen) :: ErrMsg2

ErrStat = ErrID_None
ErrMsg = ""

! Get new unit for writing file
call GetNewUnit(unit, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
if (ErrStat >= AbortErrLev) return

! Switch based on field type
select case (FF%FieldType)

case (Uniform_FieldType)
call Uniform_to_Grid3D(FF%Uniform, FF%VelInterpCubic, G3D, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
if (ErrStat < AbortErrLev) then
call Grid3D_WriteVTKsliceXY(G3D, FileRootName, XYslice_height, unit, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
end if

case (Grid3D_FieldType)
call Grid3D_WriteVTKsliceXY(FF%Grid3D, FileRootName, XYslice_height, unit, ErrStat, ErrMsg)

case default
ErrStat = ErrID_Warn
ErrMsg = RoutineName//': Field type '//TRIM(Num2LStr(FF%FieldType))// &
' cannot be converted to VTK format.'
end select
end subroutine IfW_WriteXYslice


!> This routine exists only to support the development of the module. It will not be needed after the module is complete.
SUBROUTINE printSettings( DvrFlags, DvrSettings )
! The arguments
Expand Down Expand Up @@ -2671,6 +2748,10 @@ SUBROUTINE printSettings( DvrFlags, DvrSettings )
CALL WrScr(' FFTOutputInit: '//FLAG(DvrSettings%FFTOutput%Initialized)// ' Unit #: '//TRIM(Num2LStr(DvrSettings%FFTOutput%Unit)))
CALL WrScr(' PointsVelOutputInit: '//FLAG(DvrSettings%PointsVelOutput%Initialized)// ' Unit #: '//TRIM(Num2LStr(DvrSettings%PointsVelOutput%Unit)))
CALL WrScr(' PointsAccOutputInit: '//FLAG(DvrSettings%PointsVelOutput%Initialized)// ' Unit #: '//TRIM(Num2LStr(DvrSettings%PointsVelOutput%Unit)))
CALL WrScr(' XYslice: '//FLAG(DvrFlags%XYslice))
if (DvrFlags%XYslice) then
CALL WrScr(' XYslice_height '//TRIM(Num2LStr(DvrSettings%XYslice_height)))
end if
END SUBROUTINE printSettings


Expand Down
3 changes: 3 additions & 0 deletions modules/inflowwind/src/InflowWind_Driver_Types.f90
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,8 @@ MODULE InflowWind_Driver_Types
LOGICAL :: WrBladed = .FALSE. !< Requested file conversion to Bladed format?
LOGICAL :: WrVTK = .FALSE. !< Requested file output as VTK?
LOGICAL :: WrUniform = .FALSE. !< Requested file output as Uniform wind format?

LOGICAL :: XYslice = .FALSE. !< Take XY slice at one elevation
END TYPE IfWDriver_Flags


Expand Down Expand Up @@ -105,6 +107,7 @@ MODULE InflowWind_Driver_Types
TYPE(OutputFile) :: FFTOutput
TYPE(OutputFile) :: PointsVelOutput

REAL(ReKi) :: XYslice_height = 0.0_ReKi !< height to take XY slice
END TYPE IfWDriver_Settings


Expand Down
78 changes: 77 additions & 1 deletion modules/inflowwind/src/InflowWind_IO.f90
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,8 @@ module InflowWind_IO
public :: Uniform_WriteHH, &
Grid3D_WriteBladed, &
Grid3D_WriteHAWC, &
Grid3D_WriteVTK
Grid3D_WriteVTK, &
Grid3D_WriteVTKsliceXY

type(ProgDesc), parameter :: InflowWind_IO_Ver = ProgDesc('InflowWind_IO', '', '')

Expand Down Expand Up @@ -2924,4 +2925,79 @@ subroutine Grid3D_WriteHAWC(G3D, FileRootName, unit, ErrStat, ErrMsg)

end subroutine Grid3D_WriteHAWC


subroutine Grid3D_WriteVTKsliceXY(G3D, FileRootName, XYslice_height, unit, ErrStat, ErrMsg)

type(Grid3DFieldType), intent(in) :: G3D !< Parameters
character(*), intent(in) :: FileRootName !< RootName for output files
real(ReKi), intent(in) :: XYslice_height
integer(IntKi), intent(in) :: unit !< Error status of the operation
integer(IntKi), intent(out) :: ErrStat !< Error status of the operation
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are ErrStat and ErrMsg initialized in this function? Maybe I missed it.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

They aren't explicitly initialized, but get set either by OpenFOutFile or by the check on the slice height. For completeness I can add them.

I also just noticed that there is no upper bounds check on the height, which could lead to an invalid index on iz. I'll fix that also.

character(*), intent(out) :: ErrMsg !< Error message if ErrStat /= ErrID_None

character(*), parameter :: RoutineName = 'Grid3D_WriteVTKsliceXY'
character(1024) :: RootPathName
character(1024) :: FileName
character(3) :: ht_str
integer :: it !< time index for slice
integer :: ix
integer :: iy
integer :: iz
real(ReKi) :: time !< time for this slice
real(ReKi) :: ht !< nearest grid slice elevation
integer(IntKi) :: ErrStat2
character(ErrMsgLen) :: ErrMsg2

call GetPath(FileRootName, RootPathName)
RootPathName = trim(RootPathName)//PathSep//"vtk"
call MkDir(trim(RootPathName)) ! make this directory if it doesn't already exist

! get indices for this slice
iz = nint((G3D%GridBase + XYslice_height)*G3D%InvDZ)
ht = real(iz,ReKi) / G3D%InvDZ + G3D%GridBase ! nearest height index
write(ht_str,'(i3)') nint(ht)

! check for errors in slice height
if (iz <= 0_IntKi) then
call SetErrStat(ErrID_Warn,"No grid points near XY slice height of "//trim(num2lstr(XYslice_height))//". Skipping writing slice file.",ErrStat,ErrMsg,RoutineName)
return
endif

! Loop through time steps
do it = 1, G3D%NSteps

time = real(it - 1, ReKi)*G3D%DTime

! Create the output vtk file with naming <WindFilePath>/vtk/DisYZ.t<i>.vtk
FileName = trim(RootPathName)//PathSep//"DisXY.Z"//ht_str//".t"//trim(num2lstr(it))//".vtp"

! see WrVTK_SP_header
call OpenFOutFile(unit, TRIM(FileName), ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
if (ErrStat >= AbortErrLev) return

write (unit, '(A)') '# vtk DataFile Version 3.0'
write (unit, '(A)') "InflowWind XY Slice at T= "//trim(num2lstr(time))//" s"
write (unit, '(A)') 'ASCII'
write (unit, '(A)') 'DATASET STRUCTURED_POINTS'

! Note: gridVals must be stored such that the left-most dimension is X
! and the right-most dimension is Z (see WrVTK_SP_vectors3D)
write (unit, '(A,3(i5,1X))') 'DIMENSIONS ', G3D%NSteps, G3D%NYGrids, 1
write (unit, '(A,3(f10.2,1X))') 'ORIGIN ', G3D%InitXPosition+time*G3D%MeanWS, -G3D%YHWid, ht
write (unit, '(A,3(f10.2,1X))') 'SPACING ', -G3D%Dtime*G3D%MeanWS, 1.0_ReKi/G3D%InvDY, 0.0_ReKi
write (unit, '(A,i9)') 'POINT_DATA ', G3D%NSteps*G3D%NYGrids
write (unit, '(A)') 'VECTORS DisXY float'

do iy = 1, G3D%NYGrids
do ix = 1, G3D%NSteps ! time and X are interchangeable
write (unit, '(3(f10.2,1X))') G3D%Vel(:, iy, iz, ix)
end do
end do

enddo

close (unit)
andrew-platt marked this conversation as resolved.
Show resolved Hide resolved
end subroutine Grid3D_WriteVTKsliceXY

end module InflowWind_IO