Fix InflowWind Grid3D cubic velocity interpolation NaN bug

[andrew-platt]

Ready to merge

Feature or improvement description

This is a bug fix targeting rc-5.0.1, not a feature. It fixes long-standing defects in InflowWind's Grid3D flow field (TurbSim/Bladed/HAWC full-field wind files, WindType=3 and similar) that produce NaN wind velocities whenever VelInterpCubic = True is set in the InflowWind input file and the calling code queries velocity without also requesting acceleration output. FAST.Farm's AWAE module is exactly such a caller: it samples ambient wind velocity for its high-resolution grid without asking for acceleration, so any FAST.Farm case with VelInterpCubic = True returns a 100%-NaN ambient wind field and fails immediately at t=0 with:

FARM_InitialCO: ... the rotor plane has left the low-resolution domain

Root cause

IfW_FlowField_GetVelAcc (in modules/inflowwind/src/IfW_FlowField.f90) determines whether the caller wants acceleration output from whether it allocated the optional AccelUVW output array:

OutputAccel = allocated(AccelUVW)

For the Grid3D_FieldType case, this flag (combined with VelInterpCubic) selects one of four interpolation code paths. When VelInterpCubic = True and OutputAccel = False ("cubic velocity, no acceleration"), Grid3DField_GetCell is called with CalcAccel = OutputAccel = .false., which causes it to skip populating the local AccCell(8,3) array entirely — it is only filled if (CalcAccel) then ....

The velocity is then computed by Grid3DField_GetVelAccCubic, which implements a cubic Hermite spline in time. Critically, the Hermite value formula (not just the optional acceleration output) depends on AccCell, since a proper cubic-in-time interpolant needs derivative (tangent) information at both time endpoints to compute the interpolated value itself:

Velocity = C1*P(:, 1) + C2*PP(:, 1) + C3*P(:, 2) + C4*PP(:, 2)

where PP is derived from AccCell. Since AccCell is a local, stack-allocated array that is never initialized in this code path, PP is computed from garbage memory, corrupting Velocity for every single query — confirmed via debugger to contain Inf in every element for this test case, while the companion VelCell (actual wind data) was fully valid. This explains why the reported failure showed a 100%-NaN velocity array with a 100%-valid position array: the defect is purely a function of which interpolation-flag combination is selected, not of any particular spatial location.

Two further, related defects were found in the same area, both stemming from IfW_Grid3DField_CalcAccel (which precomputes the time-derivative arrays G3D%Acc/G3D%AccTower needed for cubic interpolation) confusing a spatial grid-point count with the temporal step count needed by the cubic-spline-in-time derivative calculation:

1. The main-grid guard checked G3D%NTGrids < 3 (NTGrids = number of tower grid points) instead of G3D%NSteps < 3 (number of time steps):

   if (G3D%NTGrids < 3) then
      G3D%Acc = 0.0_SiKi
      return
   end if

   For the very common case of no tower file (NTGrids = 0), this unconditionally forced G3D%Acc to zero and skipped the real cubic-spline derivative calculation — regardless of how many time steps were actually available (4004 in the reported case, 300 in the ifw_BoxExceed reg test). This doesn't itself produce NaN (zero is well-defined), but it silently defeats cubic-in-time interpolation whenever no tower file is present, which is likely most cases that use VelInterpCubic = True.
2. The tower-grid branch further down had the identical mix-up, checking G3D%NTGrids < 3 to decide whether to compute G3D%AccTower, even though each tower height's time-derivative is computed independently and has no minimum spatial-count requirement (the temporal NSteps < 3 case is already handled by the earlier guard). This incorrectly zeroed valid tower acceleration whenever a wind file had only 1 or 2 tower grid points.

Both original defects are pre-existing in rc-5.0.1 (and earlier) — reproduced and confirmed on an unmodified rc-5.0.1 checkout with identical input files, so this is not a regression from any other branch. It is essentially untested in the existing test suite because every other reg-test input file in the repository sets VelInterpCubic = False; the reported case is the only one using True.

Fix

