From 117e1d4bab3874a8bbe73584a2d0c8718dcd88c6 Mon Sep 17 00:00:00 2001 From: Emmanuel Branlard Date: Wed, 9 Aug 2023 19:10:56 -0600 Subject: [PATCH 1/2] FF: fix plane output for wakedynamics --- glue-codes/fast-farm/src/FAST_Farm_Subs.f90 | 11 +++ modules/wakedynamics/src/WakeDynamics.f90 | 76 +++++++++++++-------- 2 files changed, 59 insertions(+), 28 deletions(-) diff --git a/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 b/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 index 8fa232a625..6450e4e05a 100644 --- a/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 +++ b/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 @@ -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 diff --git a/modules/wakedynamics/src/WakeDynamics.f90 b/modules/wakedynamics/src/WakeDynamics.f90 index afa31e4304..9360723a0d 100644 --- a/modules/wakedynamics/src/WakeDynamics.f90 +++ b/modules/wakedynamics/src/WakeDynamics.f90 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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) @@ -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) @@ -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 = "" @@ -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 @@ -1616,13 +1609,39 @@ 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 + 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' @@ -1647,11 +1666,11 @@ 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) + call vtk_dataset_structured_points((/xd%p_plane(1,i),-p%dr*p%NumRadii,-p%dr*p%NumRadii/),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) @@ -1664,11 +1683,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 @@ -1691,10 +1713,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' @@ -1779,8 +1801,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 From 2a743d0ca593027b53a26ea34ff13b2135c1d1b8 Mon Sep 17 00:00:00 2001 From: Emmanuel Branlard Date: Wed, 9 Aug 2023 19:51:40 -0600 Subject: [PATCH 2/2] FF: plane outputs based on p_hub --- modules/wakedynamics/src/WakeDynamics.f90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/modules/wakedynamics/src/WakeDynamics.f90 b/modules/wakedynamics/src/WakeDynamics.f90 index 9360723a0d..017d2e00e7 100644 --- a/modules/wakedynamics/src/WakeDynamics.f90 +++ b/modules/wakedynamics/src/WakeDynamics.f90 @@ -1634,6 +1634,7 @@ subroutine WD_WritePlaneOutputs( t, u, p, x, xd, z, OtherState, y, m, errStat, e character(1024) :: Filename type(VTK_Misc) :: mvtk real(ReKi), dimension(3) :: dx + real(ReKi), dimension(3) :: x0 errStat = ErrID_None errMsg = "" @@ -1666,7 +1667,10 @@ subroutine WD_WritePlaneOutputs( t, u, p, x, xd, z, OtherState, y, m, errStat, e dx(1) = 0.0 dx(2) = p%dr dx(3) = p%dr - call vtk_dataset_structured_points((/xd%p_plane(1,i),-p%dr*p%NumRadii,-p%dr*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(y%Vx_wake2(:,:,i),'Vx',mvtk) call vtk_point_data_scalar(y%Vy_wake2(:,:,i),'Vy',mvtk)