diff --git a/glue-codes/fast-farm/src/FASTWrapper.f90 b/glue-codes/fast-farm/src/FASTWrapper.f90 index 3207a139b2..4ff6b281f2 100644 --- a/glue-codes/fast-farm/src/FASTWrapper.f90 +++ b/glue-codes/fast-farm/src/FASTWrapper.f90 @@ -44,7 +44,7 @@ MODULE FASTWrapper PUBLIC :: FWrap_t0 ! call to compute outputs at t0 [and initialize some more variables] PUBLIC :: FWrap_Increment ! call to update states to n+1 and compute outputs at n+1 - PUBLIC :: FWrap_SetInputs + PUBLIC :: FWrap_SetWindTStart PUBLIC :: FWrap_CalcOutput @@ -119,7 +119,7 @@ SUBROUTINE FWrap_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, Init ExternInitData%windGrid_n(1) = InitInp%nX_high ExternInitData%windGrid_n(2) = InitInp%nY_high ExternInitData%windGrid_n(3) = InitInp%nZ_high - ExternInitData%windGrid_n(4) = InitInp%n_high_low + ExternInitData%windGrid_n(4) = InitInp%n_high_low+1 ! include a step at t-dt_high ExternInitData%windGrid_delta(1) = InitInp%dX_high ExternInitData%windGrid_delta(2) = InitInp%dY_high @@ -392,7 +392,7 @@ SUBROUTINE FWrap_Increment( t, n, u, p, x, xd, z, OtherState, y, m, ErrStat, Err !ELSE ! ! set the inputs needed for FAST - !call FWrap_SetInputs(u, m, t) <<< moved up into FAST.Farm FARM_UpdateStates + !call FWrap_SetWindTStart(u, m, t) <<< moved up into FAST.Farm FARM_UpdateStates ! call FAST p%n_FAST_low times (p%n_FAST_low is simply the number of steps to make per wrapper call. It is affected by MooringMod) do n_ss = 1, p%n_FAST_low @@ -434,7 +434,7 @@ SUBROUTINE FWrap_t0( u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg ) ErrMsg = '' ! set the inputs needed for FAST: - call FWrap_SetInputs(u, m, 0.0_DbKi) + call FWrap_SetWindTStart(u, m, 0.0_DbKi) ! compute the FAST t0 solution: call FAST_Solution0_T(m%Turbine, ErrStat2, ErrMsg2 ) @@ -685,16 +685,18 @@ SUBROUTINE FWrap_CalcOutput(p, u, y, m, ErrStat, ErrMsg) END SUBROUTINE FWrap_CalcOutput !---------------------------------------------------------------------------------------------------------------------------------- !> This subroutine sets the inputs needed before calling an instance of FAST -SUBROUTINE FWrap_SetInputs(u, m, t) +SUBROUTINE FWrap_SetWindTStart(u, m, t) TYPE(FWrap_InputType), INTENT(INOUT) :: u !< Inputs at t TYPE(FWrap_MiscVarType), INTENT(INOUT) :: m !< Misc variables for optimization (not copied in glue code) REAL(DbKi), INTENT(IN ) :: t !< current simulation time ! set the 4d-wind-inflow input array (a bit of a hack [simplification] so that we don't have large amounts of data copied in multiple data structures): - m%Turbine%IfW%p%FlowField%Grid4D%TimeStart = t + ! NOTE: the wind data starts at `t - DT_high` as one extra slice of wind data is added at start. If AeroDyn is updated to not require the `t-DT_high` + ! timestep, this can be changed + m%Turbine%IfW%p%FlowField%Grid4D%TimeStart = t - m%Turbine%IfW%p%FlowField%Grid4D%delta(4) -END SUBROUTINE FWrap_SetInputs +END SUBROUTINE FWrap_SetWindTStart !---------------------------------------------------------------------------------------------------------------------------------- END MODULE FASTWrapper !********************************************************************************************************************************** diff --git a/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 b/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 index 344241616b..5cd2381d1b 100644 --- a/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 +++ b/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 @@ -1204,7 +1204,7 @@ subroutine FARM_UpdateStates(t, n, farm, ErrStat, ErrMsg) ! set the inputs needed for FAST (these are slow-varying so can just be done once per farm time step) do nt = 1,farm%p%NumTurbines - call FWrap_SetInputs(farm%FWrap(nt)%u, farm%FWrap(nt)%m, t) + call FWrap_SetWindTStart(farm%FWrap(nt)%u, farm%FWrap(nt)%m, t) end do diff --git a/modules/awae/src/AWAE.f90 b/modules/awae/src/AWAE.f90 index d2a77cd576..2eb80fff69 100644 --- a/modules/awae/src/AWAE.f90 +++ b/modules/awae/src/AWAE.f90 @@ -751,7 +751,7 @@ subroutine HighResGridCalcOutput(n, u, p, xd, y, m, errStat, errMsg) real(ReKi) :: WAT_k ! WAT scaling factor (averaged from overlapping wakes) real(SiKi) :: WAT_V(3) ! WAT velocity contribution real(ReKi) :: Pos_global(3) ! global position - real(ReKi), allocatable :: WAT_B_BoxHi(:,:) ! position of WAT box (global) for each intermediate steps, shape: 3 x n_high_low + real(ReKi), allocatable :: WAT_B_BoxHi(:,:) ! position of WAT box (global) for each intermediate steps, shape: 3 x n_high_low_p1 real(ReKi), allocatable :: wk_R_p2i(:,:,:)!< Orientations from plane to inertial for each wake, shape: 3x3xnWake real(ReKi), allocatable :: wk_V(:,:) !< Wake velocity from each overlapping wake, shape: 3xnWake real(ReKi), allocatable :: wk_WAT_k(:) !< WAT scaling factors for all wakes (for overlap) @@ -760,7 +760,7 @@ subroutine HighResGridCalcOutput(n, u, p, xd, y, m, errStat, errMsg) integer(IntKi) :: maxPln integer(IntKi) :: maxN_wake integer(IntKi) :: NumGrid_high !< number of points in high res grid grid - integer(IntKi) :: n_high_low + integer(IntKi) :: n_high_low_p1 integer(IntKi) :: WAT_iT,WAT_iY,WAT_iZ !< indexes for WAT point (Time interchangeable with X) integer(IntKi) :: errStat2 character(*), parameter :: RoutineName = 'HighResGridCalcOutput' @@ -771,9 +771,9 @@ subroutine HighResGridCalcOutput(n, u, p, xd, y, m, errStat, errMsg) ! We only need one high res file for that last simulation time if ( (n/p%n_high_low) == (p%NumDT-1) ) then - n_high_low = 0 + n_high_low_p1 = 1 else - n_high_low = p%n_high_low + n_high_low_p1 = p%n_high_low_p1 end if maxN_wake = p%NumTurbines*( p%NumPlanes-1 ) @@ -786,10 +786,10 @@ subroutine HighResGridCalcOutput(n, u, p, xd, y, m, errStat, errMsg) ! Convect WAT Box tracer for each intermediate step ! Note: we substract because the high-res points are "before" current low res point if (p%WAT_Enabled) then - allocate ( WAT_B_BoxHi ( 3, 0:n_high_low), STAT=errStat2 ); if (errStat2 /= 0) call SetErrStat ( ErrID_Fatal, 'Could not allocate memory for WAT_B_BoxHi.', errStat, errMsg, RoutineName ) + allocate ( WAT_B_BoxHi ( 3, 0:n_high_low_p1), STAT=errStat2 ); if (errStat2 /= 0) call SetErrStat ( ErrID_Fatal, 'Could not allocate memory for WAT_B_BoxHi.', errStat, errMsg, RoutineName ) if (ErrStat >= AbortErrLev) return - do i_hl=0, n_high_low - WAT_B_BoxHi(1:3, i_hl) = xd%WAT_B_Box(1:3) - (n_high_low-i_hl) * xd%Ufarm(1:3) * real(p%DT_high,ReKi) + do i_hl=0, n_high_low_p1 + WAT_B_BoxHi(1:3, i_hl) = xd%WAT_B_Box(1:3) - (n_high_low_p1-i_hl) * xd%Ufarm(1:3) * real(p%DT_high,ReKi) enddo endif @@ -817,7 +817,7 @@ subroutine HighResGridCalcOutput(n, u, p, xd, y, m, errStat, errMsg) !$OMP& V_qs, & !$OMP& i_hl, Pos_global,& !$OMP& wk_WAT_k, WAT_k, WAT_iT, WAT_iY, WAT_iZ, WAT_V)& - !$OMP SHARED(NumGrid_High, m, u, p, y, xd, nt, maxPln, n_high_low, WAT_B_BoxHi, errStat, errMsg) + !$OMP SHARED(NumGrid_High, m, u, p, y, xd, nt, maxPln, n_high_low_p1, WAT_B_BoxHi, errStat, errMsg) ! Loop over all points of the high resolution ambient wind do iXYZ=1, NumGrid_high ! From flat index iXYZ to grid indices @@ -849,7 +849,7 @@ subroutine HighResGridCalcOutput(n, u, p, xd, y, m, errStat, errMsg) endif ! --- Store full velocity (Ambient + Wake QS + WAT) in grid - do i_hl=0, n_high_low + do i_hl=0, n_high_low_p1 ! Compute WAT velocity if (p%WAT_Enabled) then ! find location of grid point in the turbulent box, accounting for the convection of the box in between high res and low res @@ -929,7 +929,8 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO p%NumRadii = InitInp%InputFileData%NumRadii p%NumTurbines = InitInp%InputFileData%NumTurbines p%WindFilePath = InitInp%InputFileData%WindFilePath ! TODO: Make sure this wasn't specified with the trailing folder separator. Note: on Windows a trailing / or \ causes no problem! GJH - p%n_high_low = InitInp%n_high_low + p%n_high_low = InitInp%n_high_low ! number of timesteps between low res steps + p%n_high_low_p1 = InitInp%n_high_low + 1 ! include a time slice at t_low-DT_high (for interpolation in AeroDyn -- this is a hack) p%NumDT = InitInp%NumDT p%NOutDisWindXY = InitInp%InputFileData%NOutDisWindXY p%NOutDisWindYZ = InitInp%InputFileData%NOutDisWindYZ @@ -991,7 +992,7 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO ! This will establish certain parameters as well as all of the initialization outputs ! Sets: ! Parameters: nX_low, nY_low, nZ_low, nX_high, nY_high, nZ_high, Grid_low, - ! Grid_high, n_high_low, n_rp_max + ! Grid_high, n_high_low_p1, n_rp_max ! InitOutput: X0_high, Y0_high, Z0_high, dX_high, dY_high, dZ_high, nX_high, nY_high, nZ_high call AWAE_IO_InitGridInfo(InitInp, p, InitOut, errStat2, errMsg2); if(Failed()) return; @@ -1156,7 +1157,7 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO allocate ( y%V_plane(3,0:p%NumPlanes-1,1:p%NumTurbines), STAT=ErrStat2 ); if (Failed0('y%V_plane.' )) return; allocate ( y%Vdist_High(1:p%NumTurbines), STAT=ErrStat2 ); if (Failed0('y%Vdist_High.')) return; do i = 1, p%NumTurbines - allocate ( y%Vdist_High(i)%data(3,0:p%nX_high-1,0:p%nY_high-1,0:p%nZ_high-1,0:p%n_high_low), STAT=ErrStat2 ); if (Failed0('y%Vdist_High%data.')) return; + allocate ( y%Vdist_High(i)%data(3,0:p%nX_high-1,0:p%nY_high-1,0:p%nZ_high-1,0:p%n_high_low_p1), STAT=ErrStat2 ); if (Failed0('y%Vdist_High%data.')) return; y%Vdist_High(i)%data = 0.0_Siki end do @@ -1200,7 +1201,7 @@ subroutine AWAE_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitO allocate ( m%Vamb_high(1:p%NumTurbines), STAT=ErrStat2 ); if (Failed0('Could not allocate memory for m%Vamb_high.')) return; do i = 1, p%NumTurbines - allocate ( m%Vamb_high(i)%data(3,0:p%nX_high-1,0:p%nY_high-1,0:p%nZ_high-1,0:p%n_high_low), STAT=ErrStat2 ); if (Failed0('m%Vamb_high%data.')) return; + allocate ( m%Vamb_high(i)%data(3,0:p%nX_high-1,0:p%nY_high-1,0:p%nZ_high-1,0:p%n_high_low_p1), STAT=ErrStat2 ); if (Failed0('m%Vamb_high%data.')) return; end do allocate ( m%parallelFlag( 0:p%NumPlanes-2,1:p%NumTurbines ), STAT=errStat2 ); if (Failed0('m%parallelFlag.')) return; @@ -1388,128 +1389,225 @@ subroutine AWAE_UpdateStates( t, n, u, p, x, xd, z, OtherState, m, errStat, errM !! Output: Constraint states at t+dt type(AWAE_OtherStateType), intent(inout) :: OtherState !< Input: Other states at t; !! Output: Other states at t+dt - type(AWAE_MiscVarType), intent(inout) :: m !< Misc/optimization variables + type(AWAE_MiscVarType), target, intent(inout) :: m !< Misc/optimization variables integer(IntKi), intent( out) :: errStat !< Error status of the operation character(*), intent( out) :: errMsg !< Error message if errStat /= ErrID_None - type(AWAE_InputType) :: uInterp ! Interpolated/Extrapolated input + character(*), parameter :: RoutineName = 'AWAE_UpdateStates' integer(intKi) :: errStat2 ! temporary Error status character(ErrMsgLen) :: errMsg2 ! temporary Error message - character(*), parameter :: RoutineName = 'AWAE_UpdateStates' - integer(IntKi) :: n_high_low, nt, i_hl, i,j,k,c + type(AWAE_InputType) :: uInterp ! Interpolated/Extrapolated input + integer(IntKi) :: n_high_low, nt, i_hl + integer(IntKi) :: i,j,k,c + real(ReKi), pointer :: V_Grid(:,:,:,:) errStat = ErrID_None errMsg = "" - - ! Read the ambient wind data that is needed for t+dt, i.e., n+1 - + + ! If last time step, don't populate high-resolution grid if ( (n+1) == (p%NumDT-1) ) then n_high_low = 0 else n_high_low = p%n_high_low end if - if ( p%Mod_AmbWind == 1 ) then - ! read from file the ambient flow for the n+1 time step + !---------------------------------------------------------------------------- + ! Populate low resolution grids based on ambient wind source + !---------------------------------------------------------------------------- + + select case (p%Mod_AmbWind) + + ! File-based ambient wind + case (1) + + ! Read from file the ambient flow for the n+1 time step call ReadLowResWindFile(n+1, p, m%Vamb_Low, errStat2, errMsg2); if (Failed()) return; - !$OMP PARALLEL DO DEFAULT(Shared) PRIVATE(nt, i_hl, errStat2, errMsg2) !Private(nt,tm2,tm3) + ! InflowWind-based ambient wind (single or multiple instances) + case (2, 3) + +!FIXME:merge5.0 remove next 3 lines + ! Set the hub position and orientation to pass to IfW (IfW always calculates hub and disk avg vel) -- note that this is garbage data. + m%u_IfW_Low%HubPosition = (/ p%X0_low + 0.5*p%nX_low*p%dX_low, p%Y0_low + 0.5*p%nY_low*p%dY_low, p%Z0_low + 0.5*p%nZ_low*p%dZ_low /) + call Eye(m%u_IfW_Low%HubOrientation,ErrStat2,ErrMsg2); if (Failed()) return; + + ! Calculate the low-resolution grid inflow velocities + call InflowWind_CalcOutput(t+p%dt_low, m%u_IfW_Low, p%IfW(0), x%IfW(0), xd%IfW(0), z%IfW(0), OtherState%IfW(0), m%y_IfW_Low, m%IfW(0), errStat2, errMsg2) + if (Failed()) return + + ! Transfer velocities to low resolution grid + V_Grid(lbound(m%Vamb_low,1):ubound(m%Vamb_low,1),& + lbound(m%Vamb_low,2):ubound(m%Vamb_low,2),& + lbound(m%Vamb_low,3):ubound(m%Vamb_low,3),& + lbound(m%Vamb_low,4):ubound(m%Vamb_low,4)) => m%y_IfW_Low%VelocityUVW + m%Vamb_Low = V_Grid + + end select + + !---------------------------------------------------------------------------- + ! Populate high-resolution grid based on ambient wind source + !---------------------------------------------------------------------------- + + select case (p%Mod_AmbWind) + + ! File-based ambient wind + case (1) + + !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(nt, i_hl, errStat2, errMsg2) & + !$OMP SHARED(p, n_high_low, n, m, errStat, errMsg, AbortErrLev) do nt = 1,p%NumTurbines + + ! Copy T=T_low_previous-DT_high (end-1 index in Vamb_high) into T=T_low_now-DT_high (0 index in Vamb_high). Note that n starts at -1 + if (n/=-1_IntKi) m%Vamb_high(nt)%data(:,:,:,:,0) = m%Vamb_high(nt)%data(:,:,:,:,ubound(m%Vamb_high(nt)%data,5)-1) + do i_hl=0, n_high_low - ! read from file the ambient flow for the current time step - call ReadHighResWindFile(nt, (n+1)*p%n_high_low + i_hl, p, m%Vamb_high(nt)%data(:,:,:,:,i_hl), errStat2, errMsg2) + + ! read from file the ambient flow for the current time step +!FIXME:merge5.0 replace next line with the following + call ReadHighResWindFile(nt, (n+1)*p%n_high_low + i_hl, p, m%Vamb_high(nt)%data(:,:,:,:,i_hl+1), errStat2, errMsg2) +! call ReadHighResWindVTK(nt, n*p%n_high_low + i_hl, p, m%Vamb_high(nt)%data(:,:,:,:,i_hl+1), errStat2, errMsg2) if (ErrStat2 >= AbortErrLev) then !$OMP CRITICAL ! Needed to avoid data race on ErrStat and ErrMsg call SetErrStat( ErrStat2, ErrMsg2, errStat, errMsg, RoutineName ) !$OMP END CRITICAL endif end do + + ! Special handling at T=0 for time slice at -DT_high (0 index in Vamb_high). Note that n starts at -1 + ! -> Copy T=0 data into T=-DT_high for AD extrap/interp + if (n==-1_IntKi) m%Vamb_high(nt)%data(:,:,:,:,0) = m%Vamb_high(nt)%data(:,:,:,:,1) + end do !$OMP END PARALLEL DO + if (errStat >= AbortErrLev) return - else ! p%Mod_AmbWind == 2 .or. 3 + ! Single InflowWind instance + case (2) + + ! Loop through turbines + do nt = 1, p%NumTurbines +!FIXME:merge5.0 remove next 4 lines ! Set the hub position and orientation to pass to IfW (IfW always calculates hub and disk avg vel) - m%u_IfW_Low%HubPosition = (/ p%X0_low + 0.5*p%nX_low*p%dX_low, p%Y0_low + 0.5*p%nY_low*p%dY_low, p%Z0_low + 0.5*p%nZ_low*p%dZ_low /) - call Eye(m%u_IfW_Low%HubOrientation,ErrStat2,ErrMsg2); if (Failed()) return; - ! Set low-resolution inflow wind velocities - call InflowWind_CalcOutput(t+p%dt_low, m%u_IfW_Low, p%IfW(0), x%IfW(0), xd%IfW(0), z%IfW(0), OtherState%IfW(0), m%y_IfW_Low, m%IfW(0), errStat2, errMsg2); if (Failed()) return; - c = 1 - do k = 0,p%nZ_low-1 - do j = 0,p%nY_low-1 - do i = 0,p%nX_low-1 - m%Vamb_Low(:,i,j,k) = m%y_IfW_Low%VelocityUVW(:,c) - c = c+1 - end do + m%u_IfW_High%HubPosition = (/ p%X0_high(nt) + 0.5*p%nX_high*p%dX_high(nt), p%Y0_high(nt) + 0.5*p%nY_high*p%dY_high(nt), p%Z0_high(nt) + 0.5*p%nZ_high*p%dZ_high(nt) /) + call Eye(m%u_IfW_High%HubOrientation,ErrStat2,ErrMsg2) + +!FIXME:merge5.0 remove next 3 lines + ! Set input position + m%u_IfW_High%PositionXYZ = p%Grid_high(:,:,nt) + + ! Copy T=T_low_previous-DT_high (end-1 index in Vamb_high) into T=T_low_now-DT_high (0 index in Vamb_high). Note that n starts at -1 + if (n/=-1_IntKi) m%Vamb_high(nt)%data(:,:,:,:,0) = m%Vamb_high(nt)%data(:,:,:,:,ubound(m%Vamb_high(nt)%data,5)-1) + + ! Loop through high resolution grids + do i_hl = 0, n_high_low + + ! Calculate wind velocities at grid locations from InflowWind +!FIXME:merge5.0 replace next line with the next + call InflowWind_CalcOutput(t+p%dt_low+i_hl*p%DT_high, m%u_IfW_High, p%IfW(0), x%IfW(0), xd%IfW(0), z%IfW(0), OtherState%IfW(0), m%y_IfW_High, m%IfW(0), errStat2, errMsg2) +! call IfW_FlowField_GetVelAcc(p%IfW(0)%FlowField, 1, t + i_hl*p%DT_high, & +! m%u_IfW_High(nt)%PositionXYZ, & +! m%y_IfW_High(nt)%VelocityUVW, AccUVW, errStat2, errMsg2) + if (Failed()) return + + ! Transfer velocities to high resolution grid +!FIXME:merge5.0 remove following 4 lines and uncomment block after + V_Grid(lbound(m%Vamb_high(nt)%data,1):ubound(m%Vamb_high(nt)%data,1),& + lbound(m%Vamb_high(nt)%data,2):ubound(m%Vamb_high(nt)%data,2),& + lbound(m%Vamb_high(nt)%data,3):ubound(m%Vamb_high(nt)%data,3),& + lbound(m%Vamb_high(nt)%data,4):ubound(m%Vamb_high(nt)%data,4)) => m%y_IfW_High%VelocityUVW +! V_Grid(lbound(m%Vamb_high(nt)%data,1):ubound(m%Vamb_high(nt)%data,1),& +! lbound(m%Vamb_high(nt)%data,2):ubound(m%Vamb_high(nt)%data,2),& +! lbound(m%Vamb_high(nt)%data,3):ubound(m%Vamb_high(nt)%data,3),& +! lbound(m%Vamb_high(nt)%data,4):ubound(m%Vamb_high(nt)%data,4)) => m%y_IfW_High(nt)%VelocityUVW + m%Vamb_high(nt)%data(:,:,:,:,i_hl+1) = V_Grid end do + + ! Special handling at T=0 for time slice at -DT_high (0 index in Vamb_high). Note that n starts at -1 + ! -> Copy T=0 data into T=-DT_high for AD extrap/interp + if (n==-1_IntKi) m%Vamb_high(nt)%data(:,:,:,:,0) = m%Vamb_high(nt)%data(:,:,:,:,1) + end do - ! Set the high-resoultion inflow wind velocities for each turbine - if ( p%Mod_AmbWind == 2 ) then - do nt = 1,p%NumTurbines - m%u_IfW_High%PositionXYZ = p%Grid_high(:,:,nt) - ! Set the hub position and orientation to pass to IfW (IfW always calculates hub and disk avg vel) - m%u_IfW_High%HubPosition = (/ p%X0_high(nt) + 0.5*p%nX_high*p%dX_high(nt), p%Y0_high(nt) + 0.5*p%nY_high*p%dY_high(nt), p%Z0_high(nt) + 0.5*p%nZ_high*p%dZ_high(nt) /) - call Eye(m%u_IfW_High%HubOrientation,ErrStat2,ErrMsg2) - do i_hl=0, n_high_low - call InflowWind_CalcOutput(t+p%dt_low+i_hl*p%DT_high, m%u_IfW_High, p%IfW(0), x%IfW(0), xd%IfW(0), z%IfW(0), OtherState%IfW(0), m%y_IfW_High, m%IfW(0), errStat2, errMsg2); if (Failed()) return; - c = 1 - do k = 0,p%nZ_high-1 - do j = 0,p%nY_high-1 - do i = 0,p%nX_high-1 - m%Vamb_high(nt)%data(:,i,j,k,i_hl) = m%y_IfW_High%VelocityUVW(:,c) - c = c+1 - end do - end do + + ! Multiple InflowWind instances (one per turbine) + case (3) + + ! Loop through turbines + do nt = 1, p%NumTurbines + +!FIXME:merge5.0 remove next 10 lines + ! Set input velocity + c = 1 + do k = 0,p%nZ_high-1 + do j = 0,p%nY_high-1 + do i = 0,p%nX_high-1 + m%u_IfW_High%PositionXYZ(:,c) = p%Grid_high(:,c,nt) - p%WT_Position(:,nt) + c = c+1 end do end do end do - else !p%Mod_AmbWind == 3 - do nt = 1,p%NumTurbines - c = 1 - do k = 0,p%nZ_high-1 - do j = 0,p%nY_high-1 - do i = 0,p%nX_high-1 - m%u_IfW_High%PositionXYZ(:,c) = p%Grid_high(:,c,nt) - p%WT_Position(:,nt) - c = c+1 - end do - end do - end do - do i_hl=0, n_high_low - ! Set the hub position and orientation to pass to IfW (IfW always calculates hub and disk avg vel) - m%u_IfW_High%HubPosition = (/ p%X0_high(nt) + 0.5*p%nX_high*p%dX_high(nt), p%Y0_high(nt) + 0.5*p%nY_high*p%dY_high(nt), p%Z0_high(nt) + 0.5*p%nZ_high*p%dZ_high(nt) /) - p%WT_Position(:,nt) - call Eye(m%u_IfW_High%HubOrientation,ErrStat2,ErrMsg2) - call InflowWind_CalcOutput(t+p%dt_low+i_hl*p%DT_high, m%u_IfW_High, p%IfW(nt), x%IfW(nt), xd%IfW(nt), z%IfW(nt), OtherState%IfW(nt), m%y_IfW_High, m%IfW(nt), errStat2, errMsg2); if (Failed()) return; - c = 1 - do k = 0,p%nZ_high-1 - do j = 0,p%nY_high-1 - do i = 0,p%nX_high-1 - m%Vamb_high(nt)%data(:,i,j,k,i_hl) = m%y_IfW_High%VelocityUVW(:,c) - c = c+1 - end do - end do - end do - end do + ! Copy T=T_low_previous-DT_high (end-1 index in Vamb_high) into T=T_low_now-DT_high (0 index in Vamb_high). Note that n starts at -1 + if (n/=-1_IntKi) m%Vamb_high(nt)%data(:,:,:,:,0) = m%Vamb_high(nt)%data(:,:,:,:,ubound(m%Vamb_high(nt)%data,5)-1) + + ! Loop through high resolution grids + do i_hl = 0, n_high_low + +!FIXME:merge5.0 remove next 4 lines + ! Set the hub position and orientation to pass to IfW (IfW always calculates hub and disk avg vel) + m%u_IfW_High%HubPosition = (/ p%X0_high(nt) + 0.5*p%nX_high*p%dX_high(nt), p%Y0_high(nt) + 0.5*p%nY_high*p%dY_high(nt), p%Z0_high(nt) + 0.5*p%nZ_high*p%dZ_high(nt) /) - p%WT_Position(:,nt) + call Eye(m%u_IfW_High%HubOrientation,ErrStat2,ErrMsg2) + + ! Calculate wind velocities at grid locations from InflowWind +!FIXME:merge5.0 replace next line with the next + call InflowWind_CalcOutput(t+p%dt_low+i_hl*p%DT_high, m%u_IfW_High, p%IfW(nt), x%IfW(nt), xd%IfW(nt), z%IfW(nt), OtherState%IfW(nt), m%y_IfW_High, m%IfW(nt), errStat2, errMsg2) +! call IfW_FlowField_GetVelAcc(p%IfW(nt)%FlowField, 1, t + i_hl*p%DT_high, & +! m%u_IfW_High(nt)%PositionXYZ, & +! m%y_IfW_High(nt)%VelocityUVW, AccUVW, errStat2, errMsg2) + if (Failed()) return + + ! Transfer velocities to high resolution grid +!FIXME:merge5.0 remove following 4 lines and uncomment block after + V_Grid(lbound(m%Vamb_high(nt)%data,1):ubound(m%Vamb_high(nt)%data,1),& + lbound(m%Vamb_high(nt)%data,2):ubound(m%Vamb_high(nt)%data,2),& + lbound(m%Vamb_high(nt)%data,3):ubound(m%Vamb_high(nt)%data,3),& + lbound(m%Vamb_high(nt)%data,4):ubound(m%Vamb_high(nt)%data,4)) => m%y_IfW_High%VelocityUVW +! V_Grid(lbound(m%Vamb_high(nt)%data,1):ubound(m%Vamb_high(nt)%data,1),& +! lbound(m%Vamb_high(nt)%data,2):ubound(m%Vamb_high(nt)%data,2),& +! lbound(m%Vamb_high(nt)%data,3):ubound(m%Vamb_high(nt)%data,3),& +! lbound(m%Vamb_high(nt)%data,4):ubound(m%Vamb_high(nt)%data,4)) => m%y_IfW_High(nt)%VelocityUVW + m%Vamb_high(nt)%data(:,:,:,:,i_hl+1) = V_Grid end do - end if - end if - ! WAT tracer propagation + ! Special handling at T=0 for time slice at -DT_high. Note that n starts at -1 + ! -> Copy T=0 data into T=-DT_high for AD extrap/interp + if (n==-1_IntKi) m%Vamb_high(nt)%data(:,:,:,:,0) = m%Vamb_high(nt)%data(:,:,:,:,1) + + end do + + end select + + !---------------------------------------------------------------------------- + ! Propagate WAT tracer + !---------------------------------------------------------------------------- + if (p%WAT_Enabled) then - ! find mean velocity of all turbine disks + + ! Find mean velocity of all turbine disks xd%Ufarm = 0.0_ReKi do nt=1,p%NumTurbines xd%Ufarm(1:3) = xd%Ufarm(1:3) + m%V_amb_low_disk(1:3,nt) enddo xd%Ufarm(1:3) = xd%Ufarm(1:3) / real(p%NumTurbines,ReKi) + ! add mean velocity * dt to the tracer for the position of the WAT box xd%WAT_B_Box(1:3) = xd%WAT_B_Box(1:3) + xd%Ufarm(1:3)*real(p%dt_low,ReKi) endif contains logical function Failed() - call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) Failed = ErrStat >= AbortErrLev end function Failed end subroutine AWAE_UpdateStates diff --git a/modules/awae/src/AWAE_IO.f90 b/modules/awae/src/AWAE_IO.f90 index 86ee0b77c5..17e79eed66 100644 --- a/modules/awae/src/AWAE_IO.f90 +++ b/modules/awae/src/AWAE_IO.f90 @@ -110,13 +110,14 @@ subroutine WriteDisWindFiles( n, WrDisSkp1, p, y, m, errStat, errMsg ) if (ErrStat >= AbortErrLev) return do nt= 1,p%NumTurbines - ! We are only writing out the first of the high res data for a given low res time step + ! 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" 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 - call WrVTK_SP_vectors3D( Un, "Velocity", (/p%nX_high,p%nY_high,p%nZ_high/), (/p%X0_high(nt),p%Y0_high(nt),p%Z0_high(nt)/), (/p%dX_high(nt),p%dY_high(nt),p%dZ_high(nt)/), y%Vdist_high(nt)%data(:,:,:,:,0), errStat2, errMsg2 ) + call WrVTK_SP_vectors3D( Un, "Velocity", (/p%nX_high,p%nY_high,p%nZ_high/), (/p%X0_high(nt),p%Y0_high(nt),p%Z0_high(nt)/), (/p%dX_high(nt),p%dY_high(nt),p%dZ_high(nt)/), y%Vdist_high(nt)%data(:,:,:,:,1), 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 e6afddb9d0..5627f31bf3 100644 --- a/modules/awae/src/AWAE_Registry.txt +++ b/modules/awae/src/AWAE_Registry.txt @@ -188,6 +188,7 @@ typedef ^ ParameterType ReKi Grid_low {:}{:} - - "XYZ typedef ^ ParameterType ReKi Grid_high {:}{:}{:} - - "XYZ components (global positions) of the spatial discretization of the high-resolution spatial domain for each turbine " m typedef ^ ParameterType ReKi WT_Position {:}{:} - - "X-Y-Z position of each wind turbine; index 1 = XYZ; index 2 = turbine number" meters typedef ^ ParameterType IntKi n_high_low - - - "Number of high-resolution time steps per low" - +typedef ^ ParameterType IntKi n_high_low_p1 - - - "Number of high-resolution time steps per low, plus one at t_low-dt_high" - typedef ^ ParameterType DbKi dt_low - - - "Low-resolution (FAST.Farm driver/glue code) time step" s typedef ^ ParameterType DbKi dt_high - - - "High-resolution (FAST) time step" s typedef ^ ParameterType IntKi NumDT - - - "Number of low-resolution (FAST.Farm driver/glue code) time steps" - diff --git a/modules/awae/src/AWAE_Types.f90 b/modules/awae/src/AWAE_Types.f90 index fbd5bde33a..3cbe84a334 100644 --- a/modules/awae/src/AWAE_Types.f90 +++ b/modules/awae/src/AWAE_Types.f90 @@ -211,6 +211,7 @@ MODULE AWAE_Types REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: Grid_high !< XYZ components (global positions) of the spatial discretization of the high-resolution spatial domain for each turbine [m] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: WT_Position !< X-Y-Z position of each wind turbine; index 1 = XYZ; index 2 = turbine number [meters] INTEGER(IntKi) :: n_high_low = 0_IntKi !< Number of high-resolution time steps per low [-] + INTEGER(IntKi) :: n_high_low_p1 = 0_IntKi !< Number of high-resolution time steps per low, plus one at t_low-dt_high [-] REAL(DbKi) :: dt_low = 0.0_R8Ki !< Low-resolution (FAST.Farm driver/glue code) time step [s] REAL(DbKi) :: dt_high = 0.0_R8Ki !< High-resolution (FAST) time step [s] INTEGER(IntKi) :: NumDT = 0_IntKi !< Number of low-resolution (FAST.Farm driver/glue code) time steps [-] @@ -1999,6 +2000,7 @@ subroutine AWAE_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%WT_Position = SrcParamData%WT_Position end if DstParamData%n_high_low = SrcParamData%n_high_low + DstParamData%n_high_low_p1 = SrcParamData%n_high_low_p1 DstParamData%dt_low = SrcParamData%dt_low DstParamData%dt_high = SrcParamData%dt_high DstParamData%NumDT = SrcParamData%NumDT @@ -2221,6 +2223,7 @@ subroutine AWAE_PackParam(RF, Indata) call RegPackAlloc(RF, InData%Grid_high) call RegPackAlloc(RF, InData%WT_Position) call RegPack(RF, InData%n_high_low) + call RegPack(RF, InData%n_high_low_p1) call RegPack(RF, InData%dt_low) call RegPack(RF, InData%dt_high) call RegPack(RF, InData%NumDT) @@ -2306,6 +2309,7 @@ subroutine AWAE_UnPackParam(RF, OutData) call RegUnpackAlloc(RF, OutData%Grid_high); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%WT_Position); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%n_high_low); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%n_high_low_p1); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%dt_low); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%dt_high); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%NumDT); if (RegCheckErr(RF, RoutineName)) return