Skip to content

FF: fix plane output for wakedynamics #1729

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ebranlard merged 2 commits into from
Aug 11, 2023
Merged

FF: fix plane output for wakedynamics #1729

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
117e1d4
FF: fix plane output for wakedynamics
ebranlard Aug 10, 2023
2a743d0
FF: plane outputs based on p_hub
ebranlard Aug 10, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 11 additions & 0 deletions glue-codes/fast-farm/src/FAST_Farm_Subs.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2710,6 +2710,17 @@ subroutine FARM_CalcOutput(t, farm, ErrStat, ErrMsg)
!$OMP END PARALLEL DO
if (ErrStat >= AbortErrLev) return

! IO operation, not done using OpenMP
DO nt = 1,farm%p%NumTurbines
call WD_WritePlaneOutputs( t, farm%WD(nt)%u, farm%WD(nt)%p, farm%WD(nt)%x, farm%WD(nt)%xd, farm%WD(nt)%z, &
farm%WD(nt)%OtherSt, farm%WD(nt)%y, farm%WD(nt)%m, ErrStat2, ErrMsg2 )
if (ErrStat2 >= AbortErrLev) then
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'T'//trim(num2lstr(nt))//':'//RoutineName)
endif
END DO
if (ErrStat >= AbortErrLev) return


call Transfer_WD_to_AWAE(farm)

if ( farm%p%UseSC ) then
Expand Down
80 changes: 52 additions & 28 deletions modules/wakedynamics/src/WakeDynamics.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ module WakeDynamics
public :: WD_UpdateStates ! Loose coupling routine for solving for constraint states, integrating
! continuous states, and updating discrete states
public :: WD_CalcOutput ! Routine for computing outputs
public :: WD_WritePlaneOutputs ! Routine for IO Operation
public :: WD_CalcConstrStateResidual ! Tight coupling routine for returning the constraint state residual

public :: WD_TEST_Axi2Cart
Expand Down Expand Up @@ -427,7 +428,6 @@ subroutine WD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut
! Define parameters
!............................................................................................
p%TurbNum = InitInp%TurbNum
p%OutFileRoot = InitInp%OutFileRoot
p%DT_low = interval
! Parameters from input file
p%Mod_Wake = InitInp%InputFileData%Mod_Wake
Expand Down Expand Up @@ -479,10 +479,9 @@ subroutine WD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut
end do

! Path for VTK outputs
call GetPath( p%OutFileRoot, rootDir, baseName )
call GetPath( InitInp%OutFileRoot, rootDir, baseName )
p%OutFileRoot = baseName
p%OutFileVTKDir = trim(rootDir) // 'vtk_ff_planes'



p%filtParam = exp(-2.0_ReKi*pi*p%dt_low*InitInp%InputFileData%f_c)
p%oneMinusFiltParam = 1.0_ReKi - p%filtParam
Expand Down Expand Up @@ -1224,7 +1223,6 @@ subroutine AddSwirl(r, Vt_wake, y, z, Vy_curl, Vz_curl)
real(ReKi), dimension(:,:), intent(inout) :: Vy_curl !< Curl velocity in the y direction (m/s)
real(ReKi), dimension(:,:), intent(inout) :: Vz_curl !< Curl velocity in the z direction (m/s)

real(ReKi) :: alpha
integer(IntKi) :: iz, iy, iLow, nr
real(ReKi) :: r_tmp, r_max
real(ReKi) :: Vt, S, C ! Sine and cosine
Expand Down Expand Up @@ -1257,9 +1255,7 @@ end subroutine AddSwirl
subroutine WD_TEST_AddVelocityCurl()

