diff --git a/.gitignore b/.gitignore index 09166cd156..74be88339b 100644 --- a/.gitignore +++ b/.gitignore @@ -64,3 +64,8 @@ varcache # Python cache files openfast_io/dist/ openfast_io/openfast_io/_version.py + +# docker with VScode config +*.code-workspace +docker-compose.yml +.devcontainer/ diff --git a/docs/OtherSupporting/FAST.Farm/FAST.Farm_Plan_Rev26.docx b/docs/OtherSupporting/FAST.Farm/FAST.Farm_Plan_Rev26.docx new file mode 100644 index 0000000000..291484f998 Binary files /dev/null and b/docs/OtherSupporting/FAST.Farm/FAST.Farm_Plan_Rev26.docx differ diff --git a/docs/OtherSupporting/FAST.Farm/FAST.Farm_Plan_WakeExtentAndBuffering_Rev3.docx b/docs/OtherSupporting/FAST.Farm/FAST.Farm_Plan_WakeExtentAndBuffering_Rev3.docx new file mode 100644 index 0000000000..5bc7691131 Binary files /dev/null and b/docs/OtherSupporting/FAST.Farm/FAST.Farm_Plan_WakeExtentAndBuffering_Rev3.docx differ diff --git a/docs/OtherSupporting/FAST.Farm_Plan_Rev25.doc b/docs/OtherSupporting/FAST.Farm_Plan_Rev25.doc deleted file mode 100644 index 58ad2aad71..0000000000 Binary files a/docs/OtherSupporting/FAST.Farm_Plan_Rev25.doc and /dev/null differ diff --git a/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 b/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 index 5a2eabde5c..862e9ff6df 100644 --- a/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 +++ b/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 @@ -211,7 +211,11 @@ SUBROUTINE Farm_Initialize( farm, InputFile, ErrStat, ErrMsg ) call AllocAry( farm%p%MaxNumPlanes, farm%p%NumTurbines, 'farm%p%MaxNumPlanes', ErrStat2, ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName); if (Failed()) return do i=1,farm%p%NumTurbines ! Eventually, we will have different settings for different rotors - farm%p%MaxNumPlanes(i) = ceiling( 15.0 * ( WD_InitInput%InputFileData%NumDFull + WD_InitInput%InputFileData%NumDBuff ) / AWAE_InitInput%InputFileData%C_Meander ) + if (WD_InitInput%InputFileData%Mod_Wake == Mod_Wake_Polar) then + farm%p%MaxNumPlanes(i) = ceiling( 18.0 * ( WD_InitInput%InputFileData%NumDFull + WD_InitInput%InputFileData%NumDBuff ) / AWAE_InitInput%InputFileData%C_Meander ) + else + farm%p%MaxNumPlanes(i) = ceiling( 54.0 * ( WD_InitInput%InputFileData%NumDFull + WD_InitInput%InputFileData%NumDBuff ) ) + endif farm%p%MaxNumPlanes(i) = max( 2, min( farm%p%MaxNumPlanes(i) , farm%p%n_TMax + 2 ) ) end do diff --git a/modules/awae/CMakeLists.txt b/modules/awae/CMakeLists.txt index fdb66b91b6..289484c3ce 100644 --- a/modules/awae/CMakeLists.txt +++ b/modules/awae/CMakeLists.txt @@ -26,6 +26,7 @@ add_library(awaelib STATIC src/AWAE.f90 src/AWAE_IO.f90 src/AWAE_Types.f90 + src/AWAE_vtk.f90 src/amrex_utils.F90 ) target_link_libraries(awaelib awaelib_c ifwlib nwtclibs) diff --git a/modules/awae/src/AWAE.f90 b/modules/awae/src/AWAE.f90 index 0f8360cc54..e9f386c440 100644 --- a/modules/awae/src/AWAE.f90 +++ b/modules/awae/src/AWAE.f90 @@ -27,6 +27,7 @@ module AWAE use NWTC_Library use AWAE_Types use AWAE_IO + use AWAE_vtk use InflowWind use IfW_FlowField use KdTree @@ -59,6 +60,12 @@ module AWAE contains +!---------------------------------------------------------------------------------------------------------------------------------- +!> Extract a 2D slice from a 3D vector field `V` at the requested physical coordinate `s` along the slice-normal axis. +!! The slice orientation is selected by `sliceType` (XYSlice, YZSlice, or XZSlice). The location `s` (in meters) is +!! converted to grid units using the grid origin `s0` and spacing `ds`, the two bracketing grid planes are identified, +!! and the output `slice` is filled by linear interpolation between them. If `s` coincides with the last grid index, +!! the upper bracketing index is clamped so no out-of-bounds access occurs. subroutine ExtractSlice( sliceType, s, s0, szs, sz1, sz2, ds, V, slice) integer(IntKi), intent(in ) :: sliceType !< Type of slice: XYSlice, YZSlice, XZSlice @@ -106,8 +113,16 @@ subroutine ExtractSlice( sliceType, s, s0, szs, sz1, sz2, ds, V, slice) end subroutine ExtractSlice !---------------------------------------------------------------------------------------------------------------------------------- -!> This subroutine -!! +!> Precompute, for every pair of adjacent wake planes (np, np+1) of every turbine, the geometric quantities that +!! describe the relative orientation of the two planes. For each pair, this routine evaluates the cosine and sine of +!! the angle between the plane normals `u%xhat_plane(:,np,nt)` and `u%xhat_plane(:,np+1,nt)` and uses them, together +!! with the offset between the plane centers `u%p_plane`, to determine whether the planes are (numerically) parallel. +!! When they are not parallel, the routine computes and caches in the misc-var struct `m` the perpendicular distances +!! from each plane center to the line of intersection of the two planes (`r_s`, `r_e`), the in-plane unit vectors +!! pointing from that intersection line toward each plane center (`rhat_s`, `rhat_e`), and the closest points on the +!! intersection line to each plane center (`pvec_cs`, `pvec_ce`). The boolean `m%parallelFlag(np,nt)` records the +!! parallel/non-parallel decision. These cached quantities are reused downstream (e.g., in `interp_planes_2_point`) to +!! interpolate the wake-plane center and orientation between adjacent skewed wake planes. subroutine ComputeLocals(n, u, p, y, m, errStat, errMsg) integer(IntKi), intent(in ) :: n !< Current simulation time increment (zero-based) type(AWAE_InputType), intent(in ) :: u !< Inputs at Time t @@ -1172,7 +1187,7 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO integer(IntKi), intent( out) :: errStat !< Error status of the operation character(*), intent( out) :: errMsg !< Error message if errStat /= ErrID_None - character(1024) :: rootDir, baseName, OutFileVTKDir ! Simulation root dir, basename for outputs + character(1024) :: rootDir, baseName, OutFileVTKDir, OutFileVTKwakeDir ! Simulation root dir, basename for outputs integer(IntKi) :: i,j,nt,c ! loop counter real(ReKi) :: gridLoc ! Location of requested output slice in grid coordinates [0,sz-1] integer(IntKi) :: errStat2 ! temporary error status of the operation @@ -1253,13 +1268,33 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO ! --- Vtk Outputs call GetPath( p%OutFileRoot, rootDir, baseName ) - OutFileVTKDir = trim(rootDir) // 'vtk_ff' ! Directory for VTK outputs - p%OutFileVTKRoot = trim(rootDir) // 'vtk_ff' // PathSep // trim(baseName) ! Basename for VTK files + OutFileVTKDir = trim(rootDir) // 'vtk_ff' ! Directory for VTK outputs + p%OutFileFFvtkRoot = trim(OutFileVTKDir) // PathSep // trim(baseName) ! Basename for VTK files p%VTK_tWidth = CEILING( log10( real(p%NumDT, ReKi)/real(p%WrDisSkp1, ReKi) ) + 1) ! Length for time stamp if (p%WrDisWind .or. p%NOutDisWindXY>0 .or. p%NOutDisWindYZ>0 .or. p%NOutDisWindXZ>0) then - call MKDIR(OutFileVTKDir) ! creating output directory + call MKDIR(OutFileVTKDir) + ! placeholder for writing planes -- this will eventually be an input (revise logic here then) + p%WrPlanes = .true. end if + ! Setup wake plane writing + if (p%WrPlanes) then + OutFileVTKwakeDir = trim(OutFileVTKDir) // PathSep // 'wakes' ! Directory for VTK wake outputs + call MKDIR(OutFileVTKDir) ! we may not be writing out any other vtk, so create dir if doesn't exist + call MKDIR(OutFileVTKwakeDir) + p%OutFileFFvtkWakeRoot = trim(OutFileVTKwakeDir) // PathSep // trim(baseName) ! Basename for VTK wake files + p%VTK_tWidthPlanes = CEILING( log10(real(max(p%MaxPlanes, 1), ReKi)) + 1) ! length for the number of planes + ! Since a huge number of wake planes will be empty initially, we don't want to write all those out. + p%OutFileFFvtkWakeNullData = "FF.WakePlane_Null.vtk" + + ! track when a plane is first written out + allocate( m%WakeVTK_StartN(0:p%MaxPlanes-1,p%NumTurbines), stat=ErrStat2); if (Failed0('Could not allocate memory for m%WakeVTK_StartN')) return; + m%WakeVTK_StartN = huge(1_IntKi) + + ! write out the null wake plane + call Write_NullPlane(OutFileVTKwakeDir, p) + endif + ! Plane grids allocate( p%y(-p%Numradii+1:p%NumRadii-1), stat=errStat2); if (Failed0('Could not allocate memory for p%y.')) return; allocate( p%z(-p%Numradii+1:p%NumRadii-1), stat=errStat2); if (Failed0('Could not allocate memory for p%z.')) return; @@ -1428,10 +1463,14 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO allocate ( u%D_wake ( 0:p%MaxPlanes-1,1:p%NumTurbines), STAT=ErrStat2 ); if (Failed0('u%D_wake.' )) return; allocate ( u%WAT_k (1-p%NumRadii:p%NumRadii-1, 1-p%NumRadii:p%NumRadii-1, 0:p%MaxPlanes-1,1:p%NumTurbines), STAT=ErrStat2 ); if (Failed0('u%WAT_k.' )) return; - u%NumPlanes = 2.0_ReKi - u%Vx_wake=0.0_ReKi - u%Vy_wake=0.0_ReKi - u%Vz_wake=0.0_ReKi + u%NumPlanes = 2.0_ReKi + u%xhat_plane = 0.0_ReKi + u%p_plane = 0.0_ReKi + u%Vx_wake = 0.0_ReKi + u%Vy_wake = 0.0_ReKi + u%Vz_wake = 0.0_ReKi + u%D_wake = 0.0_ReKi + u%WAT_k = 0.0_ReKi !---------------------------------------------------------------------------- @@ -1471,12 +1510,15 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO if ( p%NOutDisWindXY > 0 ) then ALLOCATE ( m%OutVizXYPlane(3,p%LowRes%nXYZ(1), p%LowRes%nXYZ(2),1) , STAT=ErrStat2 ); if (Failed0('the Fast.Farm OutVizXYPlane arrays.')) return; + m%OutVizXYPlane = 0.0_SiKi end if if ( p%NOutDisWindYZ > 0 ) then ALLOCATE ( m%OutVizYZPlane(3,p%LowRes%nXYZ(2), p%LowRes%nXYZ(3),1) , STAT=ErrStat2 ); if (Failed0('the Fast.Farm OutVizYZPlane arrays.')) return; + m%OutVizYZPlane = 0.0_SiKi end if if ( p%NOutDisWindXZ > 0 ) then ALLOCATE ( m%OutVizXZPlane(3,p%LowRes%nXYZ(1), p%LowRes%nXYZ(3),1) , STAT=ErrStat2 ); if (Failed0('the Fast.Farm OutVizXZPlane arrays.')) return; + m%OutVizXZPlane = 0.0_SiKi end if ! miscvars to avoid the allocation per timestep @@ -1484,20 +1526,25 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO allocate(m%Vamb_low( 3, 0:p%LowRes%nXYZ(1)-1 , 0:p%LowRes%nXYZ(2)-1 , 0:p%LowRes%nXYZ(3)-1 ), STAT=errStat2); if (Failed0('m%Vamb_low.' )) return; allocate(m%Vdist_low( 3, 0:p%LowRes%nXYZ(1)-1 , 0:p%LowRes%nXYZ(2)-1 , 0:p%LowRes%nXYZ(3)-1 ), STAT=errStat2); if (Failed0('m%Vdist_low.' )) return; allocate(m%Vdist_low_full( 3, 0:p%LowRes%nXYZ(1)-1 , 0:p%LowRes%nXYZ(2)-1 , 0:p%LowRes%nXYZ(3)-1 ), STAT=errStat2); if (Failed0('m%Vdist_low_full')) return; + m%Vamb_lowpol = 0.0_ReKi + m%Vamb_low = 0.0_SiKi + m%Vdist_low = 0.0_SiKi + m%Vdist_low_full = 0.0_SiKi allocate(m%Vamb_high(1:p%NumTurbines), STAT=ErrStat2); if (Failed0('Could not allocate memory for m%Vamb_high.')) return; do nt = 1, p%NumTurbines allocate(m%Vamb_high(nt)%data(3,0:p%HighRes(nt)%nXYZ(1)-1, 0:p%HighRes(nt)%nXYZ(2)-1, 0:p%HighRes(nt)%nXYZ(3)-1, 0:p%n_high_low_p1), STAT=ErrStat2) if (Failed0('m%Vamb_high%data.')) return; + m%Vamb_high(nt)%data = 0.0_SiKi end do - allocate(m%parallelFlag( 0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%parallelFlag.')) return; - allocate(m%r_s( 0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%r_s.' )) return; - allocate(m%r_e( 0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%r_e.' )) return; - allocate(m%rhat_s( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%rhat_s.' )) return; - allocate(m%rhat_e( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%rhat_e.' )) return; - allocate(m%pvec_cs( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%pvec_cs.' )) return; - allocate(m%pvec_ce( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%pvec_ce.' )) return; + allocate(m%parallelFlag( 0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%parallelFlag.')) return; m%parallelFlag = .false. + allocate(m%r_s( 0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%r_s.' )) return; m%r_s = 0.0_ReKi + allocate(m%r_e( 0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%r_e.' )) return; m%r_e = 0.0_ReKi + allocate(m%rhat_s( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%rhat_s.' )) return; m%rhat_s = 0.0_ReKi + allocate(m%rhat_e( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%rhat_e.' )) return; m%rhat_e = 0.0_ReKi + allocate(m%pvec_cs( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%pvec_cs.' )) return; m%pvec_cs = 0.0_ReKi + allocate(m%pvec_ce( 3,0:p%MaxPlanes-2,1:p%NumTurbines ), STAT=errStat2); if (Failed0('m%pvec_ce.' )) return; m%pvec_ce = 0.0_ReKi ! WAT - store array of disk average velocities for all turbines call AllocAry(m%V_amb_low_disk,3,p%NumTurbines,'m%V_amb_low_disk', ErrStat2, ErrMsg2); if(Failed()) return; @@ -1660,6 +1707,26 @@ subroutine AWAE_End( u, p, x, xd, z, OtherState, y, m, errStat, errMsg ) errStat = ErrID_None errMsg = "" + ! Write .vtk.series files for disturbed-wind output slices + do nt = 1, p%NOutDisWindXY + if (.not. p%OutDisWindZvalid(nt)) cycle + call Write_DisWind_Series(p, "DisXY", nt) + end do + do nt = 1, p%NOutDisWindYZ + if (.not. p%OutDisWindXvalid(nt)) cycle + call Write_DisWind_Series(p, "DisYZ", nt) + end do + do nt = 1, p%NOutDisWindXZ + if (.not. p%OutDisWindYvalid(nt)) cycle + call Write_DisWind_Series(p, "DisXZ", nt) + end do + + if (p%WrPlanes) then + ! Write final ParaView .vtk.series files for wake planes + call Write_WakePlane_Series(p, m) + call Write_WireFrame_Series(p) + endif + ! Destroy InflowWind data select case(p%Mod_AmbWind) case (2) @@ -2004,7 +2071,7 @@ subroutine AWAE_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg call ExtractSlice(XYSlice, p%OutDisWindZ(k), p%LowRes%oXYZ(3), p%LowRes%nXYZ(3), p%LowRes%nXYZ(1), p%LowRes%nXYZ(2), p%LowRes%dXYZ(3), m%Vdist_low_full, m%outVizXYPlane(:,:,:,1)) ! Create the output vtk file with naming /Low/DisXY.t.vtk - FileName = trim(p%OutFileVTKRoot)//".Low.DisXY"//PlaneNumStr//"."//trim(Tstr)//".vtk" + FileName = trim(p%OutFileFFvtkRoot)//".Low.DisXY"//PlaneNumStr//"."//trim(Tstr)//".vtk" call WrVTK_SP_header(FileName, "Low resolution, disturbed wind of XY Slice at time = "//trim(num2lstr(t))//" seconds.", Un, ErrStat2, ErrMsg2 ); if (Failed()) return; call WrVTK_SP_vectors3D(Un, "Velocity", & [p%LowRes%nXYZ(1), p%LowRes%nXYZ(2), 1_IntKi], & @@ -2021,7 +2088,7 @@ subroutine AWAE_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg call ExtractSlice(YZSlice, p%OutDisWindX(k), p%LowRes%oXYZ(1), p%LowRes%nXYZ(1), p%LowRes%nXYZ(2), p%LowRes%nXYZ(3), p%LowRes%dXYZ(1), m%Vdist_low_full, m%outVizYZPlane(:,:,:,1)) ! Create the output vtk file with naming /Low/DisYZ.t.vtk - FileName = trim(p%OutFileVTKRoot)//".Low.DisYZ"//PlaneNumStr//"."//trim(Tstr)//".vtk" + FileName = trim(p%OutFileFFvtkRoot)//".Low.DisYZ"//PlaneNumStr//"."//trim(Tstr)//".vtk" call WrVTK_SP_header(FileName, "Low resolution, disturbed wind of YZ Slice at time = "//trim(num2lstr(t))//" seconds.", Un, ErrStat2, ErrMsg2 ); if (Failed()) return; call WrVTK_SP_vectors3D(Un, "Velocity", & [1, p%LowRes%nXYZ(2), p%LowRes%nXYZ(3)], & @@ -2038,7 +2105,7 @@ subroutine AWAE_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg call ExtractSlice(XZSlice, p%OutDisWindY(k), p%LowRes%oXYZ(2), p%LowRes%nXYZ(2), p%LowRes%nXYZ(1), p%LowRes%nXYZ(3), p%LowRes%dXYZ(2), m%Vdist_low_full, m%outVizXZPlane(:,:,:,1)) ! Create the output vtk file with naming /Low/DisXZ.t.vtk - FileName = trim(p%OutFileVTKRoot)//".Low.DisXZ"//PlaneNumStr//"."//trim(Tstr)//".vtk" + FileName = trim(p%OutFileFFvtkRoot)//".Low.DisXZ"//PlaneNumStr//"."//trim(Tstr)//".vtk" call WrVTK_SP_header(FileName, "Low resolution, disturbed wind of XZ Slice at time = "//trim(num2lstr(t))//" seconds.", Un, ErrStat2, ErrMsg2); if (Failed()) return; call WrVTK_SP_vectors3D(Un, "Velocity", & [p%LowRes%nXYZ(1), 1, p%LowRes%nXYZ(3)], & @@ -2046,10 +2113,27 @@ subroutine AWAE_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg p%LowRes%dXYZ, m%outVizXZPlane, ErrStat2, ErrMsg2) if (Failed()) return end do + + if (p%WrPlanes) then + !------------------------------------------------------------------------- + ! Write a VTK polydata file containing the four corners of every active + ! wake plane for every turbine as a set of quads. + !------------------------------------------------------------------------- + call Write_Planes_WireFrame(p, u, t, Tstr) + + !------------------------------------------------------------------------- + ! Write one VTK STRUCTURED_GRID file per wake plane (per turbine) with + ! the wake velocity sampled on the plane's structured Y-Z grid, plus a + ! ParaView .vtk.series JSON index of all wake-plane files written so far. + !------------------------------------------------------------------------- + call Write_Planes_Data(p, u, m, n, t, Tstr) + endif + end if contains + logical function Failed() call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) Failed = ErrStat >= AbortErrLev @@ -2581,5 +2665,4 @@ real(ReKi) function testf(y,z) end function end subroutine - end module AWAE diff --git a/modules/awae/src/AWAE_IO.f90 b/modules/awae/src/AWAE_IO.f90 index 663c9e2dc5..577943aaaa 100644 --- a/modules/awae/src/AWAE_IO.f90 +++ b/modules/awae/src/AWAE_IO.f90 @@ -80,7 +80,7 @@ subroutine WriteDisWindFiles( n, WrDisSkp1, p, y, m, errStat, errMsg ) ! TimeStamp write(Tstr, '(i' // trim(Num2LStr(p%VTK_tWidth)) //'.'// trim(Num2LStr(p%VTK_tWidth)) // ')') n_out ! TODO use n instead.. - FileName = trim(p%OutFileVTKRoot)//".Low.Dis."//trim(Tstr)//".vtk" + FileName = trim(p%OutFileFFvtkRoot)//".Low.Dis."//trim(Tstr)//".vtk" call WrVTK_SP_header( FileName, "Low resolution disturbed wind for time = "//trim(num2lstr(t_out))//" seconds.", Un, errStat2, errMsg2 ) call SetErrStat(errStat2, errMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return @@ -92,7 +92,7 @@ subroutine WriteDisWindFiles( n, WrDisSkp1, p, y, m, errStat, errMsg ) ! We are only writing out the first of the high res data for a given low res time step ! NOTE: y%Vdist_high(nt)%data(:,:,:,:,1) is at T=t_low, and index 0 is at T=t_low-DT_high - FileName = trim(p%OutFileVTKRoot)//".HighT"//trim(num2lstr(nt))//".Dis."//trim(Tstr)//".vtk" + FileName = trim(p%OutFileFFvtkRoot)//".HighT"//trim(num2lstr(nt))//".Dis."//trim(Tstr)//".vtk" call WrVTK_SP_header( FileName, "High resolution disturbed wind for time = "//trim(num2lstr(t_out))//" seconds.", Un, errStat2, errMsg2 ) call SetErrStat(errStat2, errMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return diff --git a/modules/awae/src/AWAE_Registry.txt b/modules/awae/src/AWAE_Registry.txt index 7a52d7e940..f484e120cc 100644 --- a/modules/awae/src/AWAE_Registry.txt +++ b/modules/awae/src/AWAE_Registry.txt @@ -127,13 +127,13 @@ typedef ^ MiscVarType IntKi iPlaneTurbTurb {:}{:}{:} - - "Fir typedef ^ MiscVarType IntKi iPlaneTurbChunk {:}{:}{:} - - "First and Last plane index by source turbine and destination chunk index" - typedef ^ MiscVarType Logical LowResChunkHasWake {:} - - "Low-res gridFirst and Last plane index by source turbine and destination chunk index" - typedef ^ MiscVarType ReKi MaxWakePointSep - - - "Maximum separation between wake points" - -typedef ^ MiscVarType Logical parallelFlag {:}{:} - - "" - -typedef ^ MiscVarType ReKi r_s {:}{:} - - "" - -typedef ^ MiscVarType ReKi r_e {:}{:} - - "" - -typedef ^ MiscVarType ReKi rhat_s {:}{:}{:} - - "" - -typedef ^ MiscVarType ReKi rhat_e {:}{:}{:} - - "" - -typedef ^ MiscVarType ReKi pvec_cs {:}{:}{:} - - "" - -typedef ^ MiscVarType ReKi pvec_ce {:}{:}{:} - - "" - +typedef ^ MiscVarType Logical parallelFlag {:}{:} - - "Flag indicating whether adjacent wake planes np and np+1 are (numerically) parallel for each turbine; dims: (plane-pair index np, turbine index nt)" - +typedef ^ MiscVarType ReKi r_s {:}{:} - - "Perpendicular distance from the start wake-plane center p_plane(:,np,nt) to the intersection line of adjacent (non-parallel) wake planes np and np+1; dims: (plane-pair index np, turbine index nt)" m +typedef ^ MiscVarType ReKi r_e {:}{:} - - "Perpendicular distance from the end wake-plane center p_plane(:,np+1,nt) to the intersection line of adjacent (non-parallel) wake planes np and np+1; dims: (plane-pair index np, turbine index nt)" m +typedef ^ MiscVarType ReKi rhat_s {:}{:}{:} - - "Unit vector lying in the start wake plane (np) pointing from the plane-plane intersection line toward the start plane center; dims: (XYZ component, plane-pair index np, turbine index nt)" - +typedef ^ MiscVarType ReKi rhat_e {:}{:}{:} - - "Unit vector lying in the end wake plane (np+1) pointing from the plane-plane intersection line toward the end plane center; dims: (XYZ component, plane-pair index np, turbine index nt)" - +typedef ^ MiscVarType ReKi pvec_cs {:}{:}{:} - - "Closest point on the plane-plane intersection line to the start wake-plane center, p_plane(:,np,nt) - r_s*rhat_s; dims: (XYZ component, plane-pair index np, turbine index nt)" m +typedef ^ MiscVarType ReKi pvec_ce {:}{:}{:} - - "Closest point on the plane-plane intersection line to the end wake-plane center, p_plane(:,np+1,nt) - r_e*rhat_e; dims: (XYZ component, plane-pair index np, turbine index nt)" m # typedef ^ MiscVarType SiKi outVizXYPlane {:}{:}{:}{:} - - "An array holding the output data for a 2D visualization slice" - typedef ^ MiscVarType SiKi outVizYZPlane {:}{:}{:}{:} - - "An array holding the output data for a 2D visualization slice" - @@ -147,6 +147,9 @@ typedef ^ MiscVarType InflowWind_OutputType y_IfW_High {:} - - "InflowWi typedef ^ MiscVarType ReKi V_amb_low_disk {:}{:} - - "Rotor averaged ambiend wind speed for each wind turbine (3 x nWT)" m/s typedef ^ MiscVarType IntKi planeDomainExit {:}{:} 0 - "Value indicates edge number (0: still in domain, +/-1: +/-X, +/-2: +/-Y, +/-3: +/-Z) the plane crossed" - +#visualization +typedef ^ MiscVarType IntKi WakeVTK_StartN {:}{:} - - "Time step when wake plane starts - counted by N_dtLow. Indices [wakenum,turbnum]" - + # Low-resolution grid chunk data typedef ^ LRGChunkType IntKi iChunk {3} - - "XYZ index of chunk" - typedef ^ LRGChunkType IntKi iSubGridX {2} - - "start and end grid indices in X direction" - @@ -225,8 +228,12 @@ typedef ^ ParameterType IntKi NOutDisWindXZ - - - "Number typedef ^ ParameterType ReKi OutDisWindY {:} - - "Y coordinates of XZ planes for output of disturbed wind data across the low-resolution domain [1 to NOutDisWindXZ]" meters typedef ^ ParameterType LOGICAL OutDisWindYvalid {:} - - "Valid XZ planes for output of disturbed wind data across the low-resolution domain [1 to NOutDisWindXZ]" - typedef ^ ParameterType CHARACTER(1024) OutFileRoot - - - "The root name derived from the primary FAST.Farm input file" - -typedef ^ ParameterType CHARACTER(1024) OutFileVTKRoot - - - "The root name for VTK outputs" - +typedef ^ ParameterType CHARACTER(1024) OutFileFFvtkRoot - - - "The root name for VTK outputs" - +typedef ^ ParameterType CHARACTER(1024) OutFileFFvtkWakeRoot - - - "The root name for VTK outputs for wake planes" - +typedef ^ ParameterType CHARACTER(1024) OutFileFFvtkWakeNullData - - - "Null data for an unpopulated wake plane" - typedef ^ ParameterType IntKi VTK_tWidth - - - "Number of characters for VTK timestamp outputs" - +typedef ^ ParameterType IntKi VTK_tWidthPlanes - 0 - "Number of charactes for the VTK plane numbers" - +typedef ^ ParameterType LOGICAL WrPlanes - .false. - "Write plane data out" - #typedef ^ ParameterType OutParmType OutParam {:} - - "Names and units (and other characteristics) of all requested output parameters" - #wake added turbulence typedef ^ ParameterType Logical WAT_Enabled - - - "Switch for turning on and off wake-added turbulence" - diff --git a/modules/awae/src/AWAE_Types.f90 b/modules/awae/src/AWAE_Types.f90 index ad08ebfd35..b17d3b1a0c 100644 --- a/modules/awae/src/AWAE_Types.f90 +++ b/modules/awae/src/AWAE_Types.f90 @@ -154,13 +154,13 @@ MODULE AWAE_Types INTEGER(IntKi) , DIMENSION(:,:,:), ALLOCATABLE :: iPlaneTurbChunk !< First and Last plane index by source turbine and destination chunk index [-] LOGICAL , DIMENSION(:), ALLOCATABLE :: LowResChunkHasWake !< Low-res gridFirst and Last plane index by source turbine and destination chunk index [-] REAL(ReKi) :: MaxWakePointSep = 0.0_ReKi !< Maximum separation between wake points [-] - LOGICAL , DIMENSION(:,:), ALLOCATABLE :: parallelFlag !< [-] - REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: r_s !< [-] - REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: r_e !< [-] - REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: rhat_s !< [-] - REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: rhat_e !< [-] - REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: pvec_cs !< [-] - REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: pvec_ce !< [-] + LOGICAL , DIMENSION(:,:), ALLOCATABLE :: parallelFlag !< Flag indicating whether adjacent wake planes np and np+1 are (numerically) parallel for each turbine; dims: (plane-pair index np, turbine index nt) [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: r_s !< Perpendicular distance from the start wake-plane center p_plane(:,np,nt) to the intersection line of adjacent (non-parallel) wake planes np and np+1; dims: (plane-pair index np, turbine index nt) [m] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: r_e !< Perpendicular distance from the end wake-plane center p_plane(:,np+1,nt) to the intersection line of adjacent (non-parallel) wake planes np and np+1; dims: (plane-pair index np, turbine index nt) [m] + REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: rhat_s !< Unit vector lying in the start wake plane (np) pointing from the plane-plane intersection line toward the start plane center; dims: (XYZ component, plane-pair index np, turbine index nt) [-] + REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: rhat_e !< Unit vector lying in the end wake plane (np+1) pointing from the plane-plane intersection line toward the end plane center; dims: (XYZ component, plane-pair index np, turbine index nt) [-] + REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: pvec_cs !< Closest point on the plane-plane intersection line to the start wake-plane center, p_plane(:,np,nt) - r_s*rhat_s; dims: (XYZ component, plane-pair index np, turbine index nt) [m] + REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: pvec_ce !< Closest point on the plane-plane intersection line to the end wake-plane center, p_plane(:,np+1,nt) - r_e*rhat_e; dims: (XYZ component, plane-pair index np, turbine index nt) [m] REAL(SiKi) , DIMENSION(:,:,:,:), ALLOCATABLE :: outVizXYPlane !< An array holding the output data for a 2D visualization slice [-] REAL(SiKi) , DIMENSION(:,:,:,:), ALLOCATABLE :: outVizYZPlane !< An array holding the output data for a 2D visualization slice [-] REAL(SiKi) , DIMENSION(:,:,:,:), ALLOCATABLE :: outVizXZPlane !< An array holding the output data for a 2D visualization slice [-] @@ -170,6 +170,7 @@ MODULE AWAE_Types TYPE(InflowWind_OutputType) , DIMENSION(:), ALLOCATABLE :: y_IfW_High !< InflowWind module outputs for the high-resolution grid [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: V_amb_low_disk !< Rotor averaged ambiend wind speed for each wind turbine (3 x nWT) [m/s] INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: planeDomainExit !< Value indicates edge number (0: still in domain, +/-1: +/-X, +/-2: +/-Y, +/-3: +/-Z) the plane crossed [-] + INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: WakeVTK_StartN !< Time step when wake plane starts - counted by N_dtLow. Indices [wakenum,turbnum] [-] END TYPE AWAE_MiscVarType ! ======================= ! ========= LRGChunkType ======= @@ -251,8 +252,12 @@ MODULE AWAE_Types REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: OutDisWindY !< Y coordinates of XZ planes for output of disturbed wind data across the low-resolution domain [1 to NOutDisWindXZ] [meters] LOGICAL , DIMENSION(:), ALLOCATABLE :: OutDisWindYvalid !< Valid XZ planes for output of disturbed wind data across the low-resolution domain [1 to NOutDisWindXZ] [-] CHARACTER(1024) :: OutFileRoot !< The root name derived from the primary FAST.Farm input file [-] - CHARACTER(1024) :: OutFileVTKRoot !< The root name for VTK outputs [-] + CHARACTER(1024) :: OutFileFFvtkRoot !< The root name for VTK outputs [-] + CHARACTER(1024) :: OutFileFFvtkWakeRoot !< The root name for VTK outputs for wake planes [-] + CHARACTER(1024) :: OutFileFFvtkWakeNullData !< Null data for an unpopulated wake plane [-] INTEGER(IntKi) :: VTK_tWidth = 0_IntKi !< Number of characters for VTK timestamp outputs [-] + INTEGER(IntKi) :: VTK_tWidthPlanes = 0 !< Number of charactes for the VTK plane numbers [-] + LOGICAL :: WrPlanes = .false. !< Write plane data out [-] LOGICAL :: WAT_Enabled = .false. !< Switch for turning on and off wake-added turbulence [-] TYPE(FlowFieldType) , POINTER :: WAT_FlowField => NULL() !< Pointer to the InflowWinds flow field data type [-] END TYPE AWAE_ParameterType @@ -1452,6 +1457,18 @@ subroutine AWAE_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) end if DstMiscData%planeDomainExit = SrcMiscData%planeDomainExit end if + if (allocated(SrcMiscData%WakeVTK_StartN)) then + LB(1:2) = lbound(SrcMiscData%WakeVTK_StartN) + UB(1:2) = ubound(SrcMiscData%WakeVTK_StartN) + if (.not. allocated(DstMiscData%WakeVTK_StartN)) then + allocate(DstMiscData%WakeVTK_StartN(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstMiscData%WakeVTK_StartN.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstMiscData%WakeVTK_StartN = SrcMiscData%WakeVTK_StartN + end if end subroutine subroutine AWAE_DestroyMisc(MiscData, ErrStat, ErrMsg) @@ -1564,6 +1581,9 @@ subroutine AWAE_DestroyMisc(MiscData, ErrStat, ErrMsg) if (allocated(MiscData%planeDomainExit)) then deallocate(MiscData%planeDomainExit) end if + if (allocated(MiscData%WakeVTK_StartN)) then + deallocate(MiscData%WakeVTK_StartN) + end if end subroutine subroutine AWAE_PackMisc(RF, Indata) @@ -1626,6 +1646,7 @@ subroutine AWAE_PackMisc(RF, Indata) end if call RegPackAlloc(RF, InData%V_amb_low_disk) call RegPackAlloc(RF, InData%planeDomainExit) + call RegPackAlloc(RF, InData%WakeVTK_StartN) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -1703,6 +1724,7 @@ subroutine AWAE_UnPackMisc(RF, OutData) end if call RegUnpackAlloc(RF, OutData%V_amb_low_disk); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%planeDomainExit); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%WakeVTK_StartN); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine AWAE_CopyLRGChunkType(SrcLRGChunkTypeData, DstLRGChunkTypeData, CtrlCode, ErrStat, ErrMsg) @@ -2172,8 +2194,12 @@ subroutine AWAE_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%OutDisWindYvalid = SrcParamData%OutDisWindYvalid end if DstParamData%OutFileRoot = SrcParamData%OutFileRoot - DstParamData%OutFileVTKRoot = SrcParamData%OutFileVTKRoot + DstParamData%OutFileFFvtkRoot = SrcParamData%OutFileFFvtkRoot + DstParamData%OutFileFFvtkWakeRoot = SrcParamData%OutFileFFvtkWakeRoot + DstParamData%OutFileFFvtkWakeNullData = SrcParamData%OutFileFFvtkWakeNullData DstParamData%VTK_tWidth = SrcParamData%VTK_tWidth + DstParamData%VTK_tWidthPlanes = SrcParamData%VTK_tWidthPlanes + DstParamData%WrPlanes = SrcParamData%WrPlanes DstParamData%WAT_Enabled = SrcParamData%WAT_Enabled DstParamData%WAT_FlowField => SrcParamData%WAT_FlowField end subroutine @@ -2298,8 +2324,12 @@ subroutine AWAE_PackParam(RF, Indata) call RegPackAlloc(RF, InData%OutDisWindY) call RegPackAlloc(RF, InData%OutDisWindYvalid) call RegPack(RF, InData%OutFileRoot) - call RegPack(RF, InData%OutFileVTKRoot) + call RegPack(RF, InData%OutFileFFvtkRoot) + call RegPack(RF, InData%OutFileFFvtkWakeRoot) + call RegPack(RF, InData%OutFileFFvtkWakeNullData) call RegPack(RF, InData%VTK_tWidth) + call RegPack(RF, InData%VTK_tWidthPlanes) + call RegPack(RF, InData%WrPlanes) call RegPack(RF, InData%WAT_Enabled) call RegPack(RF, associated(InData%WAT_FlowField)) if (associated(InData%WAT_FlowField)) then @@ -2384,8 +2414,12 @@ subroutine AWAE_UnPackParam(RF, OutData) call RegUnpackAlloc(RF, OutData%OutDisWindY); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%OutDisWindYvalid); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%OutFileRoot); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%OutFileVTKRoot); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%OutFileFFvtkRoot); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%OutFileFFvtkWakeRoot); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%OutFileFFvtkWakeNullData); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%VTK_tWidth); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%VTK_tWidthPlanes); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%WrPlanes); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WAT_Enabled); if (RegCheckErr(RF, RoutineName)) return if (associated(OutData%WAT_FlowField)) deallocate(OutData%WAT_FlowField) call RegUnpack(RF, IsAllocAssoc); if (RegCheckErr(RF, RoutineName)) return diff --git a/modules/awae/src/AWAE_vtk.f90 b/modules/awae/src/AWAE_vtk.f90 new file mode 100644 index 0000000000..ca8df71867 --- /dev/null +++ b/modules/awae/src/AWAE_vtk.f90 @@ -0,0 +1,508 @@ +!********************************************************************************************************************************** +! LICENSING +! Copyright (C) 2015-2016 National Renewable Energy Laboratory +! +! This file is part of Ambient Wind and Array Effects model for FAST.Farm. +! +! Licensed under the Apache License, Version 2.0 (the "License"); +! you may not use this file except in compliance with the License. +! You may obtain a copy of the License at +! +! http://www.apache.org/licenses/LICENSE-2.0 +! +! Unless required by applicable law or agreed to in writing, software +! distributed under the License is distributed on an "AS IS" BASIS, +! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +! See the License for the specific language governing permissions and +! limitations under the License. +! +!********************************************************************************************************************************** +!> VTK output routines for the AWAE module (wake-plane wireframes, structured +!! grid data files, and ParaView .vtk.series index files). +module AWAE_vtk + + use NWTC_Library + use AWAE_Types + use VTK + +#ifdef _OPENMP + use OMP_LIB +#endif + + implicit none + + private + + public :: PlaneAxes + public :: PlaneCorners + public :: Write_Planes_WireFrame + public :: Write_Planes_Data + public :: Write_WakePlane_Data_File + public :: Write_WakePlane_Series + public :: Write_WireFrame_Series + public :: Write_DisWind_Series + public :: Write_NullPlane + +contains + +!> Construct the orthonormal in-plane basis (yhat, zhat) in the global +!! inertial frame for a wake plane with unit normal `xhat`. yhat is the +!! horizontal in-plane direction (no Z component); zhat = xhat x yhat. +!! If xhat is purely vertical, yhat falls back to the global Y axis. +subroutine PlaneAxes(xhat, yhat, zhat) + real(ReKi), intent(in ) :: xhat(3) + real(ReKi), intent( out) :: yhat(3) + real(ReKi), intent( out) :: zhat(3) + real(ReKi) :: ynorm + + yhat = (/ -xhat(2), xhat(1), 0.0_ReKi /) + ynorm = TwoNorm(yhat) + if (ynorm > 0.0_ReKi) then + yhat = yhat / ynorm + else + yhat = (/ 0.0_ReKi, 1.0_ReKi, 0.0_ReKi /) + end if + + zhat(1) = xhat(2)*yhat(3) - xhat(3)*yhat(2) + zhat(2) = xhat(3)*yhat(1) - xhat(1)*yhat(3) + zhat(3) = xhat(1)*yhat(2) - xhat(2)*yhat(1) +end subroutine PlaneAxes + + +!> Compute the four corners of a single wake plane in the global inertial +!! frame from the plane normal `xhat`, the plane center `pc`, and the +!! in-plane half-extents `hy` (horizontal) and `hz` (vertical-ish). +!! Corners are returned counter-clockwise about +xhat. +subroutine PlaneCorners(xhat, pc, hy, hz, corners) + real(ReKi), intent(in ) :: xhat(3) !< Plane normal (unit vector) + real(ReKi), intent(in ) :: pc(3) !< Plane center, global frame + real(ReKi), intent(in ) :: hy !< In-plane horizontal half-extent + real(ReKi), intent(in ) :: hz !< In-plane vertical-ish half-extent + real(ReKi), intent( out) :: corners(3,4) !< Four corner positions, global frame + + real(ReKi) :: yhat(3), zhat(3) + + call PlaneAxes(xhat, yhat, zhat) + + ! Four corners, ordered counter-clockwise about +xhat + corners(:,1) = pc - hy*yhat - hz*zhat + corners(:,2) = pc + hy*yhat - hz*zhat + corners(:,3) = pc + hy*yhat + hz*zhat + corners(:,4) = pc - hy*yhat + hz*zhat +end subroutine PlaneCorners + + +!> Write the four edges of every active wake plane for each turbine as a +!! wireframe to a per-turbine VTK polydata file. One file is produced per +!! turbine per output step, named .T.WakePlanesWireFrame..vtk. +subroutine Write_Planes_WireFrame(p, u, t, Tstr) + type(AWAE_ParameterType), intent(in) :: p !< AWAE parameters + type(AWAE_InputType), intent(in) :: u !< AWAE inputs + real(DbKi), intent(in) :: t !< Current simulation time in seconds + character(*), intent(in) :: Tstr !< Zero-padded time-step string for filenames + + integer(IntKi) :: nt_wp, np_wp, nActive, ipt + real(ReKi) :: pc(3), corners(3,4) + real(ReKi) :: hy, hz + real(ReKi), allocatable :: WPPoints(:,:) + integer(IntKi), allocatable :: WPLines(:,:) + type(VTK_Misc) :: mvtk + character(1024) :: WPFileName + integer(IntKi) :: ntWidth + character(16) :: TurbNum, FmtStrT + + ! Plane half-extents in the local Y and Z (in-plane) directions + hy = p%y(p%NumRadii-1) + hz = p%z(p%NumRadii-1) + + ! Zero-padded width sufficient to represent all turbine indices + ntWidth = max(1, int(floor(log10(real(max(p%NumTurbines, 1), ReKi))) + 1, IntKi)) + write(FmtStrT, '(A,I0,A,I0,A)') '(I', ntWidth, '.', ntWidth, ')' + + do nt_wp = 1, p%NumTurbines + nActive = NINT(u%NumPlanes(nt_wp)) + if (nActive <= 0) cycle + + allocate(WPPoints(3, 4*nActive)) + allocate(WPLines(2, 4*nActive)) + + do np_wp = 0, nActive - 1 + pc = u%p_plane(:, np_wp, nt_wp) + call PlaneCorners(u%xhat_plane(:, np_wp, nt_wp), pc, hy, hz, corners) + + ipt = 4*np_wp + WPPoints(:, ipt+1) = corners(:,1) + WPPoints(:, ipt+2) = corners(:,2) + WPPoints(:, ipt+3) = corners(:,3) + WPPoints(:, ipt+4) = corners(:,4) + + ! Four edges of the closed quad outline (0-based point indices for VTK) + WPLines(:, ipt+1) = (/ ipt, ipt+1 /) + WPLines(:, ipt+2) = (/ ipt+1, ipt+2 /) + WPLines(:, ipt+3) = (/ ipt+2, ipt+3 /) + WPLines(:, ipt+4) = (/ ipt+3, ipt /) + end do + + write(TurbNum, FmtStrT) nt_wp + WPFileName = trim(p%OutFileFFvtkWakeRoot)//".T"//trim(TurbNum)// & + ".WakePlanesWireFrame."//trim(Tstr)//".vtk" + + call vtk_misc_init(mvtk) + if (vtk_new_ascii_file(WPFileName, & + "Wake plane wireframes for turbine "//trim(TurbNum)// & + " at time = "//trim(num2lstr(t))//" seconds.", mvtk)) then + call vtk_dataset_polydata(WPPoints, mvtk, .false.) + call vtk_lines(WPLines, mvtk) + call vtk_close_file(mvtk) + end if + + deallocate(WPPoints, WPLines) + + end do +end subroutine Write_Planes_WireFrame + + +!> Write one VTK STRUCTURED_GRID file per active wake plane containing the +!! wake velocity sampled on the plane's structured Y-Z grid. Point +!! coordinates and velocity vectors are written in the global inertial frame. +subroutine Write_Planes_Data(p, u, m, n, t, Tstr) + type(AWAE_ParameterType), intent(in) :: p !< AWAE parameters + type(AWAE_InputType), intent(in) :: u !< AWAE inputs + type(AWAE_MiscVarType), intent(inout) :: m !< AWAE misc variables + integer(IntKi), intent(in) :: n !< Current low-resolution time step index + real(DbKi), intent(in) :: t !< Current simulation time in seconds + character(*), intent(in) :: Tstr !< Zero-padded time-step string for filenames + + integer(IntKi) :: nt_wp, np_wp, nActive, ntWidth + integer(IntKi) :: jy, kz, idx, nY, nZ + real(ReKi) :: pc(3), xhat(3), yhat(3), zhat(3) + real(ReKi), allocatable :: Pts(:,:), Vel(:,:) + character(1024) :: WPFileName + character(16) :: FmtStrWk, FmtStrT + character(p%VTK_tWidthPlanes) :: PlaneNum + character(6) :: TurbNumStr + + ! Plane grid dimensions (p%y and p%z run from -NumRadii+1 to NumRadii-1) + nY = 2*p%NumRadii - 1 + nZ = 2*p%NumRadii - 1 + + ! Field width for plane index, based on MaxPlanes (1-based count) + write(FmtStrWk, '(A,I0,A,I0,A)') '(I', p%VTK_tWidthPlanes, '.', p%VTK_tWidthPlanes, ')' + ! Zero-padded width sufficient to represent all turbine indices + ntWidth = max(1, int(floor(log10(real(max(p%NumTurbines, 1), ReKi))) + 1, IntKi)) + write(FmtStrT, '(A,I0,A,I0,A)') '(A1,I', ntWidth, '.', ntWidth, ')' + + allocate(Pts(3, nY*nZ)) + allocate(Vel(3, nY*nZ)) + + do nt_wp = 1, p%NumTurbines + ! Turbine number + write(TurbNumStr, FmtStrT) "T", nt_wp + + ! Number of planes active + nActive = NINT(u%NumPlanes(nt_wp)) + + do np_wp = 0, nActive - 1 + pc = u%p_plane(:, np_wp, nt_wp) + xhat = u%xhat_plane(:, np_wp, nt_wp) + call PlaneAxes(xhat, yhat, zhat) + + ! Fill points and velocity in VTK natural order (jy fastest, then kz) + do kz = 1, nZ + do jy = 1, nY + idx = jy + (kz-1)*nY + Pts(:, idx) = pc + p%y(jy-p%NumRadii) * yhat & + + p%z(kz-p%NumRadii) * zhat + Vel(:, idx) = u%Vx_wake(jy-p%NumRadii, kz-p%NumRadii, np_wp, nt_wp) * xhat & + + u%Vy_wake(jy-p%NumRadii, kz-p%NumRadii, np_wp, nt_wp) * yhat & + + u%Vz_wake(jy-p%NumRadii, kz-p%NumRadii, np_wp, nt_wp) * zhat + end do + end do + + ! Per-plane index string with consistent zero-padded width + write(PlaneNum, FmtStrWk) np_wp + + WPFileName = trim(p%OutFileFFvtkWakeRoot)//"."//trim(TurbNumStr)//".WakePlane_"//trim(PlaneNum)// & + "."//trim(Tstr)//".vtk" + + call Write_WakePlane_Data_File(WPFileName, & + "Wake plane "//trim(PlaneNum)//" at time = "// & + trim(num2lstr(t))//" seconds.", nY, nZ, Pts, Vel) + + ! track what plane this was first written out at (initialized to huge, so this will catch only the first) + if (m%WakeVTK_StartN(np_wp,nt_wp) > n) m%WakeVTK_StartN(np_wp,nt_wp) = n + end do + end do + + deallocate(Pts, Vel) +end subroutine Write_Planes_Data + + +!> Helper: write a single 2D VTK STRUCTURED_GRID file for one wake plane +!! containing point coordinates and a "WakeVelocity" point-data vector. +subroutine Write_WakePlane_Data_File(WPFileName, label, n1, n2, Pts, Vel) + character(*), intent(in) :: WPFileName + character(*), intent(in) :: label + integer(IntKi), intent(in) :: n1, n2 + real(ReKi), intent(in) :: Pts(:,:) !< 3 x (n1*n2) + real(ReKi), intent(in) :: Vel(:,:) !< 3 x (n1*n2) + type(VTK_Misc) :: mvtk + + call vtk_misc_init(mvtk) + if (vtk_new_ascii_file(WPFileName, label, mvtk)) then + call vtk_dataset_structured_grid(Pts, n1, n2, 1, mvtk) + call vtk_point_data_init(mvtk) + call vtk_point_data_vector(Vel, "WakeVelocityDeficit", mvtk) + call vtk_close_file(mvtk) + end if +end subroutine Write_WakePlane_Data_File + + +!> Write the final ParaView .vtk.series index files for all wake planes +!! that were output during the simulation. Called once from AWAE_End. +!! Each series file lists every VTK output step, referencing either the +!! actual per-plane VTK file or a null-data placeholder for steps before +!! the plane first appeared. +subroutine Write_WakePlane_Series(p, m) + type(AWAE_ParameterType), intent(in) :: p !< AWAE parameters + type(AWAE_MiscVarType), intent(in) :: m !< AWAE misc variables + + integer(IntKi) :: nt_wp, np_wp, ntWidth, n_final, n_out + integer(IntKi) :: UnSer, out_idx, SerErrStat + character(ErrMsgLen) :: SerErrMsg + character(16) :: FmtStrWk, FmtStrT + character(p%VTK_tWidthPlanes) :: PlaneNum + character(6) :: TurbNumStr + character(1024) :: SeriesFile, baseName, EntryName, VTKprefix + character(p%VTK_tWidth) :: TstrOut + character(32) :: TimeStr + real(DbKi) :: t_out + logical :: firstEntry + + if (.not. p%WrPlanes) return + + ! Final low-res time step and total number of VTK output steps + n_final = p%NumDT - 1 + n_out = n_final / p%WrDisSkp1 + + ! Compute basename once (shared across all series files) + call GetPath(p%OutFileFFvtkWakeRoot, EntryName, baseName) + + ! Format strings for zero-padded plane and turbine numbers + write(FmtStrWk, '(A,I0,A,I0,A)') '(I', p%VTK_tWidthPlanes, '.', p%VTK_tWidthPlanes, ')' + ntWidth = max(1, int(floor(log10(real(max(p%NumTurbines, 1), ReKi))) + 1, IntKi)) + write(FmtStrT, '(A,I0,A,I0,A)') '(A1,I', ntWidth, '.', ntWidth, ')' + + do nt_wp = 1, p%NumTurbines + write(TurbNumStr, FmtStrT) "T", nt_wp + do np_wp = 0, p%MaxPlanes - 1 + ! Skip planes that were never written during the simulation + if (m%WakeVTK_StartN(np_wp, nt_wp) > n_final) cycle + + write(PlaneNum, FmtStrWk) np_wp + VTKprefix = trim(TurbNumStr)//".WakePlane_"//trim(PlaneNum) + SeriesFile = trim(p%OutFileFFvtkWakeRoot)//"."//trim(VTKprefix)//".vtk.series" + + !$OMP critical(fileopen_critical) + call GetNewUnit(UnSer, SerErrStat, SerErrMsg) + call OpenFOutFile(UnSer, SeriesFile, SerErrStat, SerErrMsg) + !$OMP end critical(fileopen_critical) + if (SerErrStat >= AbortErrLev) cycle + + write(UnSer, '(A)') '{' + write(UnSer, '(A)') ' "file-series-version" : "1.0",' + write(UnSer, '(A)') ' "files" : [' + + firstEntry = .true. + do out_idx = 0, n_out + t_out = real(out_idx, DbKi) * real(p%WrDisSkp1, DbKi) * p%dt_low + write(TstrOut, '(i'//trim(Num2LStr(p%VTK_tWidth))//'.'// & + trim(Num2LStr(p%VTK_tWidth))//')') out_idx + + if (m%WakeVTK_StartN(np_wp, nt_wp) > out_idx) then + EntryName = trim(p%OutFileFFvtkWakeNullData) + else + EntryName = trim(baseName)//"."//trim(VTKprefix)//"."//trim(TstrOut)//".vtk" + endif + + write(TimeStr, '(F14.5)') t_out + if (firstEntry) then + write(UnSer, '(A,A,A,A,A)') ' { "name" : "', trim(EntryName), & + '", "time" : ', trim(TimeStr), ' }' + firstEntry = .false. + else + write(UnSer, '(A,A,A,A,A)') ' ,{ "name" : "', trim(EntryName), & + '", "time" : ', trim(TimeStr), ' }' + end if + end do + + write(UnSer, '(A)') ' ]' + write(UnSer, '(A)') '}' + + !$OMP critical(fileopen_critical) + close(UnSer) + !$OMP end critical(fileopen_critical) + end do + end do +end subroutine Write_WakePlane_Series + +!> Write a zeroed-out "null" wake-plane VTK file used as a placeholder +!! in the ParaView series for time steps before a plane first appears. +subroutine Write_NullPlane(OutFileVTKwakeDir, p) + character(*), intent(in) :: OutFileVTKwakeDir !< Directory for VTK wake output files + type(AWAE_ParameterType), intent(in) :: p !< AWAE parameters + + real(ReKi), allocatable :: PtsVel(:,:) + integer(IntKi) :: nY, nZ + + nY = 2*p%NumRadii - 1 + nZ = 2*p%NumRadii - 1 + allocate(PtsVel(3, nY*nZ)) + PtsVel = 0.0_ReKi + call Write_WakePlane_Data_File(trim(OutFileVTKwakeDir)//PathSep//p%OutFileFFvtkWakeNullData, & + "Wake plane NULL at time = NULL seconds.", nY, nZ, PtsVel, PtsVel) + deallocate(PtsVel) +end subroutine Write_NullPlane + + +!> Write one ParaView .vtk.series index file per turbine for the wireframe +!! VTK files produced by Write_Planes_WireFrame. Called once from AWAE_End. +subroutine Write_WireFrame_Series(p) + type(AWAE_ParameterType), intent(in) :: p !< AWAE parameters + + integer(IntKi) :: nt_wp, ntWidth, n_out + integer(IntKi) :: UnSer, out_idx, SerErrStat + character(ErrMsgLen) :: SerErrMsg + character(16) :: FmtStrT, TurbNum + character(1024) :: SeriesFile, baseName, EntryName + character(p%VTK_tWidth) :: TstrOut + character(32) :: TimeStr + real(DbKi) :: t_out + logical :: firstEntry + + if (.not. p%WrPlanes) return + + ! Total number of VTK output steps + n_out = (p%NumDT - 1) / p%WrDisSkp1 + + ! Compute basename once (filename portion of OutFileFFvtkWakeRoot) + call GetPath(p%OutFileFFvtkWakeRoot, EntryName, baseName) + + ! Format string for zero-padded turbine number + ntWidth = max(1, int(floor(log10(real(max(p%NumTurbines, 1), ReKi))) + 1, IntKi)) + write(FmtStrT, '(A,I0,A,I0,A)') '(I', ntWidth, '.', ntWidth, ')' + + do nt_wp = 1, p%NumTurbines + write(TurbNum, FmtStrT) nt_wp + + SeriesFile = trim(p%OutFileFFvtkWakeRoot)//".T"//trim(TurbNum)// & + ".WakePlanesWireFrame.vtk.series" + + !$OMP critical(fileopen_critical) + call GetNewUnit(UnSer, SerErrStat, SerErrMsg) + call OpenFOutFile(UnSer, SeriesFile, SerErrStat, SerErrMsg) + !$OMP end critical(fileopen_critical) + if (SerErrStat >= AbortErrLev) cycle + + write(UnSer, '(A)') '{' + write(UnSer, '(A)') ' "file-series-version" : "1.0",' + write(UnSer, '(A)') ' "files" : [' + + firstEntry = .true. + do out_idx = 0, n_out + t_out = real(out_idx, DbKi) * real(p%WrDisSkp1, DbKi) * p%dt_low + write(TstrOut, '(i'//trim(Num2LStr(p%VTK_tWidth))//'.'// & + trim(Num2LStr(p%VTK_tWidth))//')') out_idx + + EntryName = trim(baseName)//".T"//trim(TurbNum)// & + ".WakePlanesWireFrame."//trim(TstrOut)//".vtk" + + write(TimeStr, '(F14.5)') t_out + if (firstEntry) then + write(UnSer, '(A,A,A,A,A)') ' { "name" : "', trim(EntryName), & + '", "time" : ', trim(TimeStr), ' }' + firstEntry = .false. + else + write(UnSer, '(A,A,A,A,A)') ' ,{ "name" : "', trim(EntryName), & + '", "time" : ', trim(TimeStr), ' }' + end if + end do + + write(UnSer, '(A)') ' ]' + write(UnSer, '(A)') '}' + + !$OMP critical(fileopen_critical) + close(UnSer) + !$OMP end critical(fileopen_critical) + end do +end subroutine Write_WireFrame_Series + + +!> Write a ParaView .vtk.series index file for one disturbed-wind output +!! slice (XY, YZ, or XZ). The slice label and 1-based index select which +!! series of per-timestep VTK files to reference. +subroutine Write_DisWind_Series(p, sliceLabel, iSlice) + type(AWAE_ParameterType), intent(in) :: p !< AWAE parameters + character(*), intent(in) :: sliceLabel !< "DisXY", "DisYZ", or "DisXZ" + integer(IntKi), intent(in) :: iSlice !< 1-based slice index + + integer(IntKi) :: n_out + integer(IntKi) :: UnSer, out_idx, SerErrStat + character(ErrMsgLen) :: SerErrMsg + character(1024) :: SeriesFile, baseName, EntryName + character(p%VTK_tWidth) :: TstrOut + character(32) :: TimeStr + character(3) :: PlaneNumStr + real(DbKi) :: t_out + logical :: firstEntry + + ! Total number of VTK output steps + n_out = (p%NumDT - 1) / p%WrDisSkp1 + + ! 3-digit zero-padded slice index (matches CalcOutput naming) + write(PlaneNumStr, '(i3.3)') iSlice + + ! Compute basename (filename portion of OutFileFFvtkRoot) + call GetPath(p%OutFileFFvtkRoot, EntryName, baseName) + + ! Series filename + SeriesFile = trim(p%OutFileFFvtkRoot)//".Low."//trim(sliceLabel)//PlaneNumStr//".vtk.series" + + !$OMP critical(fileopen_critical) + call GetNewUnit(UnSer, SerErrStat, SerErrMsg) + call OpenFOutFile(UnSer, SeriesFile, SerErrStat, SerErrMsg) + !$OMP end critical(fileopen_critical) + if (SerErrStat >= AbortErrLev) return + + write(UnSer, '(A)') '{' + write(UnSer, '(A)') ' "file-series-version" : "1.0",' + write(UnSer, '(A)') ' "files" : [' + + firstEntry = .true. + do out_idx = 0, n_out + t_out = real(out_idx, DbKi) * real(p%WrDisSkp1, DbKi) * p%dt_low + write(TstrOut, '(i'//trim(Num2LStr(p%VTK_tWidth))//'.'// & + trim(Num2LStr(p%VTK_tWidth))//')') out_idx + + EntryName = trim(baseName)//".Low."//trim(sliceLabel)//PlaneNumStr//"."//trim(TstrOut)//".vtk" + + write(TimeStr, '(F14.5)') t_out + if (firstEntry) then + write(UnSer, '(A,A,A,A,A)') ' { "name" : "', trim(EntryName), & + '", "time" : ', trim(TimeStr), ' }' + firstEntry = .false. + else + write(UnSer, '(A,A,A,A,A)') ' ,{ "name" : "', trim(EntryName), & + '", "time" : ', trim(TimeStr), ' }' + end if + end do + + write(UnSer, '(A)') ' ]' + write(UnSer, '(A)') '}' + + !$OMP critical(fileopen_critical) + close(UnSer) + !$OMP end critical(fileopen_critical) +end subroutine Write_DisWind_Series + +end module AWAE_vtk diff --git a/modules/wakedynamics/src/WakeDynamics.f90 b/modules/wakedynamics/src/WakeDynamics.f90 index 05b357088e..cdeaf256e6 100644 --- a/modules/wakedynamics/src/WakeDynamics.f90 +++ b/modules/wakedynamics/src/WakeDynamics.f90 @@ -500,7 +500,9 @@ subroutine WD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut allocate( u%Ct_azavg ( 0:p%NumRadii-1 ),stat=errStat2); if (Failed0('u%Ct_azavg.')) return; allocate( u%Cq_azavg ( 0:p%NumRadii-1 ),stat=errStat2); if (Failed0('u%Cq_azavg.')) return; if (errStat /= ErrID_None) return - + u%V_plane = 0.0_ReKi + u%Ct_azavg = 0.0_ReKi + u%Cq_azavg = 0.0_ReKi @@ -548,6 +550,8 @@ subroutine WD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut xd%Vx_wake = 0.0_ReKi xd%Vr_wake = 0.0_ReKi xd%Vx_wake2 = 0.0_ReKi + xd%Vy_wake2 = 0.0_ReKi + xd%Vz_wake2 = 0.0_ReKi xd%V_plane_filt = 0.0_ReKi xd%Vx_wind_disk_filt = 0.0_ReKi xd%TI_amb_filt = 0.0_ReKi @@ -574,12 +578,15 @@ subroutine WD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut allocate ( m%vt_tot (0:p%NumRadii-1,0:p%MaxNumPlanes-1 ) , STAT=ErrStat2 ); if (Failed0('m%vt_tot.')) return; allocate ( m%vt_amb (0:p%NumRadii-1,0:p%MaxNumPlanes-1 ) , STAT=ErrStat2 ); if (Failed0('m%vt_amb.')) return; allocate ( m%vt_shr (0:p%NumRadii-1,0:p%MaxNumPlanes-1 ) , STAT=ErrStat2 ); if (Failed0('m%vt_shr.')) return; + m%dvtdr = 0.0_ReKi + m%vt_tot = 0.0_ReKi + m%vt_amb = 0.0_ReKi + m%vt_shr = 0.0_ReKi else if (p%Mod_Wake == Mod_Wake_Cartesian .or. p%Mod_Wake == Mod_Wake_Curl) then allocate ( m%nu_dvx_dy(-p%NumRadii+1:p%NumRadii-1,-p%NumRadii+1:p%NumRadii-1), STAT=ErrStat2 ); if (Failed0('m%nu_dvx_dy.')) return; allocate ( m%nu_dvx_dz(-p%NumRadii+1:p%NumRadii-1,-p%NumRadii+1:p%NumRadii-1), STAT=ErrStat2 ); if (Failed0('m%nu_dvx_dz.')) return; allocate ( m%dnuvx_dy (-p%NumRadii+1:p%NumRadii-1,-p%NumRadii+1:p%NumRadii-1), STAT=ErrStat2 ); if (Failed0('m%dnuvx_dy.' )) return; allocate ( m%dnuvx_dz (-p%NumRadii+1:p%NumRadii-1,-p%NumRadii+1:p%NumRadii-1), STAT=ErrStat2 ); if (Failed0('m%dnuvx_dz.' )) return; - if (errStat /= ErrID_None) return m%nu_dvx_dy = 0.0_ReKi m%nu_dvx_dz = 0.0_ReKi m%dnuvx_dy = 0.0_ReKi @@ -595,10 +602,18 @@ subroutine WD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut allocate ( m%d(0:p%NumRadii-1 ), STAT=ErrStat2 ); if (Failed0('m%d.')) return; allocate ( m%r_wake(0:p%NumRadii-1 ), STAT=ErrStat2 ); if (Failed0('m%r_wake.' )) return; allocate ( m%Vx_high(0:p%NumRadii-1 ), STAT=ErrStat2 ); if (Failed0('m%Vx_high.' )) return; - allocate ( m%Vt_wake(0:p%NumRadii-1 ), STAT=ErrStat2 ); if (Failed0('m%Vx_high.' )) return; + allocate ( m%Vt_wake(0:p%NumRadii-1 ), STAT=ErrStat2 ); if (Failed0('m%Vt_wake.' )) return; allocate ( m%Vx_polar(0:p%NumRadii-1 ), STAT=ErrStat2 ); if (Failed0('m%Vx_polar.')) return; + + m%a = 0.0_ReKi + m%b = 0.0_ReKi + m%c = 0.0_ReKi + m%d = 0.0_ReKi + m%r_wake = 0.0_ReKi + m%Vx_high = 0.0_ReKi m%Vx_polar = 0.0_ReKi - m%Vt_wake = 0.0_ReKi + m%Vt_wake = 0.0_ReKi + !............................................................................................ ! Define initialization output here !............................................................................................ @@ -952,6 +967,11 @@ subroutine WD_UpdateStates( t, n, u, p, x, xd, z, OtherState, m, errStat, errMsg !Used for debugging: write(51,'(I5,100(1x,ES10.2E2))') n, xd%x_plane(n), xd%x_plane(n)/xd%D_rotor_filt(n), xd%Vx_wind_disk_filt(n) + xd%Vx_wake(:,n), xd%Vr_wake(:,n) + + ! -------------------------------------------------------------------------------- + ! Drop planes that exit the buffer, and merge planes that collide. + ! -------------------------------------------------------------------------------- + xd%NumPlanes = xd%NumPlanes + 1.0 if ( NINT(xd%NumPlanes) > p%MaxNumPlanes ) then xd%NumPlanes = real(p%MaxNumPlanes,ReKi) @@ -1123,8 +1143,8 @@ subroutine updateVelocityCartesian() call gradient_z(m%nu_dvx_dz, p%dr, m%dnuvx_dz ) ! Loop through all the points on the plane (y, z) - do iz = -p%NumRadii+2, p%NumRadii-2 - do iy = -p%NumRadii+2, p%NumRadii-2 + do iz = -p%NumRadii+1, p%NumRadii-1 + do iy = -p%NumRadii+1, p%NumRadii-1 ! Eddy viscosity term divTau = m%dnuvx_dy(iy,iz) + m%dnuvx_dz(iy,iz) @@ -1616,7 +1636,17 @@ subroutine WD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg ) if ( p%OutAllPlanes ) then call mkdir(p%OutFileVTKDir) endif - + + ! set initial outputs for Cartesian/Curl wake models from the polar info. + if (p%Mod_Wake == Mod_Wake_Cartesian .or. p%Mod_Wake == Mod_Wake_Curl) then + call Axisymmetric2CartesianVx(y%Vx_wake(:,0), p%r, p%y, p%z, y%Vx_wake2(:,:,0)) + y%Vy_wake2(:,:,0) = 0.0_ReKi + y%Vz_wake2(:,:,0) = 0.0_ReKi + y%Vx_wake2(:,:,1) = y%Vx_wake2(:,:,0) + y%Vy_wake2(:,:,1) = 0.0_ReKi + y%Vz_wake2(:,:,1) = 0.0_ReKi + end if + else y%x_plane = xd%x_plane y%p_plane = xd%p_plane @@ -1631,13 +1661,14 @@ subroutine WD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg ) y%NumPlanes = xd%NumPlanes end if - ! --- Linearly decay wake deficits in the buffer region based on distance + ! --- Linearly decay wake deficits in the buffer region based on distance. + ! Polar arrays are tapered first; for the Polar wake model these tapered values + ! are then converted to the Cartesian output below, so the taper carries through. do i = 0,maxPln - if ( xd%x_plane(i) > p%x_full ) then - ! Note: Clamp to zero just in case, but all wake planes that propagated past x_buff should have been removed. - ScBuff = max( ( p%x_buff - xd%x_plane(i) ) / p%d_buff , 0.0 ) - y%Vx_wake(:,i) = y%Vx_wake(:,i) * ScBuff - y%Vr_wake(:,i) = y%Vr_wake(:,i) * ScBuff + ScBuff = BufferScale(xd%x_plane(i)) + if ( ScBuff < 1.0_ReKi ) then + y%Vx_wake(:,i) = y%Vx_wake(:,i) * ScBuff + y%Vr_wake(:,i) = y%Vr_wake(:,i) * ScBuff end if end do @@ -1662,11 +1693,29 @@ subroutine WD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg ) enddo endif else if (p%Mod_Wake == Mod_Wake_Cartesian .or. p%Mod_Wake == Mod_Wake_Curl) then - do i = 0, maxPln - y%Vx_wake2(:,:,i) = xd%Vx_wake2(:,:,i) - y%Vy_wake2(:,:,i) = xd%Vy_wake2(:,:,i) - y%Vz_wake2(:,:,i) = xd%Vz_wake2(:,:,i) - enddo + ! For Cartesian/Curl, the Cartesian output arrays are copied from the (untapered) + ! discrete state. Apply the buffer-region taper here so that the wake deficit + ! handed to AWAE fades linearly from x_Full to x_Buff. Do NOT modify xd%*_wake2: + ! tapering the state would compound the scaling at every time step. + ! Skip on firstPass: y%Vx_wake2 was already set from NearWakeCorrection above, + ! and xd%Vx_wake2 has not yet been initialized (still zero). + if (.not. OtherState%firstPass) then + do i = 0, maxPln + ScBuff = BufferScale(xd%x_plane(i)) + y%Vx_wake2(:,:,i) = xd%Vx_wake2(:,:,i) * ScBuff + y%Vy_wake2(:,:,i) = xd%Vy_wake2(:,:,i) * ScBuff + y%Vz_wake2(:,:,i) = xd%Vz_wake2(:,:,i) * ScBuff + enddo + ! Recompute Cartesian gradients from the tapered output field so that the WAT + ! gradient term in Calc_k_WAT and the dvx_dy/dz VTK diagnostics fade consistently + ! with the velocity deficit in the buffer region. + if ( p%WAT .or. p%OutAllPlanes ) then + do i = 0, maxPln + call gradient_y(y%Vx_wake2(:,:,i), p%dr, m%dvx_dy(:,:,i)) + call gradient_z(y%Vx_wake2(:,:,i), p%dr, m%dvx_dz(:,:,i)) + end do + endif + end if endif ! Curl or Polar ! --- WAT - Compute k_mt and add turbulence @@ -1675,6 +1724,22 @@ subroutine WD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg ) end if contains + !> Linear buffer-region taper applied to the output wake fields. + !! Returns 1 for x_plane <= x_Full, fades linearly to 0 at x_Buff, and is + !! clamped to 0 beyond x_Buff (or whenever d_Buff <= 0, which avoids a + !! divide-by-zero when the user sets NumDBuff = 0). + pure function BufferScale(x_plane) result(ScBuff_loc) + real(ReKi), intent(in) :: x_plane + real(ReKi) :: ScBuff_loc + if ( x_plane <= p%x_Full ) then + ScBuff_loc = 1.0_ReKi + else if ( p%d_Buff > 0.0_ReKi ) then + ScBuff_loc = max( ( p%x_Buff - x_plane ) / p%d_Buff , 0.0_ReKi ) + else + ScBuff_loc = 0.0_ReKi + end if + end function BufferScale + subroutine Calc_k_WAT() integer(intKi) :: i, iy, iz real(ReKi) :: C, S, dvdr, dvdtheta_r, R, r_tmp @@ -1850,47 +1915,58 @@ subroutine InitStatesWithInputs(numPlanes, numRadii, u, p, xd, m, errStat, errMs character(ErrMsgLen) :: ErrMsg2 real(ReKi) :: correction(3) ! Note, all of these states will have been set to zero in the WD_Init routine - - + ErrStat = ErrID_None ErrMsg = "" - - + correction = 0.0_ReKi do i = 0, 1 xd%x_plane (i) = u%Vx_rel_disk*real(i,ReKi)*real(p%DT_low,ReKi) xd%YawErr_filt (i) = u%YawErr xd%psi_skew_filt = u%psi_skew xd%chi_skew_filt = u%chi_skew - + correction = correction + GetYawCorrection(u%YawErr, u%xhat_disk, xd%x_plane(i), p, errStat2, errMsg2) call SetErrStat(errStat2, errMsg2, errStat, errMsg, RoutineName) if (errStat >= AbortErrLev) then ! TEST: E3 return end if - + !correction = ( p%C_HWkDfl_x + p%C_HWkDfl_xY*u%YawErr )*xd%x_plane(i) + correctionA - + xd%p_plane (:,i) = u%p_hub(:) + xd%x_plane(i)*u%xhat_disk(:) + correction xd%xhat_plane(:,i) = u%xhat_disk(:) xd%V_plane_filt(:,i) = u%V_plane(:,i) xd%Vx_wind_disk_filt(i) = u%Vx_wind_disk xd%TI_amb_filt (i) = u%TI_amb xd%D_rotor_filt (i) = u%D_rotor - - end do - + xd%Vx_rel_disk_filt = u%Vx_rel_disk - + ! Initialze Ct_azavg_filt, Cq_azavg_filt, and Vx_wake; Vr_wake is already initialized to zero, so, we don't need to do that here. xd%Ct_azavg_filt (:) = u%Ct_azavg(:) xd%Cq_azavg_filt (:) = u%Cq_azavg(:) - + call NearWakeCorrection( xd%Ct_azavg_filt, xd%Cq_azavg_filt, xd%Vx_rel_disk_filt, p, m, xd%Vx_wake(:,0), m%Vt_wake, xd%D_rotor_filt(0), errStat, errMsg ) xd%Vx_wake(:,1) = xd%Vx_wake(:,0) - + + + ! Initialize states for cartesian and curled wake formulations + if (p%Mod_Wake == Mod_Wake_Cartesian .or. p%Mod_Wake == Mod_Wake_Curl) then + ! Compute Vx(r) + call NearWakeCorrection( xd%Ct_azavg_filt, xd%Cq_azavg_filt, xd%Vx_rel_disk_filt, p, m, m%Vx_polar(:), m%Vt_wake, xd%D_rotor_filt(0), errStat2, errMsg2 ) + call SetErrStat(ErrStat2, ErrMsg2, errStat, errMsg, RoutineName) + if (errStat >= AbortErrLev) return + call Axisymmetric2CartesianVx(m%Vx_polar, p%r, p%y, p%z, xd%Vx_wake2(:,:,0)) + xd%Vy_wake2(:,:,0) = 0.0_ReKi + xd%Vz_wake2(:,:,0) = 0.0_ReKi + xd%Vx_wake2(:,:,1) = xd%Vx_wake2(:,:,0) + xd%Vy_wake2(:,:,1) = xd%Vy_wake2(:,:,0) + xd%Vz_wake2(:,:,1) = xd%Vz_wake2(:,:,0) + endif + end subroutine InitStatesWithInputs !---------------------------------------------------------------------------------------------------------------------------------- diff --git a/vs-build/modules/AWAE.vfproj b/vs-build/modules/AWAE.vfproj index 95ceae90ea..7e6227167b 100644 --- a/vs-build/modules/AWAE.vfproj +++ b/vs-build/modules/AWAE.vfproj @@ -94,6 +94,7 @@ +