Skip to content

Fix InflowWind Grid3DField fatal error for round-off just below t = 0 - #3466

Open
RBergua wants to merge 3 commits into
OpenFAST:rc-5.0.1from
RBergua:InflowWind_Grid3DField
Open

RBergua wants to merge 3 commits into
OpenFAST:rc-5.0.1from
RBergua:InflowWind_Grid3DField

Conversation

@RBergua

@RBergua RBergua commented Sep 15, 2026

Copy link
Copy Markdown
Contributor

Feature or improvement description
In IfW_FlowField.f90, the subroutine GetBoundsT can throw a fatal error if we are trying to access a time close to 0 but slightly negative due to a floating-point rounding.

For example, running an AeroDyn standalone model with 2 wind turbines using OLAF and periodic turbulent wind, I got:

Time: 3 of 40.002 seconds.  Estimated final completion at 02:46:03 (in 1.329 days).
Dvr_TimeStep:ADI_UpdateStates:AD_UpdateStates:AD_CalcWind:IfW_FlowField_GetVelAcc:Grid3DField_GetCell: Error: GF wind array was exhausted at 2.784 seconds (trying to access data
at -2.62260E-06 seconds). IT_Lo=0, IT_HI=1

The value that makes the solver crash is very small (-2.62260E-06 s). The flag IT_LO=0 confirms that a point is asking for data before t = 0. Interestingly, this error only occurs when compiling in single precision (e.g., using GitHub Actions default). When building in double precision, there is no problem.

Looking at GetBoundsT, we can see that MODULO is applied twice to try to avoid this condition:

! In distance, X: InputInfo%PosX - p%InitXPosition - TIME*p%MeanWS
TimeShifted = real(Time, ReKi) + (G3D%InitXPosition - PosX)*G3D%InvMWS
! Get position on T grid
T_GRID = TimeShifted*G3D%Rate
! If field is periodic
if (G3D%Periodic) then
     ! Take modulus of negative grid to get positive value between 0 and NSteps
     T_GRID = MODULO(T_GRID, real(G3D%NSteps, ReKi))
     ! For very small negative numbers, the above modulus will return exactly NSteps
     ! so take modulus again to ensure that T_GRID is less than NSteps
     T_GRID = MODULO(T_GRID, real(G3D%NSteps, ReKi))
end if

To confirm this was the actual issue, I added a debug print before the fatal error:
print *, 'DEBUG: Periodic=', G3D%Periodic, ' T_GRID=', T_GRID, ' NSteps=', G3D%NSteps
which returned:
DEBUG: Periodic= T T_GRID= -2.4414062E-04 NSteps= 4000

So even with Periodic = .true., T_GRID could come out negative after both MODULO calls when using single precision.

Looking at the code, GetBoundsY and GetBoundsZ follow a similar logic to GetBoundsT. But those subroutines have an additional condition for the lower bound with a grid tolerance (e.g., IY_LO == 0 .and. DY >= 1.0_ReKi - GridTol). I followed a similar approach for IT_LO == 0, so a point within GridTol (e.g., 0.001) of t = 0 s at the negative side is clamped onto the first time sample instead of raising a fatal error:

! Adjust indices and interpolant 
if (IT_LO >= 1 .and. IT_HI <= G3D%NSteps) then 
	! Point is within grid 
else if (IT_LO == 0 .and. T_GRID >= -GridTol) then
	! Within tolerance of the first time step 
	IT_LO = 1
	IT_HI = 1 
	DT = -1.0_ReKi 
else if (IT_LO == G3D%NSteps) then
    ...

@stefanoc1396

Copy link
Copy Markdown

I found the same issue when running the standalone aerodyn version compiled in 32bits. This change solved the error.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants