IfW: add VTK slice in XY output to driver

This allows visualizing a slice through hub-height (or other height) when compairing with a simulation.

modules/inflowwind/src/InflowWind_Driver.f90

@@ -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 -=-=-
    !--------------------------------------------------------------------------------------------------------------------------------

modules/inflowwind/src/InflowWind_Driver_Subs.f90

@@ -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 )
@@ -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
@@ -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
 
 

modules/inflowwind/src/InflowWind_Driver_Types.f90

@@ -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
 
 
@@ -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
 
 

modules/inflowwind/src/InflowWind_IO.f90

@@ -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', '', '')
 
@@ -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
+   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)
+end subroutine Grid3D_WriteVTKsliceXY
+
 end module InflowWind_IO