1. IfW_FlowField_GetVelAcc: pass OutputAccel .or. FF%VelInterpCubic as the CalcAccel argument to Grid3DField_GetCell, so AccCell is always populated whenever cubic-in-time interpolation is active, independent of whether the caller also wants acceleration returned. Also widened the existing AccFieldValid guard to fire whenever OutputAccel .or. FF%VelInterpCubic is true but the acceleration field hasn't actually been computed, so a misconfigured caller gets a clear fatal error instead of dereferencing an unallocated array.
2. IfW_Grid3DField_CalcAccel: changed the main-grid early-return guard from G3D%NTGrids < 3 to G3D%NSteps < 3, so real time derivatives are computed whenever enough time samples exist, regardless of whether a tower file is present. This early-return path now also allocates and zeroes G3D%AccTower (when a tower grid is present) so downstream tower interpolation always has a valid array to reference.
3. IfW_Grid3DField_CalcAccel: removed the erroneous G3D%NTGrids < 3 guard in the tower-grid branch entirely — G3D%NSteps >= 3 is already guaranteed by fix DWM Driver #2's early return, and there's no minimum tower-height-count requirement for the per-height time derivative.

Related issue, if one exists

None filed yet; found during investigation of a reported FAST.Farm FARM_InitialCO fatal error ("rotor plane has left the low-resolution domain") on a 2-turbine wind-tunnel test case.

Impacted areas of the software

• modules/inflowwind/src/IfW_FlowField.f90 (IfW_FlowField_GetVelAcc, IfW_Grid3DField_CalcAccel)
• Any glue code or driver using InflowWind's Grid3D flow field type (TurbSim, Bladed, HAWC full-field wind, WindType = 3, 4, 5, 7) with VelInterpCubic = True: FAST.Farm (AWAE ambient-wind sampling and per-turbine InflowWind instances), OpenFAST, InflowWind driver, AeroDyn driver.
• No effect on the (much more common) VelInterpCubic = False linear-interpolation path.

Additional supporting information

Diagnosed via GDB by tracing a NaN observed in AeroDyn's RotInflow%Blade%InflowVel backwards through FASTWrapper, AWAE's ambient-wind fill, and into InflowWind's Grid3D cubic interpolation, where AccCell was directly observed to contain Inf in all elements at the point of use. Confirmed the same failure reproduces on an unmodified rc-5.0.1 build with AMReX disabled, ruling out any connection to unrelated work-in-progress changes on another branch (including a separately-identified, unrelated WakeDynamics array-bounds issue, which is not part of this PR).

Generative AI usage

Root-cause investigation (GDB tracing, source review), the code fixes, and the new unit tests were developed with substantial assistance from AI coding agents (GitHub Copilot chat agent, running Claude Sonnet 5), including responding to GitHub Copilot's automated PR review comments across two review rounds.

Co-authored-by: GitHub Copilot 175728472+Copilot@users.noreply.github.com
Co-authored-by: Claude Sonnet 5 noreply@anthropic.com
Reviewed by: Google Gemini gemini@google.com

Testing

• Added modules/inflowwind/tests/test_grid3d_field.F90 with three new unit tests, registered in inflowwind_utest:
  • test_grid3d_cubic_vel_only — builds a minimal synthetic Grid3D field with a spatially-uniform, linear-in-time velocity ramp and VelInterpCubic = True; queries velocity both with and without requesting acceleration output. Fails today with "not finite (NaN/Inf)" without fix AeroDyn14 Driver #1; passes with it.
  • test_grid3d_calcaccel_no_tower — calls IfW_Grid3DField_CalcAccel directly on a field with no tower grid and checks that the computed derivative recovers the known true value. Fails today (derivative forced to zero) without fix DWM Driver #2; passes with it.
  • test_grid3d_calcaccel_few_tower_points — calls IfW_Grid3DField_CalcAccel on a field with only 2 tower grid points and checks that the tower acceleration is actually computed rather than zeroed. Fails today without fix DLLEXPORT directives in FAST_Library #3; passes with it.