real(ReKi) :: Vy_curl(2,2)=0.0_ReKi
real(ReKi) :: Vy_curl_ref(2,2)
real(ReKi) :: Vz_curl(2,2)=0.0_ReKi
real(ReKi) :: Vz_curl_ref(2,2)
real(ReKi) :: y(2)=(/ 0., 2./)
real(ReKi) :: z(2)=(/-1.,1./)
real(ReKi) :: Gamma0
Expand Down Expand Up @@ -1287,7 +1283,6 @@ subroutine filter_angles2(psi_filt, chi_filt, psi, chi, alpha, alpha_bar)
real(ReKi), intent(in) :: chi !< skew angle
real(ReKi), intent(in) :: alpha !< filter weight
real(ReKi), intent(in) :: alpha_bar !< 1-alpha
real(ReKi) :: t_filt !< output
real(ReKi) :: x,y
real(ReKi) :: lambda(3,2)
real(ReKi) :: theta_out(3)
Expand Down Expand Up @@ -1391,7 +1386,6 @@ subroutine Axisymmetric2CartesianVx(Vx_axi, r, y, z, Vx)
real(ReKi), dimension(:,:), intent(inout) :: Vx !< Axial velocity, distributed across the plane (m/s)
integer(IntKi) :: iz, iy, nr, iLow
real(ReKi) :: r_tmp, r_max
real(ReKi) :: Vr, S, C ! Sine and cosine
nr = size(r)
r_max = r(nr)
do iz = 1,size(z)
Expand Down Expand Up @@ -1498,16 +1492,12 @@ subroutine WD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg )
CHARACTER(*), INTENT( OUT) :: errMsg !< Error message if errStat /= ErrID_None


integer, parameter :: indx = 1
integer(intKi) :: n, i, iy, iz, maxPln
integer(intKi) :: ErrStat2
character(ErrMsgLen) :: ErrMsg2
character(*), parameter :: RoutineName = 'WD_CalcOutput'
real(ReKi) :: correction(3)
real(ReKi) :: C, S, dvdr, dvdtheta_r, R, r_tmp
character(1024) :: Filename
type(VTK_Misc) :: mvtk
real(ReKi), dimension(3) :: dx
errStat = ErrID_None
errMsg = ""

Expand Down Expand Up @@ -1550,7 +1540,10 @@ subroutine WD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg )
! Misc approx
m%Ct_avg = get_Ctavg(p%r, u%Ct_azavg, u%D_rotor)
m%GammaCurl = u%D_Rotor/2. * u%Vx_wind_disk * m%Ct_avg * sin(u%chi_skew) * cos(u%chi_skew)


if ( p%OutAllPlanes ) then
call mkdir(p%OutFileVTKDir)
endif

else
y%x_plane = xd%x_plane
Expand Down Expand Up @@ -1616,13 +1609,40 @@ subroutine WD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg )

end if

end subroutine WD_CalcOutput

