From 83c7e89062abc9cf76d40fcd1b4d2826fbd747b9 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Wed, 21 Jun 2023 14:17:32 -0600 Subject: [PATCH] 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 | 17 ++++ .../inflowwind/src/InflowWind_Driver_Subs.f90 | 81 +++++++++++++++++++ .../src/InflowWind_Driver_Types.f90 | 3 + modules/inflowwind/src/InflowWind_IO.f90 | 78 +++++++++++++++++- reg_tests/r-test | 2 +- 5 files changed, 179 insertions(+), 2 deletions(-) diff --git a/modules/inflowwind/src/InflowWind_Driver.f90 b/modules/inflowwind/src/InflowWind_Driver.f90 index 0a20b2e333..560cea6df5 100644 --- a/modules/inflowwind/src/InflowWind_Driver.f90 +++ b/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 -=-=- !-------------------------------------------------------------------------------------------------------------------------------- diff --git a/modules/inflowwind/src/InflowWind_Driver_Subs.f90 b/modules/inflowwind/src/InflowWind_Driver_Subs.f90 index d7033445b2..a89e13fe9e 100644 --- a/modules/inflowwind/src/InflowWind_Driver_Subs.f90 +++ b/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 diff --git a/modules/inflowwind/src/InflowWind_Driver_Types.f90 b/modules/inflowwind/src/InflowWind_Driver_Types.f90 index 669303cf60..c16a2ceb65 100644 --- a/modules/inflowwind/src/InflowWind_Driver_Types.f90 +++ b/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 diff --git a/modules/inflowwind/src/InflowWind_IO.f90 b/modules/inflowwind/src/InflowWind_IO.f90 index 96c2d7e21c..eae5f741b4 100644 --- a/modules/inflowwind/src/InflowWind_IO.f90 +++ b/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 /vtk/DisYZ.t.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 diff --git a/reg_tests/r-test b/reg_tests/r-test index 3e92b0bdc7..bca4292bb2 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 3e92b0bdc7b1a641649d94b92f99397c5669b9b3 +Subproject commit bca4292bb2c9e04e63f36d02fb5a8c88928c62f8