From 8516a85679b2934ccdc5dc7849bfe76a79288274 Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Thu, 30 Jul 2026 18:19:14 -0600 Subject: [PATCH 01/35] Optimize GridInterp tensor-product cubic interpolation Speed up the 3D/4D grid interpolation routines in the NWTC Library while preserving the public interface and functionality. Results are numerically equivalent, differing only at floating-point truncation level from a strictly serial evaluation. Changes: - Setup routines (3D/4D/N): build the weight tensors with array syntax over the contiguous first dimension so the inner writes vectorize, and hoist shared subexpressions (e.g. Nkl = N1D(k,3)*N1D(l,4)). - Scalar and vector interpolation routines: accumulate over the innermost dimension into four independent partial sums to break the serial add dependency chain (latency hiding) and expose ILP, then combine at the end. Index lookups are hoisted out of the inner loop. - Vector routines: move the component (vi) loop outermost for unit-stride access and per-component accumulation. - Slope routines (N/S): accumulate both derivatives with dual accumulators and reuse a single data load per grid point instead of loading the same value twice. Validation: all 10 SeaState regression tests pass, and the HydroDyn WAMIT interpolation paths (including hd_NBodyMod1) pass. gfortran does not reassociate reduction sums by default, so the multi-accumulator regrouping is the primary win. Co-authored-by: GitHub Copilot <175728472+Copilot@users.noreply.github.com> Co-authored-by: Claude Opus --- modules/nwtc-library/src/GridInterp.f90 | 314 +++++++++++++++--------- 1 file changed, 202 insertions(+), 112 deletions(-) diff --git a/modules/nwtc-library/src/GridInterp.f90 b/modules/nwtc-library/src/GridInterp.f90 index 769b905fca..aa73489137 100644 --- a/modules/nwtc-library/src/GridInterp.f90 +++ b/modules/nwtc-library/src/GridInterp.f90 @@ -320,9 +320,7 @@ Subroutine GridInterpSetup3D( position, p, m, ErrStat, ErrMsg ) do k = 1,4 do j = 1,4 - do i = 1,4 - m%N3D(i,j,k) = N1D(i,1)*N1D(j,2)*N1D(k,3) - end do + m%N3D(:,j,k) = N1D(:,1) * (N1D(j,2)*N1D(k,3)) end do end do @@ -346,6 +344,7 @@ Subroutine GridInterpSetup4D( position, p, m, ErrStat, ErrMsg ) integer(IntKi) :: dim,i,j,k,l integer(IntKi) :: support real(ReKi) :: N1D(4,4) + real(ReKi) :: Nkl real(ReKi) :: isopc ! isoparametric coordinates integer(IntKi) :: ErrStat2 character(ErrMsgLen) :: ErrMsg2 @@ -361,10 +360,9 @@ Subroutine GridInterpSetup4D( position, p, m, ErrStat, ErrMsg ) do l = 1,4 do k = 1,4 + Nkl = N1D(k,3)*N1D(l,4) do j = 1,4 - do i = 1,4 - m%N4D(i,j,k,l) = N1D(i,1)*N1D(j,2)*N1D(k,3)*N1D(l,4) - end do + m%N4D(:,j,k,l) = N1D(:,1) * (N1D(j,2)*Nkl) end do end do end do @@ -409,10 +407,8 @@ Subroutine GridInterpSetupN( position, p, m, ErrStat, ErrMsg ) ! Need two sets of weights for d(.)/dx and d(.)/dy. Borrow m%N4D for this. do k = 1,4 do j = 1,4 - do i = 1,4 - m%N4D(i,j,k,1) = N1D(i,1)*N1Ddx(j,2)*N1D (k,3) - m%N4D(i,j,k,2) = N1D(i,1)*N1D (j,2)*N1Ddx(k,3) - end do + m%N4D(:,j,k,1) = N1D(:,1) * (N1Ddx(j,2)*N1D (k,3)) + m%N4D(:,j,k,2) = N1D(:,1) * (N1D (j,2)*N1Ddx(k,3)) end do end do @@ -436,17 +432,21 @@ function GridInterp3DR4( data, m ) character(*), parameter :: RoutineName = 'GridInterp3DR4' real(SiKi) :: GridInterp3DR4 - integer(IntKi) :: i,j,k + real(SiKi) :: acc(4) + integer(IntKi) :: i,j,k,jj,kk ! interpolate - GridInterp3DR4 = 0.0_SiKi + acc = 0.0_SiKi do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - GridInterp3DR4 = GridInterp3DR4 + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + acc(i) = acc(i) + m%N3D(i,j,k) * data( m%Indx(i,1), jj, kk ) end do end do end do + GridInterp3DR4 = acc(1) + acc(2) + acc(3) + acc(4) end function GridInterp3DR4 @@ -456,17 +456,21 @@ function GridInterp3DR8( data, m ) character(*), parameter :: RoutineName = 'GridInterp3DR8' real(DbKi) :: GridInterp3DR8 - integer(IntKi) :: i,j,k + real(DbKi) :: acc(4) + integer(IntKi) :: i,j,k,jj,kk ! interpolate - GridInterp3DR8 = 0.0_DbKi + acc = 0.0_DbKi do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - GridInterp3DR8 = GridInterp3DR8 + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + acc(i) = acc(i) + m%N3D(i,j,k) * data( m%Indx(i,1), jj, kk ) end do end do end do + GridInterp3DR8 = acc(1) + acc(2) + acc(3) + acc(4) end function GridInterp3DR8 @@ -482,18 +486,23 @@ function GridInterp3DVecR4( data, m ) character(*), parameter :: RoutineName = 'GridInterp3DVecR4' integer(IntKi), parameter :: vDim = 3 integer(IntKi) :: i,j,k,vi + integer(IntKi) :: jj,kk + real(SiKi) :: acc(4) real(SiKi) :: GridInterp3DVecR4(vDim) ! interpolate - GridInterp3DVecR4 = 0.0_SiKi - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp3DVecR4(vi) = GridInterp3DVecR4(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + do vi = 1,vDim + acc = 0.0_SiKi + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N3D(i,j,k) * data( m%Indx(i,1), jj, kk, vi ) end do end do end do + GridInterp3DVecR4(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp3DVecR4 @@ -505,18 +514,23 @@ function GridInterp3DVecR8( data, m ) character(*), parameter :: RoutineName = 'GridInterp3DVecR8' integer(IntKi), parameter :: vDim = 3 integer(IntKi) :: i,j,k,vi + integer(IntKi) :: jj,kk + real(DbKi) :: acc(4) real(DbKi) :: GridInterp3DVecR8(vDim) ! interpolate - GridInterp3DVecR8 = 0.0_DbKi - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp3DVecR8(vi) = GridInterp3DVecR8(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + do vi = 1,vDim + acc = 0.0_DbKi + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N3D(i,j,k) * data( m%Indx(i,1), jj, kk, vi ) end do end do end do + GridInterp3DVecR8(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp3DVecR8 @@ -533,18 +547,23 @@ function GridInterp3DVec6R4( data, m ) character(*), parameter :: RoutineName = 'GridInterp3DVec6R4' integer(IntKi), parameter :: vDim = 6 integer(IntKi) :: i,j,k,vi + integer(IntKi) :: jj,kk + real(SiKi) :: acc(4) real(SiKi) :: GridInterp3DVec6R4(vDim) ! interpolate - GridInterp3DVec6R4 = 0.0_SiKi - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp3DVec6R4(vi) = GridInterp3DVec6R4(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + do vi = 1,vDim + acc = 0.0_SiKi + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N3D(i,j,k) * data( m%Indx(i,1), jj, kk, vi ) end do end do end do + GridInterp3DVec6R4(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp3DVec6R4 @@ -556,18 +575,23 @@ function GridInterp3DVec6R8( data, m ) character(*), parameter :: RoutineName = 'GridInterp3DVec6R8' integer(IntKi), parameter :: vDim = 6 integer(IntKi) :: i,j,k,vi + integer(IntKi) :: jj,kk + real(DbKi) :: acc(4) real(DbKi) :: GridInterp3DVec6R8(vDim) ! interpolate - GridInterp3DVec6R8 = 0.0_DbKi - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp3DVec6R8(vi) = GridInterp3DVec6R8(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + do vi = 1,vDim + acc = 0.0_DbKi + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N3D(i,j,k) * data( m%Indx(i,1), jj, kk, vi ) end do end do end do + GridInterp3DVec6R8(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp3DVec6R8 @@ -583,19 +607,24 @@ function GridInterp4DR4( data, m ) character(*), parameter :: RoutineName = 'GridInterp4DR4' real(SiKi) :: GridInterp4DR4 - integer(IntKi) :: i,j,k,l + real(SiKi) :: acc(4) + integer(IntKi) :: i,j,k,l,jj,kk,ll ! interpolate - GridInterp4DR4 = 0.0_SiKi + acc = 0.0_SiKi do l = 1,4 + ll = m%Indx(l,4) do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - GridInterp4DR4 = GridInterp4DR4 + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4) ) + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll ) end do end do end do end do + GridInterp4DR4 = acc(1) + acc(2) + acc(3) + acc(4) end function GridInterp4DR4 @@ -605,19 +634,24 @@ function GridInterp4DR8( data, m ) character(*), parameter :: RoutineName = 'GridInterp4DR8' real(DbKi) :: GridInterp4DR8 - integer(IntKi) :: i,j,k,l + real(DbKi) :: acc(4) + integer(IntKi) :: i,j,k,l,jj,kk,ll ! interpolate - GridInterp4DR8 = 0.0_DbKi + acc = 0.0_DbKi do l = 1,4 + ll = m%Indx(l,4) do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - GridInterp4DR8 = GridInterp4DR8 + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4) ) + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll ) end do end do end do end do + GridInterp4DR8 = acc(1) + acc(2) + acc(3) + acc(4) end function GridInterp4DR8 @@ -633,20 +667,26 @@ function GridInterp4DVecR4( data, m ) character(*), parameter :: RoutineName = 'GridInterp4DVecR4' integer(IntKi), parameter :: vDim = 3 integer(IntKi) :: i,j,k,l,vi + integer(IntKi) :: jj,kk,ll + real(SiKi) :: acc(4) real(SiKi) :: GridInterp4DVecR4(vDim) ! interpolate - GridInterp4DVecR4 = 0.0_SiKi - do l = 1,4 - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp4DVecR4(vi) = GridInterp4DVecR4(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + do vi = 1,vDim + acc = 0.0_SiKi + do l = 1,4 + ll = m%Indx(l,4) + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll, vi ) end do end do end do end do + GridInterp4DVecR4(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp4DVecR4 @@ -658,20 +698,26 @@ function GridInterp4DVecR8( data, m ) character(*), parameter :: RoutineName = 'GridInterp4DVecR8' integer(IntKi), parameter :: vDim = 3 integer(IntKi) :: i,j,k,l,vi + integer(IntKi) :: jj,kk,ll + real(DbKi) :: acc(4) real(DbKi) :: GridInterp4DVecR8(vDim) ! interpolate - GridInterp4DVecR8 = 0.0_DbKi - do l = 1,4 - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp4DVecR8(vi) = GridInterp4DVecR8(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + do vi = 1,vDim + acc = 0.0_DbKi + do l = 1,4 + ll = m%Indx(l,4) + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll, vi ) end do end do end do end do + GridInterp4DVecR8(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp4DVecR8 @@ -688,20 +734,26 @@ function GridInterp4DVec6R4( data, m ) character(*), parameter :: RoutineName = 'GridInterp4DVec6R4' integer(IntKi), parameter :: vDim = 6 integer(IntKi) :: i,j,k,l,vi + integer(IntKi) :: jj,kk,ll + real(SiKi) :: acc(4) real(SiKi) :: GridInterp4DVec6R4(vDim) ! interpolate - GridInterp4DVec6R4 = 0.0_SiKi - do l = 1,4 - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp4DVec6R4(vi) = GridInterp4DVec6R4(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + do vi = 1,vDim + acc = 0.0_SiKi + do l = 1,4 + ll = m%Indx(l,4) + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll, vi ) end do end do end do end do + GridInterp4DVec6R4(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp4DVec6R4 @@ -713,20 +765,26 @@ function GridInterp4DVec6R8( data, m ) character(*), parameter :: RoutineName = 'GridInterp4DVec6R8' integer(IntKi), parameter :: vDim = 6 integer(IntKi) :: i,j,k,l,vi + integer(IntKi) :: jj,kk,ll + real(DbKi) :: acc(4) real(DbKi) :: GridInterp4DVec6R8(vDim) ! interpolate - GridInterp4DVec6R8 = 0.0_DbKi - do l = 1,4 - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp4DVec6R8(vi) = GridInterp4DVec6R8(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + do vi = 1,vDim + acc = 0.0_DbKi + do l = 1,4 + ll = m%Indx(l,4) + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll, vi ) end do end do end do end do + GridInterp4DVec6R8(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp4DVec6R8 @@ -743,20 +801,26 @@ function GridInterp4DVecNR4( vDim, data, m ) character(*), parameter :: RoutineName = 'GridInterp4DVecNR4' integer(IntKi) :: i,j,k,l,vi + integer(IntKi) :: jj,kk,ll + real(SiKi) :: acc(4) real(SiKi) :: GridInterp4DVecNR4(vDim) ! interpolate - GridInterp4DVecNR4 = 0.0_SiKi - do l = 1,4 - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp4DVecNR4(vi) = GridInterp4DVecNR4(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + do vi = 1,vDim + acc = 0.0_SiKi + do l = 1,4 + ll = m%Indx(l,4) + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll, vi ) end do end do end do end do + GridInterp4DVecNR4(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp4DVecNR4 @@ -768,20 +832,26 @@ function GridInterp4DVecNR8( vDim, data, m ) character(*), parameter :: RoutineName = 'GridInterp4DVecNR8' integer(IntKi) :: i,j,k,l,vi + integer(IntKi) :: jj,kk,ll + real(DbKi) :: acc(4) real(DbKi) :: GridInterp4DVecNR8(vDim) ! interpolate - GridInterp4DVecNR8 = 0.0_DbKi - do l = 1,4 - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp4DVecNR8(vi) = GridInterp4DVecNR8(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + do vi = 1,vDim + acc = 0.0_DbKi + do l = 1,4 + ll = m%Indx(l,4) + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll, vi ) end do end do end do end do + GridInterp4DVecNR8(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp4DVecNR8 @@ -799,21 +869,25 @@ function GridInterpNR4( data, p, m ) character(*), parameter :: RoutineName = 'GridInterpNR4' real(SiKi) :: GridInterpNR4(3) real(SiKi) :: dZetadx, dZetady - integer(IntKi) :: i,j,k + real(SiKi) :: accx(4), accy(4), d + integer(IntKi) :: i,j,k,jj,kk ! interpolate slope - dZetadx = 0.0_SiKi - dZetady = 0.0_SiKi + accx = 0.0_SiKi + accy = 0.0_SiKi do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - dZetadx = dZetadx + m%N4D(i,j,k,1) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - dZetady = dZetady + m%N4D(i,j,k,2) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + d = data( m%Indx(i,1), jj, kk ) + accx(i) = accx(i) + m%N4D(i,j,k,1) * d + accy(i) = accy(i) + m%N4D(i,j,k,2) * d end do end do end do - dZetadx = dZetadx / p%delta(2) - dZetady = dZetady / p%delta(3) + dZetadx = ( accx(1) + accx(2) + accx(3) + accx(4) ) / p%delta(2) + dZetady = ( accy(1) + accy(2) + accy(3) + accy(4) ) / p%delta(3) GridInterpNR4 = [-dZetadx,-dZetady,1.0_SiKi] GridInterpNR4 = GridInterpNR4 / TwoNorm(GridInterpNR4) @@ -828,21 +902,25 @@ function GridInterpNR8( data, p, m ) character(*), parameter :: RoutineName = 'GridInterpNR8' real(DbKi) :: GridInterpNR8(3) real(DbKi) :: dZetadx, dZetady - integer(IntKi) :: i,j,k + real(DbKi) :: accx(4), accy(4), d + integer(IntKi) :: i,j,k,jj,kk ! interpolate slope - dZetadx = 0.0_DbKi - dZetady = 0.0_DbKi + accx = 0.0_DbKi + accy = 0.0_DbKi do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - dZetadx = dZetadx + m%N4D(i,j,k,1) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - dZetady = dZetady + m%N4D(i,j,k,2) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + d = data( m%Indx(i,1), jj, kk ) + accx(i) = accx(i) + m%N4D(i,j,k,1) * d + accy(i) = accy(i) + m%N4D(i,j,k,2) * d end do end do end do - dZetadx = dZetadx / p%delta(2) - dZetady = dZetady / p%delta(3) + dZetadx = ( accx(1) + accx(2) + accx(3) + accx(4) ) / p%delta(2) + dZetady = ( accy(1) + accy(2) + accy(3) + accy(4) ) / p%delta(3) GridInterpNR8 = (/-dZetadx,-dZetady,1.0_DbKi/) GridInterpNR8 = GridInterpNR8 / TwoNorm(GridInterpNR8) @@ -861,19 +939,25 @@ function GridInterpSR4( data, p, m ) character(*), parameter :: RoutineName = 'GridInterpSR4' real(SiKi) :: GridInterpSR4(2) - integer(IntKi) :: i,j,k,dir + real(SiKi) :: accx(4), accy(4), d + integer(IntKi) :: i,j,k,jj,kk ! interpolate slope - GridInterpSR4 = 0.0_SiKi + accx = 0.0_SiKi + accy = 0.0_SiKi do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - do dir = 1,2 - GridInterpSR4(dir) = GridInterpSR4(dir) + m%N4D(i,j,k,dir) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - end do + d = data( m%Indx(i,1), jj, kk ) + accx(i) = accx(i) + m%N4D(i,j,k,1) * d + accy(i) = accy(i) + m%N4D(i,j,k,2) * d end do end do end do + GridInterpSR4(1) = accx(1) + accx(2) + accx(3) + accx(4) + GridInterpSR4(2) = accy(1) + accy(2) + accy(3) + accy(4) GridInterpSR4 = GridInterpSR4 / p%delta(2:3) end function GridInterpSR4 @@ -885,19 +969,25 @@ function GridInterpSR8( data, p, m ) character(*), parameter :: RoutineName = 'GridInterpSR8' real(DbKi) :: GridInterpSR8(2) - integer(IntKi) :: i,j,k,dir + real(DbKi) :: accx(4), accy(4), d + integer(IntKi) :: i,j,k,jj,kk ! interpolate slope - GridInterpSR8 = 0.0_DbKi + accx = 0.0_DbKi + accy = 0.0_DbKi do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - do dir = 1,2 - GridInterpSR8(dir) = GridInterpSR8(dir) + m%N4D(i,j,k,dir) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - end do + d = data( m%Indx(i,1), jj, kk ) + accx(i) = accx(i) + m%N4D(i,j,k,1) * d + accy(i) = accy(i) + m%N4D(i,j,k,2) * d end do end do end do + GridInterpSR8(1) = accx(1) + accx(2) + accx(3) + accx(4) + GridInterpSR8(2) = accy(1) + accy(2) + accy(3) + accy(4) GridInterpSR8 = GridInterpSR8 / p%delta(2:3) end function GridInterpSR8 From 02d93217803adc82625f681c6929b6452fb004c5 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 08:12:48 -0600 Subject: [PATCH 02/35] AeroDyn/OLAF Registry: add xPts/yPts/zPts and list-file fields to GridOutType --- modules/aerodyn/src/FVW_Registry.txt | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/modules/aerodyn/src/FVW_Registry.txt b/modules/aerodyn/src/FVW_Registry.txt index 45a7aff474..2d9ce90812 100644 --- a/modules/aerodyn/src/FVW_Registry.txt +++ b/modules/aerodyn/src/FVW_Registry.txt @@ -28,6 +28,12 @@ typedef ^ ^ IntKi typedef ^ ^ ReKi uGrid {:}{:}{:}{:} - - "Grid velocity 3 x nz x ny x nx" - typedef ^ ^ ReKi omGrid {:}{:}{:}{:} - - "Grid vorticity 3 x nz x ny x nx" - typedef ^ ^ DbKi tLastOutput - - - "Last output time" - +typedef ^ ^ ReKi xPts {:} - - "Explicit x coordinates (non-equidistant grid, if used)" m +typedef ^ ^ ReKi yPts {:} - - "Explicit y coordinates (non-equidistant grid, if used)" m +typedef ^ ^ ReKi zPts {:} - - "Explicit z coordinates (non-equidistant grid, if used)" m +typedef ^ ^ CHARACTER(1024) xListFile - - - "File with explicit x coordinates (empty if equidistant)" - +typedef ^ ^ CHARACTER(1024) yListFile - - - "File with explicit y coordinates (empty if equidistant)" - +typedef ^ ^ CHARACTER(1024) zListFile - - - "File with explicit z coordinates (empty if equidistant)" - ##################### Segments ############### typedef FVW/FVW T_Sgmt ReKi Points :: - - "Points delimiting the segments" - From 1463748db34cc0c1b737fab0604643847ce48815 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 08:30:55 -0600 Subject: [PATCH 03/35] AeroDyn/OLAF: parse a number or quoted filename in ReadGridOut --- modules/aerodyn/src/FVW_IO.f90 | 43 ++++++++++++++++++++++++++-------- 1 file changed, 33 insertions(+), 10 deletions(-) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 4cde5064d4..014d36f607 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -327,19 +327,42 @@ subroutine ReadGridOut(sLine, GridOut) return endif endif - ! x,y,z + ! x ErrMsg2='Error reading OLAF "x" inputs for grid outputs line: '//trim(sLine) - if (.not. is_numeric(StrArray( 6), GridOut%xStart) ) return - if (.not. is_numeric(StrArray( 7), GridOut%xEnd ) ) return - if (.not. is_integer(StrArray( 8), GridOut%nx ) ) return + GridOut%xListFile = '' + if ( is_numeric(StrArray(6), GridOut%xStart) ) then + if (.not. is_numeric(StrArray(7), GridOut%xEnd) ) return + if (.not. is_integer(StrArray(8), GridOut%nx ) ) return + else + GridOut%xListFile = StrArray(6) + GridOut%xStart = 0.0_ReKi + GridOut%xEnd = 0.0_ReKi + GridOut%nx = -1 + endif + ! y ErrMsg2='Error reading OLAF "y" inputs for grid outputs line: '//trim(sLine) - if (.not. is_numeric(StrArray( 9), GridOut%yStart) ) return - if (.not. is_numeric(StrArray(10), GridOut%yEnd ) ) return - if (.not. is_integer(StrArray(11), GridOut%ny ) ) return + GridOut%yListFile = '' + if ( is_numeric(StrArray(9), GridOut%yStart) ) then + if (.not. is_numeric(StrArray(10), GridOut%yEnd) ) return + if (.not. is_integer(StrArray(11), GridOut%ny ) ) return + else + GridOut%yListFile = StrArray(9) + GridOut%yStart = 0.0_ReKi + GridOut%yEnd = 0.0_ReKi + GridOut%ny = -1 + endif + ! z ErrMsg2='Error reading OLAF "z" inputs for grid outputs line: '//trim(sLine) - if (.not. is_numeric(StrArray(12), GridOut%zStart) ) return - if (.not. is_numeric(StrArray(13), GridOut%zEnd ) ) return - if (.not. is_integer(StrArray(14), GridOut%nz ) ) return + GridOut%zListFile = '' + if ( is_numeric(StrArray(12), GridOut%zStart) ) then + if (.not. is_numeric(StrArray(13), GridOut%zEnd) ) return + if (.not. is_integer(StrArray(14), GridOut%nz ) ) return + else + GridOut%zListFile = StrArray(12) + GridOut%zStart = 0.0_ReKi + GridOut%zEnd = 0.0_ReKi + GridOut%nz = -1 + endif ! Success ErrStat2=ErrID_None ErrMsg2='' From 1ae78909c2d969110bfcca9c2ac94aed21185cd4 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 08:55:38 -0600 Subject: [PATCH 04/35] AeroDyn/OLAF: call ResolveGridAxis in FVW_ReadInputFile and require GridType=1 for list-based grids --- modules/aerodyn/src/FVW_IO.f90 | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 014d36f607..5a16dea07e 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -100,13 +100,28 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) allocate(m%GridOutputs(p%nGridOut), stat=ErrStat2); CALL ReadCom (UnIn,FileName, 'GridOutHeaders', ErrStat2,ErrMsg2); if(Failed()) return CALL ReadCom (UnIn,FileName, 'GridOutUnits', ErrStat2,ErrMsg2); if(Failed()) return - do i =1, p%nGridOut + do i =1, p%nGridOut ErrMsg2='Error reading OLAF grid outputs line '//trim(num2lstr(i)) read(UnIn, fmt='(A)', iostat=ErrStat2) sLine ; if(Failed()) return call ReadGridOut(sLine, m%GridOutputs(i)); if(Failed()) return + ! Resolve each axis to an explicit coordinate array (regardless of whether it is a list file or a range) + call ResolveGridAxis(m%GridOutputs(i)%xStart, m%GridOutputs(i)%xEnd, m%GridOutputs(i)%nx, & + m%GridOutputs(i)%xListFile, PriPath, m%GridOutputs(i)%xPts, ErrStat2, ErrMsg2); if(Failed()) return + call ResolveGridAxis(m%GridOutputs(i)%yStart, m%GridOutputs(i)%yEnd, m%GridOutputs(i)%ny, & + m%GridOutputs(i)%yListFile, PriPath, m%GridOutputs(i)%yPts, ErrStat2, ErrMsg2); if(Failed()) return + call ResolveGridAxis(m%GridOutputs(i)%zStart, m%GridOutputs(i)%zEnd, m%GridOutputs(i)%nz, & + m%GridOutputs(i)%zListFile, PriPath, m%GridOutputs(i)%zPts, ErrStat2, ErrMsg2); if(Failed()) return + ! Error checking if (Check(m%GridOutputs(i)%nx<1, 'Grid output nx needs to be >=1')) return if (Check(m%GridOutputs(i)%ny<1, 'Grid output ny needs to be >=1')) return if (Check(m%GridOutputs(i)%nz<1, 'Grid output nz needs to be >=1')) return + ! Vorticity requires equidistant spacing. GridType must be 1 when using list files. + if (m%GridOutputs(i)%type==idGridVelVorticity) then + if (Check( len_trim(m%GridOutputs(i)%xListFile)>0 .or. & + len_trim(m%GridOutputs(i)%yListFile)>0 .or. & + len_trim(m%GridOutputs(i)%zListFile)>0, & + 'Grid "'//trim(m%GridOutputs(i)%name)//'": vorticity requires equidistant spacing. nGridOut must be 1 when using list files.')) return + endif enddo endif From aa525a00b2e84894a134691fb613076bb9a78aa3 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 09:19:43 -0600 Subject: [PATCH 05/35] AeroDyn/OLAF: New ResolveGridAxis for list-based or range --- modules/aerodyn/src/FVW_IO.f90 | 67 +++++++++++++++++++++++++++++++++- 1 file changed, 66 insertions(+), 1 deletion(-) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 5a16dea07e..756db5686d 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -120,7 +120,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) if (Check( len_trim(m%GridOutputs(i)%xListFile)>0 .or. & len_trim(m%GridOutputs(i)%yListFile)>0 .or. & len_trim(m%GridOutputs(i)%zListFile)>0, & - 'Grid "'//trim(m%GridOutputs(i)%name)//'": vorticity requires equidistant spacing. nGridOut must be 1 when using list files.')) return + 'Grid "'//trim(m%GridOutputs(i)%name)//'": vorticity requires equidistant spacing. GridType must be 1 when using list files.')) return endif enddo endif @@ -383,6 +383,71 @@ subroutine ReadGridOut(sLine, GridOut) ErrMsg2='' end subroutine ReadGridOut + ! Resolve one grid axis to an explicit, strictly increasing coordinate array, + ! either read from ListFile (if given) or built from an equidistant Start/End/n range. + subroutine ResolveGridAxis(AStart, AEnd, n, ListFile, Pts) + real(ReKi), intent(in) :: AStart, AEnd + integer(IntKi), intent(inout) :: n + character(*), intent(in) :: ListFile + real(ReKi), allocatable, intent(out) :: Pts(:) + character(1024) :: FullFile + integer(IntKi) :: UnList, j, IOS + real(ReKi) :: val + + if (len_trim(ListFile) == 0) then + ! Equidistant points, expressed as an explicit array + allocate(Pts(max(n,1))) + do j = 1, n + Pts(j) = AStart + (AEnd - AStart) * real(j-1,ReKi) / real(max(n-1,1),ReKi) + enddo + ErrStat2 = ErrID_None; ErrMsg2 = '' + return + endif + + ! Explicit list from file + FullFile = ListFile + if (PathIsRelative(FullFile)) FullFile = trim(PriPath)//trim(FullFile) + call GetNewUnit(UnList) + call OpenFInpFile(UnList, trim(FullFile), ErrStat2, ErrMsg2) + if (ErrStat2 /= ErrID_None) return + + ! Count valid lines + n = 0 + do + read(UnList, *, iostat=IOS) val + if (IOS /= 0) exit + n = n + 1 + enddo + if (n < 1) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = 'Grid point list file "'//trim(FullFile)//'" contains no valid values.' + close(UnList) + return + endif + + ! Read values + rewind(UnList) + allocate(Pts(n)) + do j = 1, n + read(UnList, *) Pts(j) + enddo + close(UnList) + + ! Validate strictly increasing, no duplicates + do j = 2, n + if (Pts(j) <= Pts(j-1)) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = 'Grid point list file "'//trim(FullFile)//'" must be strictly increasing (no duplicates); '// & + 'value at line '//trim(Num2LStr(j))//' ('//trim(Num2LStr(Pts(j)))//') is not greater than '// & + 'the previous value ('//trim(Num2LStr(Pts(j-1)))//').' + return + endif + enddo + + ErrStat2 = ErrID_None + ErrMsg2 = '' + end subroutine ResolveGridAxis + END SUBROUTINE FVW_ReadInputFile From 9b7ff6d9b5b05c501abf2c4222141f1855270cf4 Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Fri, 31 Jul 2026 09:30:50 -0600 Subject: [PATCH 06/35] GridInterp.f90: removed now unused iteration counter i in some subroutines --- modules/nwtc-library/src/GridInterp.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/modules/nwtc-library/src/GridInterp.f90 b/modules/nwtc-library/src/GridInterp.f90 index aa73489137..22ba5a9de1 100644 --- a/modules/nwtc-library/src/GridInterp.f90 +++ b/modules/nwtc-library/src/GridInterp.f90 @@ -302,7 +302,7 @@ Subroutine GridInterpSetup3D( position, p, m, ErrStat, ErrMsg ) character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None character(*), parameter :: RoutineName = 'GridInterpSetup3D' - integer(IntKi) :: dim,i,j,k + integer(IntKi) :: dim,j,k integer(IntKi) :: support real(ReKi) :: N1D(4,3) real(ReKi) :: isopc ! isoparametric coordinates @@ -341,7 +341,7 @@ Subroutine GridInterpSetup4D( position, p, m, ErrStat, ErrMsg ) character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None character(*), parameter :: RoutineName = 'GridInterpSetup4D' - integer(IntKi) :: dim,i,j,k,l + integer(IntKi) :: dim,j,k,l integer(IntKi) :: support real(ReKi) :: N1D(4,4) real(ReKi) :: Nkl @@ -384,7 +384,7 @@ Subroutine GridInterpSetupN( position, p, m, ErrStat, ErrMsg ) character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None character(*), parameter :: RoutineName = 'GridInterpSetupN' - integer(IntKi) :: dim,i,j,k + integer(IntKi) :: dim,j,k integer(IntKi) :: support real(ReKi) :: N1D(4,3) real(ReKi) :: N1Ddx(4,2:3) From 11a716246fc28e97a1f3b309c261fb9d669a006f Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 09:34:58 -0600 Subject: [PATCH 07/35] AeroDyn/OLAF: export grid velocity as RectilinearGrid (list-based) or StructuredPoints (range) --- modules/aerodyn/src/FVW_IO.f90 | 33 ++++++++++++++++++++++----------- 1 file changed, 22 insertions(+), 11 deletions(-) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 756db5686d..b05f5007a0 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -600,6 +600,7 @@ subroutine WrVTK_FVW_Grid(p, m, iGrid, FileRootName, VTKcount, Twidth, HubOrient character(255) :: Label character(Twidth) :: Tstr ! string for current VTK write-out step (padded with zeros) real(ReKi), dimension(3) :: dx + logical :: bListBased type(GridOutType), pointer :: g type(VTK_Misc) :: mvtk @@ -613,19 +614,29 @@ subroutine WrVTK_FVW_Grid(p, m, iGrid, FileRootName, VTKcount, Twidth, HubOrient g => m%GridOutputs(iGrid) Label=trim(g%name) Filename = TRIM(FileRootName)//'.'//trim(Label)//'.'//Tstr//'.vtk' + bListBased = len_trim(g%xListFile)>0 .or. len_trim(g%yListFile)>0 .or. len_trim(g%zListFile)>0 if ( vtk_new_ascii_file(trim(filename),Label,mvtk) ) then - dx(1) = (g%xEnd- g%xStart)/max(g%nx-1,1) - dx(2) = (g%yEnd- g%yStart)/max(g%ny-1,1) - dx(3) = (g%zEnd- g%zStart)/max(g%nz-1,1) - call vtk_dataset_structured_points((/g%xStart, g%yStart, g%zStart/),dx,(/g%nx,g%ny,g%nz/),mvtk) - call vtk_point_data_init(mvtk) - call vtk_point_data_vector(g%uGrid(1:3,:,:,:),'Velocity',mvtk) - ! Compute vorticity on the fly - if (g%type==idGridVelVorticity) then - call curl_regular_grid(g%uGrid, g%omgrid, 1,1,1, g%nx,g%ny,g%nz, dx(1),dx(2),dx(3)) - call vtk_point_data_vector(g%omGrid(1:3,:,:,:),'Vorticity',mvtk) + if (bListBased) then + ! List-based grid + ! At least one axis is non-equidistant: RectilinearGrid (explicit coordinates) + ! NOTE: vorticity is blocked for this case at read time. GridType = 1 (velocity only). + call vtk_dataset_rectilinear(g%xPts, g%yPts, g%zPts, mvtk) + call vtk_point_data_init(mvtk) + call vtk_point_data_vector(g%uGrid(1:3,:,:,:),'Velocity',mvtk) + else + ! Equidistant grid: StructuredPoints (implicit coordinates) + dx(1) = (g%xEnd- g%xStart)/max(g%nx-1,1) + dx(2) = (g%yEnd- g%yStart)/max(g%ny-1,1) + dx(3) = (g%zEnd- g%zStart)/max(g%nz-1,1) + call vtk_dataset_structured_points((/g%xStart, g%yStart, g%zStart/),dx,(/g%nx,g%ny,g%nz/),mvtk) + call vtk_point_data_init(mvtk) + call vtk_point_data_vector(g%uGrid(1:3,:,:,:),'Velocity',mvtk) + ! Compute vorticity on the fly + if (g%type==idGridVelVorticity) then + call curl_regular_grid(g%uGrid, g%omgrid, 1,1,1, g%nx,g%ny,g%nz, dx(1),dx(2),dx(3)) + call vtk_point_data_vector(g%omGrid(1:3,:,:,:),'Vorticity',mvtk) + endif endif - ! call vtk_close_file(mvtk) endif From 231c6eab9b94f7bbeff53702b638f931fa024291 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 09:38:54 -0600 Subject: [PATCH 08/35] FVW_IO.f90: trim trailing whitespace --- modules/aerodyn/src/FVW_IO.f90 | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index b05f5007a0..1ea8dad016 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -18,7 +18,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None ! Local variables character(1024) :: PriPath ! the path to the primary input file - character(1024) :: sDummy, sLine ! string to temporarially hold value of read line + character(1024) :: sDummy, sLine ! string to temporarially hold value of read line integer(IntKi) :: UnIn, i integer(IntKi) :: ErrStat2 character(ErrMsgLen) :: ErrMsg2 @@ -26,7 +26,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) ErrMsg = "" Inp%SrcPnlFile = '' ! TODO registry init for empty strings ! Open file - CALL GetNewUnit( UnIn ) + CALL GetNewUnit( UnIn ) CALL OpenFInpfile(UnIn, TRIM(FileName), ErrStat2, ErrMsg2) if (Check( ErrStat2 /= ErrID_None , 'Could not open input file')) return CALL GetPath( FileName, PriPath ) ! Input files will be relative to the path where the primary input file is located. @@ -128,7 +128,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) ! --- Advanced Options ! NOTE: no error handling since this is for debug ! Default options are typically "true" - CALL ReadCom(UnIn,FileName, '=== Separator' ,ErrStat2,ErrMsg2); + CALL ReadCom(UnIn,FileName, '=== Separator' ,ErrStat2,ErrMsg2); CALL ReadCom(UnIn,FileName, '--- Advanced options header' ,ErrStat2,ErrMsg2); if(ErrStat2==ErrID_None) then call WrScr(' - Reading advanced options for OLAF:') @@ -218,7 +218,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) if (Check(Inp%WingRegParam<0 , 'Wing regularization parameter (WakeRegParam) should be positive')) return if (Check(Inp%CoreSpreadEddyVisc<0 , 'Core spreading eddy viscosity (CoreSpreadEddyVisc) should be positive')) return - ! Removing the shed vorticity is a dangerous option if this is done too close to the blades. + ! Removing the shed vorticity is a dangerous option if this is done too close to the blades. ! To be safe, we will no matter what ensure that the last segments of NW are 0 if FWShedVorticity is False (see PackPanelsToSegments) ! Still we force the user to be responsible. if (Check((.not.(Inp%FWShedVorticity)) .and. Inp%nNWPanels<30, '`FWShedVorticity` should be true if `nNWPanels`<30. Alternatively, use a larger number of NWPanels ')) return @@ -231,7 +231,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) CONTAINS logical function Failed() - call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'FVW_ReadInputFile') + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'FVW_ReadInputFile') Failed = ErrStat >= AbortErrLev if (Failed) call CleanUp() end function Failed @@ -304,7 +304,7 @@ subroutine ReadGridOut(sLine, GridOut) ErrStat2=ErrID_Fatal ErrMsg2='Error reading OLAF grid outputs line: '//trim(sLine) ! Name - GridOut%name =StrArray(1) + GridOut%name =StrArray(1) ! Type if (.not. is_integer(StrArray(2), GridOut%type ) ) then ErrMsg2=trim(ErrMsg2)//NewLine//'GridType needs to be an integer.' @@ -315,7 +315,7 @@ subroutine ReadGridOut(sLine, GridOut) if ( index(StrArray(3), "DEFAULT" ) == 1 ) then GridOut%tStart = 0.0_ReKi else - if (.not. is_numeric(StrArray(3), GridOut%tStart) ) then + if (.not. is_numeric(StrArray(3), GridOut%tStart) ) then ErrMsg2=trim(ErrMsg2)//NewLine//'TStart needs to be numeric or "default".' return endif @@ -490,7 +490,7 @@ subroutine WrVTK_FVW(p, x, z, m, FileRootName, VTKcount, Twidth, bladeFrame, Hub Call ProgAbort('Programming error in WrVTK_FVW call: Cannot use the WrVTK_FVW with bladeFrame==TRUE without the optional arguments of HubOrientation and HubPosition') endif endif - + if (DEV_VERSION) then print*,'------------------------------------------------------------------------------' print'(A,L1,A,I0,A,I0,A,I0)','VTK Output - First call ',m%FirstCall, ' nNW:',m%nNW,' nFW:',m%nFW,' i:',VTKCount @@ -502,7 +502,7 @@ subroutine WrVTK_FVW(p, x, z, m, FileRootName, VTKcount, Twidth, bladeFrame, Hub write(Tstr, '(i' // trim(Num2LStr(Twidth)) //'.'// trim(Num2LStr(Twidth)) // ')') VTKcount ! --------------------------------------------------------------------------------} - ! --- Blade + ! --- Blade ! --------------------------------------------------------------------------------{ ! --- Blade Quarter chord points (AC) do iW=1,p%VTKBlades @@ -530,7 +530,7 @@ subroutine WrVTK_FVW(p, x, z, m, FileRootName, VTKcount, Twidth, bladeFrame, Hub ! call WrVTK_Lattice(FileName, mvtk, m%W(iW)%r_LL(1:3,:,:), m%W(iW)%Gamma_LL(:), bladeFrame=bladeFrame) ! enddo ! --------------------------------------------------------------------------------} - ! --- Near wake + ! --- Near wake ! --------------------------------------------------------------------------------{ ! --- Near wake panels do iW=1,p%VTKBlades @@ -548,7 +548,7 @@ subroutine WrVTK_FVW(p, x, z, m, FileRootName, VTKcount, Twidth, bladeFrame, Hub endif enddo ! --------------------------------------------------------------------------------} - ! --- Far wake + ! --- Far wake ! --------------------------------------------------------------------------------{ ! --- Far wake panels do iW=1,p%VTKBlades @@ -572,7 +572,7 @@ subroutine WrVTK_FVW(p, x, z, m, FileRootName, VTKcount, Twidth, bladeFrame, Hub endif if (nSeg>0) then Filename = TRIM(FileRootName)//'.AllSeg.'//Tstr//'.vtk' - CALL WrVTK_Segments(Filename, mvtk, m%Sgmt%Points(:,1:nSegP), m%Sgmt%Connct(:,1:nSeg), m%Sgmt%Gamma(1:nSeg), m%Sgmt%Epsilon(1:nSeg), bladeFrame) + CALL WrVTK_Segments(Filename, mvtk, m%Sgmt%Points(:,1:nSegP), m%Sgmt%Connct(:,1:nSeg), m%Sgmt%Gamma(1:nSeg), m%Sgmt%Epsilon(1:nSeg), bladeFrame) endif if (p%SrcPnl%n>0) then @@ -683,14 +683,14 @@ subroutine WrVTK_Panels(filename, mvtk, p, m, z) endsubroutine WrVTK_Panels -subroutine WrVTK_Segments(filename, mvtk, SegPoints, SegConnct, SegGamma, SegEpsilon, bladeFrame) +subroutine WrVTK_Segments(filename, mvtk, SegPoints, SegConnct, SegGamma, SegEpsilon, bladeFrame) use VTK character(len=*),intent(in) :: filename type(VTK_Misc), intent(inout) :: mvtk !< miscvars for VTK output - real(ReKi), dimension(:,:), intent(in) :: SegPoints !< - integer(IntKi), dimension(:,:), intent(in) :: SegConnct !< - real(ReKi), dimension(:) , intent(in) :: SegGamma !< - real(ReKi), dimension(:) , intent(in) :: SegEpsilon !< + real(ReKi), dimension(:,:), intent(in) :: SegPoints !< + integer(IntKi), dimension(:,:), intent(in) :: SegConnct !< + real(ReKi), dimension(:) , intent(in) :: SegGamma !< + real(ReKi), dimension(:) , intent(in) :: SegEpsilon !< logical, intent(in ) :: bladeFrame !< Output in blade coordinate frame if ( vtk_new_ascii_file(filename,'Sgmt',mvtk) ) then call vtk_dataset_polydata(SegPoints(1:3,:),mvtk,bladeFrame) From b204b29e3de18c2415be0c075b209de4e2264784 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 09:47:56 -0600 Subject: [PATCH 09/35] AeroDyn/OLAF: use resolved xPts/yPts/zPts in SetRequestedWindPoints --- modules/aerodyn/src/FVW_Subs.f90 | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/modules/aerodyn/src/FVW_Subs.f90 b/modules/aerodyn/src/FVW_Subs.f90 index 29b310b274..904fff13cd 100644 --- a/modules/aerodyn/src/FVW_Subs.f90 +++ b/modules/aerodyn/src/FVW_Subs.f90 @@ -572,7 +572,7 @@ subroutine SetRequestedWindPoints(r_wind, x, p, m) type(FVW_MiscVarType), intent(in ), target :: m !< Initial misc/optimization variables integer(IntKi) :: iP_start,iP_end ! Current index of point, start and end of range integer(IntKi) :: iGrid,i,j,k,iW - real(ReKi) :: xP,yP,zP,dx,dy,dz + real(ReKi) :: xP,yP,zP type(GridOutType), pointer :: g ! Using array reshaping to ensure a given near or far wake point is always at the same location in the array. @@ -608,15 +608,12 @@ subroutine SetRequestedWindPoints(r_wind, x, p, m) iP_start=iP_end+1 do iGrid=1,p%nGridOut g => m%GridOutputs(iGrid) - dx = (g%xEnd- g%xStart)/max(g%nx-1,1) - dy = (g%yEnd- g%yStart)/max(g%ny-1,1) - dz = (g%zEnd- g%zStart)/max(g%nz-1,1) do k=1,g%nz - zP = g%zStart + (k-1)*dz + zP = g%zPts(k) do j=1,g%ny - yP = g%yStart + (j-1)*dy + yP = g%yPts(j) do i=1,g%nx - xP = g%xStart + (i-1)*dx + xP = g%xPts(i) r_wind(1:3,iP_start) = (/xP,yP,zP/) iP_start=iP_start+1 enddo From ea53472ffdeab0cfd8da295aa338985e4953d43e Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 09:50:48 -0600 Subject: [PATCH 10/35] AeroDyn/OLAF: use resolved xPts/yPts/zPts in InducedVelocitiesAll_OnGrid --- modules/aerodyn/src/FVW_Subs.f90 | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/modules/aerodyn/src/FVW_Subs.f90 b/modules/aerodyn/src/FVW_Subs.f90 index 904fff13cd..e9b60f81d0 100644 --- a/modules/aerodyn/src/FVW_Subs.f90 +++ b/modules/aerodyn/src/FVW_Subs.f90 @@ -1135,7 +1135,7 @@ subroutine InducedVelocitiesAll_OnGrid(g, p, x, m, ErrStat, ErrMsg) ! Local variables integer(IntKi) :: nCPs, iHeadP integer(IntKi) :: i,j,k - real(ReKi) :: xP,yP,zP,dx,dy,dz + real(ReKi) :: xP,yP,zP ! TODO new options type(T_Tree) :: Tree type(T_Panl) :: Panl @@ -1149,15 +1149,12 @@ subroutine InducedVelocitiesAll_OnGrid(g, p, x, m, ErrStat, ErrMsg) nCPs = g%nx * g%ny * g%nz allocate(CPs(3, nCPs), stat=ErrStat) iHeadP=1 - dx = (g%xEnd- g%xStart)/max(g%nx-1,1) - dy = (g%yEnd- g%yStart)/max(g%ny-1,1) - dz = (g%zEnd- g%zStart)/max(g%nz-1,1) do k=1,g%nz - zP = g%zStart + (k-1)*dz + zP = g%zPts(k) do j=1,g%ny - yP = g%yStart + (j-1)*dy + yP = g%yPts(j) do i=1,g%nx - xP = g%xStart + (i-1)*dx + xP = g%xPts(i) CPs(1:3,iHeadP) = (/xP,yP,zP/) iHeadP=iHeadP+1 enddo From edb2fc0ae3aeb1ffb2f42c336fadae69d80b8e72 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 10:01:17 -0600 Subject: [PATCH 11/35] FVW_Subs.f90: trim trailing whitespace --- modules/aerodyn/src/FVW_Subs.f90 | 120 +++++++++++++++---------------- 1 file changed, 60 insertions(+), 60 deletions(-) diff --git a/modules/aerodyn/src/FVW_Subs.f90 b/modules/aerodyn/src/FVW_Subs.f90 index e9b60f81d0..3c9061e537 100644 --- a/modules/aerodyn/src/FVW_Subs.f90 +++ b/modules/aerodyn/src/FVW_Subs.f90 @@ -504,7 +504,7 @@ subroutine find_nan_1D(array, varname) tot=tot+1 endif enddo - if (found) then + if (found) then print*,'OLAF NAN ',trim(varname),tot,n STOP endif @@ -526,7 +526,7 @@ subroutine find_nan_2D(array, varname) tot=tot+1 endif enddo - if (found) then + if (found) then print*,'OLAF NAN ',trim(varname),tot,n STOP endif @@ -784,8 +784,8 @@ subroutine FVW_InitStates( x, p, ErrStat, ErrMsg ) allocate(x%W(p%nWings)) do iW=1,p%nWings - call AllocAry( x%W(iW)%Gamma_NW, p%W(iW)%nSpan , p%nNWMax , 'NW Panels Circulation', ErrStat2, ErrMsg2 );call SetErrStat ( ErrStat2, ErrMsg2, ErrStat,ErrMsg,'FVW_InitStates' ); - call AllocAry( x%W(iW)%Gamma_FW, FWnSpan , p%nFWMax , 'FW Panels Circulation', ErrStat2, ErrMsg2 );call SetErrStat ( ErrStat2, ErrMsg2, ErrStat,ErrMsg,'FVW_InitStates' ); + call AllocAry( x%W(iW)%Gamma_NW, p%W(iW)%nSpan , p%nNWMax , 'NW Panels Circulation', ErrStat2, ErrMsg2 );call SetErrStat ( ErrStat2, ErrMsg2, ErrStat,ErrMsg,'FVW_InitStates' ); + call AllocAry( x%W(iW)%Gamma_FW, FWnSpan , p%nFWMax , 'FW Panels Circulation', ErrStat2, ErrMsg2 );call SetErrStat ( ErrStat2, ErrMsg2, ErrStat,ErrMsg,'FVW_InitStates' ); call AllocAry( x%W(iW)%Eps_NW , 3, p%W(iW)%nSpan , p%nNWMax , 'NW Panels Reg Param' , ErrStat2, ErrMsg2 );call SetErrStat ( ErrStat2, ErrMsg2, ErrStat,ErrMsg,'FVW_InitStates' ); call AllocAry( x%W(iW)%Eps_FW , 3, FWnSpan , p%nFWMax , 'FW Panels Reg Param' , ErrStat2, ErrMsg2 );call SetErrStat ( ErrStat2, ErrMsg2, ErrStat,ErrMsg,'FVW_InitStates' ); ! set x%W(iW)%r_NW and x%W(iW)%r_FW to (0,0,0) so that InflowWind can shortcut the calculations @@ -794,10 +794,10 @@ subroutine FVW_InitStates( x, p, ErrStat, ErrMsg ) if (ErrStat >= AbortErrLev) return x%W(iW)%r_NW = 0.0_ReKi x%W(iW)%r_FW = 0.0_ReKi - x%W(iW)%Gamma_NW = 0.0_ReKi ! First call of calcoutput, states might not be set + x%W(iW)%Gamma_NW = 0.0_ReKi ! First call of calcoutput, states might not be set x%W(iW)%Gamma_FW = 0.0_ReKi ! NOTE, these values might be mapped from z%W(iW)%Gamma_LL at init - x%W(iW)%Eps_NW = 0.001_ReKi - x%W(iW)%Eps_FW = 0.001_ReKi + x%W(iW)%Eps_NW = 0.001_ReKi + x%W(iW)%Eps_FW = 0.001_ReKi enddo end subroutine FVW_InitStates @@ -836,7 +836,7 @@ subroutine FVW_InitMiscVarsPostParam( p, m, ErrStat, ErrMsg ) bWakeNeedsPart = p%VelocityMethod(1)==idVelocityPart .or. p%VelocityMethod(1)==idVelocityTreePart bLLNeedsPart = p%VelocityMethod(2)==idVelocityPart .or. p%VelocityMethod(2)==idVelocityTreePart if (bLLNeedsPart .or. bWakeNeedsPart) then - nPart = 0 + nPart = 0 if (bWakeNeedsPart) nPart = max(nPart, nSeg * p%PartPerSegment(1)) if (bLLNeedsPart) nPart = max(nPart, nSeg * p%PartPerSegment(2)) call AllocAry( m%Part%P , 3, nPart, 'PartP' , ErrStat2, ErrMsg2 ); if(Failed())return; m%Part%P = -999999_ReKi; @@ -851,7 +851,7 @@ subroutine FVW_InitMiscVarsPostParam( p, m, ErrStat, ErrMsg ) call AllocAry( m%Uind , 3, nCPs, 'Uind' , ErrStat2, ErrMsg2 ); if(Failed())return; m%Uind= -999999_ReKi; contains logical function Failed() - call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'FVW_InitMiscVarsPostParam') + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'FVW_InitMiscVarsPostParam') Failed = ErrStat >= AbortErrLev end function Failed end subroutine FVW_InitMiscVarsPostParam @@ -1205,7 +1205,7 @@ subroutine SegmentsToPartWrap(Sgmt, nSeg, PartPerSegment, RegFunction, Part, all else ! check that we have enough space if (.not. allocated(Part%P)) then - print*,'>>> PartP not allocated'; + print*,'>>> PartP not allocated'; STOP endif if (size(Part%P,2) Initialize / allocated main variables for source panels. +!> Initialize / allocated main variables for source panels. !! If an non empty input file is provided, the panels points and connectivity are read !! Otherwise, Points and IDs should be provided in "p" -!! -!! Acknowledgements: +!! +!! Acknowledgements: !! The original implementation of the source panel method was funded by Accelerate Wind, !! and implemented by E. Branlard. !! For more acknowledgements, visit: https://openfast.readthedocs.io/en/main/source/acknowledgements.html -!! +!! subroutine srcPnl_init(p, m, z, errStat, errMsg, filename) use VTK !, only: ReadVTK_PD_info, ReadVTK_PD_fields type(T_SrcPanlParam), intent(inout) :: p @@ -1894,7 +1894,7 @@ subroutine srcPnl_init(p, m, z, errStat, errMsg, filename) ! --- Compute influence matrix ! For now, panels don't move so we compute this only once, otherwise, put this in FVW_CalcConstrStateResidual - call srcPnl_build_mat(p, m%AI, m%UUI) + call srcPnl_build_mat(p, m%AI, m%UUI) ! --- Factorization call linalg_factor(m%AI, m%IPIV, errStat2, errMsg2); if(Failed()) return @@ -1912,27 +1912,27 @@ subroutine srcPnl_geometry(Panl, errStat, errMsg) type(T_SrcPanlParam), intent(inout) :: Panl integer(IntKi) , intent(out) :: errStat !< Error status of the operation character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None - real(ReKi) :: alpha !< - real(ReKi) :: d1 !< - real(ReKi) :: DLastRingTE !< - real(ReKi) :: eta0 !< - integer(IntKi), dimension(4) :: IDs !< - real(ReKi), dimension(3) :: P1 !< - real(ReKi), dimension(3) :: P1e !< - real(ReKi), dimension(3) :: P1es !< - real(ReKi), dimension(3) :: P1p !< - real(ReKi), dimension(3) :: P2 !< - real(ReKi), dimension(3) :: P2e !< - real(ReKi), dimension(3) :: P2es !< - real(ReKi), dimension(3) :: P2p !< - real(ReKi), dimension(3) :: P3 !< - real(ReKi), dimension(3) :: P3e !< - real(ReKi), dimension(3) :: P3es !< - real(ReKi), dimension(3) :: P3p !< - real(ReKi), dimension(3) :: P4 !< - real(ReKi), dimension(3) :: P4e !< - real(ReKi), dimension(3) :: P4es !< - real(ReKi), dimension(3) :: P4p !< + real(ReKi) :: alpha !< + real(ReKi) :: d1 !< + real(ReKi) :: DLastRingTE !< + real(ReKi) :: eta0 !< + integer(IntKi), dimension(4) :: IDs !< + real(ReKi), dimension(3) :: P1 !< + real(ReKi), dimension(3) :: P1e !< + real(ReKi), dimension(3) :: P1es !< + real(ReKi), dimension(3) :: P1p !< + real(ReKi), dimension(3) :: P2 !< + real(ReKi), dimension(3) :: P2e !< + real(ReKi), dimension(3) :: P2es !< + real(ReKi), dimension(3) :: P2p !< + real(ReKi), dimension(3) :: P3 !< + real(ReKi), dimension(3) :: P3e !< + real(ReKi), dimension(3) :: P3es !< + real(ReKi), dimension(3) :: P3p !< + real(ReKi), dimension(3) :: P4 !< + real(ReKi), dimension(3) :: P4e !< + real(ReKi), dimension(3) :: P4es !< + real(ReKi), dimension(3) :: P4p !< real(ReKi), dimension(3) :: T1 real(ReKi), dimension(3) :: T2 real(ReKi), dimension(3) :: Ptmp !< temp for computing norm of delta points @@ -1941,7 +1941,7 @@ subroutine srcPnl_geometry(Panl, errStat, errMsg) real(ReKi) :: norm_T1 real(ReKi) :: norm_T2 integer :: ip - real(ReKi) :: xi0 !< + real(ReKi) :: xi0 !< integer(IntKi) :: errStat2 !< temporary Error status character(ErrMsgLen) :: errMsg2 !< temporary Error message errStat = ErrID_None @@ -1964,7 +1964,7 @@ subroutine srcPnl_geometry(Panl, errStat, errMsg) call AllocAry(Panl%eta , 4, Panl%n,'eta ' ,errStat2,errMsg2); if(Failed())return call AllocAry(Panl%xi , 4, Panl%n,'xi ' ,errStat2,errMsg2); if(Failed())return call AllocAry(Panl%Area , Panl%n,'Area ' ,errStat2,errMsg2); if(Failed())return - + do ip = 1, Panl%n IDs = Panl%IDs(:,ip) P1 = Panl%P(:,IDs(1)) @@ -1974,10 +1974,10 @@ subroutine srcPnl_geometry(Panl, errStat, errMsg) Panl%Pmid(1:3,ip) = (P1+P2+P3+p4)/4 ! that's hess's barred coordinates T1 = P3-P1 T2 = P4-P2 - ! maximum diagonal + ! maximum diagonal norm_T1 = sqrt(T1(1)**2+ T1(2)**2+ T1(3)**2) norm_T2 = sqrt(T2(1)**2+ T2(2)**2+ T2(3)**2) - ! flat panel coordinate system + ! flat panel coordinate system Ptmp(1) = T2(2) * T1(3) - T2(3) * T1(2) Ptmp(2) = T2(3) * T1(1) - T2(1) * T1(3) Ptmp(3) = T2(1) * T1(2) - T2(2) * T1(1) @@ -1995,22 +1995,22 @@ subroutine srcPnl_geometry(Panl, errStat, errMsg) Panl%R_g2p(2, 1:3, ip) = T2 Panl%R_g2p(3, 1:3, ip) = Panl%Normal(:,ip) Mat = Panl%R_g2p(:,:,ip) - ! Projection of the surface into a flat panel - Hess primed coordinates + ! Projection of the surface into a flat panel - Hess primed coordinates d1 = dot_product(Panl%Normal(:,ip),Panl%Pmid(:,ip)-P1) P1p = P1+Panl%Normal(:,ip)*(-1)**(1-1)*d1 P2p = P2+Panl%Normal(:,ip)*(-1)**(2-1)*d1 P3p = P3+Panl%Normal(:,ip)*(-1)**(3-1)*d1 P4p = P4+Panl%Normal(:,ip)*(-1)**(4-1)*d1 - !Coordinates of flat panel points in panel coordinate system - Hess starred coordinates with greek letters - ! the transformation is such that the zeta coordinate will always be zero + !Coordinates of flat panel points in panel coordinate system - Hess starred coordinates with greek letters + ! the transformation is such that the zeta coordinate will always be zero P1es = matmul(Mat,(P1p-Panl%Pmid(:,ip))) P2es = matmul(Mat,(P2p-Panl%Pmid(:,ip))) P3es = matmul(Mat,(P3p-Panl%Pmid(:,ip))) P4es = matmul(Mat,(P4p-Panl%Pmid(:,ip))) - ! Coordinates of the centroid + ! Coordinates of the centroid xi0 = 1._ReKi/3._ReKi*1.0_ReKi/(P2es(2)-P4es(2)) * (P4es(1)*(P1es(2)-P2es(2))+P2es(1)*(P4es(2)-P1es(2) )) eta0 = -1._ReKi/3._ReKi * P1es(2) - ! Coordinates based on centroid - Hess greek letters coordinates + ! Coordinates based on centroid - Hess greek letters coordinates P1e = P1es-(/ xi0,eta0,0.0_ReKi /) P2e = P2es-(/ xi0,eta0,0.0_ReKi /) P3e = P3es-(/ xi0,eta0,0.0_ReKi /) @@ -2019,7 +2019,7 @@ subroutine srcPnl_geometry(Panl, errStat, errMsg) Panl%eta(:,ip) = (/ P1e(2), P2e(2), P3e(2), P4e(2)/) ! Centroid in reference frame Panl%Pcent(:,ip) = Panl%Pmid(:,ip) + matmul( (/ xi0,eta0,0.0_ReKi /),Mat) - ! Area + ! Area Panl%Area(ip) = 0.5_ReKi*(Panl%xi(3,ip)-Panl%xi(1,ip))*(Panl%eta(2,ip)-Panl%eta(4,ip)) end do ! Loop on panels contains @@ -2035,8 +2035,8 @@ subroutine srcPnl_build_mat(Panl, AI, UUI) real(ReKi), dimension(:,:), intent(out) :: AI !< (nCPs x nPanels) Self Induced Velocities matrix along normal real(ReKi), dimension(:,:,:), intent(out) :: UUI !< (3 x nCPs x nPanels) Unit induced velocity ! Variables - real(ReKi), parameter :: UnitIntensity=1.0_ReKi !< - real(ReKi), dimension(3) :: Uind_tmp !< + real(ReKi), parameter :: UnitIntensity=1.0_ReKi !< + real(ReKi), dimension(3) :: Uind_tmp !< integer :: icp, ip !< loop variables if (OLAF_PROFILING) call tic('SrcPanel build matrix') !$OMP PARALLEL DEFAULT(shared) @@ -2047,12 +2047,12 @@ subroutine srcPnl_build_mat(Panl, AI, UUI) ! ---- loop on all panls do ip = 1, Panl%n call ui_quad_src_11(Panl%Pcent(:,icp), UnitIntensity, Panl%xi(1:4,ip), Panl%eta(1:4,ip), Panl%Pcent(1:3,ip), Panl%R_g2p(1:3,1:3,ip), Uind_tmp) - ! AI= Vi . N + ! AI= Vi . N AI(icp, ip) = dot_product(Uind_tmp, Panl%Normal(1:3,icp)) UUI(:, icp, ip) = Uind_tmp - end do - end do - !$OMP END DO + end do + end do + !$OMP END DO !$OMP END PARALLEL if (OLAF_PROFILING) call toc() end subroutine srcPnl_build_mat @@ -2079,7 +2079,7 @@ subroutine srcPnl_ExtVelocities_OnPanels(u, p, x, m, errStat, errMsg) m%SrcPnl%Uext(1:3,:) = 0.0_ReKi ! Due to side effects of ui_ functions ! Convert Panels to segments, segments to particles, particles to tree call InducedVelocitiesAll_Init(p, x, m, m%Sgmt, m%Part, Tree, Panl, errStat, errMsg, allocPart=.false.) - ! We don't want the influence of panels so we nullify + ! We don't want the influence of panels so we nullify nullify(Panl%p_Src); nullify(Panl%m_Src) call InducedVelocitiesAll_Calc(p%SrcPnl%Pcent(1:3,:), p%SrcPnl%n, m%SrcPnl%Uext, p, m%Sgmt, m%Part, Tree, Panl, errStat, errMsg) call InducedVelocitiesAll_End(p, Tree, m%Part, Panl, errStat, errMsg, deallocPart=.false.) @@ -2139,7 +2139,7 @@ subroutine srcPnl_calcOutput(p, m, z, rho) !, errStat, errMsg ! Static pressure ps = 1/2 rho Utot**2 ! Reference pressure qinf = 1/2 rho Uwnd**2 ! Pressure Coefficient Cp = ps-pinf/qinf (by convention) - ! Pressure force F = (ps-p0) A e_n + ! Pressure force F = (ps-p0) A e_n if (Uwnd_norm2>0) then Cp = max( 1-(norm2(m%Utot(1:3,ip))**2)/Uwnd_norm2, -10._ReKi) ! Bernoulli. else @@ -2167,7 +2167,7 @@ subroutine linalg_factor(AA, IPIV, errStat, errMsg) n = size(AA,1) m = n call LAPACK_GETRF(m, n, AA, IPIV, errStat, errMsg) -endsubroutine +endsubroutine subroutine linalg_solve(AFact, RHS, IPIV, errStat, errMsg) use NWTC_LAPACK, only : LAPACK_GETRS @@ -2179,7 +2179,7 @@ subroutine linalg_solve(AFact, RHS, IPIV, errStat, errMsg) integer :: n n = size(AFact,1) call LAPACK_GETRS('N', n, AFact, IPIV, RHS, errStat, errMsg) -end subroutine +end subroutine !> Solve A x = B subroutine linalg_solveWrap(AA, B, X, errStat, errMsg) From fdab91ba29cc71e8bd22e1fd82f4e246daf59fa7 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 12:40:20 -0600 Subject: [PATCH 12/35] FVW_Types.f90: updated FVW registry outputs --- modules/aerodyn/src/FVW_Types.f90 | 66 +++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/modules/aerodyn/src/FVW_Types.f90 b/modules/aerodyn/src/FVW_Types.f90 index 53ff476c2f..9f0b0b1aa7 100644 --- a/modules/aerodyn/src/FVW_Types.f90 +++ b/modules/aerodyn/src/FVW_Types.f90 @@ -56,6 +56,12 @@ MODULE FVW_Types REAL(ReKi) , DIMENSION(:,:,:,:), ALLOCATABLE :: uGrid !< Grid velocity 3 x nz x ny x nx [-] REAL(ReKi) , DIMENSION(:,:,:,:), ALLOCATABLE :: omGrid !< Grid vorticity 3 x nz x ny x nx [-] REAL(DbKi) :: tLastOutput = 0.0_R8Ki !< Last output time [-] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: xPts !< Explicit x coordinates (non-equidistant grid, if used) [m] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: yPts !< Explicit y coordinates (non-equidistant grid, if used) [m] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: zPts !< Explicit z coordinates (non-equidistant grid, if used) [m] + CHARACTER(1024) :: xListFile !< File with explicit x coordinates (empty if equidistant) [-] + CHARACTER(1024) :: yListFile !< File with explicit y coordinates (empty if equidistant) [-] + CHARACTER(1024) :: zListFile !< File with explicit z coordinates (empty if equidistant) [-] END TYPE GridOutType ! ======================= ! ========= T_Sgmt ======= @@ -466,6 +472,45 @@ subroutine FVW_CopyGridOutType(SrcGridOutTypeData, DstGridOutTypeData, CtrlCode, DstGridOutTypeData%omGrid = SrcGridOutTypeData%omGrid end if DstGridOutTypeData%tLastOutput = SrcGridOutTypeData%tLastOutput + if (allocated(SrcGridOutTypeData%xPts)) then + LB(1:1) = lbound(SrcGridOutTypeData%xPts) + UB(1:1) = ubound(SrcGridOutTypeData%xPts) + if (.not. allocated(DstGridOutTypeData%xPts)) then + allocate(DstGridOutTypeData%xPts(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstGridOutTypeData%xPts.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstGridOutTypeData%xPts = SrcGridOutTypeData%xPts + end if + if (allocated(SrcGridOutTypeData%yPts)) then + LB(1:1) = lbound(SrcGridOutTypeData%yPts) + UB(1:1) = ubound(SrcGridOutTypeData%yPts) + if (.not. allocated(DstGridOutTypeData%yPts)) then + allocate(DstGridOutTypeData%yPts(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstGridOutTypeData%yPts.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstGridOutTypeData%yPts = SrcGridOutTypeData%yPts + end if + if (allocated(SrcGridOutTypeData%zPts)) then + LB(1:1) = lbound(SrcGridOutTypeData%zPts) + UB(1:1) = ubound(SrcGridOutTypeData%zPts) + if (.not. allocated(DstGridOutTypeData%zPts)) then + allocate(DstGridOutTypeData%zPts(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstGridOutTypeData%zPts.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstGridOutTypeData%zPts = SrcGridOutTypeData%zPts + end if + DstGridOutTypeData%xListFile = SrcGridOutTypeData%xListFile + DstGridOutTypeData%yListFile = SrcGridOutTypeData%yListFile + DstGridOutTypeData%zListFile = SrcGridOutTypeData%zListFile end subroutine subroutine FVW_DestroyGridOutType(GridOutTypeData, ErrStat, ErrMsg) @@ -481,6 +526,15 @@ subroutine FVW_DestroyGridOutType(GridOutTypeData, ErrStat, ErrMsg) if (allocated(GridOutTypeData%omGrid)) then deallocate(GridOutTypeData%omGrid) end if + if (allocated(GridOutTypeData%xPts)) then + deallocate(GridOutTypeData%xPts) + end if + if (allocated(GridOutTypeData%yPts)) then + deallocate(GridOutTypeData%yPts) + end if + if (allocated(GridOutTypeData%zPts)) then + deallocate(GridOutTypeData%zPts) + end if end subroutine subroutine FVW_PackGridOutType(RF, Indata) @@ -505,6 +559,12 @@ subroutine FVW_PackGridOutType(RF, Indata) call RegPackAlloc(RF, InData%uGrid) call RegPackAlloc(RF, InData%omGrid) call RegPack(RF, InData%tLastOutput) + call RegPackAlloc(RF, InData%xPts) + call RegPackAlloc(RF, InData%yPts) + call RegPackAlloc(RF, InData%zPts) + call RegPack(RF, InData%xListFile) + call RegPack(RF, InData%yListFile) + call RegPack(RF, InData%zListFile) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -533,6 +593,12 @@ subroutine FVW_UnPackGridOutType(RF, OutData) call RegUnpackAlloc(RF, OutData%uGrid); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%omGrid); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%tLastOutput); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%xPts); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%yPts); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%zPts); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%xListFile); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%yListFile); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%zListFile); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine FVW_CopyT_Sgmt(SrcT_SgmtData, DstT_SgmtData, CtrlCode, ErrStat, ErrMsg) From 242241e824a8a2a7c2bac5f28718140fed30ee72 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 13:07:20 -0600 Subject: [PATCH 13/35] AeroDyn/OLAF: fix ResolveGridAxis call (5 arguments) --- modules/aerodyn/src/FVW_IO.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 1ea8dad016..cc43223ba7 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -106,11 +106,11 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) call ReadGridOut(sLine, m%GridOutputs(i)); if(Failed()) return ! Resolve each axis to an explicit coordinate array (regardless of whether it is a list file or a range) call ResolveGridAxis(m%GridOutputs(i)%xStart, m%GridOutputs(i)%xEnd, m%GridOutputs(i)%nx, & - m%GridOutputs(i)%xListFile, PriPath, m%GridOutputs(i)%xPts, ErrStat2, ErrMsg2); if(Failed()) return + m%GridOutputs(i)%xListFile, m%GridOutputs(i)%xPts); if(Failed()) return call ResolveGridAxis(m%GridOutputs(i)%yStart, m%GridOutputs(i)%yEnd, m%GridOutputs(i)%ny, & - m%GridOutputs(i)%yListFile, PriPath, m%GridOutputs(i)%yPts, ErrStat2, ErrMsg2); if(Failed()) return + m%GridOutputs(i)%yListFile, m%GridOutputs(i)%yPts); if(Failed()) return call ResolveGridAxis(m%GridOutputs(i)%zStart, m%GridOutputs(i)%zEnd, m%GridOutputs(i)%nz, & - m%GridOutputs(i)%zListFile, PriPath, m%GridOutputs(i)%zPts, ErrStat2, ErrMsg2); if(Failed()) return + m%GridOutputs(i)%zListFile, m%GridOutputs(i)%zPts); if(Failed()) return ! Error checking if (Check(m%GridOutputs(i)%nx<1, 'Grid output nx needs to be >=1')) return if (Check(m%GridOutputs(i)%ny<1, 'Grid output ny needs to be >=1')) return From 7244c06479cb90881ae1bf2b39922f1f6b071a1b Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 14:21:05 -0600 Subject: [PATCH 14/35] AeroDyn/OLAF docs: document non-equidistant grid output (list-file) --- docs/source/user/aerodyn-olaf/InputFiles.rst | 53 ++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/docs/source/user/aerodyn-olaf/InputFiles.rst b/docs/source/user/aerodyn-olaf/InputFiles.rst index 3c6fabd8ad..afc81f7b85 100644 --- a/docs/source/user/aerodyn-olaf/InputFiles.rst +++ b/docs/source/user/aerodyn-olaf/InputFiles.rst @@ -312,6 +312,59 @@ of a box of shape 5x20x30 and dimension 1200x300x295. The grid contains both th The two other grids are vertical and horizontal planes containing only the velocity. +Non-equidistant grid points +~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +By default, grid points along each axis (X, Y, Z) are equidistant, defined by +the ``Start``, ``End``, and ``n`` columns of the grid output table. + +Alternatively, a ``Start`` cell can be given as a quoted filename instead of a +number. In that case, the corresponding axis uses an explicit, user-defined +list of coordinates read from that file, and the ``End``/``n`` columns for +that axis are ignored. + +This choice is made independently for each axis, so a single grid can mix +equidistant and list-defined axes (e.g., equidistant in X and Z, list-defined +in Y), and a filename can be given for ``XStart``, ``YStart``, ``ZStart``, or +any combination of the three. + +The referenced file must: + +- contain exactly one coordinate value per line, +- list values in strictly ascending order (no duplicates), +- contain no comments or blank lines, +- be located in the same directory as the main OLAF input file. + +Example, requesting a non-equidistant Y-axis:: + + GridName GridType TStart TEnd DTOut XStart XEnd nX YStart YEnd nY ZStart ZEnd nZ + (-) (-) (s) (s) (s) (m) (m) (-) (m) (m) (-) (m) (m) (-) + "Yline" 1 default default default 90 90 1 "Ypoints.dat" - - 100 100 1 + +with ``Ypoints.dat`` (located next to the OLAF input file) containing:: + + -125 + -100 + -75 + -50 + -37.5 + -25 + -12.5 + 0 + 12.5 + 25 + 37.5 + 50 + 75 + 100 + 125 + +.. note:: + Vorticity output (``GridType=2``) requires equidistant spacing + and is not available for a grid that uses a list-defined axis. + Use ``GridType=1`` (velocity only) for such grids. + + Advanced Options ~~~~~~~~~~~~~~~~ From 4fcd5a1da046f0f8b85b58a3461b124f290a012f Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 21:55:12 -0600 Subject: [PATCH 15/35] AeroDyn/OLAF: make error reporting explicit in ResolveGridAxis and guard all allocations --- modules/aerodyn/src/FVW_IO.f90 | 68 ++++++++++++++++++++++++---------- 1 file changed, 48 insertions(+), 20 deletions(-) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index cc43223ba7..a7c252bf04 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -106,16 +106,16 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) call ReadGridOut(sLine, m%GridOutputs(i)); if(Failed()) return ! Resolve each axis to an explicit coordinate array (regardless of whether it is a list file or a range) call ResolveGridAxis(m%GridOutputs(i)%xStart, m%GridOutputs(i)%xEnd, m%GridOutputs(i)%nx, & - m%GridOutputs(i)%xListFile, m%GridOutputs(i)%xPts); if(Failed()) return + m%GridOutputs(i)%xListFile, m%GridOutputs(i)%xPts, ErrStat2, ErrMsg2); if(Failed()) return call ResolveGridAxis(m%GridOutputs(i)%yStart, m%GridOutputs(i)%yEnd, m%GridOutputs(i)%ny, & - m%GridOutputs(i)%yListFile, m%GridOutputs(i)%yPts); if(Failed()) return + m%GridOutputs(i)%yListFile, m%GridOutputs(i)%yPts, ErrStat2, ErrMsg2); if(Failed()) return call ResolveGridAxis(m%GridOutputs(i)%zStart, m%GridOutputs(i)%zEnd, m%GridOutputs(i)%nz, & - m%GridOutputs(i)%zListFile, m%GridOutputs(i)%zPts); if(Failed()) return + m%GridOutputs(i)%zListFile, m%GridOutputs(i)%zPts, ErrStat2, ErrMsg2); if(Failed()) return ! Error checking if (Check(m%GridOutputs(i)%nx<1, 'Grid output nx needs to be >=1')) return if (Check(m%GridOutputs(i)%ny<1, 'Grid output ny needs to be >=1')) return if (Check(m%GridOutputs(i)%nz<1, 'Grid output nz needs to be >=1')) return - ! Vorticity requires equidistant spacing. GridType must be 1 when using list files. + ! Vorticity (GridType=2) requires equidistant spacing. GridType must be 1 when using list files. if (m%GridOutputs(i)%type==idGridVelVorticity) then if (Check( len_trim(m%GridOutputs(i)%xListFile)>0 .or. & len_trim(m%GridOutputs(i)%yListFile)>0 .or. & @@ -385,31 +385,53 @@ end subroutine ReadGridOut ! Resolve one grid axis to an explicit, strictly increasing coordinate array, ! either read from ListFile (if given) or built from an equidistant Start/End/n range. - subroutine ResolveGridAxis(AStart, AEnd, n, ListFile, Pts) + subroutine ResolveGridAxis(AStart, AEnd, n, ListFile, Pts, ErrStat, ErrMsg) real(ReKi), intent(in) :: AStart, AEnd integer(IntKi), intent(inout) :: n character(*), intent(in) :: ListFile real(ReKi), allocatable, intent(out) :: Pts(:) - character(1024) :: FullFile - integer(IntKi) :: UnList, j, IOS - real(ReKi) :: val + integer(IntKi), intent(out) :: ErrStat + character(*), intent(out) :: ErrMsg + ! Locals + character(1024) :: FullFile + integer(IntKi) :: UnList, j, IOS + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + real(ReKi) :: val + + ErrStat = ErrID_None + ErrMsg = '' if (len_trim(ListFile) == 0) then ! Equidistant points, expressed as an explicit array - allocate(Pts(max(n,1))) + if (n < 1) then + call SetErrStat(ErrID_Fatal, 'ResolveGridAxis: grid axis definition requires n >= 1 (got '//trim(Num2LStr(n))//').', & + ErrStat, ErrMsg, 'ResolveGridAxis') + return + endif + allocate(Pts(n), stat=IOS) + if (IOS /= 0) then + call SetErrStat(ErrID_Fatal, 'ResolveGridAxis: error allocating grid axis points array (n='//trim(Num2LStr(n))//').', & + ErrStat, ErrMsg, 'ResolveGridAxis') + return + endif do j = 1, n Pts(j) = AStart + (AEnd - AStart) * real(j-1,ReKi) / real(max(n-1,1),ReKi) enddo - ErrStat2 = ErrID_None; ErrMsg2 = '' return endif ! Explicit list from file FullFile = ListFile if (PathIsRelative(FullFile)) FullFile = trim(PriPath)//trim(FullFile) + call GetNewUnit(UnList) call OpenFInpFile(UnList, trim(FullFile), ErrStat2, ErrMsg2) - if (ErrStat2 /= ErrID_None) return + if (ErrStat2 /= ErrID_None) then + call SetErrStat(ErrID_Fatal, 'ResolveGridAxis: could not open grid point list file "'//trim(FullFile)//'": '//trim(ErrMsg2), & + ErrStat, ErrMsg, 'ResolveGridAxis') + return + endif ! Count valid lines n = 0 @@ -418,16 +440,23 @@ subroutine ResolveGridAxis(AStart, AEnd, n, ListFile, Pts) if (IOS /= 0) exit n = n + 1 enddo + if (n < 1) then - ErrStat2 = ErrID_Fatal - ErrMsg2 = 'Grid point list file "'//trim(FullFile)//'" contains no valid values.' + call SetErrStat(ErrID_Fatal, 'ResolveGridAxis: grid point list file "'//trim(FullFile)//'" contains no valid values.', & + ErrStat, ErrMsg, 'ResolveGridAxis') close(UnList) return endif ! Read values rewind(UnList) - allocate(Pts(n)) + allocate(Pts(n), stat=IOS) + if (IOS /= 0) then + call SetErrStat(ErrID_Fatal, 'ResolveGridAxis: error allocating grid axis points array (n='//trim(Num2LStr(n))//').', & + ErrStat, ErrMsg, 'ResolveGridAxis') + close(UnList) + return + endif do j = 1, n read(UnList, *) Pts(j) enddo @@ -436,16 +465,15 @@ subroutine ResolveGridAxis(AStart, AEnd, n, ListFile, Pts) ! Validate strictly increasing, no duplicates do j = 2, n if (Pts(j) <= Pts(j-1)) then - ErrStat2 = ErrID_Fatal - ErrMsg2 = 'Grid point list file "'//trim(FullFile)//'" must be strictly increasing (no duplicates); '// & - 'value at line '//trim(Num2LStr(j))//' ('//trim(Num2LStr(Pts(j)))//') is not greater than '// & - 'the previous value ('//trim(Num2LStr(Pts(j-1)))//').' + call SetErrStat(ErrID_Fatal, & + 'ResolveGridAxis: grid point list file "'//trim(FullFile)//'" must be strictly increasing (no duplicates); '// & + 'value at line '//trim(Num2LStr(j))//' ('//trim(Num2LStr(Pts(j)))//') is not greater than '// & + 'the previous value ('//trim(Num2LStr(Pts(j-1)))//').', & + ErrStat, ErrMsg, 'ResolveGridAxis') return endif enddo - ErrStat2 = ErrID_None - ErrMsg2 = '' end subroutine ResolveGridAxis END SUBROUTINE FVW_ReadInputFile From 37ea967ac498f0b1a774ec989de95fee954a988a Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Mon, 3 Aug 2026 14:28:55 -0600 Subject: [PATCH 16/35] Increased length for console output --- modules/nwtc-library/src/SysIVF.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/nwtc-library/src/SysIVF.f90 b/modules/nwtc-library/src/SysIVF.f90 index 41af20386a..f1fc8e067f 100644 --- a/modules/nwtc-library/src/SysIVF.f90 +++ b/modules/nwtc-library/src/SysIVF.f90 @@ -55,7 +55,7 @@ MODULE SysSubs MODULE PROCEDURE NWTC_gammaR8 END INTERFACE - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 300 ! The record length for console output. INTEGER, PUBLIC :: CU = 7 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr LOGICAL, PARAMETER :: KBInputOK = .TRUE. ! A flag to tell the program that keyboard input is allowed in the environment. From 16a5cd40ae0251f0d856382e49fbd2e2495e4edb Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Mon, 3 Aug 2026 14:30:58 -0600 Subject: [PATCH 17/35] Delete trailing whitespaces --- modules/nwtc-library/src/SysIVF.f90 | 70 ++++++++++++++--------------- 1 file changed, 35 insertions(+), 35 deletions(-) diff --git a/modules/nwtc-library/src/SysIVF.f90 b/modules/nwtc-library/src/SysIVF.f90 index f1fc8e067f..4ecc6d0c04 100644 --- a/modules/nwtc-library/src/SysIVF.f90 +++ b/modules/nwtc-library/src/SysIVF.f90 @@ -20,19 +20,19 @@ MODULE SysSubs ! This module contains routines with system-specific logic and references, including all references to the console unit, CU. ! It also contains standard (but not system-specific) routines it uses. - ! SysIVF.f90 is specifically for the Intel Visual Fortran for Windows compiler. + ! SysIVF.f90 is specifically for the Intel Visual Fortran for Windows compiler. ! It contains the following routines: ! FUNCTION FileSize( Unit ) ! Returns the size (in bytes) of an open file. ! FUNCTION Is_NaN( DblNum ) ! Please use IEEE_IS_NAN() instead ! FUNCTION NWTC_ERF( x ) - ! FUNCTION NWTC_gamma( x ) ! Returns the gamma value of its argument. + ! FUNCTION NWTC_gamma( x ) ! Returns the gamma value of its argument. ! SUBROUTINE FlushOut ( Unit ) ! SUBROUTINE GET_CWD( DirName, Status ) ! SUBROUTINE MKDIR( new_directory_path ) ! SUBROUTINE OpenCon ! SUBROUTINE OpenUnfInpBEFile ( Un, InFile, RecLen, Error ) ! SUBROUTINE ProgExit ( StatCode ) - ! SUBROUTINE Set_IEEE_Constants( NaN_D, Inf_D, NaN, Inf, NaN_S, Inf_S ) + ! SUBROUTINE Set_IEEE_Constants( NaN_D, Inf_D, NaN, Inf, NaN_S, Inf_S ) ! SUBROUTINE UsrAlarm ! SUBROUTINE WrNR ( Str ) ! SUBROUTINE WrOver ( Str ) @@ -43,11 +43,11 @@ MODULE SysSubs USE NWTC_Base IMPLICIT NONE - + INTERFACE NWTC_ERF ! Returns the ERF value of its argument MODULE PROCEDURE NWTC_ERFR4 MODULE PROCEDURE NWTC_ERFR8 - END INTERFACE + END INTERFACE INTERFACE NWTC_gamma ! Returns the gamma value of its argument ! note: gamma is part of the F08 standard, but may not be implemented everywhere... @@ -115,20 +115,20 @@ FUNCTION Is_NaN( DblNum ) RETURN END FUNCTION Is_NaN ! ( DblNum ) !======================================================================= -!> This function returns the ERF value of its argument. The result has a value equal -!! to the Gauss error function: +!> This function returns the ERF value of its argument. The result has a value equal +!! to the Gauss error function: !! \f{equation}{ !! \mathrm{erf}(x)=\frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}\,dt !! \f} \n !! Use NWTC_ERF (syssubs::nwtc_erf) instead of directly calling a specific routine in the generic interface. FUNCTION NWTC_ERFR4( x ) - ! Returns the ERF value of its argument. The result has a value equal - ! to the error function: 2/pi * integral_from_0_to_x of e^(-t^2) dt. + ! Returns the ERF value of its argument. The result has a value equal + ! to the error function: 2/pi * integral_from_0_to_x of e^(-t^2) dt. - REAL(SiKi), INTENT(IN) :: x !< input + REAL(SiKi), INTENT(IN) :: x !< input REAL(SiKi) :: NWTC_ERFR4 !< \f$\mathrm{erf}(x)\f$ - + NWTC_ERFR4 = ERF( x ) END FUNCTION NWTC_ERFR4 @@ -136,14 +136,14 @@ END FUNCTION NWTC_ERFR4 !> \copydoc syssubs::nwtc_erfr4 FUNCTION NWTC_ERFR8( x ) - REAL(R8Ki), INTENT(IN) :: x ! input + REAL(R8Ki), INTENT(IN) :: x ! input REAL(R8Ki) :: NWTC_ERFR8 ! this function - + NWTC_ERFR8 = ERF( x ) END FUNCTION NWTC_ERFR8 !======================================================================= -!> Returns the gamma value of its argument. The result has a value equal +!> Returns the gamma value of its argument. The result has a value equal !! to a processor-dependent approximation to the gamma function of x: !! \f{equation}{ !! \mathrm{ \Gamma }(x) = \int_0^\infty t^{x-1}e^{-t}\,dt @@ -153,9 +153,9 @@ FUNCTION NWTC_GammaR4( x ) ! \mathrm{ \Gamma }(x) = \int_0^\infinity t^{x-1}e^{-t}\,dt - REAL(SiKi), INTENT(IN) :: x !< input + REAL(SiKi), INTENT(IN) :: x !< input REAL(SiKi) :: NWTC_GammaR4 !< \f$\mathrm{\Gamma}(x)\f$ - + NWTC_GammaR4 = gamma( x ) END FUNCTION NWTC_GammaR4 @@ -163,9 +163,9 @@ END FUNCTION NWTC_GammaR4 !> \copydoc syssubs::nwtc_gammar4 FUNCTION NWTC_GammaR8( x ) - REAL(R8Ki), INTENT(IN) :: x ! input + REAL(R8Ki), INTENT(IN) :: x ! input REAL(R8Ki) :: NWTC_GammaR8 ! result - + NWTC_GammaR8 = gamma( x ) END FUNCTION NWTC_GammaR8 @@ -230,7 +230,7 @@ SUBROUTINE OpenCon() END SUBROUTINE OpenCon !======================================================================= !> This routine opens a binary input file with data stored in Big Endian format (created on a UNIX machine.) -!! Data are stored in RecLen-byte records. (This routine is used for +!! Data are stored in RecLen-byte records. (This routine is used for SUBROUTINE OpenUnfInpBEFile ( Un, InFile, RecLen, Error ) IMPLICIT NONE @@ -265,7 +265,7 @@ SUBROUTINE ProgExit ( StatCode ) INTEGER, INTENT(IN) :: StatCode ! The status code to pass to the OS. CALL EXIT ( StatCode ) - + ! IF ( StatCode == 0 ) THEN ! STOP 0 ! ELSE @@ -276,12 +276,12 @@ SUBROUTINE ProgExit ( StatCode ) ! END IF END SUBROUTINE ProgExit ! ( StatCode ) !======================================================================= -!> This routine sets the values of NaN_D, Inf_D, NaN, Inf (IEEE -!! values for not-a-number and infinity in sindle and double -!! precision) This uses standard F03 intrinsic routines, +!> This routine sets the values of NaN_D, Inf_D, NaN, Inf (IEEE +!! values for not-a-number and infinity in sindle and double +!! precision) This uses standard F03 intrinsic routines, !! however Gnu has not yet implemented it, so we've placed this !! routine in the system-specific code. -SUBROUTINE Set_IEEE_Constants( NaN_D, Inf_D, NaN, Inf, NaN_S, Inf_S ) +SUBROUTINE Set_IEEE_Constants( NaN_D, Inf_D, NaN, Inf, NaN_S, Inf_S ) USE, INTRINSIC :: ieee_arithmetic ! use this for compilers that have implemented ieee_arithmetic from F03 standard (otherwise see logic in SysGnu*.f90) @@ -290,7 +290,7 @@ SUBROUTINE Set_IEEE_Constants( NaN_D, Inf_D, NaN, Inf, NaN_S, Inf_S ) REAL(ReKi), INTENT(inout) :: Inf !< IEEE value for NaN (not-a-number) REAL(ReKi), INTENT(inout) :: NaN !< IEEE value for Inf (infinity) - + REAL(SiKi), INTENT(inout) :: Inf_S !< IEEE value for NaN (not-a-number) in single precision REAL(SiKi), INTENT(inout) :: NaN_S !< IEEE value for Inf (infinity) in single precision @@ -298,12 +298,12 @@ SUBROUTINE Set_IEEE_Constants( NaN_D, Inf_D, NaN, Inf, NaN_S, Inf_S ) Inf_D = ieee_value(0.0_DbKi, ieee_positive_inf) NaN = ieee_value(0.0_ReKi, ieee_quiet_nan) - Inf = ieee_value(0.0_ReKi, ieee_positive_inf) + Inf = ieee_value(0.0_ReKi, ieee_positive_inf) NaN_S = ieee_value(0.0_SiKi, ieee_quiet_nan) Inf_S = ieee_value(0.0_SiKi, ieee_positive_inf) -END SUBROUTINE Set_IEEE_Constants +END SUBROUTINE Set_IEEE_Constants !======================================================================= !> This routine generates an alarm to warn the user that something went wrong. SUBROUTINE UsrAlarm() @@ -327,11 +327,11 @@ END SUBROUTINE WrNR ! ( Str ) SUBROUTINE WrOver ( Str ) CHARACTER(*), INTENT(IN) :: Str !< The string to write to the screen. - + INTEGER :: MaxLen !< maximum length of string to be written to the screen (cannot exceed ConRecL here) ! When the file is opened using CARRIAGECONTROL='FORTRAN', the "+" character allows writing over the previous line. However, the Fortran carriage control has been deleted from the Fortran standard. - + MaxLen = min(ConRecL-1, len(Str)) if (MaxLen > 0) then WRITE (CU,'("+",A)') Str(1:MaxLen) @@ -404,11 +404,11 @@ SUBROUTINE LoadDynamicLibProc ( DLL, ErrStat, ErrMsg ) ErrStat = ErrID_None ErrMsg = '' - + ! Get the procedure addresses: do i=1,NWTC_MAX_DLL_PROC if ( len_trim( DLL%ProcName(i) ) > 0 ) then - + ProcAddr = GetProcAddress( DLL%FileAddr, TRIM(DLL%ProcName(i))//C_NULL_CHAR ) !the "C_NULL_CHAR" converts the Fortran string to a C-type string (i.e., adds //CHAR(0) to the end) DLL%ProcAddr(i) = TRANSFER(ProcAddr, DLL%ProcAddr(i)) !convert INTEGER(LPVOID) to INTEGER(C_FUNPTR) [used only for compatibility with gfortran] @@ -417,7 +417,7 @@ SUBROUTINE LoadDynamicLibProc ( DLL, ErrStat, ErrMsg ) ErrMsg = 'The procedure '//TRIM(DLL%ProcName(i))//' in file '//TRIM(DLL%FileName)//' could not be loaded.' RETURN END IF - + end if end do @@ -435,12 +435,12 @@ SUBROUTINE FreeDynamicLib ( DLL, ErrStat, ErrMsg ) CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None INTEGER(HANDLE) :: FileAddr ! The address of file FileName. (RETURN value from LoadLibrary in kernel32.f90) INTEGER(BOOL) :: Success ! Whether or not the call to FreeLibrary was successful - + ErrStat = ErrID_None ErrMsg = '' - + IF ( DLL%FileAddr == INT(0,C_INTPTR_T) ) RETURN - + FileAddr = TRANSFER(DLL%FileAddr, FileAddr) !convert INTEGER(C_INTPTR_T) to INTEGER(HANDLE) [used only for compatibility with gfortran] ! Free the DLL: From 1a2c30fd1ccf49fdfca4f12f0e817a211da82636 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Tue, 4 Aug 2026 10:26:54 -0600 Subject: [PATCH 18/35] Refactor console output length limits --- modules/nwtc-library/src/SysFlangLinux.f90 | 4 ++-- modules/nwtc-library/src/SysGnuLinux.f90 | 4 ++-- modules/nwtc-library/src/SysIVF.f90 | 4 ++-- modules/nwtc-library/src/SysIVF_Labview.f90 | 4 ++-- modules/nwtc-library/src/SysMatlabLinuxGnu.f90 | 4 ++-- modules/nwtc-library/src/SysMatlabLinuxIntel.f90 | 4 ++-- 6 files changed, 12 insertions(+), 12 deletions(-) diff --git a/modules/nwtc-library/src/SysFlangLinux.f90 b/modules/nwtc-library/src/SysFlangLinux.f90 index 8805951090..ee97ce4702 100644 --- a/modules/nwtc-library/src/SysFlangLinux.f90 +++ b/modules/nwtc-library/src/SysFlangLinux.f90 @@ -54,9 +54,9 @@ MODULE SysSubs MODULE PROCEDURE NWTC_gammaR8 END INTERFACE - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 6 ! The I/O unit for the console. Unit 6 causes ADAMS to crash. - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .TRUE. ! A flag to tell the program that keyboard input is allowed in the environment. CHARACTER(*), PARAMETER :: NewLine = ACHAR(10) ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}] CHARACTER(*), PARAMETER :: OS_Desc = 'GNU Fortran for Linux' ! Description of the language/OS diff --git a/modules/nwtc-library/src/SysGnuLinux.f90 b/modules/nwtc-library/src/SysGnuLinux.f90 index 6eddb9b410..44d7fe73ea 100644 --- a/modules/nwtc-library/src/SysGnuLinux.f90 +++ b/modules/nwtc-library/src/SysGnuLinux.f90 @@ -55,9 +55,9 @@ MODULE SysSubs MODULE PROCEDURE NWTC_gammaR8 END INTERFACE - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 6 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .TRUE. ! A flag to tell the program that keyboard input is allowed in the environment. CHARACTER(*), PARAMETER :: NewLine = ACHAR(10) ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}] CHARACTER(*), PARAMETER :: OS_Desc = 'GNU Fortran for Linux' ! Description of the language/OS diff --git a/modules/nwtc-library/src/SysIVF.f90 b/modules/nwtc-library/src/SysIVF.f90 index 4ecc6d0c04..71bfcf7b35 100644 --- a/modules/nwtc-library/src/SysIVF.f90 +++ b/modules/nwtc-library/src/SysIVF.f90 @@ -55,9 +55,9 @@ MODULE SysSubs MODULE PROCEDURE NWTC_gammaR8 END INTERFACE - INTEGER, PARAMETER :: ConRecL = 300 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 7 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .TRUE. ! A flag to tell the program that keyboard input is allowed in the environment. CHARACTER(*), PARAMETER :: NewLine = ACHAR(10) ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}]; Note: NewLine change to ACHAR(10) here on Windows to fix issues with C/Fortran interoperability using WrScr CHARACTER(*), PARAMETER :: OS_Desc = 'Intel Visual Fortran for Windows'! Description of the language/OS diff --git a/modules/nwtc-library/src/SysIVF_Labview.f90 b/modules/nwtc-library/src/SysIVF_Labview.f90 index 51a57144f5..4b2d433723 100644 --- a/modules/nwtc-library/src/SysIVF_Labview.f90 +++ b/modules/nwtc-library/src/SysIVF_Labview.f90 @@ -72,9 +72,9 @@ MODULE SysSubs !======================================================================= - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 7 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .FALSE. ! A flag to tell the program that keyboard input is allowed in the environment. diff --git a/modules/nwtc-library/src/SysMatlabLinuxGnu.f90 b/modules/nwtc-library/src/SysMatlabLinuxGnu.f90 index 1072ca8e3b..e80a4ae599 100644 --- a/modules/nwtc-library/src/SysMatlabLinuxGnu.f90 +++ b/modules/nwtc-library/src/SysMatlabLinuxGnu.f90 @@ -58,9 +58,9 @@ MODULE SysSubs MODULE PROCEDURE NWTC_gammaR8 END INTERFACE - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 6 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .FALSE. ! A flag to tell the program that keyboard input is allowed in the environment. CHARACTER(*), PARAMETER :: NewLine = ACHAR(10) ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}] CHARACTER(*), PARAMETER :: OS_Desc = 'GNU Fortran for Linux with Matlab' ! Description of the language/OS diff --git a/modules/nwtc-library/src/SysMatlabLinuxIntel.f90 b/modules/nwtc-library/src/SysMatlabLinuxIntel.f90 index 303206a9a4..2d55265fcc 100644 --- a/modules/nwtc-library/src/SysMatlabLinuxIntel.f90 +++ b/modules/nwtc-library/src/SysMatlabLinuxIntel.f90 @@ -58,9 +58,9 @@ MODULE SysSubs MODULE PROCEDURE NWTC_gammaR8 END INTERFACE - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 6 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .FALSE. ! A flag to tell the program that keyboard input is allowed in the environment. CHARACTER(*), PARAMETER :: NewLine = ACHAR(10) ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}] CHARACTER(*), PARAMETER :: OS_Desc = 'Intel Fortran for Linux with Matlab' ! Description of the language/OS From 0718c221f68d3627f8861cdd47e50bc5f757bea7 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Fri, 31 Jul 2026 11:53:55 -0600 Subject: [PATCH 19/35] Updated the .gitignore for vscode containers / workspaces --- .gitignore | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.gitignore b/.gitignore index 09166cd156..74be88339b 100644 --- a/.gitignore +++ b/.gitignore @@ -64,3 +64,8 @@ varcache # Python cache files openfast_io/dist/ openfast_io/openfast_io/_version.py + +# docker with VScode config +*.code-workspace +docker-compose.yml +.devcontainer/ From 1db891bc3c25757513329072895a8eb91513562f Mon Sep 17 00:00:00 2001 From: Andy Platt Date: Wed, 5 Aug 2026 14:06:38 -0600 Subject: [PATCH 20/35] typo in description Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- modules/nwtc-library/src/SysIVF.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/nwtc-library/src/SysIVF.f90 b/modules/nwtc-library/src/SysIVF.f90 index 71bfcf7b35..0286de427d 100644 --- a/modules/nwtc-library/src/SysIVF.f90 +++ b/modules/nwtc-library/src/SysIVF.f90 @@ -277,7 +277,7 @@ SUBROUTINE ProgExit ( StatCode ) END SUBROUTINE ProgExit ! ( StatCode ) !======================================================================= !> This routine sets the values of NaN_D, Inf_D, NaN, Inf (IEEE -!! values for not-a-number and infinity in sindle and double +!! values for not-a-number and infinity in single and double !! precision) This uses standard F03 intrinsic routines, !! however Gnu has not yet implemented it, so we've placed this !! routine in the system-specific code. From 588a0233ba2c13d7bbd74e6840a73beb55eab3e1 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Wed, 5 Aug 2026 14:17:26 -0600 Subject: [PATCH 21/35] Missing files for console output length limits --- modules/nwtc-library/src/SysGnuWin.f90 | 4 ++-- modules/nwtc-library/src/SysIFL.f90 | 4 ++-- modules/nwtc-library/src/SysMatlabWindows.f90 | 3 ++- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/modules/nwtc-library/src/SysGnuWin.f90 b/modules/nwtc-library/src/SysGnuWin.f90 index 1422989cb2..cad5fda7d4 100644 --- a/modules/nwtc-library/src/SysGnuWin.f90 +++ b/modules/nwtc-library/src/SysGnuWin.f90 @@ -55,9 +55,9 @@ MODULE SysSubs MODULE PROCEDURE NWTC_gammaR8 END INTERFACE - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 6 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .TRUE. ! A flag to tell the program that keyboard input is allowed in the environment. CHARACTER(*), PARAMETER :: NewLine = ACHAR(10) ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}] CHARACTER(*), PARAMETER :: OS_Desc = 'GNU Fortran for Windows' ! Description of the language/OS diff --git a/modules/nwtc-library/src/SysIFL.f90 b/modules/nwtc-library/src/SysIFL.f90 index bf427e8fb6..17e7db0e28 100644 --- a/modules/nwtc-library/src/SysIFL.f90 +++ b/modules/nwtc-library/src/SysIFL.f90 @@ -55,9 +55,9 @@ MODULE SysSubs MODULE PROCEDURE NWTC_gammaR8 END INTERFACE - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 6 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .TRUE. ! A flag to tell the program that keyboard input is allowed in the environment. CHARACTER(*), PARAMETER :: NewLine = ACHAR(10) ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}] CHARACTER(*), PARAMETER :: OS_Desc = 'Intel Fortran for Linux' ! Description of the language/OS diff --git a/modules/nwtc-library/src/SysMatlabWindows.f90 b/modules/nwtc-library/src/SysMatlabWindows.f90 index 417a0e9ef6..c082043b8d 100644 --- a/modules/nwtc-library/src/SysMatlabWindows.f90 +++ b/modules/nwtc-library/src/SysMatlabWindows.f90 @@ -47,8 +47,9 @@ MODULE SysSubs !======================================================================= + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 6 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .FALSE. ! A flag to tell the program that keyboard input is allowed in the environment. From 00d469986446a465025e24033bee9b1f25307d1c Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Thu, 6 Aug 2026 17:11:21 -0600 Subject: [PATCH 22/35] Fix #3421: Guard DT_mooring validation with mooring enabled check - Wrap DT_mooring validation in 'if (p%MooringMod == 3)' conditional - Add file existence check for MD_FileName when mooring is enabled - Prevents fatal error when Mod_SharedMooring=0 and mooring is not used - Follows established OpenFAST pattern for module-specific validations Co-authored-by: GitHub Copilot Co-authored-by: Claude Sonnet 4.5 --- glue-codes/fast-farm/src/FAST_Farm_IO.f90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/glue-codes/fast-farm/src/FAST_Farm_IO.f90 b/glue-codes/fast-farm/src/FAST_Farm_IO.f90 index 32fce6d53c..d954a4ea39 100644 --- a/glue-codes/fast-farm/src/FAST_Farm_IO.f90 +++ b/glue-codes/fast-farm/src/FAST_Farm_IO.f90 @@ -1047,6 +1047,7 @@ SUBROUTINE Farm_ValidateInput( p, WD_InitInp, AWAE_InitInp, ErrStat, ErrMsg ) CHARACTER(*), PARAMETER :: RoutineName = 'Farm_ValidateInput' INTEGER(IntKi) :: n_disDT_dt character(60) :: tmpStr + logical :: file_exists ! Flag for file existence check ErrStat = ErrID_None ErrMsg = "" @@ -1066,8 +1067,11 @@ SUBROUTINE Farm_ValidateInput( p, WD_InitInp, AWAE_InitInp, ErrStat, ErrMsg ) ! TODO : Verify that the DLL file exists ! --- SHARED MOORING SYSTEM --- - ! TODO : Verify that p%MD_FileName file exists - if ((p%DT_mooring <= 0.0_ReKi) .or. (p%DT_mooring > p%DT_high)) CALL SetErrStat(ErrID_Fatal,'DT_mooring must be greater than zero and no greater than dt_high.',ErrStat,ErrMsg,RoutineName) + if (p%MooringMod == 3) then + inquire(file=trim(p%MD_FileName), exist=file_exists) + if (.not. file_exists) call SetErrStat(ErrID_Fatal,'Cannot find MoorDyn input file '//trim(p%MD_FileName),ErrStat,ErrMsg,RoutineName) + if ((p%DT_mooring <= 0.0_ReKi) .or. (p%DT_mooring > p%DT_high)) CALL SetErrStat(ErrID_Fatal,'DT_mooring must be greater than zero and no greater than dt_high.',ErrStat,ErrMsg,RoutineName) + end if ! --- AMBIENT WIND: INFLOWWIND MODULE --- [used only for Mod_AmbWind=2 or 3] --- ! FIXME: this really should be checked with the turbine specific size diameter -- maybe relocate this check to AWAE or in FF after initializing all turbines? From 92c1ac34ccf62190cde6f4a28a6ea85b1006860f Mon Sep 17 00:00:00 2001 From: ptrbortolotti Date: Mon, 10 Aug 2026 16:23:15 -0600 Subject: [PATCH 23/35] fix bug in bldprcnt following issue #3416 --- modules/aerodyn/src/AeroAcoustics.f90 | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/modules/aerodyn/src/AeroAcoustics.f90 b/modules/aerodyn/src/AeroAcoustics.f90 index 9142a8eea4..3727f0e33c 100644 --- a/modules/aerodyn/src/AeroAcoustics.f90 +++ b/modules/aerodyn/src/AeroAcoustics.f90 @@ -271,7 +271,15 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) p%BlSpn = InitInp%BlSpn p%BlChord = InitInp%BlChord - p%startnode = max(1, p%NumBlNds - 1) + IF (InputFileData%AA_Bl_Prcntge .lt. (100.*(p%BlSpn(p%NumBlNds,1) - p%BlSpn(p%NumBlNds-1,1))/p%BlSpn(p%NumBlNds,1)) THEN + CALL SetErrStat(ErrID_Warn, 'AA_Bl_Prcntge is smaller than the last blade element size. '// & + 'OpenFAST will move on assuming the last blade element, which is the minimum blade span used for noise calculations. '// & + 'Either increase BldPrcnt in your aeroacoustic input file '// & + 'or refine the aerodynamic mesh in your blade AeroDyn input file.', ErrStat2, ErrMsg2, RoutineName ) + if(Failed()) return + endif + + p%startnode = min(p%NumBlNds, 2) BladeSpanUsedForNoise = p%BlSpn(p%NumBlNds,1)*(1.0 - InputFileData%AA_Bl_Prcntge/100.0) do j=p%NumBlNds-1,2,-1 IF ( p%BlSpn(j,1) .lt. BladeSpanUsedForNoise )THEN @@ -279,7 +287,7 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) exit ! exit the loop endif enddo - p%startnode = max(min(p%NumBlNds,2),p%startnode) + p%BlElemSpn = 0; DO I = 1,p%numBlades From 970d5c54f0189f2bb6ad4a04114416444255e0eb Mon Sep 17 00:00:00 2001 From: ptrbortolotti Date: Mon, 10 Aug 2026 16:31:27 -0600 Subject: [PATCH 24/35] missed a parenthesis --- modules/aerodyn/src/AeroAcoustics.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/aerodyn/src/AeroAcoustics.f90 b/modules/aerodyn/src/AeroAcoustics.f90 index 3727f0e33c..f7a5eb112e 100644 --- a/modules/aerodyn/src/AeroAcoustics.f90 +++ b/modules/aerodyn/src/AeroAcoustics.f90 @@ -271,7 +271,7 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) p%BlSpn = InitInp%BlSpn p%BlChord = InitInp%BlChord - IF (InputFileData%AA_Bl_Prcntge .lt. (100.*(p%BlSpn(p%NumBlNds,1) - p%BlSpn(p%NumBlNds-1,1))/p%BlSpn(p%NumBlNds,1)) THEN + IF (InputFileData%AA_Bl_Prcntge .lt. (100.*(p%BlSpn(p%NumBlNds,1) - p%BlSpn(p%NumBlNds-1,1))/p%BlSpn(p%NumBlNds,1))) THEN CALL SetErrStat(ErrID_Warn, 'AA_Bl_Prcntge is smaller than the last blade element size. '// & 'OpenFAST will move on assuming the last blade element, which is the minimum blade span used for noise calculations. '// & 'Either increase BldPrcnt in your aeroacoustic input file '// & From 37dcc8a059eabbb37846c7660c22b879f707155c Mon Sep 17 00:00:00 2001 From: Pietro Bortolotti Date: Tue, 11 Aug 2026 08:42:47 -0600 Subject: [PATCH 25/35] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- modules/aerodyn/src/AeroAcoustics.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/aerodyn/src/AeroAcoustics.f90 b/modules/aerodyn/src/AeroAcoustics.f90 index f7a5eb112e..c22f4c3d2b 100644 --- a/modules/aerodyn/src/AeroAcoustics.f90 +++ b/modules/aerodyn/src/AeroAcoustics.f90 @@ -272,7 +272,7 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) p%BlChord = InitInp%BlChord IF (InputFileData%AA_Bl_Prcntge .lt. (100.*(p%BlSpn(p%NumBlNds,1) - p%BlSpn(p%NumBlNds-1,1))/p%BlSpn(p%NumBlNds,1))) THEN - CALL SetErrStat(ErrID_Warn, 'AA_Bl_Prcntge is smaller than the last blade element size. '// & + CALL SetErrStat(ErrID_Warn, 'BldPrcnt is smaller than the last blade element size. '// & 'OpenFAST will move on assuming the last blade element, which is the minimum blade span used for noise calculations. '// & 'Either increase BldPrcnt in your aeroacoustic input file '// & 'or refine the aerodynamic mesh in your blade AeroDyn input file.', ErrStat2, ErrMsg2, RoutineName ) From add13f51b693584c222b137bafbcc4b1b6a88297 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Wed, 12 Aug 2026 14:02:07 -0600 Subject: [PATCH 26/35] Update r-test pointer after merge there --- reg_tests/r-test | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reg_tests/r-test b/reg_tests/r-test index 65babdd4da..d376069e24 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 65babdd4da8419ab9cf25c5662a9e1ee45c22049 +Subproject commit d376069e24765fc3e6456dbf16ae39e6fd8e78ba From 2942295fa78f79f2ee74430b94c5b6a8684023ee Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Wed, 12 Aug 2026 14:31:30 -0600 Subject: [PATCH 27/35] Fix InflowWind Grid3D cubic velocity interpolation NaN bug IfW_FlowField_GetVelAcc only populated the local AccCell array when the caller requested acceleration output (OutputAccel). But the cubic Hermite velocity formula in Grid3DField_GetVelAccCubic depends on AccCell for its tangent terms regardless of whether acceleration is separately requested. Callers that use VelInterpCubic=True but only want velocity (e.g. AWAE's ambient-wind sampling) therefore computed velocity from uninitialized memory, corrupting every returned value. Also fixes IfW_Grid3DField_CalcAccel, which checked G3D%NTGrids (tower grid point count) instead of G3D%NSteps (time step count) to decide whether to compute real cubic-spline time derivatives, forcing G3D%Acc to zero for every case without a tower file regardless of how many time steps were available. Adds unit tests (test_grid3d_field.F90) that reproduce both defects in isolation and verify the fix. Co-authored-by: Claude Sonnet 5 --- modules/inflowwind/src/IfW_FlowField.f90 | 10 +- modules/inflowwind/tests/inflowwind_utest.F90 | 2 + .../inflowwind/tests/test_grid3d_field.F90 | 163 ++++++++++++++++++ unit_tests/CMakeLists.txt | 1 + 4 files changed, 172 insertions(+), 4 deletions(-) create mode 100644 modules/inflowwind/tests/test_grid3d_field.F90 diff --git a/modules/inflowwind/src/IfW_FlowField.f90 b/modules/inflowwind/src/IfW_FlowField.f90 index 922b831843..7c058c0e4d 100644 --- a/modules/inflowwind/src/IfW_FlowField.f90 +++ b/modules/inflowwind/src/IfW_FlowField.f90 @@ -224,8 +224,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,8 +1444,8 @@ 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 return end if 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..0a36d383af --- /dev/null +++ b/modules/inflowwind/tests/test_grid3d_field.F90 @@ -0,0 +1,163 @@ +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) & + ] +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 + +end module 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 From c49fdbfe59be43ece0506bf725cd2456cb370650 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Wed, 12 Aug 2026 15:26:29 -0600 Subject: [PATCH 28/35] Update ifw_boxexceed results for acceleration --- reg_tests/r-test | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reg_tests/r-test b/reg_tests/r-test index 65babdd4da..5d3c808893 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 65babdd4da8419ab9cf25c5662a9e1ee45c22049 +Subproject commit 5d3c808893544cefeff424790089e787aaa1fdae From 4567015194b4180d3448d51b24dfa0e87400d1d1 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Wed, 12 Aug 2026 14:31:30 -0600 Subject: [PATCH 29/35] Fix InflowWind Grid3D cubic velocity interpolation NaN bug IfW_FlowField_GetVelAcc only populated the local AccCell array when the caller requested acceleration output (OutputAccel). But the cubic Hermite velocity formula in Grid3DField_GetVelAccCubic depends on AccCell for its tangent terms regardless of whether acceleration is separately requested. Callers that use VelInterpCubic=True but only want velocity (e.g. AWAE's ambient-wind sampling) therefore computed velocity from uninitialized memory, corrupting every returned value. Also fixes IfW_Grid3DField_CalcAccel, which checked G3D%NTGrids (tower grid point count) instead of G3D%NSteps (time step count) to decide whether to compute real cubic-spline time derivatives, forcing G3D%Acc to zero for every case without a tower file regardless of how many time steps were available. Adds unit tests (test_grid3d_field.F90) that reproduce both defects in isolation and verify the fix. Co-authored-by: Claude Sonnet 5 --- modules/inflowwind/src/IfW_FlowField.f90 | 23 ++- modules/inflowwind/tests/inflowwind_utest.F90 | 2 + .../inflowwind/tests/test_grid3d_field.F90 | 163 ++++++++++++++++++ unit_tests/CMakeLists.txt | 1 + 4 files changed, 184 insertions(+), 5 deletions(-) create mode 100644 modules/inflowwind/tests/test_grid3d_field.F90 diff --git a/modules/inflowwind/src/IfW_FlowField.f90 b/modules/inflowwind/src/IfW_FlowField.f90 index 922b831843..332879e6f1 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 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..0a36d383af --- /dev/null +++ b/modules/inflowwind/tests/test_grid3d_field.F90 @@ -0,0 +1,163 @@ +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) & + ] +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 + +end module 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 From 56bd27bd623866b384db617dc3fd9dfd1ed5e76b Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Thu, 13 Aug 2026 07:05:04 -0600 Subject: [PATCH 30/35] Fix Grid3D tower acceleration guard using wrong variable (NTGrids vs NSteps) IfW_Grid3DField_CalcAccel's tower-grid branch checked G3D%NTGrids < 3 to decide whether to zero out tower acceleration, but the actual time-derivative calculation needs G3D%NSteps >= 3 (already guaranteed by the earlier NSteps<3 early return). This mirrors the NTGrids/NSteps mix-up fixed elsewhere in this PR, and incorrectly zeroed valid tower acceleration whenever a wind file had only 1 or 2 tower grid points. Removed the erroneous guard; the tower acceleration loop needs no minimum tower height count. Added test_grid3d_calcaccel_few_tower_points to cover this case. Co-authored-by: GitHub Copilot <175728472+Copilot@users.noreply.github.com> Co-authored-by: Claude Sonnet 5 --- modules/inflowwind/src/IfW_FlowField.f90 | 16 +++--- .../inflowwind/tests/test_grid3d_field.F90 | 57 ++++++++++++++++++- 2 files changed, 63 insertions(+), 10 deletions(-) diff --git a/modules/inflowwind/src/IfW_FlowField.f90 b/modules/inflowwind/src/IfW_FlowField.f90 index 332879e6f1..a8212790e1 100644 --- a/modules/inflowwind/src/IfW_FlowField.f90 +++ b/modules/inflowwind/src/IfW_FlowField.f90 @@ -1490,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/test_grid3d_field.F90 b/modules/inflowwind/tests/test_grid3d_field.F90 index 0a36d383af..16a0441e43 100644 --- a/modules/inflowwind/tests/test_grid3d_field.F90 +++ b/modules/inflowwind/tests/test_grid3d_field.F90 @@ -21,7 +21,8 @@ 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_no_tower", test_grid3d_calcaccel_no_tower), & + new_unittest("test_grid3d_calcaccel_few_tower_points", test_grid3d_calcaccel_few_tower_points) & ] end subroutine @@ -160,4 +161,58 @@ subroutine test_grid3d_calcaccel_no_tower(error) 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 From 7f2d1ea942bd5bbf6098326a3a1e22fa7ae26182 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Thu, 13 Aug 2026 08:49:14 -0600 Subject: [PATCH 31/35] Update several MHK regression tests that had incorrect accelerations calculated before --- reg_tests/r-test | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reg_tests/r-test b/reg_tests/r-test index 04b5f4faf5..f8d1d131de 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 04b5f4faf5efa7e284f7db27d4839974e6f13627 +Subproject commit f8d1d131debad662dffa749bf37b440f8320095f From 53f11d85d750901776af320d5577d53f77813d9c Mon Sep 17 00:00:00 2001 From: jcvem Date: Mon, 17 Aug 2026 12:55:36 +0800 Subject: [PATCH 32/35] Found and fixed a type cast missing in FAST_writer.py At L:620 of FAST_writer.py contains a str + int type error. Patched with the missing str() casting function --- openfast_io/openfast_io/FAST_writer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openfast_io/openfast_io/FAST_writer.py b/openfast_io/openfast_io/FAST_writer.py index e9998d6fbf..2db63a18ee 100644 --- a/openfast_io/openfast_io/FAST_writer.py +++ b/openfast_io/openfast_io/FAST_writer.py @@ -617,7 +617,7 @@ def write_ElastoDynBlade(self, bldInd = None): blade_file = os.path.join(self.FAST_runDirectory,self.fst_vt['ElastoDyn']['BldFile1']) else: EDbld_dict = self.fst_vt['ElastoDynBlade'][bldInd] - blade_file = os.path.join(self.FAST_runDirectory,self.fst_vt['ElastoDyn']['BldFile'+(bldInd+1)]) + blade_file = os.path.join(self.FAST_runDirectory,self.fst_vt['ElastoDyn']['BldFile'+ str(bldInd+1)]) f = open(blade_file, 'w') From 36e61151bdcd475b64ce5a03ca6859b62dc625ee Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Mon, 17 Aug 2026 14:50:43 -0600 Subject: [PATCH 33/35] Fix GitHub Copilot identified issues in rc-5.0.1 cross-merge Fix three code quality issues identified by GitHub Copilot in PR #3432: 1. AeroAcoustics.f90: Prevent out-of-bounds array access - Add guard for single-node case (NumBlNds==1) - Calculate LastElemPct safely without accessing NumBlNds-1 index - When NumBlNds==1, treat element as spanning entire blade (100%) 2. IfW_FlowField.f90: Improve error messaging for cubic interpolation - Clarify why acceleration field is required for cubic velocity interpolation - Provide context-specific error messages for three scenarios: * Both acceleration output and cubic interpolation enabled * Only acceleration output requested * Only cubic velocity interpolation enabled 3. FVW_IO.f90: Add I/O error handling in ResolveGridAxis - Add iostat check to grid point list file read loop - Return controlled error with line number and filename on read failure - Prevent unhandled runtime crashes from malformed/truncated files Co-authored-by: GitHub Copilot Co-authored-by: Claude Sonnet 4.5 --- modules/inflowwind/src/IfW_FlowField.f90 | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/modules/inflowwind/src/IfW_FlowField.f90 b/modules/inflowwind/src/IfW_FlowField.f90 index a8212790e1..1985647fb5 100644 --- a/modules/inflowwind/src/IfW_FlowField.f90 +++ b/modules/inflowwind/src/IfW_FlowField.f90 @@ -78,11 +78,19 @@ subroutine IfW_FlowField_GetVelAcc(FF, IStart, Time, PositionXYZ, VelocityUVW, A ! Determine if acceleration should be calculated and returned OutputAccel = allocated(AccelUVW) - ! Cubic velocity interpolation also requires a valid acceleration field, since its - ! formula uses the derivative data even when acceleration output is not requested. + ! 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) + if (OutputAccel .and. FF%VelInterpCubic) then + call SetErrStat(ErrID_Fatal, "Acceleration output and cubic velocity interpolation both require a valid acceleration field, but accel field is not valid", & + ErrStat, ErrMsg, RoutineName) + else if (OutputAccel) then + call SetErrStat(ErrID_Fatal, "Acceleration output requested, but accel field is not valid", & + ErrStat, ErrMsg, RoutineName) + else ! FF%VelInterpCubic + call SetErrStat(ErrID_Fatal, "Cubic velocity interpolation requires a valid acceleration field, but accel field is not valid", & + ErrStat, ErrMsg, RoutineName) + end if return end if From fedc284f55186920a41ef5d393b287e600be4ad7 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Mon, 17 Aug 2026 14:38:12 -0600 Subject: [PATCH 34/35] Fix GitHub Copilot identified issues in rc-5.0.1 cross-merge Fix three code quality issues identified by GitHub Copilot in PR #3432: 1. AeroAcoustics.f90: Prevent out-of-bounds array access - Add guard for single-node case (NumBlNds==1) - Calculate LastElemPct safely without accessing NumBlNds-1 index - When NumBlNds==1, treat element as spanning entire blade (100%) 2. IfW_FlowField.f90: Improve error messaging for cubic interpolation - Clarify why acceleration field is required for cubic velocity interpolation - Provide context-specific error messages for three scenarios: * Both acceleration output and cubic interpolation enabled * Only acceleration output requested * Only cubic velocity interpolation enabled 3. FVW_IO.f90: Add I/O error handling in ResolveGridAxis - Add iostat check to grid point list file read loop - Return controlled error with line number and filename on read failure - Prevent unhandled runtime crashes from malformed/truncated files Co-authored-by: GitHub Copilot Co-authored-by: Claude Sonnet 4.5 --- modules/aerodyn/src/AeroAcoustics.f90 | 10 +++++++++- modules/aerodyn/src/FVW_IO.f90 | 8 +++++++- 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/modules/aerodyn/src/AeroAcoustics.f90 b/modules/aerodyn/src/AeroAcoustics.f90 index c22f4c3d2b..708ab8fe41 100644 --- a/modules/aerodyn/src/AeroAcoustics.f90 +++ b/modules/aerodyn/src/AeroAcoustics.f90 @@ -165,6 +165,7 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) character(*), parameter :: RoutineName = 'SetParameters' REAL(ReKi) :: val1,val10,f2,f4, dist1, dist10 REAL(ReKi) :: BladeSpanUsedForNoise + REAL(ReKi) :: LastElemPct ! Initialize variables for this routine ErrStat = ErrID_None @@ -271,7 +272,14 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) p%BlSpn = InitInp%BlSpn p%BlChord = InitInp%BlChord - IF (InputFileData%AA_Bl_Prcntge .lt. (100.*(p%BlSpn(p%NumBlNds,1) - p%BlSpn(p%NumBlNds-1,1))/p%BlSpn(p%NumBlNds,1))) THEN + ! Calculate last element size as percentage of blade span + IF (p%NumBlNds > 1) THEN + LastElemPct = 100.0 * (p%BlSpn(p%NumBlNds,1) - p%BlSpn(p%NumBlNds-1,1)) / p%BlSpn(p%NumBlNds,1) + ELSE + LastElemPct = 100.0 ! Single node means element spans entire blade + ENDIF + + IF (InputFileData%AA_Bl_Prcntge .lt. LastElemPct) THEN CALL SetErrStat(ErrID_Warn, 'BldPrcnt is smaller than the last blade element size. '// & 'OpenFAST will move on assuming the last blade element, which is the minimum blade span used for noise calculations. '// & 'Either increase BldPrcnt in your aeroacoustic input file '// & diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index a7c252bf04..6e56e84547 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -458,7 +458,13 @@ subroutine ResolveGridAxis(AStart, AEnd, n, ListFile, Pts, ErrStat, ErrMsg) return endif do j = 1, n - read(UnList, *) Pts(j) + read(UnList, *, iostat=IOS) Pts(j) + if (IOS /= 0) then + call SetErrStat(ErrID_Fatal, 'ResolveGridAxis: error reading value at line '//trim(Num2LStr(j))//' from grid point list file "'//trim(FullFile)//'".', & + ErrStat, ErrMsg, 'ResolveGridAxis') + close(UnList) + return + end if enddo close(UnList) From fb960d80748a5e4b9ad7feb2212addbec259ef78 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Mon, 17 Aug 2026 14:48:47 -0600 Subject: [PATCH 35/35] Apply GitHub Copilot PR review suggestions Address two wording improvements suggested by Copilot in PR #3433: 1. IfW_FlowField.f90: Use consistent terminology - Changed 'accel field' to 'the acceleration field' in all three error messages for consistency within the same conditional block 2. FVW_IO.f90: Improve read error diagnostics - Changed 'line j' to 'grid point #j' (j is entry index, not line) - Added iostat code to error message for better troubleshooting - List-directed READ can consume multiple values per line, so reporting it as a line number was misleading Co-authored-by: GitHub Copilot --- modules/aerodyn/src/FVW_IO.f90 | 2 +- modules/inflowwind/src/IfW_FlowField.f90 | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 6e56e84547..a037c85dc3 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -460,7 +460,7 @@ subroutine ResolveGridAxis(AStart, AEnd, n, ListFile, Pts, ErrStat, ErrMsg) do j = 1, n read(UnList, *, iostat=IOS) Pts(j) if (IOS /= 0) then - call SetErrStat(ErrID_Fatal, 'ResolveGridAxis: error reading value at line '//trim(Num2LStr(j))//' from grid point list file "'//trim(FullFile)//'".', & + call SetErrStat(ErrID_Fatal, 'ResolveGridAxis: error reading grid point #'//trim(Num2LStr(j))//' from grid point list file "'//trim(FullFile)//'" (iostat='//trim(Num2LStr(IOS))//').', & ErrStat, ErrMsg, 'ResolveGridAxis') close(UnList) return diff --git a/modules/inflowwind/src/IfW_FlowField.f90 b/modules/inflowwind/src/IfW_FlowField.f90 index 1985647fb5..43f21947e0 100644 --- a/modules/inflowwind/src/IfW_FlowField.f90 +++ b/modules/inflowwind/src/IfW_FlowField.f90 @@ -82,13 +82,13 @@ subroutine IfW_FlowField_GetVelAcc(FF, IStart, Time, PositionXYZ, VelocityUVW, A ! since its formula uses the derivative data even when acceleration output is not requested. if ((OutputAccel .or. FF%VelInterpCubic) .and. .not. FF%AccFieldValid) then if (OutputAccel .and. FF%VelInterpCubic) then - call SetErrStat(ErrID_Fatal, "Acceleration output and cubic velocity interpolation both require a valid acceleration field, but accel field is not valid", & + call SetErrStat(ErrID_Fatal, "Acceleration output and cubic velocity interpolation both require a valid acceleration field, but the acceleration field is not valid", & ErrStat, ErrMsg, RoutineName) else if (OutputAccel) then - call SetErrStat(ErrID_Fatal, "Acceleration output requested, but accel field is not valid", & + call SetErrStat(ErrID_Fatal, "Acceleration output requested, but the acceleration field is not valid", & ErrStat, ErrMsg, RoutineName) else ! FF%VelInterpCubic - call SetErrStat(ErrID_Fatal, "Cubic velocity interpolation requires a valid acceleration field, but accel field is not valid", & + call SetErrStat(ErrID_Fatal, "Cubic velocity interpolation requires a valid acceleration field, but the acceleration field is not valid", & ErrStat, ErrMsg, RoutineName) end if return