FF: Curl/Cartesian wake fixes and AWAE VTK wake-plane output

.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/

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
 

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)

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,33 +1510,41 @@ 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
    allocate(m%Vamb_lowpol(    3, 0:p%n_rp_max*8 ),                             STAT=errStat2);  if (Failed0('m%Vamb_lowpol.'  )) return;
    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 <WindFilePath>/Low/DisXY<k>.t<n/p%WrDisSkp1>.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 <WindFilePath>/Low/DisYZ<k>.t<n/p%WrDisSkp1>.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,18 +2105,35 @@ 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 <WindFilePath>/Low/DisXZ<k>.t<n/p%WrDisSkp1>.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)], &
                                  [p%LowRes%oXYZ(1), p%OutDisWindY(k), p%LowRes%oXYZ(3)], &
                                  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

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