diff --git a/modules/inflowwind/src/IfW_FlowField.f90 b/modules/inflowwind/src/IfW_FlowField.f90 index 922b831843..a8212790e1 100644 --- a/modules/inflowwind/src/IfW_FlowField.f90 +++ b/modules/inflowwind/src/IfW_FlowField.f90 @@ -78,7 +78,9 @@ subroutine IfW_FlowField_GetVelAcc(FF, IStart, Time, PositionXYZ, VelocityUVW, A ! Determine if acceleration should be calculated and returned OutputAccel = allocated(AccelUVW) - if (OutputAccel .and. .not. FF%AccFieldValid) then + ! Cubic velocity interpolation also requires a valid acceleration field, since its + ! formula uses the derivative data even when acceleration output is not requested. + if ((OutputAccel .or. FF%VelInterpCubic) .and. .not. FF%AccFieldValid) then call SetErrStat(ErrID_Fatal, "Accel output requested, but accel field is not valid", & ErrStat, ErrMsg, RoutineName) return @@ -224,8 +226,10 @@ subroutine IfW_FlowField_GetVelAcc(FF, IStart, Time, PositionXYZ, VelocityUVW, A ! Calculate grid cells for interpolation, returns velocity and acceleration ! components at corners of grid cell containing time and position. Also - ! returns interpolation values Xi. - call Grid3DField_GetCell(FF%Grid3D, Time, Position(:, i), OutputAccel, GridExceedAllow, & + ! returns interpolation values Xi. AccCell is required whenever cubic + ! velocity interpolation is used (its formula uses the derivative data), + ! not only when acceleration is explicitly requested as an output. + call Grid3DField_GetCell(FF%Grid3D, Time, Position(:, i), OutputAccel .or. FF%VelInterpCubic, GridExceedAllow, & VelCell, AccCell, Xi, Is3D, TmpErrStat, TmpErrMsg) if (TmpErrStat >= AbortErrLev) then call SetErrStat(TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName) @@ -1442,9 +1446,18 @@ subroutine IfW_Grid3DField_CalcAccel(G3D, ErrStat, ErrMsg) call SetErrStat(TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return - ! If number of time grids is 1 or 2, set all accelerations to zero, return - if (G3D%NTGrids < 3) then + ! If number of time steps is 1 or 2, set all accelerations to zero, return + if (G3D%NSteps < 3) then G3D%Acc = 0.0_SiKi + ! Also allocate/zero tower acceleration so GetCellInTower has a valid array to reference + if (G3D%NTGrids > 0) then + call AllocAry(G3D%AccTower, size(G3D%VelTower, dim=1), & + size(G3D%VelTower, dim=2), size(G3D%VelTower, dim=3), & + 'tower wind acceleration data.', TmpErrStat, TmpErrMsg) + call SetErrStat(TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return + G3D%AccTower = 0.0_SiKi + end if return end if @@ -1477,16 +1490,14 @@ subroutine IfW_Grid3DField_CalcAccel(G3D, ErrStat, ErrMsg) call SetErrStat(TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return - ! If number of time grids is 1 or 2, set all accelerations to zero - if (G3D%NTGrids < 3) then - G3D%Acc = 0.0_SiKi - else ! Otherwise, calculate acceleration at each grid point - do iz = 1, G3D%NTGrids - do ic = 1, G3D%NComp - call CalcCubicSplineDeriv(G3D%NSteps, G3D%DTime, G3D%VelTower(ic, iz, :), G3D%AccTower(ic, iz, :)) - end do + ! Calculate acceleration at each tower grid point. NSteps < 3 is already handled by the + ! early return above, and each tower height is interpolated independently in time, so no + ! minimum tower height count is required here. + do iz = 1, G3D%NTGrids + do ic = 1, G3D%NComp + call CalcCubicSplineDeriv(G3D%NSteps, G3D%DTime, G3D%VelTower(ic, iz, :), G3D%AccTower(ic, iz, :)) end do - end if + end do contains diff --git a/modules/inflowwind/tests/inflowwind_utest.F90 b/modules/inflowwind/tests/inflowwind_utest.F90 index 758c25de13..30c2a5e258 100644 --- a/modules/inflowwind/tests/inflowwind_utest.F90 +++ b/modules/inflowwind/tests/inflowwind_utest.F90 @@ -3,6 +3,7 @@ program inflowwind_utest use testdrive, only: run_testsuite, new_testsuite, testsuite_type use test_bladed_wind, only: test_bladed_wind_suite +use test_grid3d_field, only: test_grid3d_field_suite use test_hawc_wind, only: test_hawc_wind_suite use test_outputs, only: test_outputs_suite use test_steady_wind, only: test_steady_wind_suite @@ -21,6 +22,7 @@ program inflowwind_utest testsuites = [ & new_testsuite("Bladed Wind", test_bladed_wind_suite), & + new_testsuite("Grid3D Field", test_grid3d_field_suite), & new_testsuite("HAWC Wind", test_hawc_wind_suite), & new_testsuite("Outputs", test_outputs_suite), & new_testsuite("Steady Wind", test_steady_wind_suite), & diff --git a/modules/inflowwind/tests/test_grid3d_field.F90 b/modules/inflowwind/tests/test_grid3d_field.F90 new file mode 100644 index 0000000000..16a0441e43 --- /dev/null +++ b/modules/inflowwind/tests/test_grid3d_field.F90 @@ -0,0 +1,218 @@ +module test_grid3d_field + +use, intrinsic :: ieee_arithmetic, only: ieee_is_finite +use testdrive, only: new_unittest, unittest_type, error_type, check +use ifw_test_tools +use IfW_FlowField +use IfW_FlowField_Types +use NWTC_Library + +implicit none +private +public :: test_grid3d_field_suite + +integer(IntKi), parameter :: NY = 4, NZ = 4, NT = 10 +real(ReKi), parameter :: DTIME = 0.1_ReKi + +contains + +!> Collect all exported unit tests +subroutine test_grid3d_field_suite(testsuite) + type(unittest_type), allocatable, intent(out) :: testsuite(:) + testsuite = [ & + new_unittest("test_grid3d_cubic_vel_only", test_grid3d_cubic_vel_only), & + new_unittest("test_grid3d_calcaccel_no_tower", test_grid3d_calcaccel_no_tower), & + new_unittest("test_grid3d_calcaccel_few_tower_points", test_grid3d_calcaccel_few_tower_points) & + ] +end subroutine + +!> Reproduces a bug in IfW_FlowField_GetVelAcc: when VelInterpCubic is enabled and the +!! caller does not request acceleration output (AccelUVW left unallocated), the local +!! AccCell array used by the cubic Hermite velocity formula is never populated, so the +!! returned velocity is computed from uninitialized memory. This test builds a minimal +!! Grid3D flow field with a velocity that is a known linear ramp in time (spatially +!! uniform), so the exact correct answer at any query time is known, and checks that +!! querying velocity without requesting acceleration gives the same (finite, correct) +!! answer as querying with acceleration requested. +subroutine test_grid3d_cubic_vel_only(error) + type(error_type), allocatable, intent(out) :: error + + type(FlowFieldType) :: FF + real(ReKi), allocatable :: Position(:, :), VelNoAcc(:, :), VelWithAcc(:, :) + real(ReKi), allocatable :: AccelUVW(:, :) + integer(IntKi) :: it, iy, iz + integer(IntKi) :: TmpErrStat + character(ErrMsgLen) :: TmpErrMsg + real(DbKi) :: QueryTime + real(ReKi) :: Expected + + ! Build minimal Grid3D field: NY x NZ spatial points, NT time steps, no tower grid. + FF%FieldType = Grid3D_FieldType + FF%VelInterpCubic = .true. + FF%RotateWindBox = .false. + + FF%Grid3D%NComp = 3 + FF%Grid3D%NYGrids = NY + FF%Grid3D%NZGrids = NZ + FF%Grid3D%NTGrids = 0 + FF%Grid3D%NSteps = NT + FF%Grid3D%DTime = DTIME + FF%Grid3D%Rate = 1.0_ReKi/DTIME + FF%Grid3D%YHWid = 5.0_ReKi + FF%Grid3D%ZHWid = 5.0_ReKi + FF%Grid3D%GridBase = 5.0_ReKi + FF%Grid3D%InvDY = real(NY - 1, ReKi)/(2.0_ReKi*FF%Grid3D%YHWid) + FF%Grid3D%InvDZ = real(NZ - 1, ReKi)/(2.0_ReKi*FF%Grid3D%ZHWid) + FF%Grid3D%MeanWS = 8.0_ReKi + FF%Grid3D%InvMWS = 1.0_ReKi/FF%Grid3D%MeanWS + FF%Grid3D%InitXPosition = 0.0_ReKi + FF%Grid3D%TotalTime = real(NT - 1, ReKi)*DTIME + FF%Grid3D%Periodic = .false. + FF%Grid3D%InterpTower = .false. + + allocate (FF%Grid3D%Vel(3, NY, NZ, NT)) + FF%Grid3D%Vel = 0.0_SiKi + ! U component is a spatially-uniform linear ramp in time: U(t) = t (seconds -> m/s) + do it = 1, NT + do iz = 1, NZ + do iy = 1, NY + FF%Grid3D%Vel(1, iy, iz, it) = real((it - 1), SiKi)*real(DTIME, SiKi) + end do + end do + end do + + ! Set the true time-derivative (dU/dt = 1, dV/dt = dW/dt = 0) directly rather than via + ! IfW_Grid3DField_CalcAccel so this test isolates the IfW_FlowField_GetVelAcc behavior. + ! (test_grid3d_calcaccel_no_tower below covers IfW_Grid3DField_CalcAccel itself.) + allocate (FF%Grid3D%Acc(3, NY, NZ, NT)) + FF%Grid3D%Acc = 0.0_SiKi + FF%Grid3D%Acc(1, :, :, :) = 1.0_SiKi + FF%AccFieldValid = .true. + + ! Query at a position centered in the grid, at a time between samples + allocate (Position(3, 1)) + Position(:, 1) = [0.0_ReKi, 0.0_ReKi, FF%Grid3D%GridBase + FF%Grid3D%ZHWid] + QueryTime = 0.25_DbKi + Expected = real(QueryTime, ReKi) + + ! Case 1: acceleration NOT requested (AccelUVW left unallocated) - buggy path + allocate (VelNoAcc(3, 1)) + call IfW_FlowField_GetVelAcc(FF, 1, QueryTime, Position, VelNoAcc, AccelUVW, TmpErrStat, TmpErrMsg) + call check(error, TmpErrStat, ErrID_None, message='GetVelAcc (no accel) error: '//trim(TmpErrMsg)); if (allocated(error)) return + + call check(error, ieee_is_finite(VelNoAcc(1, 1)), message='Velocity (no accel requested) is not finite (NaN/Inf)'); if (allocated(error)) return + call check(error, VelNoAcc(1, 1), Expected, thr=1.0e-3_ReKi); if (allocated(error)) return + + ! Case 2: acceleration requested (AccelUVW allocated) - reference path + allocate (VelWithAcc(3, 1)) + allocate (AccelUVW(3, 1)) + call IfW_FlowField_GetVelAcc(FF, 1, QueryTime, Position, VelWithAcc, AccelUVW, TmpErrStat, TmpErrMsg) + call check(error, TmpErrStat, ErrID_None, message='GetVelAcc (with accel) error: '//trim(TmpErrMsg)); if (allocated(error)) return + + call check(error, ieee_is_finite(VelWithAcc(1, 1)), message='Velocity (accel requested) is not finite (NaN/Inf)'); if (allocated(error)) return + call check(error, VelWithAcc(1, 1), Expected, thr=1.0e-3_ReKi); if (allocated(error)) return + + ! The two calls must agree: requesting acceleration must not change the velocity result + call check(error, VelNoAcc(1, 1), VelWithAcc(1, 1), thr=1.0e-6_ReKi); if (allocated(error)) return + +end subroutine + +!> Reproduces a bug in IfW_Grid3DField_CalcAccel: it checks G3D%NTGrids (the tower-grid +!! point count) instead of G3D%NSteps (the number of time samples) to decide whether to +!! compute real cubic-spline time derivatives. For the common case of no tower file +!! (NTGrids=0), this unconditionally zeroes G3D%Acc regardless of how many time steps +!! are actually available, silently defeating cubic-in-time interpolation. This test +!! uses a spatially-uniform linear velocity ramp in time (dU/dt = 1 exactly), so the +!! correct computed derivative is known, and checks that it is recovered when there is +!! no tower grid but plenty of time steps. +subroutine test_grid3d_calcaccel_no_tower(error) + type(error_type), allocatable, intent(out) :: error + + type(Grid3DFieldType) :: G3D + integer(IntKi) :: it, iy, iz + integer(IntKi) :: TmpErrStat + character(ErrMsgLen) :: TmpErrMsg + + G3D%NComp = 3 + G3D%NYGrids = NY + G3D%NZGrids = NZ + G3D%NTGrids = 0 + G3D%NSteps = NT + G3D%DTime = DTIME + G3D%Periodic = .false. + + allocate (G3D%Vel(3, NY, NZ, NT)) + G3D%Vel = 0.0_SiKi + do it = 1, NT + do iz = 1, NZ + do iy = 1, NY + G3D%Vel(1, iy, iz, it) = real((it - 1), SiKi)*real(DTIME, SiKi) + end do + end do + end do + + call IfW_Grid3DField_CalcAccel(G3D, TmpErrStat, TmpErrMsg) + call check(error, TmpErrStat, ErrID_None, message='CalcAccel error: '//trim(TmpErrMsg)); if (allocated(error)) return + + ! Interior time points should recover close to the true derivative (dU/dt = 1); a + ! natural cubic spline has some boundary-driven error, but zero here would indicate + ! the NTGrids-vs-NSteps bug has regressed (computation skipped entirely). + call check(error, real(G3D%Acc(1, 2, 2, 5), ReKi), 1.0_ReKi, thr=0.01_ReKi); if (allocated(error)) return + +end subroutine + +!> Reproduces the mirror-image bug in IfW_Grid3DField_CalcAccel's tower branch: it checked +!! G3D%NTGrids (tower height count) against the same <3 threshold meant for time-step counts, +!! zeroing out tower acceleration whenever there were only 1 or 2 tower grid points, even though +!! each tower height's time derivative is computed independently and only needs G3D%NSteps >= 3 +!! (already guaranteed by the earlier NSteps<3 early return). This test uses NTGrids=2 with a +!! spatially-uniform linear velocity ramp in time (dU/dt = 1 exactly) and checks that the tower +!! acceleration is actually computed rather than silently zeroed. +subroutine test_grid3d_calcaccel_few_tower_points(error) + type(error_type), allocatable, intent(out) :: error + + integer(IntKi), parameter :: NTWR = 2 + type(Grid3DFieldType) :: G3D + integer(IntKi) :: it, iy, iz + integer(IntKi) :: TmpErrStat + character(ErrMsgLen) :: TmpErrMsg + + G3D%NComp = 3 + G3D%NYGrids = NY + G3D%NZGrids = NZ + G3D%NTGrids = NTWR + G3D%NSteps = NT + G3D%DTime = DTIME + G3D%Periodic = .false. + + allocate (G3D%Vel(3, NY, NZ, NT)) + G3D%Vel = 0.0_SiKi + do it = 1, NT + do iz = 1, NZ + do iy = 1, NY + G3D%Vel(1, iy, iz, it) = real((it - 1), SiKi)*real(DTIME, SiKi) + end do + end do + end do + + allocate (G3D%VelTower(3, NTWR, NT)) + G3D%VelTower = 0.0_SiKi + do it = 1, NT + do iz = 1, NTWR + G3D%VelTower(1, iz, it) = real((it - 1), SiKi)*real(DTIME, SiKi) + end do + end do + + call IfW_Grid3DField_CalcAccel(G3D, TmpErrStat, TmpErrMsg) + call check(error, TmpErrStat, ErrID_None, message='CalcAccel error: '//trim(TmpErrMsg)); if (allocated(error)) return + + call check(error, allocated(G3D%AccTower), message='AccTower was not allocated'); if (allocated(error)) return + + ! Interior time point at either tower height should recover close to the true derivative + ! (dU/dt = 1); zero here would indicate the NTGrids<3 tower guard bug has regressed. + call check(error, real(G3D%AccTower(1, 1, 5), ReKi), 1.0_ReKi, thr=0.01_ReKi); if (allocated(error)) return + call check(error, real(G3D%AccTower(1, 2, 5), ReKi), 1.0_ReKi, thr=0.01_ReKi); if (allocated(error)) return + +end subroutine + +end module diff --git a/reg_tests/r-test b/reg_tests/r-test index d376069e24..f8d1d131de 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit d376069e24765fc3e6456dbf16ae39e6fd8e78ba +Subproject commit f8d1d131debad662dffa749bf37b440f8320095f diff --git a/unit_tests/CMakeLists.txt b/unit_tests/CMakeLists.txt index c441a90a20..c8847a0957 100644 --- a/unit_tests/CMakeLists.txt +++ b/unit_tests/CMakeLists.txt @@ -51,6 +51,7 @@ add_test(NAME beamdyn_utest COMMAND beamdyn_utest) add_executable(inflowwind_utest ${PROJECT_SOURCE_DIR}/modules/inflowwind/tests/inflowwind_utest.F90 ${PROJECT_SOURCE_DIR}/modules/inflowwind/tests/test_bladed_wind.F90 + ${PROJECT_SOURCE_DIR}/modules/inflowwind/tests/test_grid3d_field.F90 ${PROJECT_SOURCE_DIR}/modules/inflowwind/tests/test_hawc_wind.F90 ${PROJECT_SOURCE_DIR}/modules/inflowwind/tests/test_outputs.F90 ${PROJECT_SOURCE_DIR}/modules/inflowwind/tests/test_steady_wind.F90