• Full inflowwind_utest suite passes with no regressions in any pre-existing test.
• Rebuilt FAST.Farm with the fix and reran the original, unmodified reported test case (VelInterpCubic = True): the instant t=0 FARM_InitialCO fatal error no longer occurs; the simulation progresses normally through simulated time.
• ifw_BoxExceed regression baseline updated to reflect corrected (previously erroneously-zeroed) acceleration values; other InflowWind regression tests unaffected.
• Three MHK regression cases also need baseline updates: ad_MHK_RM1_Fixed (RtFldMyg), MHK_RM1_Floating (B2FldMx, RtFldFzg, RtFldMxh), and MHK_RM1_Floating_MR (R2HbMbx, R2HbMby). All three set VelInterpCubic = False but are MHK cases, where OutputAccel is always forced True — so IfW_Grid3DField_CalcAccel has always run for them, independent of fix AeroDyn14 Driver #1. All three also use a TurbSim wind file with no tower extension (WrADTWR = False, i.e. NTGrids = 0), which is exactly the condition fix DWM Driver #2 addresses: previously G3D%Acc was forced to zero for the whole run regardless of the ~700 available time steps, silently feeding zero fluid acceleration into AeroDyn's MHK added-mass/inertial force terms (InflowAcc → aFBTemp/aFTTemp). No other regression cases are affected.
• r-test branch merge required (pointer updated to include corrected ifw_BoxExceed baseline)

[Copilot]

Pull request overview

Fixes InflowWind Grid3D cubic-in-time velocity interpolation behavior so callers that request velocity-only (no acceleration output) no longer receive NaN/Inf wind velocities, and corrects the Grid3D acceleration precompute guard to use the number of time steps rather than the tower-grid count.

Changes:

• Ensure Grid3DField_GetCell populates AccCell whenever VelInterpCubic is enabled, independent of whether AccelUVW is allocated by the caller.
• Fix IfW_Grid3DField_CalcAccel early-return condition to check G3D%NSteps < 3 (time steps) instead of G3D%NTGrids < 3 (tower grid points).
• Add and register new InflowWind unit-test coverage for the cubic-velocity-only bug path.

Reviewed changes

Copilot reviewed 4 out of 4 changed files in this pull request and generated 3 comments.

File	Description
unit_tests/CMakeLists.txt	Adds the new Grid3D unit test source to the InflowWind unit test target.
modules/inflowwind/tests/test_grid3d_field.F90	Introduces a new unit test module exercising the cubic velocity interpolation path.
modules/inflowwind/tests/inflowwind_utest.F90	Registers the new Grid3D test suite in the InflowWind unit test runner.
modules/inflowwind/src/IfW_FlowField.f90	Fixes the Grid3D cubic interpolation path selection and corrects the CalcAccel guard condition.

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[Copilot]

Pull request overview

Copilot reviewed 4 out of 4 changed files in this pull request and generated 2 comments.

[andrew-platt]

ifw_BoxExceed regression baseline comparison

The build-all-test-modules-debug CI failure on ifw_BoxExceed is an expected consequence of this fix, not a regression.

ifw_BoxExceed uses a TurbSim .bts wind file with no tower grid (NTGrids=0) and requests acceleration output (CalcAccel=true). Acceleration is always computed via the cubic-Hermite-in-time formula, which depends on IfW_Grid3DField_CalcAccel. Under the old NTGrids<3 guard, this routine zeroed the time-derivative input for any no-tower case, so the reported acceleration was computed as if dU/dt = 0 at the endpoints. With the fix (NSteps-based guard), the real cubic-spline-derived derivative is used instead, which is why the accelerations reported here differ from the current committed baseline.

Position and velocity channels are byte-for-byte identical (0.0 diff across all 300 points) — only the acceleration channels (UA, VA, WA) change:

Channel	Max |diff|	Mean |diff|	Old max magnitude	New max magnitude
T, X, Y, Z, U, V, W	0.0	0.0	—	—
UA	10.16	1.30	5.37	13.57
VA	11.97	1.39	5.57	15.50
WA	9.26	1.28	7.97	12.63

The new accelerations are ~2-3x larger in typical magnitude, consistent with the old code systematically biasing derivatives toward zero rather than a small numeric drift.

[Copilot]

Pull request overview

Copilot reviewed 5 out of 5 changed files in this pull request and generated no new comments.

Suppressed comments (1)

modules/inflowwind/src/IfW_FlowField.f90:85

• The AccFieldValid guard now also triggers when FF%VelInterpCubic is enabled even if AccelUVW is not allocated, but the fatal message still says "Accel output requested". This is misleading for the cubic-velocity/no-accel caller case (the one this PR fixes). Update the message to mention cubic velocity interpolation (or more generally that an acceleration field is required).

   ! 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)