!----------------------------------------------------------------------------------------------------------------------------------
!>
subroutine WD_WritePlaneOutputs( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg )
use VTK !
REAL(DbKi), INTENT(IN ) :: t !< Current simulation time in seconds
TYPE(WD_InputType), INTENT(IN ) :: u !< Inputs at Time t
TYPE(WD_ParameterType), INTENT(IN ) :: p !< Parameters
TYPE(WD_ContinuousStateType), INTENT(IN ) :: x !< Continuous states at t
TYPE(WD_DiscreteStateType), INTENT(IN ) :: xd !< Discrete states at t
TYPE(WD_ConstraintStateType), INTENT(IN ) :: z !< Constraint states at t
TYPE(WD_OtherStateType), INTENT(IN ) :: OtherState !< Other states at t
TYPE(WD_OutputType), INTENT(IN ) :: y !< Outputs computed at t (Input only so that mesh con-
type(WD_MiscVarType), intent(IN ) :: m !< Misc/optimization variables
INTEGER(IntKi), INTENT( OUT) :: errStat !< Error status of the operation
CHARACTER(*), INTENT( OUT) :: errMsg !< Error message if errStat /= ErrID_None
integer(intKi) :: n, i
integer(intKi) :: ErrStat2
character(ErrMsgLen) :: ErrMsg2
character(*), parameter :: RoutineName = 'WD_WritePlaneOutputs'
real(ReKi) :: correction(3)
character(1024) :: Filename
type(VTK_Misc) :: mvtk
real(ReKi), dimension(3) :: dx
real(ReKi), dimension(3) :: x0
errStat = ErrID_None
errMsg = ""

n = nint(t/p%DT_low)
! --- VTK outputs per plane
if (p%OutAllPlanes) then
call vtk_misc_init(mvtk)
call set_vtk_binary_format(.false., mvtk)
if ( OtherState%firstPass ) then
call MKDIR(p%OutFileVTKDir)
endif
do i = 0, min(n-1,p%NumPlanes-1), 1
! if (EqualRealNos(t,0.0_DbKi) ) then
! write(Filename,'(A,I4.4,A)') trim(p%OutFileVTKDir)//'/PlaneOutputsAtPlane_',i,'_Init.vtk'
Expand All @@ -1647,11 +1667,14 @@ subroutine WD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg )
dx(1) = 0.0
dx(2) = p%dr
dx(3) = p%dr
call vtk_dataset_structured_points((/xd%p_plane(1,i),-dx*p%NumRadii,-dx*p%NumRadii/),dx,(/1,p%NumRadii*2-1,p%NumRadii*2-1/),mvtk)
x0(1) = xd%p_plane(1,i)
x0(2) = xd%p_plane(2,i) - p%dr*p%NumRadii
x0(3) = xd%p_plane(3,i) - p%dr*p%NumRadii
call vtk_dataset_structured_points(x0, dx, (/1,p%NumRadii*2-1,p%NumRadii*2-1/),mvtk)
call vtk_point_data_init(mvtk)
call vtk_point_data_scalar(xd%Vx_wake2(:,:,i),'Vx',mvtk)
call vtk_point_data_scalar(xd%Vy_wake2(:,:,i),'Vy',mvtk)
call vtk_point_data_scalar(xd%Vz_wake2(:,:,i),'Vz',mvtk)
call vtk_point_data_scalar(y%Vx_wake2(:,:,i),'Vx',mvtk)
call vtk_point_data_scalar(y%Vy_wake2(:,:,i),'Vy',mvtk)
call vtk_point_data_scalar(y%Vz_wake2(:,:,i),'Vz',mvtk)
call vtk_point_data_scalar(m%vt_amb2(:,:,i),'vt_amb2', mvtk)
call vtk_point_data_scalar(m%vt_shr2(:,:,i),'vt_shr2', mvtk)
call vtk_point_data_scalar(m%vt_tot2(:,:,i),'vt_tot2', mvtk)
Expand All @@ -1664,11 +1687,14 @@ subroutine WD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg )
call vtk_point_data_scalar(y%WAT_k_mt(:,:,i),'k_mt', mvtk)
endif
call vtk_close_file(mvtk)
else
call SetErrStat(ErrID_Fatal, '[INFO] Failed to write: '//trim(filename), errStat, errMsg, RoutineName)
endif
enddo ! loop on planes
endif

end subroutine WD_CalcOutput
end subroutine WD_WritePlaneOutputs


!----------------------------------------------------------------------------------------------------------------------------------
!> Tight coupling routine for solving for the residual of the constraint state equations
Expand All @@ -1691,10 +1717,10 @@ subroutine WD_CalcConstrStateResidual( Time, u, p, x, xd, z, OtherState, m, z_re


! Local variables
integer, parameter :: indx = 1
integer(intKi) :: ErrStat2
character(ErrMsgLen) :: ErrMsg2
character(*), parameter :: RoutineName = 'WD_CalcConstrStateResidual'
!integer, parameter :: indx = 1
!integer(intKi) :: ErrStat2
!character(ErrMsgLen) :: ErrMsg2
!character(*), parameter :: RoutineName = 'WD_CalcConstrStateResidual'



Expand Down Expand Up @@ -1779,8 +1805,6 @@ SUBROUTINE ValidateInitInputData( DT_low, InitInp, InputFileData, errStat, errMs


! local variables
integer(IntKi) :: k ! Blade number
integer(IntKi) :: j ! node number
character(*), parameter :: RoutineName = 'ValidateInitInputData'

errStat = ErrID_None
Expand Down