From d0e500626b3dcd6f6c8001ed46e26404e315c3c2 Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Wed, 26 Aug 2026 17:54:38 -0600 Subject: [PATCH 01/12] NWTC Library / fftpack 5.1.f: Don't error for interfaces Debug build does not work without this. FFTPACK relies on some old Fortran method of passing the different datatypes through the same interface. Note: this also reorders some files that were added out of alphabetical order. This prevents VS from giving differences every time the file is open (though it doesn't fix the problem with the registry files in the VS project) --- vs-build/.gitignore | 4 ---- vs-build/modules/NWTC-Library.vfproj | 24 ++++++++++++------------ 2 files changed, 12 insertions(+), 16 deletions(-) delete mode 100644 vs-build/.gitignore diff --git a/vs-build/.gitignore b/vs-build/.gitignore deleted file mode 100644 index 5253e137da..0000000000 --- a/vs-build/.gitignore +++ /dev/null @@ -1,4 +0,0 @@ -.vs -*.user -*.u2d -gitVersionInfo.h \ No newline at end of file diff --git a/vs-build/modules/NWTC-Library.vfproj b/vs-build/modules/NWTC-Library.vfproj index 22bb1fce7c..bf5d4164da 100644 --- a/vs-build/modules/NWTC-Library.vfproj +++ b/vs-build/modules/NWTC-Library.vfproj @@ -123,28 +123,28 @@ - + - + - + - - + + - + - - + + - + - + @@ -218,6 +218,8 @@ + + @@ -257,8 +259,6 @@ - - From 09ed542da0afd7d97d7b0c480e8288fbf33bac57 Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Wed, 26 Aug 2026 17:59:43 -0600 Subject: [PATCH 02/12] Add AeroAcoustics and AirfoilInfo drivers to VS solution --- modules/aerodyn/src/AeroAcoustics.f90 | 1 - modules/aerodyn/src/AeroAcoustics_Driver.f90 | 11 - modules/aerodyn/src/AirfoilInfo_Driver.f90 | 210 +++++++++++++++++++ vs-build/OpenFAST.sln | 80 +++++++ vs-build/drivers/AeroAcoustics_Driver.vfproj | 82 ++++++++ vs-build/drivers/AirfoilInfo_Driver.vfproj | 59 ++++++ 6 files changed, 431 insertions(+), 12 deletions(-) create mode 100644 modules/aerodyn/src/AirfoilInfo_Driver.f90 create mode 100644 vs-build/drivers/AeroAcoustics_Driver.vfproj create mode 100644 vs-build/drivers/AirfoilInfo_Driver.vfproj diff --git a/modules/aerodyn/src/AeroAcoustics.f90 b/modules/aerodyn/src/AeroAcoustics.f90 index 708ab8fe41..93b0515e9d 100644 --- a/modules/aerodyn/src/AeroAcoustics.f90 +++ b/modules/aerodyn/src/AeroAcoustics.f90 @@ -35,7 +35,6 @@ module AeroAcoustics use AeroAcoustics_Types use AeroAcoustics_IO use NWTC_LAPACK - USE NWTC_FFTPACK implicit none private diff --git a/modules/aerodyn/src/AeroAcoustics_Driver.f90 b/modules/aerodyn/src/AeroAcoustics_Driver.f90 index ea5bdbf6bd..81040948de 100644 --- a/modules/aerodyn/src/AeroAcoustics_Driver.f90 +++ b/modules/aerodyn/src/AeroAcoustics_Driver.f90 @@ -16,17 +16,6 @@ ! See the License for the specific language governing permissions and ! limitations under the License. ! - -!there will also be various control flags... this may be updated as needed: -!TBLflag = {'BPM','TNO'} -!bluntnessFlag = {'DTU','BPM'} -!BPMBLflag = {'true','false'} -!useOrigModelAtSepOnset = {'true','false'} - - - - - !********************************************************************************************************************************** program AeroAcoustics_Driver use AeroAcoustics_Driver_Subs diff --git a/modules/aerodyn/src/AirfoilInfo_Driver.f90 b/modules/aerodyn/src/AirfoilInfo_Driver.f90 new file mode 100644 index 0000000000..3c3662d0b4 --- /dev/null +++ b/modules/aerodyn/src/AirfoilInfo_Driver.f90 @@ -0,0 +1,210 @@ +!********************************************************************************************************************************** +! AFI_Driver: This code tests a stand-alone version of the AFI module +!.................................................................................................................................. +! LICENSING +! Copyright (C) 2018 Envision Energy +! +! This file is part of AirfoilInfo. +! +! Licensed under the Apache License, Version 2.0 (the "License"); +! you may not use this file except in compliance with the License. +! You may obtain a copy of the License at +! +! http://www.apache.org/licenses/LICENSE-2.0 +! +! Unless required by applicable law or agreed to in writing, software +! distributed under the License is distributed on an "AS IS" BASIS, +! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +! See the License for the specific language governing permissions and +! limitations under the License. +! +!********************************************************************************************************************************** + + + +program AFI_Driver + + use NWTC_Library + use VersionInfo + use AirfoilInfo + use AirfoilInfo_Types + + implicit none + + + TYPE(ProgDesc), PARAMETER :: AFI_Ver = ProgDesc( 'AFI_driver', '', '' ) + + + ! Variables + type(AFI_InitInputType) :: AFI_InitInputs ! Input data for initialization + integer, parameter :: NumAFI = 1; + type(AFI_ParameterType) :: AFI_p(NumAFI) ! Parameters + type(AFI_OutputType) :: AFI_interp ! interpolated AFI output values + integer(IntKi) :: ErrStat ! Status of error message + character(1024) :: ErrMsg ! Error message if ErrStat /= ErrID_None + + character(1024) :: afName + character(1024) :: outFileName + integer :: unOutFile(3) = -1 + character(*), parameter :: NumFmt = 'ES16.9E2' + character(*), parameter :: Frmt = '(1x,'//NumFmt//')' + character(*), parameter :: Ext(3) = (/ '.cl','.cd','.cm' /) + character(*), parameter :: delim = ' ' + + integer :: i, j, k, iFile + + real(ReKi), allocatable :: Re(:) + real(ReKi) :: alpha + real(ReKi), parameter :: UserProp = 0.0_ReKi + + + ! Initialize the NWTC library + call NWTC_Init(AFI_Ver%Name) + + ! Initialize error handling variables + ErrMsg = '' + ErrStat = ErrID_None + + + CALL DispNVD(AFI_Ver) + + + ! Check for command line arguments. + afName = '' ! default name for input file + CALL CheckArgs( afName ) + + CALL GetRoot( afName, outFileName ) + outFileName = trim(outFileName)//'.interp.out' + + + ! Setup Airfoil InitInput data structure (should come from an input file): + AFI_InitInputs%AFTabMod = AFITable_1 ! AFITable_2Re ! + AFI_InitInputs%InCol_Alfa = 1 + AFI_InitInputs%InCol_Cl = 2 + AFI_InitInputs%InCol_Cd = 3 + AFI_InitInputs%InCol_Cm = 4 + AFI_InitInputs%InCol_Cpmin = 0 + AFI_InitInputs%FileName = afName + AFI_InitInputs%UAMod = UA_Gonzalez + + + ! Write UA parameters to file: + call AFI_WrHeader(delim, trim(outFileName)//'.AFI.sum', unOutFile(1), ErrStat, ErrMsg) + + do iFile=1,NumAFI + if (NumAFI > 1) then + if (iFile < 10) then + AFI_InitInputs%FileName = 'af00'//trim(num2lstr(iFile))//'.dat' + else + AFI_InitInputs%FileName = 'af0'//trim(num2lstr(iFile))//'.dat' + end if + end if + + ! Initialize the Airfoil Info Params + call AFI_Init ( AFI_InitInputs, AFI_p(iFile), ErrStat, ErrMsg ) + call checkError() + + if (ErrStat < AbortErrLev) then + call AFI_WrData(iFile, unOutFile(1), delim, AFI_p(iFile)) + + call AFI_WrTables(AFI_p(iFile), AFI_InitInputs%UAMod, trim(AFI_InitInputs%FileName) ) + + end if + end do + close(unOutFile(1)) + + iFile = 1 + ! allocate Re array, based on Re in the tables + call allocAry( Re, AFI_p(iFile)%NumTabs*2 + 1, 'Re', ErrStat, ErrMsg ) + call checkError() + + Re(1) = AFI_p(iFile)%Table(1)%Re / 2.0_ReKi + do i = 1, AFI_p(iFile)%NumTabs + Re(2*i) = AFI_p(iFile)%Table(i)%Re + end do + do i = 1, AFI_p(iFile)%NumTabs-1 + Re(2*i+1) = (AFI_p(iFile)%Table(i)%Re + AFI_p(iFile)%Table(i+1)%Re)/2.0_ReKi + end do + Re(size(Re)) = AFI_p(iFile)%Table(AFI_p(iFile)%NumTabs)%Re * 2 + + + ! ------------ + do k=1,size(Ext) + call GetNewUnit( unOutFile(k) ) + call OpenFOutFile ( unOutFile(k), trim(outFileName)//Ext(k), errStat, errMsg ) + call checkError() + + write( unOutFile(k), '('//trim(num2lstr(size(Re)+1))//Frmt//')' ) NaN, Re + end do + + + ! time marching loop + do i = -180, 180 ! alpha + alpha = Real(i, ReKi) ! degrees + + do k=1,size(Ext) + write(unOutFile(k), Frmt, ADVANCE='no') alpha + end do + + alpha = alpha*pi/180.0_ReKi ! radians + + do j=1, size(Re) + + call AFI_ComputeAirfoilCoefs( alpha, Re(j), UserProp, AFI_p(iFile), AFI_interp, ErrStat, ErrMsg) + call checkError() + + write(unOutFile(1), Frmt, ADVANCE='no') AFI_interp%cl + write(unOutFile(2), Frmt, ADVANCE='no') AFI_interp%cd + write(unOutFile(3), Frmt, ADVANCE='no') AFI_interp%cm + + end do + + do k=1,size(Ext) + write (unOutFile(k),'()', IOSTAT=ErrStat) ! write the line return + end do + + end do + + + !------------------------------------------------------------------------------------------------- + ! Close our output files + !------------------------------------------------------------------------------------------------- + + + call Cleanup() + call NormStop() + + contains + + !==================================================================================================== + subroutine Cleanup() + ! The routine closes any open files. + !---------------------------------------------------------------------------------------------------- + integer :: ie + + do ie=1,size(Ext) + if (unOutFile(ie) > 0) close( unOutFile(ie), IOSTAT = ErrStat ) + end do + + end subroutine Cleanup + + !---------------------------------------------------------------------------------------------------- + subroutine checkError() + + if (ErrStat >= AbortErrLev) then + + call Cleanup() + call ProgAbort(ErrMsg) + + elseif ( ErrStat /= ErrID_None ) then + + call WrScr( trim(ErrMsg) ) + + end if + + end subroutine checkError + !---------------------------------------------------------------------------------------------------- + + +end program AFI_Driver + diff --git a/vs-build/OpenFAST.sln b/vs-build/OpenFAST.sln index 8ec8938fcb..08412001a9 100644 --- a/vs-build/OpenFAST.sln +++ b/vs-build/OpenFAST.sln @@ -589,6 +589,20 @@ Project("{6989167D-11E4-40FE-8C1A-2192A86A7E90}") = "SoilDyn_Driver", "drivers\S {F8C26AB7-4B2F-4610-A18B-A8C860926B21} = {F8C26AB7-4B2F-4610-A18B-A8C860926B21} EndProjectSection EndProject +Project("{6989167D-11E4-40FE-8C1A-2192A86A7E90}") = "AirfoilInfo_Driver", "drivers\AirfoilInfo_Driver.vfproj", "{3BBE2741-5B28-47BC-9E7F-3E1D172838FB}" + ProjectSection(ProjectDependencies) = postProject + {12DF411B-C7DA-47BA-BB85-7714D5FD2A16} = {12DF411B-C7DA-47BA-BB85-7714D5FD2A16} + {5ADBD025-C654-42C2-BA7C-10F3C3CEEB0E} = {5ADBD025-C654-42C2-BA7C-10F3C3CEEB0E} + {EAF5E602-E6CD-4194-8CCA-0827AA4CCEC9} = {EAF5E602-E6CD-4194-8CCA-0827AA4CCEC9} + EndProjectSection +EndProject +Project("{6989167D-11E4-40FE-8C1A-2192A86A7E90}") = "AeroAcoustics_Driver", "drivers\AeroAcoustics_Driver.vfproj", "{CFB604E9-0D11-4903-BF4C-97D411400C33}" + ProjectSection(ProjectDependencies) = postProject + {12DF411B-C7DA-47BA-BB85-7714D5FD2A16} = {12DF411B-C7DA-47BA-BB85-7714D5FD2A16} + {5ADBD025-C654-42C2-BA7C-10F3C3CEEB0E} = {5ADBD025-C654-42C2-BA7C-10F3C3CEEB0E} + {EAF5E602-E6CD-4194-8CCA-0827AA4CCEC9} = {EAF5E602-E6CD-4194-8CCA-0827AA4CCEC9} + EndProjectSection +EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution Debug|x64 = Debug|x64 @@ -2390,6 +2404,70 @@ Global {9779535B-8DE4-4484-8B59-6E78344D658A}.Release|x64.Build.0 = Release|x64 {9779535B-8DE4-4484-8B59-6E78344D658A}.Release|x86.ActiveCfg = Release|x64 {9779535B-8DE4-4484-8B59-6E78344D658A}.Release|x86.Build.0 = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Debug|x64.ActiveCfg = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Debug|x64.Build.0 = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Debug|x86.ActiveCfg = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Debug|x86.Build.0 = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Debug|x64.ActiveCfg = Double_Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Debug|x64.Build.0 = Double_Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Debug|x86.ActiveCfg = Double_Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Debug|x86.Build.0 = Double_Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_OpenMP_Release|x64.ActiveCfg = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_OpenMP_Release|x64.Build.0 = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_OpenMP_Release|x86.ActiveCfg = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_OpenMP_Release|x86.Build.0 = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Release|x64.ActiveCfg = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Release|x64.Build.0 = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Release|x86.ActiveCfg = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Release|x86.Build.0 = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Debug|x64.ActiveCfg = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Debug|x64.Build.0 = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Debug|x86.ActiveCfg = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Debug|x86.Build.0 = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Release|x64.ActiveCfg = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Release|x64.Build.0 = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Release|x86.ActiveCfg = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Release|x86.Build.0 = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.OpenMP_Release|x64.ActiveCfg = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.OpenMP_Release|x64.Build.0 = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.OpenMP_Release|x86.ActiveCfg = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.OpenMP_Release|x86.Build.0 = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Release|x64.ActiveCfg = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Release|x64.Build.0 = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Release|x86.ActiveCfg = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Release|x86.Build.0 = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Debug|x64.ActiveCfg = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Debug|x64.Build.0 = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Debug|x86.ActiveCfg = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Debug|x86.Build.0 = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Debug|x64.ActiveCfg = Double_Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Debug|x64.Build.0 = Double_Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Debug|x86.ActiveCfg = Double_Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Debug|x86.Build.0 = Double_Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_OpenMP_Release|x64.ActiveCfg = Double_OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_OpenMP_Release|x64.Build.0 = Double_OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_OpenMP_Release|x86.ActiveCfg = Double_OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_OpenMP_Release|x86.Build.0 = Double_OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Release|x64.ActiveCfg = Double_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Release|x64.Build.0 = Double_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Release|x86.ActiveCfg = Double_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Release|x86.Build.0 = Double_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Debug|x64.ActiveCfg = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Debug|x64.Build.0 = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Debug|x86.ActiveCfg = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Debug|x86.Build.0 = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Release|x64.ActiveCfg = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Release|x64.Build.0 = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Release|x86.ActiveCfg = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Release|x86.Build.0 = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.OpenMP_Release|x64.ActiveCfg = OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.OpenMP_Release|x64.Build.0 = OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.OpenMP_Release|x86.ActiveCfg = OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.OpenMP_Release|x86.Build.0 = OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Release|x64.ActiveCfg = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Release|x64.Build.0 = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Release|x86.ActiveCfg = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Release|x86.Build.0 = Release|x64 EndGlobalSection GlobalSection(SolutionProperties) = preSolution HideSolutionNode = FALSE @@ -2453,6 +2531,8 @@ Global {4AB5A895-5B6B-406D-AE36-85F5F5940974} = {76DE56EB-E48F-4F1A-8961-F97A71E7937B} {F8C26AB7-4B2F-4610-A18B-A8C860926B21} = {272B8080-A022-4F4A-BDD6-835871E44C23} {9779535B-8DE4-4484-8B59-6E78344D658A} = {3517E990-350F-4471-A518-8B0BC77CFDDB} + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB} = {3517E990-350F-4471-A518-8B0BC77CFDDB} + {CFB604E9-0D11-4903-BF4C-97D411400C33} = {3517E990-350F-4471-A518-8B0BC77CFDDB} EndGlobalSection GlobalSection(ExtensibilityGlobals) = postSolution SolutionGuid = {B362252D-3254-4C68-B527-CC85CE3CCF75} diff --git a/vs-build/drivers/AeroAcoustics_Driver.vfproj b/vs-build/drivers/AeroAcoustics_Driver.vfproj new file mode 100644 index 0000000000..c408345078 --- /dev/null +++ b/vs-build/drivers/AeroAcoustics_Driver.vfproj @@ -0,0 +1,82 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/vs-build/drivers/AirfoilInfo_Driver.vfproj b/vs-build/drivers/AirfoilInfo_Driver.vfproj new file mode 100644 index 0000000000..6319d2a2ee --- /dev/null +++ b/vs-build/drivers/AirfoilInfo_Driver.vfproj @@ -0,0 +1,59 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From ace18daedeba1caa5f806b897f3bfe9f3ecca0ac Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Wed, 26 Aug 2026 18:14:53 -0600 Subject: [PATCH 03/12] AeroAcoustics: ensure initialization of variables We noticed a couple of spikes that occurred in AA output files during a simulation, but when re-run with the same executable and same input files, the spike disappeared (all non-AA outputs were the same). I asked AI to take a look at the code for uninitialized variables and other memory issues in the AeroAcoustics module, and it came back with these changes. Note the warning about re-entrant code in TNO Assisted-by: Kiro:claude-opus-5 --- modules/aerodyn/src/AeroAcoustics.f90 | 126 ++++++++++++++++++---- modules/aerodyn/src/AeroAcoustics_TNO.f90 | 40 ++++++- 2 files changed, 140 insertions(+), 26 deletions(-) diff --git a/modules/aerodyn/src/AeroAcoustics.f90 b/modules/aerodyn/src/AeroAcoustics.f90 index 93b0515e9d..8a29bda3de 100644 --- a/modules/aerodyn/src/AeroAcoustics.f90 +++ b/modules/aerodyn/src/AeroAcoustics.f90 @@ -367,7 +367,12 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) DO k=1,size(AFInfo) ! for each airfoil interpolation - ! find index where LE is found + ! NOTE: AFInfo(k)%X_Coord(1)/Y_Coord(1) is the aerodynamic center, NOT a point on the airfoil surface. + ! The airfoil surface coordinates are therefore indices 2:NumCoords, and the surface is traversed + ! from the trailing edge, around the leading edge, and back to the trailing edge. + + ! find index where LE is found (i.e., the first index where the x-coordinate starts increasing again) + iLE = -1 ! initialize to an invalid index so we can tell whether the search below succeeded DO i=3,size(AFInfo(k)%X_Coord) IF (AFInfo(k)%X_Coord(i) - AFInfo(k)%X_Coord(i-1) > 0.) THEN iLE = i @@ -375,7 +380,21 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) ENDIF ENDDO - ! From LE toward TE + ! bjj: guard against the search failing; without this, iLE (and the indices derived from it) would be + ! used uninitialized, which produces garbage airfoil thicknesses (and thus garbage inflow noise) + ! that changes from run to run. + IF ( iLE < 3 .OR. iLE >= size(AFInfo(k)%X_Coord) ) THEN + call SetErrStat( ErrID_Fatal, 'Unable to locate the leading edge in the coordinates of airfoil '// & + trim(num2lstr(k))//'. The Full Guidati inflow-noise model (InflowMod=2) requires airfoil coordinates '// & + 'that start at the trailing edge, wrap around the leading edge, and return to the trailing edge.', & + ErrStat, ErrMsg, RoutineName ) + return + ENDIF + + ! From LE toward TE: find the points closest to 1% and 10% chord. + ! Seed the search with the LE point itself so that i1_1/i10_1 are always valid indices. + i1_1 = iLE + i10_1 = iLE dist1 = ABS( AFInfo(k)%X_Coord(iLE) - 0.01) dist10 = ABS( AFInfo(k)%X_Coord(iLE) - 0.10) DO i=iLE+1,size(AFInfo(k)%X_Coord) @@ -389,10 +408,15 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) ENDIF ENDDO - ! From TE to LE - dist1 = 0.99 - dist10 = 0.90 - DO i=1,iLE-1 + ! From TE to LE: find the points closest to 1% and 10% chord. + ! Start at index 2 (index 1 is the aerodynamic center, not a surface point) and seed the search with + ! that point so that i1_2/i10_2 are always valid indices. (The previous code seeded the search with + ! the arbitrary thresholds 0.99/0.90 and could leave these indices undefined.) + i1_2 = 2 + i10_2 = 2 + dist1 = ABS( AFInfo(k)%X_Coord(2) - 0.01) + dist10 = ABS( AFInfo(k)%X_Coord(2) - 0.10) + DO i=3,iLE-1 IF (ABS(AFInfo(k)%X_Coord(i) - 0.01) < dist1) THEN i1_2 = i dist1 = ABS(AFInfo(k)%X_Coord(i) - 0.01) @@ -541,17 +565,35 @@ subroutine Init_MiscVars(m, p, errStat, errMsg) call AllocAry(m%PtotalFreq , size(p%FreqList) , p%NrObsLoc , 'm%PtotalFreq' , errStat2 , errMsg2); if(Failed()) return call AllocAry(m%OASPL , p%NrObsLoc , p%NumBlNds , p%NumBlades , 'm%OASPL' , errStat2 , errMsg2); if(Failed()) return + ! bjj: initialize every allocated array here. Most of these are overwritten on each call to AA_CalcOutput before + ! they are read, but m%ChordAngleLE/m%SpanAngleLE and the SPL* arrays were previously left with whatever happened + ! to be in the freshly allocated memory, which is a latent source of run-to-run differences. m%ChordAngleTE = 0.0_ReKi m%SpanAngleTE = 0.0_ReKi + m%ChordAngleLE = 0.0_ReKi + m%SpanAngleLE = 0.0_ReKi m%rTEtoObserve = 0.0_ReKi m%rLEtoObserve = 0.0_ReKi + m%SPLLBL = 0.0_ReKi + m%SPLP = 0.0_ReKi + m%SPLS = 0.0_ReKi + m%SPLALPH = 0.0_ReKi + m%SPLBLUNT = 0.0_ReKi + m%SPLTIP = 0.0_ReKi + m%SPLTI = 0.0_ReKi m%SPLTIGui = 0.0_ReKi + m%CfVar = 0.0_ReKi m%d99Var = 0.0_ReKi m%dstarVar = 0.0_ReKi m%EdgeVelVar = 0.0_ReKi m%LE_Location = 0.0_ReKi + + m%DirectiviOutput = 0.0_ReKi + m%SumSpecNoiseSep = 0.0_ReKi + m%PtotalFreq = 0.0_ReKi + m%OASPL = 0.0_ReKi contains logical function Failed() call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) @@ -643,6 +685,14 @@ subroutine AA_UpdateStates( t, n, m, u, p, xd, OtherState, errStat, errMsg ) OtherState%allregcounter(k_minus1,rco_minus1) = OtherState%allregcounter(k_minus1,rco_minus1) + 1 ! increase the sample amount in that specific bin + ! bjj: This counter only ever increases, so on a long enough simulation it would overflow IntKi and go + ! negative, which would make the circular-buffer index below zero or negative (an out-of-bounds access). + ! Once we are past the first full window, subtracting exactly Num_total_sampleTI keeps the counter + ! bounded without changing either the test below or the mod() phase of the buffer index. + if ( OtherState%allregcounter(k_minus1,rco_minus1) > 2*p%Num_total_sampleTI ) then + OtherState%allregcounter(k_minus1,rco_minus1) = OtherState%allregcounter(k_minus1,rco_minus1) - p%Num_total_sampleTI + end if + InflowNorm = TwoNorm( u%Inflow(:,j,i) ) !note: p%Num_total_sampleTI = size(xd%RegVxStor,1) ! with storage region dependent moving average and TI @@ -650,7 +700,10 @@ subroutine AA_UpdateStates( t, n, m, u, p, xd, OtherState, errStat, errMsg ) xd%RegVxStor(OtherState%allregcounter(k_minus1,rco_minus1),k_minus1,rco_minus1) = InflowNorm xd%TIVx(j,i) = 0 ELSE - xd%RegVxStor( mod( OtherState%allregcounter(k_minus1,rco_minus1), p%Num_total_sampleTI )+1, k_minus1, rco_minus1)=InflowNorm + ! bjj: use mod(counter-1, N)+1 so that the circular buffer advances in lock step with the counter. + ! The previous expression, mod(counter, N)+1, skipped slot 1 on the first pass, so the very first + ! sample stayed in the window for two full windows instead of one. + xd%RegVxStor( mod( OtherState%allregcounter(k_minus1,rco_minus1)-1, p%Num_total_sampleTI )+1, k_minus1, rco_minus1)=InflowNorm meanInflow = SUM( xd%RegVxStor(:,k_minus1,rco_minus1) ) /p%Num_total_sampleTI if ( EqualRealNos(meanInflow,0.0_ReKi)) then @@ -668,7 +721,12 @@ subroutine AA_UpdateStates( t, n, m, u, p, xd, OtherState, errStat, errMsg ) do j=1,p%NumBlNds ! We scale the incident turbulence intensity by the ratio of average to incident wind speed ! The scaled TI is used by the Amiet model - xd%TIVx(j,i)=p%TI * p%avgV/u%Vrel(J,I) + ! bjj: Vrel can be zero (e.g., before the rotor starts spinning, or at an inboard node), which would + ! make TIVx infinite. InflowNoise() takes 10*log10(TIVx**2), so an Inf here shows up as a huge spike in + ! the inflow-turbulence noise. Apply the same lower bound on relative velocity that AA_CalcOutput uses. + InflowNorm = u%Vrel(J,I) + IF (abs(InflowNorm) < AA_u_min) InflowNorm = SIGN(AA_u_min, InflowNorm) + xd%TIVx(j,i)=p%TI * p%avgV/InflowNorm enddo enddo endif @@ -1489,11 +1547,12 @@ SUBROUTINE TIPNOIS(ALPHTIP,ALPRAT2,C,U ,THETA,PHI, R,p,SPLTIP) UM = MM * p%SpdSound ! Eq 65 from BPM Airfoil Self-noise and Prediction paper TERM = M*M*MM**3*L**2*DBARH/R**2 ! TERM = M^2 * M_max^5 *l^2 *D / r^2 according to Semi-Empirical Aeroacoustic Noise Prediction Code for Wind Turbines paper ! Term is correct according to Eq 61 from BPM Airfoil self-noise and Prediction paper - IF (TERM .NE. 0.0) THEN - SCALE = 10.*LOG10(TERM) - ELSE - SCALE = 0.0 - ENDIF + ! bjj: Use Log10AA here instead of a "TERM /= 0" test around LOG10: + ! (1) TERM is negative when U (and therefore M and MM) is negative, and LOG10 of a negative number returns NaN. + ! (2) The old ELSE branch set SCALE = 0, which left SPLTIP near 126 dB when the directivity DBARH was zero -- + ! i.e., it reported a large tip noise level precisely where there should be none. Every other noise mechanism + ! in this module treats DBARH <= 0 as "no contribution", and Log10AA now gives that behavior here too. + SCALE = 10.*Log10AA(TERM) DO I=1,size(p%FreqList) STPP = p%FreqList(I) * L / UM ! Eq 62 from BPM Airfoil Self-noise and Prediction paper SPLTIP(I) = 126.-30.5*(LOG10AA(STPP)+.3)**2 + SCALE ! Eq 61 from BPM Airfoil Self-noise and Prediction paper @@ -1792,6 +1851,15 @@ SUBROUTINE BLUNT(ALPSTAR,C,U ,THETA,PHI,L,R,H,PSI,p,d99Var2,dstarVar1,dstarVar2, ! Compute average displacement thickness DSTRAVG = (DSTRS + DSTRP) / 2. + + ! bjj: When the boundary layer properties are read from tables (X_BLMethod_Tables), DSTRS and DSTRP can be zero + ! (e.g., for table entries that were not converged). Without this check, HDSTAR becomes Inf (or NaN if H is also + ! zero), which then propagates through G5COMP into SPLBLUNT and on into the summed noise levels. + IF (DSTRAVG <= 0. .OR. H <= 0.) THEN + SPLBLUNT = -100. ! same floor used by the other noise mechanisms; effectively no contribution + RETURN + ENDIF + HDSTAR = H / DSTRAVG DSTARH = 1. /HDSTAR ! Compute directivity function @@ -1877,13 +1945,16 @@ REAL(ReKi) FUNCTION G5COMP(HDSTAR,ETA) result(G5) ETA0 = -SQRT((M*M*MU**4)/(6.25+M*M*MU*MU)) ! eq 80 from BPM Airfoil Self-noise and Prediction paper + ! bjj: the arguments of the SQRT calls below are guarded with MAX(0,...): they are all analytically non-negative + ! over the range of each branch, but the third branch in particular has a radicand that is exactly zero at the + ! branch boundary (SQRT(1.5625/1194.99) = 0.03615995), so round-off can make it slightly negative and return NaN. IF (ETA .LE. ETA0) then - K = 2.5*SQRT(1.-(ETA0/MU)**2)-2.5-M*ETA0 ! eq 81 from BPM Airfoil Self-noise and Prediction paper + K = 2.5*SQRT(MAX(0.0_ReKi, 1.-(ETA0/MU)**2))-2.5-M*ETA0 ! eq 81 from BPM Airfoil Self-noise and Prediction paper G5 = M * ETA + K ! begin eq 76 from BPM Airfoil Self-noise and Prediction paper elseif (ETA .LE. 0.) then - G5 = 2.5*SQRT(1.-(ETA/MU)**2)-2.5 + G5 = 2.5*SQRT(MAX(0.0_ReKi, 1.-(ETA/MU)**2))-2.5 elseif (ETA .LE. 0.03615995) then - G5 = SQRT(1.5625-1194.99*ETA**2)-1.25 + G5 = SQRT(MAX(0.0_ReKi, 1.5625-1194.99*ETA**2))-1.25 else G5 = -155.543 * ETA + 4.375 end if @@ -2178,8 +2249,10 @@ SUBROUTINE TBLTE_TNO(U,THETA,PHI,D,R,Cfall,d99all,EdgeVelAll,p,SPLP,SPLS) n_freq = size(p%FreqList) freq = p%FreqList - SPLS = 0.0_ReKi ! initialize in case Cfall(1) <= 0 - SPLP = 0.0_ReKi ! initialize in case Cfall(2) <= 0 + ! bjj: -100 dB (the floor used elsewhere in this module) means "no contribution". The previous value of 0 dB + ! actually adds unit mean-square pressure at every frequency when Cfall <= 0. + SPLS = -100.0_ReKi ! initialize in case Cfall(1) <= 0 + SPLP = -100.0_ReKi ! initialize in case Cfall(2) <= 0 ! Body of TNO band_ratio = 2.**(1./3.) @@ -2193,18 +2266,25 @@ SUBROUTINE TBLTE_TNO(U,THETA,PHI,D,R,Cfall,d99all,EdgeVelAll,p,SPLP,SPLS) do i_omega = 1,n_freq omega = TwoPi*p%FreqList(i_omega) !integration limits + ! bjj: use ABS(Mach) so that the upper limit stays above the lower limit. U (and therefore Mach) carries the + ! sign of Vrel, and a negative upper limit would reverse the integration interval and return a bogus spectrum. int_limits(1) = 0.0e0 - int_limits(2) = 10*omega/(Mach*p%SpdSound) + int_limits(2) = 10*omega/(ABS(Mach)*p%SpdSound) ! Convert to third octave band_width = 2. * omega * (sqrt(band_ratio)-1./sqrt(band_ratio)) + ! bjj: The log10 calls below were unguarded. "answer" is the result of a fixed-order Gauss-Kronrod quadrature + ! of an integrand that underflows to zero at high frequency and can also return a small negative value through + ! round-off cancellation. log10(0) is -Inf and log10(negative) is NaN, and the "SPL < -100" floor below does not + ! catch NaN (any comparison with NaN is false), so a NaN would propagate silently into the summed noise levels. + ! Log10AA() floors the argument at AA_EPSILON, which maps both cases to "no contribution". IF (Cfall(1) .GT. 0.) THEN answer = SPL_integrate(omega=omega,limits=int_limits,ISSUCTION=.true., & Mach=Mach,SpdSound=p%SpdSound,AirDens=p%AirDens,KinVisc=p%KinVisc, & Cfall=Cfall,d99all=d99all,EdgeVelAll=EdgeVelAll) Spectrum = D/(4.*pi*R**2)*answer - SPL_suction = 10.*log10(Spectrum*DBARH/2.e-5/2.e-5) - SPLS(i_omega) = SPL_suction + 10.*log10(band_width) + SPL_suction = 10.*Log10AA(Spectrum*DBARH/2.e-5/2.e-5) + SPLS(i_omega) = SPL_suction + 10.*Log10AA(band_width) ENDIF IF (Cfall(2) .GT. 0.) THEN @@ -2212,8 +2292,8 @@ SUBROUTINE TBLTE_TNO(U,THETA,PHI,D,R,Cfall,d99all,EdgeVelAll,p,SPLP,SPLS) Mach=Mach,SpdSound=p%SpdSound,AirDens=p%AirDens,KinVisc=p%KinVisc, & Cfall=Cfall,d99all=d99all,EdgeVelAll=EdgeVelAll) Spectrum = D/(4.*pi*R**2)*answer - SPL_press = 10.*log10(Spectrum*DBARH/2.e-5/2.e-5) - SPLP(i_omega) = SPL_press + 10.*log10(band_width) + SPL_press = 10.*Log10AA(Spectrum*DBARH/2.e-5/2.e-5) + SPLP(i_omega) = SPL_press + 10.*Log10AA(band_width) ENDIF ! Sum the noise sources SPLALPH is BPM value diff --git a/modules/aerodyn/src/AeroAcoustics_TNO.f90 b/modules/aerodyn/src/AeroAcoustics_TNO.f90 index 00053096fb..3de7e1b940 100644 --- a/modules/aerodyn/src/AeroAcoustics_TNO.f90 +++ b/modules/aerodyn/src/AeroAcoustics_TNO.f90 @@ -10,6 +10,12 @@ MODULE TNO INTEGER, PARAMETER :: TNOKi = ReKi + ! NOTE (bjj): The variables declared at module scope below are effectively SAVEd state. They are how the integrand + ! functions f_int1/f_int2/Pressure receive their parameters, because slatec_qk61 only accepts a function of one + ! variable. SPL_integrate() sets all of them on entry, so sequential calls (including calls for different rotors or + ! different blade nodes) are safe. This module is NOT reentrant, however: it must not be called from more than one + ! thread at a time, or one thread will overwrite another thread's flow conditions mid-integration. + REAL (TNOKi), PARAMETER :: Cnuk = 5.5 REAL (TNOKi), PARAMETER :: kappa = 0.41 REAL (TNOKi), PARAMETER :: Cmu = 0.09 @@ -123,14 +129,33 @@ FUNCTION f_int1(x2) delta = d99(2) Uo=Mach_TNO*co*edgevel(2) endif - if (Cfin .le. 0.) then - write(*,*) 'Cf is less than zero, Cf = ',Cfin - stop + ! bjj: Bail out (contributing nothing to the integral) instead of producing NaN/Inf or killing the program. + ! - Cf <= 0 used to execute a bare "stop", which aborts without OpenFAST's error handling and without closing + ! output files. TBLTE_TNO already skips the side whose Cf is non-positive, so this is a backstop. + ! - delta (d99) can be zero for unconverged entries in the pre-tabulated boundary-layer files. That makes L = 0/0 + ! and pi*x2/delta infinite below. + ! - u_star is zero when the edge velocity ratio or the Mach number is zero, and log(u_star*x2/nu) is then -Inf, + ! so U evaluates to 0*(-Inf) = NaN. + ! - L is zero at x2 = 0 (the lower integration limit), which makes ke = sqrt(pi)/L infinite. + if (Cfin .le. 0. .or. delta .le. 0. .or. x2 .le. 0.) then + f_int1 = 0. + RETURN endif + u_star = Uo*sqrt(Cfin/2.) + if (u_star .le. 0.) then + f_int1 = 0. + RETURN + endif + L = 0.085*delta*tanh(kappa*x2/(0.085*delta)) + if (L .le. 0.) then + f_int1 = 0. + RETURN + endif + if (x2 .gt. delta)then U = Uo dudx = 0. @@ -165,6 +190,15 @@ END FUNCTION f_int1 FUNCTION f_int2(k1_in) ! changed name from 'int2' to avoid conflicts with intrinsic of same name REAL (TNOKi), intent(in) :: k1_in REAL (TNOKi) :: f_int2 + + ! bjj: The lower integration limit passed to slatec_qk61 is exactly zero. The 61-point Gauss-Kronrod rule does not + ! evaluate the integrand at the interval end points, so k1_in is not zero in practice, but guard the 1/k1_in here + ! (and the k1**2/(k1**2+k3**2) = 0/0 in Pressure) so this does not depend on the internals of the quadrature rule. + if (k1_in .le. 0.0_TNOKi) then + f_int2 = 0.0_TNOKi + RETURN + endif + f_int2 = Omega_TNO/co/k1_in*Pressure(k1_in) RETURN END FUNCTION f_int2 From 3e0bc9ab69eeabc2db6d481b300ff79e361be743 Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Thu, 27 Aug 2026 11:52:45 -0600 Subject: [PATCH 04/12] AeroAcoustics: avoid error about 100% BldPrcnt when there is only 1 node The AA driver allows only one node for testing purposes, but if the BldPrct isn't 100 in that case, it would ALWAYS give a warning about changing the aerodynamic mesh. This check avoids printing the message in that case. Assisted-by: Kiro:claude-opus-5 --- modules/aerodyn/src/AeroAcoustics.f90 | 41 ++++++++++++++++++--------- 1 file changed, 27 insertions(+), 14 deletions(-) diff --git a/modules/aerodyn/src/AeroAcoustics.f90 b/modules/aerodyn/src/AeroAcoustics.f90 index 8a29bda3de..c369d087d7 100644 --- a/modules/aerodyn/src/AeroAcoustics.f90 +++ b/modules/aerodyn/src/AeroAcoustics.f90 @@ -271,20 +271,32 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) p%BlSpn = InitInp%BlSpn p%BlChord = InitInp%BlChord - ! Calculate last element size as percentage of blade span + ! Calculate last element size as percentage of blade span, and warn if BldPrcnt is smaller than that. + ! bjj: This check only has meaning for a blade with more than one node. A single node is an explicitly supported + ! configuration (AA_Init requires only NumBlNds >= 1, and the standalone AeroAcoustics driver uses exactly + ! one node): that one element spans the whole blade, so there is no "last element" for BldPrcnt to be smaller + ! than. LastElemPct is 100.0 for a single node, so the warning fired on every single-node + ! run with BldPrcnt < 100 and told the user to refine an AeroDyn blade mesh that has nothing to refine. + LastElemPct = 100.0_ReKi ! a single element spans the entire blade 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 + + ! bjj: guard the divide. This would be rare, and an indication that the calling program has an error. + IF (p%BlSpn(p%NumBlNds,1) > 0.0_ReKi) THEN + LastElemPct = 100.0 * (p%BlSpn(p%NumBlNds,1) - p%BlSpn(p%NumBlNds-1,1)) / p%BlSpn(p%NumBlNds,1) + ELSE + CALL SetErrStat(ErrID_Fatal, 'The outermost blade node is at zero span, which is invalid.', ErrStat2, ErrMsg2, RoutineName ) + if(Failed()) return + 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 '// & - 'or refine the aerodynamic mesh in your blade AeroDyn input file.', ErrStat2, ErrMsg2, RoutineName ) - if(Failed()) return - endif + if (InputFileData%AA_Bl_Prcntge < 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 '// & + 'or refine the aerodynamic mesh in your blade AeroDyn input file.', ErrStat2, ErrMsg2, RoutineName ) + if(Failed()) return + endif + + ENDIF p%startnode = min(p%NumBlNds, 2) BladeSpanUsedForNoise = p%BlSpn(p%NumBlNds,1)*(1.0 - InputFileData%AA_Bl_Prcntge/100.0) @@ -300,11 +312,12 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) DO I = 1,p%numBlades DO J = p%startnode,p%NumBlNds ! starts loop from startnode. IF (J < 2) THEN + ! bjj: J == 1 is reachable only when NumBlNds == 1 (otherwise startnode >= 2). p%BlElemSpn(J,I) = p%BlSpn(J,I) !assume this is the innermost node ELSEIF (J .EQ. p%NumBlNds) THEN p%BlElemSpn(J,I) = p%BlSpn(J,I)-p%BlSpn(J-1,I) ELSE - p%BlElemSpn(J,I) = (p%BlSpn(J,I)-p%BlSpn(J-1,I))/2 + (p%BlSpn(J+1,I)-p%BlSpn(J,I))/2 ! this is the average element size around this node, equivalent to (p%BlSpn(J+1,I) - p%BlSpn(J-1,I))/2 + p%BlElemSpn(J,I) = (p%BlSpn(J+1,I) - p%BlSpn(J-1,I)) / 2 ! this is the average element size around this node ENDIF end do end do @@ -1048,7 +1061,7 @@ SUBROUTINE CalcAeroAcousticsOutput(u,p,m,xd,errStat,errMsg) !--------Tip Noise--------------------------------------------------------------! - IF ( (p%ITIP == ITIP_ON) .AND. (J .EQ. p%NumBlNds) ) THEN ! calculate m%SPLTIP(1:nFreq) + IF ( (p%ITIP == ITIP_ON) .AND. (J == p%NumBlNds) ) THEN ! calculate m%SPLTIP(1:nFreq) CALL TIPNOIS(AlphaNoise_Deg,p%ALpRAT,p%BlChord(J,I),UNoise,m%ChordAngleTE(K,J,I),m%SpanAngleTE(K,J,I), & m%rTEtoObserve(K,J,I), p, m%SPLTIP) From 749864634d32e5dd8907d9570839192ea38bab9c Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Thu, 27 Aug 2026 12:01:33 -0600 Subject: [PATCH 05/12] Ignore: don't include AI steering files for specific users/machines In the future, we can upload steering that is more general, but for now, I'm ignoring my .kiro files. --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.gitignore b/.gitignore index 74be88339b..95dfda90f4 100644 --- a/.gitignore +++ b/.gitignore @@ -69,3 +69,6 @@ openfast_io/openfast_io/_version.py *.code-workspace docker-compose.yml .devcontainer/ + +# AI steering files for specific users +.kiro/* \ No newline at end of file From 5d8b16c677071adbc01c84fb5e070666bc912d0a Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Fri, 28 Aug 2026 16:32:53 -0600 Subject: [PATCH 06/12] Aeroacoustics: harden 10^x calculations Make sure the 10^x calculations don't overflow. All 10**(x) calculations have been replaced with Power10AA(x), which makes sure the exponent doesn't get larger than new variable AA_max_exp. --- modules/aerodyn/src/AeroAcoustics.f90 | 64 ++++++++++++++++----------- 1 file changed, 37 insertions(+), 27 deletions(-) diff --git a/modules/aerodyn/src/AeroAcoustics.f90 b/modules/aerodyn/src/AeroAcoustics.f90 index c369d087d7..ab795588a7 100644 --- a/modules/aerodyn/src/AeroAcoustics.f90 +++ b/modules/aerodyn/src/AeroAcoustics.f90 @@ -48,6 +48,9 @@ module AeroAcoustics REAL(ReKi), parameter :: AA_u_min = 0.1_ReKi REAL(ReKi), parameter :: AA_EPSILON = 1.E-16 ! EPSILON(AA_EPSILON) + REAL(ReKi), parameter :: AA_max_exp = log10(huge(1.0_ReKi))-8.0 ! bjj: I picked 8 here somewhat randomly...in debugging, later calculations + ! with the results of this function multiplied by factors of 1E6 + ! so I removed a couple more orders of magnitude (1E8) REAL(ReKi), parameter :: RotorRegionAlph_delta = 60.0_ReKi ! degrees : size of bin, must be a number that evenly divides 360 degrees REAL(ReKi), parameter :: RotorRegionRad_delta = 5.0_ReKi ! meters : size of bin along blade span (rotor radius) @@ -526,8 +529,8 @@ subroutine Init_y(y, m, p, errStat, errMsg) end if call AllocAry(y%WriteOutput , p%numOutsAll(1), 'y%WriteOutput' , errStat2 , errMsg2); if(Failed()) return - call AllocAry(y%WriteOutputSep , p%numOutsAll(3), 'y%WriteOutputSep' , errStat2 , errMsg2); if(Failed()) return call AllocAry(y%WriteOutputForPE , p%numOutsAll(2), 'y%WriteOutputForPE' , errStat2 , errMsg2); if(Failed()) return + call AllocAry(y%WriteOutputSep , p%numOutsAll(3), 'y%WriteOutputSep' , errStat2 , errMsg2); if(Failed()) return call AllocAry(y%WriteOutputNodes , p%numOutsAll(4), 'y%WriteOutputSepFreq' , errStat2 , errMsg2); if(Failed()) return y%WriteOutput = 0.0_reki @@ -849,6 +852,13 @@ REAL(ReKi) FUNCTION Log10AA(X) RESULT(F) F = LOG10( MAX(AA_EPSILON, X) ) END FUNCTION Log10AA +!---------------------------------------------------------------------------------------------------------------------------------- +REAL(ReKi) FUNCTION Power10AA(X) RESULT(F) + REAL(ReKi),INTENT(IN) :: X + + F = 10.0_ReKi ** min(X, AA_max_exp) + +END FUNCTION Power10AA !----------------------------------------------------------------------------------------------------------------------------------! SUBROUTINE Calc_LE_Location_Array(p,m,u) TYPE(AA_ParameterType), intent(in ) :: p !< Parameters @@ -1161,7 +1171,7 @@ subroutine TotalContributionFromType(SPL,Ptotal,NoiseMech) do III=1,size(p%FreqList) ! Loops through each 1/3rd octave center frequency - Pt = 10.0_ReKi**(SPL(III)/10.0_ReKi) ! SPL to P Conversion for III Frequency + Pt = Power10AA(SPL(III)/10.0_ReKi) ! SPL to P Conversion for III Frequency P_SumAllFreq = P_SumAllFreq + Pt ! Sum for Running Total m%SumSpecNoiseSep(NoiseMech,III,K) = m%SumSpecNoiseSep(NoiseMech,III,K) + Pt ! Running sum of observer and frequency dependent sound pressure @@ -1235,13 +1245,13 @@ SUBROUTINE LBLVS(ALPSTAR,C,U,THETA,PHI,L,R,p,d99Var2,dstarVar1,dstarVar2,SPLLAM, else ST1PRIM = .28 end if - STPKPRM = 10.**(-.04*ALPSTAR) * ST1PRIM ! Eq 56 from BPM Airfoil Self-noise and Prediction paper + STPKPRM = Power10AA(-.04*ALPSTAR) * ST1PRIM ! Eq 56 from BPM Airfoil Self-noise and Prediction paper ! compute reference reynolds number ! Eq 59 from BPM Airfoil Self-noise and Prediction paper IF (ALPSTAR .LE. 3.0) then - RC0=10.**(.215*ALPSTAR+4.978) + RC0=Power10AA(.215*ALPSTAR+4.978) else - RC0=10.**(.120*ALPSTAR+5.263) + RC0=Power10AA(.120*ALPSTAR+5.263) end if ! compute peak scaled spectrum level @@ -1390,7 +1400,7 @@ SUBROUTINE TBLTE(ALPSTAR,C,U,THETA,PHI,L,R,p,d99Var2,dstarVar1,dstarVar2,StallVa IF (ALPSTAR .LE. 1.333) then ST2 = ST1 elseif (ALPSTAR .LE. StallVal) then - ST2 = ST1*10.**(.0054*(ALPSTAR-1.333)**2) + ST2 = ST1*Power10AA(.0054*(ALPSTAR-1.333)**2) else ST2 = 4.72 * ST1 end if @@ -1501,9 +1511,9 @@ SUBROUTINE TBLTE(ALPSTAR,C,U,THETA,PHI,L,R,p,d99Var2,dstarVar1,dstarVar2,StallVa IF (SPLS(I) .LT. -100.) SPLS(I) = -100. ! Similar to Eq 29 of BPM Airfoil Self-noise and Prediction paper IF (SPLALPH(I) .LT. -100.) SPLALPH(I) = -100. ! Eq 30 of BPM Airfoil Self-noise and Prediction paper recommends SPLALPH = 10log(stuff) + A' + K2, where A' is calculated same as A but with x3 Rc - !P1 = 10.**(SPLP(I) / 10.) ! SPL_Pressure - !P2 = 10.**(SPLS(I) / 10.) ! SPL_Suction - !P4 = 10.**(SPLALPH(I) / 10.) ! SPL_AoA + !P1 = Power10AA(SPLP(I) / 10.) ! SPL_Pressure + !P2 = Power10AA(SPLS(I) / 10.) ! SPL_Suction + !P4 = Power10AA(SPLALPH(I) / 10.) ! SPL_AoA !SPLTBL(I) = 10. * LOG10AA(P1 + P2 + P4) ! Eq 24 from BPM Airfoil Self-noise and Prediction paper @@ -1911,11 +1921,11 @@ SUBROUTINE BLUNT(ALPSTAR,C,U ,THETA,PHI,L,R,H,PSI,p,d99Var2,dstarVar1,dstarVar2, G5(I) = G50 + .0714 * PSI * (G514-G50) ! interpolate G5 from G50 and G514 IF (G5(I) .GT. 0.) G5(I) = 0. - G5Sum = 10**(G5(I)/10)+G5Sum ! to be subtracted + G5Sum = Power10AA(G5(I)/10)+G5Sum ! to be subtracted if ( G5Sum .ne. 0) then - LogVal = MAX(AA_EPSILON,1/G5Sum) + LogVal = MAX(AA_EPSILON,1/G5Sum) else - LogVal = 1 + LogVal = AA_max_exp end if SPLBLUNT(I) = G4 + G5(I) + SCALE - 10*log10(LogVal) ! equation mentioned there is plus but it is stated subtract, thus ''- 10*log10(1/G5Sum)'' end do @@ -2089,12 +2099,12 @@ SUBROUTINE THICK(C,RC,ALPSTAR,p,DELTAP,DSTRS,DSTRP,StallVal) LogRC = LOG10AA( RC ) ! Boundary layer thickness - DELTA0 = 10.**(1.6569-0.9045*LogRC+0.0596*LogRC**2)*C ! (untripped) Eq. (5) of [1] - IF (p%ITRIP /= ITRIP_None) DELTA0 = 10.**(1.892 -0.9045*LogRC+0.0596*LogRC**2)*C ! (heavily tripped) Eq. (2) of [1] + DELTA0 = Power10AA(1.6569-0.9045*LogRC+0.0596*LogRC**2)*C ! (untripped) Eq. (5) of [1] + IF (p%ITRIP /= ITRIP_None) DELTA0 = Power10AA(1.892 -0.9045*LogRC+0.0596*LogRC**2)*C ! (heavily tripped) Eq. (2) of [1] IF (p%ITRIP .EQ. ITRIP_Light) DELTA0=.6*DELTA0 ! Pressure side boundary layer thickness, Eq (8) of [1] - DELTAP = 10.**(-.04175*ALPSTAR+.00106*ALPSTAR**2)*DELTA0 + DELTAP = Power10AA(-.04175*ALPSTAR+.00106*ALPSTAR**2)*DELTA0 ! Compute zero angle of attack displacement thickness IF (p%ITRIP /= ITRIP_None) THEN @@ -2102,37 +2112,37 @@ SUBROUTINE THICK(C,RC,ALPSTAR,p,DELTAP,DSTRS,DSTRP,StallVal) IF (RC .LE. .3E+06) THEN DSTR0 = .0601 * RC **(-.114)*C ELSE - DSTR0=10.**(3.411-1.5397*LogRC+.1059*LogRC**2)*C + DSTR0=Power10AA(3.411-1.5397*LogRC+.1059*LogRC**2)*C END IF ! Lightly tripped IF (p%ITRIP .EQ. ITRIP_Light) DSTR0 = DSTR0 * .6 ELSE ! Untripped, Eq. (6) of [1] - DSTR0=10.**(3.0187-1.5397*LogRC+.1059*LogRC**2)*C + DSTR0=Power10AA(3.0187-1.5397*LogRC+.1059*LogRC**2)*C ENDIF ! Pressure side displacement thickness, Eq. (9) of [1] - DSTRP = 10.**(-.0432*ALPSTAR+.00113*ALPSTAR**2)*DSTR0 + DSTRP = Power10AA(-.0432*ALPSTAR+.00113*ALPSTAR**2)*DSTR0 ! IF (p%ITRIP .EQ. 3) DSTRP = DSTRP * 1.48 ! commented since itrip is never 3 check if meant 2.(EB_DTU) ! Suction side displacement thickness IF (p%ITRIP .EQ. ITRIP_Heavy) THEN ! Heavily tripped, Eq. (12) of [1] IF (ALPSTAR .LE. 5.) THEN - DSTRS=10.**(.0679*ALPSTAR)*DSTR0 + DSTRS=Power10AA(.0679*ALPSTAR)*DSTR0 ELSEIF (ALPSTAR .LE. StallVal) THEN - DSTRS = 0.381 * 10.**(.1516*ALPSTAR)*DSTR0 + DSTRS = 0.381 * Power10AA(.1516*ALPSTAR)*DSTR0 ELSE - DSTRS = 14.296 * 10.**(.0258*ALPSTAR)*DSTR0 + DSTRS = 14.296 * Power10AA(.0258*ALPSTAR)*DSTR0 ENDIF ELSE ! Untripped or lightly tripped, Eq. (15) of [1] IF (ALPSTAR .LE. 7.5) THEN - DSTRS =10.**(.0679*ALPSTAR)*DSTR0 + DSTRS =Power10AA(.0679*ALPSTAR)*DSTR0 ELSEIF(ALPSTAR .LE. StallVal) THEN - DSTRS = .0162*10.**(.3066*ALPSTAR)*DSTR0 + DSTRS = .0162*Power10AA(.3066*ALPSTAR)*DSTR0 ELSE - DSTRS = 52.42*10.**(.0258*ALPSTAR)*DSTR0 + DSTRS = 52.42*Power10AA(.0258*ALPSTAR)*DSTR0 ENDIF ENDIF @@ -2313,9 +2323,9 @@ SUBROUTINE TBLTE_TNO(U,THETA,PHI,D,R,Cfall,d99all,EdgeVelAll,p,SPLP,SPLS) IF (SPLP(i_omega) .LT. -100.) SPLP(i_omega) = -100. IF (SPLS(i_omega) .LT. -100.) SPLS(i_omega) = -100. - !P1 = 10.**(SPLP(i_omega) / 10.) - !P2 = 10.**(SPLS(i_omega) / 10.) - !P4 = 10.**(SPLALPH(i_omega) / 10.) + !P1 = Power10AA(SPLP(i_omega) / 10.) + !P2 = Power10AA(SPLS(i_omega) / 10.) + !P4 = Power10AA(SPLALPH(i_omega) / 10.) ! !SPLTBL(i_omega) = 10. * LOG10(P1 + P2 + P4) enddo From 890e1996af4776eccdf8bf33156e2ba077c361de Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Fri, 28 Aug 2026 16:51:53 -0600 Subject: [PATCH 07/12] AeroAcoustics: ensure Re and Mach numbers are positive in case U < 0, we want to make sure that the Reynold's Number and Mach Numbers are positive --- modules/aerodyn/src/AeroAcoustics.f90 | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/modules/aerodyn/src/AeroAcoustics.f90 b/modules/aerodyn/src/AeroAcoustics.f90 index ab795588a7..c97475ef63 100644 --- a/modules/aerodyn/src/AeroAcoustics.f90 +++ b/modules/aerodyn/src/AeroAcoustics.f90 @@ -1217,8 +1217,8 @@ SUBROUTINE LBLVS(ALPSTAR,C,U,THETA,PHI,L,R,p,d99Var2,dstarVar1,dstarVar2,SPLLAM, integer(intKi) :: I ! I A generic index for DO loops. !compute reynolds number and mach number - M = U / p%SpdSound ! MACH NUMBER - RC = U * C/p%KinVisc ! REYNOLDS NUMBER BASED ON CHORD + M = abs(U / p%SpdSound) ! MACH NUMBER + RC = abs(U * C/p%KinVisc) ! REYNOLDS NUMBER BASED ON CHORD ! compute boundary layer thicknesses IF (p%X_BLMethod .eq. X_BLMethod_Tables) THEN @@ -1368,8 +1368,8 @@ SUBROUTINE TBLTE(ALPSTAR,C,U,THETA,PHI,L,R,p,d99Var2,dstarVar1,dstarVar2,StallVa LOGICAL :: SWITCH !!LOGICAL FOR COMPUTATION OF ANGLE OF ATTACK CONTRIBUTION ! Compute reynolds number and mach number - M = U / p%SpdSound - RC = U * C/p%KinVisc + M = abs(U / p%SpdSound) + RC = abs(U * C/p%KinVisc) ! Compute boundary layer thicknesses IF (p%X_BLMethod .eq. X_BLMethod_Tables) THEN @@ -1390,8 +1390,8 @@ SUBROUTINE TBLTE(ALPSTAR,C,U,THETA,PHI,L,R,p,d99Var2,dstarVar1,dstarVar2,StallVa ! RETURN ! ENDIF ! Calculate the reynolds numbers based on pressure and suction displacement thickness - RDSTRS = DSTRS * U / p%KinVisc - RDSTRP = DSTRP * U / p%KinVisc + RDSTRS = abs(DSTRS * U / p%KinVisc) !bjj: should this be absolute value? + RDSTRP = abs(DSTRP * U / p%KinVisc) ! Determine peak strouhal numbers to be used for 'a' and 'b' curve calculations ST1 = .02 * M ** (-.6) ! Eq 32 from BPM Airfoil Self-noise and Prediction paper @@ -1554,7 +1554,7 @@ SUBROUTINE TIPNOIS(ALPHTIP,ALPRAT2,C,U ,THETA,PHI, R,p,SPLTIP) ENDIF !! used to be ALPTIPP = ALPHTIP * ALPRAT2 ALPTIPP = ABS(ALPHTIP) * ALPRAT2 - M = U / p%SpdSound ! MACH NUMBER + M = abs(U / p%SpdSound) ! MACH NUMBER ! Compute directivity function DBARH = DIRECTH_TE(M,THETA,PHI) IF (p%ROUND) THEN @@ -1617,7 +1617,7 @@ SUBROUTINE InflowNoise(AlphaNoise,Chord,U,THETA,PHI,d,RObs,TINoise,p,SPLti) INTEGER(intKi) :: I !I A generic index for DO loops. !!!--- NAF NOISE IDENTICAL - Mach = U/p%SpdSound + Mach = abs(U/p%SpdSound) ! This part is recently added for height and surface roughness dependent estimation of turbulence intensity and turbulence scales !%Lturb=300*(Z/300)^(0.46+0.074*log(p%z0_aa)); !% Gives larger length scale @@ -1861,8 +1861,8 @@ SUBROUTINE BLUNT(ALPSTAR,C,U ,THETA,PHI,L,R,H,PSI,p,d99Var2,dstarVar1,dstarVar2, real(ReKi) :: LogVal ! temp variable to help us not take log10(0) --- ! Reynolds number and mach number - M = U / p%SpdSound - RC = U * C/p%KinVisc + M = abs(U / p%SpdSound) + RC = abs(U * C/p%KinVisc) ! Compute boundary layer thicknesses IF (p%X_BLMethod .eq. X_BLMethod_Tables) THEN DELTAP = d99Var2 @@ -2281,7 +2281,7 @@ SUBROUTINE TBLTE_TNO(U,THETA,PHI,D,R,Cfall,d99all,EdgeVelAll,p,SPLP,SPLS) band_ratio = 2.**(1./3.) ! Mach number - Mach = U / p%SpdSound + Mach = abs(U / p%SpdSound) ! Directivity function DBARH = DIRECTH_TE(REAL(Mach,ReKi),THETA,PHI) @@ -2292,7 +2292,7 @@ SUBROUTINE TBLTE_TNO(U,THETA,PHI,D,R,Cfall,d99all,EdgeVelAll,p,SPLP,SPLS) ! bjj: use ABS(Mach) so that the upper limit stays above the lower limit. U (and therefore Mach) carries the ! sign of Vrel, and a negative upper limit would reverse the integration interval and return a bogus spectrum. int_limits(1) = 0.0e0 - int_limits(2) = 10*omega/(ABS(Mach)*p%SpdSound) + int_limits(2) = 10*omega/(Mach*p%SpdSound) ! Convert to third octave band_width = 2. * omega * (sqrt(band_ratio)-1./sqrt(band_ratio)) @@ -2358,7 +2358,7 @@ SUBROUTINE BL_Param_Interp(p,m,U,AlphaNoise_Deg,C,whichAirfoil) !!!! this if is not used but if necessary two sets of tables can be populated for tripped and untripped cases - RC = U * C/p%KinVisc ! REYNOLDS NUMBER BASED ON CHORD + RC = abs(U * C/p%KinVisc) ! REYNOLDS NUMBER BASED ON CHORD ! find the indices into the arrays representing coordinates of each dimension: From 2341bda4f8f1b69a00d6dc9099111b7da9d9d1b6 Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Fri, 28 Aug 2026 16:57:03 -0600 Subject: [PATCH 08/12] Aeroacoustics: fix comments --- modules/aerodyn/src/AeroAcoustics.f90 | 4 ++-- modules/aerodyn/src/AeroAcoustics_IO.f90 | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/modules/aerodyn/src/AeroAcoustics.f90 b/modules/aerodyn/src/AeroAcoustics.f90 index c97475ef63..a4f19ad0cd 100644 --- a/modules/aerodyn/src/AeroAcoustics.f90 +++ b/modules/aerodyn/src/AeroAcoustics.f90 @@ -206,8 +206,8 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) ! if tno is on or one of the guidati models is on, check if we have airfoil coordinates DO k=1,size(AFInfo) ! if any of the airfoil coordinates are missing change calculation method IF( AFInfo(k)%NumCoords .lt. 5 )then - CALL WrScr( 'Airfoil coordinates are missing: If Full or Simplified Guidati or Bl Calculation is on coordinates are needed ' ) - CALL WrScr( 'Calculation methods enforced as BPM for TBLTE and only Amiet for inflow ' ) + CALL WrScr( 'Airfoil coordinates are missing: If Full or Simplified Guidati or Bl Calculation are used, coordinates are needed.' ) + CALL WrScr( 'Calculation methods have been changed: BPM for TBLTE and only Amiet (BPM) for inflow ' ) p%ITURB = ITURB_BPM p%IInflow = IInflow_BPM exit ! stop checking do loop diff --git a/modules/aerodyn/src/AeroAcoustics_IO.f90 b/modules/aerodyn/src/AeroAcoustics_IO.f90 index 24049bb023..3fa7110495 100644 --- a/modules/aerodyn/src/AeroAcoustics_IO.f90 +++ b/modules/aerodyn/src/AeroAcoustics_IO.f90 @@ -46,7 +46,7 @@ MODULE AeroAcoustics_IO integer(intKi), parameter :: ITURB_None = 0 ! TBLTE noise is not calculated integer(intKi), parameter :: ITURB_BPM = 1 ! TBLTE noise is calculated with BPM - integer(intKi), parameter :: ITURB_TNO = 2 ! TBLTE noise is calculated with TNO + integer(intKi), parameter :: ITURB_TNO = 2 ! TBLTE noise is calculated with TNO (and part of BPM) integer(intKi), parameter :: IInflow_None = 0 ! IInflow noise is not calculated integer(intKi), parameter :: IInflow_BPM = 1 ! IInflow noise is calculated with BPM @@ -452,7 +452,7 @@ SUBROUTINE ValidateInputData( InputFileData, NumBl, ErrStat, ErrMsg ) if (InputFileData%IBLUNT /= IBLUNT_None .and. InputFileData%IBLUNT /= IBLUNT_BPM) then call SetErrStat ( ErrID_Fatal, & - 'IBLUNT must '//trim(num2lstr(IBLUNT_None))//' (none) or '//trim(num2lstr(IBLUNT_BPM))//' (Bluntness noise calculated).', ErrStat, ErrMsg, RoutineName ) + 'IBLUNT must be '//trim(num2lstr(IBLUNT_None))//' (none) or '//trim(num2lstr(IBLUNT_BPM))//' (Bluntness noise calculated).', ErrStat, ErrMsg, RoutineName ) endif if (InputFileData%ILAM /= ILAM_None .and. InputFileData%ilam /= ILAM_BPM) then call SetErrStat ( ErrID_Fatal, 'ILAM must be '//trim(num2lstr(ILAM_None))//' No calculation '//& @@ -467,7 +467,7 @@ SUBROUTINE ValidateInputData( InputFileData, NumBl, ErrStat, ErrMsg ) ' (heavily tripped BL Calculation) or '//trim(num2lstr(ITRIP_Light))//' (lightly tripped BL)' ,ErrStat, ErrMsg, RoutineName ) end if if (InputFileData%ITURB /= ITURB_None .and. InputFileData%ITURB /= ITURB_BPM .and. InputFileData%ITURB /= ITURB_TNO) then - call SetErrStat ( ErrID_Fatal, 'ITURB must be 0 (off) or 1 (BPM) or 2 (TNO) .', ErrStat, ErrMsg, RoutineName ) + call SetErrStat ( ErrID_Fatal, 'ITURB must be 0 (off) or 1 (BPM) or 2 (TNO with BPM alpha) .', ErrStat, ErrMsg, RoutineName ) end if if (InputFileData%IInflow /= IInflow_None .and. InputFileData%IInflow /= IInflow_BPM & .and. InputFileData%IInflow /= IInflow_FullGuidati .and. InputFileData%IInflow /= IInflow_SimpleGuidati ) then From 9c2587bf9fd5c04ba8860df9d46d4c93beebab21 Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Fri, 28 Aug 2026 18:28:07 -0600 Subject: [PATCH 09/12] AeroAcoustics cleanup - Added a warning message if we ever end up with overflow, which indicates the noise predictions are unrealistically large. - Made `Unoise` be a positive value. It's used for Reynolds Number, Mach Number, and Strouhal Number. Making it the absolute value avoids having to add abs() around 7 more lines of code. - fixed initialization of LogVal in the Blunt routine, which I had erroneously added in a previous commit Assisted-by: Kiro:claude-opus-5 --- modules/aerodyn/src/AeroAcoustics.f90 | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/modules/aerodyn/src/AeroAcoustics.f90 b/modules/aerodyn/src/AeroAcoustics.f90 index a4f19ad0cd..4b8c52fe91 100644 --- a/modules/aerodyn/src/AeroAcoustics.f90 +++ b/modules/aerodyn/src/AeroAcoustics.f90 @@ -856,6 +856,19 @@ END FUNCTION Log10AA REAL(ReKi) FUNCTION Power10AA(X) RESULT(F) REAL(ReKi),INTENT(IN) :: X + !note that this "SAVE" is a deviation from standard OpenFAST code. Because it is used only to throw a warning if we ever + ! trigger numerical overflow, I'm not concerned about it getting reset per instance. If we don't like this behavior, + ! we could store this FirstWarning as a misc var and thread it into every call, or just remove the warning alltogether. + LOGICAL, SAVE :: FirstWarning = .TRUE. + + if (FirstWarning) then + if (X > AA_max_exp) then + call WrScr(' WARNING: noise predictions are extremely large. Use caution in interpreting results.') + FirstWarning = .FALSE. + end if + end if + + F = 10.0_ReKi ** min(X, AA_max_exp) END FUNCTION Power10AA @@ -1006,11 +1019,7 @@ SUBROUTINE CalcAeroAcousticsOutput(u,p,m,xd,errStat,errMsg) DO I = 1,p%numBlades DO J = p%startnode,p%NumBlNds ! starts loop from startnode. !------------------------------!!------------------------------!!------------------------------!!------------------------------! - - Unoise = u%Vrel(J,I) - IF (abs(Unoise) < AA_u_min) then - Unoise = SIGN(AA_u_min, Unoise) - ENDIF + Unoise = max(abs(u%Vrel(J,I)), AA_u_min) AlphaNoise= u%AoANoise(J,I) call MPi2Pi(AlphaNoise) ! make sure this is in an appropriate range [-pi,pi] @@ -1390,8 +1399,8 @@ SUBROUTINE TBLTE(ALPSTAR,C,U,THETA,PHI,L,R,p,d99Var2,dstarVar1,dstarVar2,StallVa ! RETURN ! ENDIF ! Calculate the reynolds numbers based on pressure and suction displacement thickness - RDSTRS = abs(DSTRS * U / p%KinVisc) !bjj: should this be absolute value? - RDSTRP = abs(DSTRP * U / p%KinVisc) + RDSTRS = DSTRS * U / p%KinVisc + RDSTRP = DSTRP * U / p%KinVisc ! Determine peak strouhal numbers to be used for 'a' and 'b' curve calculations ST1 = .02 * M ** (-.6) ! Eq 32 from BPM Airfoil Self-noise and Prediction paper @@ -1925,7 +1934,7 @@ SUBROUTINE BLUNT(ALPSTAR,C,U ,THETA,PHI,L,R,H,PSI,p,d99Var2,dstarVar1,dstarVar2, if ( G5Sum .ne. 0) then LogVal = MAX(AA_EPSILON,1/G5Sum) else - LogVal = AA_max_exp + LogVal = 1.0_ReKi ! note that this is reached only when the accumulated bluntness contribution is already ~ 10^-69 end if SPLBLUNT(I) = G4 + G5(I) + SCALE - 10*log10(LogVal) ! equation mentioned there is plus but it is stated subtract, thus ''- 10*log10(1/G5Sum)'' end do From 319336acac7714729ec199f681b16acadd2ba91d Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Wed, 9 Sep 2026 14:46:09 -0600 Subject: [PATCH 10/12] Update .gitignore with vs-build/.gitignore entries the vs-build/.gitignore was removed and these entries would be lost --- .gitignore | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 95dfda90f4..56647ad7a2 100644 --- a/.gitignore +++ b/.gitignore @@ -71,4 +71,10 @@ docker-compose.yml .devcontainer/ # AI steering files for specific users -.kiro/* \ No newline at end of file +.kiro/* + +# vs-build files +.vs +*.user +*.u2d +gitVersionInfo.h From 685cd85b3c04e099e0f8302d831fad51bf586598 Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Thu, 10 Sep 2026 14:11:43 -0600 Subject: [PATCH 11/12] AeroAcoustics TNO: replace non-reentrant SLATEC integrator This replaces the NWTC_SLATEC.f90 code with an integration scheme that is reentrant. Global variables are replaced with subroutine arguments stored in a data type. Assisted-by: Kiro-claude-opus-5 --- modules/aerodyn/src/AeroAcoustics_TNO.f90 | 332 +++++++++++++++++----- 1 file changed, 257 insertions(+), 75 deletions(-) diff --git a/modules/aerodyn/src/AeroAcoustics_TNO.f90 b/modules/aerodyn/src/AeroAcoustics_TNO.f90 index 3de7e1b940..c0d3c46731 100644 --- a/modules/aerodyn/src/AeroAcoustics_TNO.f90 +++ b/modules/aerodyn/src/AeroAcoustics_TNO.f90 @@ -2,7 +2,7 @@ MODULE TNO use NWTC_Library ! ReKi, DBKi, R8Ki - use NWTC_SLATEC ! slatec_qk61 -- which is all that is in that library right now. + !use NWTC_SLATEC ! slatec_qk61: bjj: this has been replaced with a local qk61_ctx that passes a context (TNO_ContextType) to the integrand. See comments below. implicit none PRIVATE @@ -10,38 +10,111 @@ MODULE TNO INTEGER, PARAMETER :: TNOKi = ReKi - ! NOTE (bjj): The variables declared at module scope below are effectively SAVEd state. They are how the integrand - ! functions f_int1/f_int2/Pressure receive their parameters, because slatec_qk61 only accepts a function of one - ! variable. SPL_integrate() sets all of them on entry, so sequential calls (including calls for different rotors or - ! different blade nodes) are safe. This module is NOT reentrant, however: it must not be called from more than one - ! thread at a time, or one thread will overwrite another thread's flow conditions mid-integration. - REAL (TNOKi), PARAMETER :: Cnuk = 5.5 REAL (TNOKi), PARAMETER :: kappa = 0.41 REAL (TNOKi), PARAMETER :: Cmu = 0.09 ! INTEGER(IntKi),PARAMETER :: limit = 5000 - !TNO variables - REAL (TNOKi) :: Omega_TNO ! NOTE: not a constant and used by function f_int1 and f_int2 + !> Everything the TNO integrands need in order to be evaluated. + !! + !! This type replaces what used to be a set of module-level (and therefore implicitly SAVEd) variables. The reason + !! those existed is that the SLATEC integrator qk61 declares its integrand as "external f" and calls it as f(x), so + !! there was nowhere to thread the flow conditions through. Passing this context as an explicit argument instead + !! makes the module reentrant: every activation of SPL_integrate owns its own context, so concurrent calls (e.g. + !! from an OpenMP region, or two rotors evaluated in parallel) can no longer overwrite each other's state + !! mid-integration. + type :: TNO_ContextType + ! frequency + real(TNOKi) :: Omega = 0.0_TNOKi !< radian frequency + ! atmosphere + real(TNOKi) :: nu = 0.0_TNOKi !< kinematic viscosity + real(TNOKi) :: co = 0.0_TNOKi !< speed of sound + real(TNOKi) :: rho = 0.0_TNOKi !< air density + ! airfoil + real(TNOKi) :: Mach = 0.0_TNOKi !< Mach number of this blade node + logical :: IsSuction = .false. !< .true. for the suction side, .false. for the pressure side + ! blade-node boundary-layer properties; index 1 = suction side, index 2 = pressure side + real(TNOKi) :: d99(2) = 0.0_TNOKi !< boundary layer thickness + real(TNOKi) :: Cf(2) = 0.0_TNOKi !< skin friction coefficient + real(TNOKi) :: edgevel(2) = 0.0_TNOKi !< edge velocity ratio + ! Wavenumbers. Unlike the fields above (which are fixed for the whole integration), these are updated by + ! Pressure() at each point of the OUTER integration over k1 and then read by f_int1 during the INNER + ! integration over x2. That mid-integration coupling is why the context is passed as intent(inout). + real(TNOKi) :: k1 = 0.0_TNOKi + real(TNOKi) :: k3 = 0.0_TNOKi + real(TNOKi) :: k = 0.0_TNOKi + end type TNO_ContextType - !atmosphere variables - REAL (TNOKi) :: nu - REAL (TNOKi) :: co - REAL (TNOKi) :: rho + abstract interface + !> Integrand accepted by qk61_ctx: a function of the integration variable plus a context. + real(TNOKi) function TNO_Integrand(x, ctx) + import :: TNOKi, TNO_ContextType + real(TNOKi), intent(in ) :: x + type(TNO_ContextType), intent(inout) :: ctx + end function TNO_Integrand + end interface - ! Wavenumber variables - REAL (TNOKi) :: k - REAL (TNOKi) :: k1 - REAL (TNOKi) :: k3 + ! ------------------------------------------------------------------------------------------------------------------ + ! Data for the 61-point Gauss-Kronrod rule, transcribed from modules/nwtc-library/src/NetLib/slatec/dqk61.f + ! (SLATEC/QUADPACK; weights and abscissae evaluated with 80-decimal-digit arithmetic by L. W. Fullerton, Bell Labs, + ! Nov. 1981). The abscissae and weights are given for the interval (-1,1); by symmetry only the non-negative + ! abscissae and their corresponding weights are stored. + ! + ! The literals below carry dqk61's full precision and are converted to TNOKi, so a single-precision build rounds to + ! exactly the values hard-coded in qk61.f and a double-precision build matches dqk61.f. That is what lets this + ! routine reproduce the previous results bit-for-bit in both precisions. + ! + ! xgk - abscissae of the 61-point Kronrod rule. + ! xgk(2), xgk(4), ... are the abscissae of the 30-point Gauss rule. + ! xgk(1), xgk(3), ... are the optimally added abscissae. + ! wgk - weights of the 61-point Kronrod rule. + ! wg - weights of the 30-point Gauss rule. + ! ------------------------------------------------------------------------------------------------------------------ + real(TNOKi), parameter :: xgk(31) = (/ & + 0.999484410050490637571325895705811_TNOKi, 0.996893484074649540271630050918695_TNOKi, & + 0.991630996870404594858628366109486_TNOKi, 0.983668123279747209970032581605663_TNOKi, & + 0.973116322501126268374693868423707_TNOKi, 0.960021864968307512216871025581798_TNOKi, & + 0.944374444748559979415831324037439_TNOKi, 0.926200047429274325879324277080474_TNOKi, & + 0.905573307699907798546522558925958_TNOKi, 0.882560535792052681543116462530226_TNOKi, & + 0.857205233546061098958658510658944_TNOKi, 0.829565762382768397442898119732502_TNOKi, & + 0.799727835821839083013668942322683_TNOKi, 0.767777432104826194917977340974503_TNOKi, & + 0.733790062453226804726171131369528_TNOKi, 0.697850494793315796932292388026640_TNOKi, & + 0.660061064126626961370053668149271_TNOKi, 0.620526182989242861140477556431189_TNOKi, & + 0.579345235826361691756024932172540_TNOKi, 0.536624148142019899264169793311073_TNOKi, & + 0.492480467861778574993693061207709_TNOKi, 0.447033769538089176780609900322854_TNOKi, & + 0.400401254830394392535476211542661_TNOKi, 0.352704725530878113471037207089374_TNOKi, & + 0.304073202273625077372677107199257_TNOKi, 0.254636926167889846439805129817805_TNOKi, & + 0.204525116682309891438957671002025_TNOKi, 0.153869913608583546963794672743256_TNOKi, & + 0.102806937966737030147096751318001_TNOKi, 0.051471842555317695833025213166723_TNOKi, & + 0.000000000000000000000000000000000_TNOKi /) - ! Blade params - REAL (TNOKi) :: d99(2) - REAL (TNOKi) :: Cf(2) - REAL (TNOKi) :: edgevel(2) + real(TNOKi), parameter :: wgk(31) = (/ & + 0.001389013698677007624551591226760_TNOKi, 0.003890461127099884051267201844516_TNOKi, & + 0.006630703915931292173319826369750_TNOKi, 0.009273279659517763428441146892024_TNOKi, & + 0.011823015253496341742232898853251_TNOKi, 0.014369729507045804812451432443580_TNOKi, & + 0.016920889189053272627572289420322_TNOKi, 0.019414141193942381173408951050128_TNOKi, & + 0.021828035821609192297167485738339_TNOKi, 0.024191162078080601365686370725232_TNOKi, & + 0.026509954882333101610601709335075_TNOKi, 0.028754048765041292843978785354334_TNOKi, & + 0.030907257562387762472884252943092_TNOKi, 0.032981447057483726031814191016854_TNOKi, & + 0.034979338028060024137499670731468_TNOKi, 0.036882364651821229223911065617136_TNOKi, & + 0.038678945624727592950348651532281_TNOKi, 0.040374538951535959111995279752468_TNOKi, & + 0.041969810215164246147147541285970_TNOKi, 0.043452539701356069316831728117073_TNOKi, & + 0.044814800133162663192355551616723_TNOKi, 0.046059238271006988116271735559374_TNOKi, & + 0.047185546569299153945261478181099_TNOKi, 0.048185861757087129140779492298305_TNOKi, & + 0.049055434555029778887528165367238_TNOKi, 0.049795683427074206357811569379942_TNOKi, & + 0.050405921402782346840893085653585_TNOKi, 0.050881795898749606492297473049805_TNOKi, & + 0.051221547849258772170656282604944_TNOKi, 0.051426128537459025933862879215781_TNOKi, & + 0.051494729429451567558340433647099_TNOKi /) - ! Airfoil - REAL(TNOKi) :: Mach_TNO - LOGICAL :: ISSUCTION_TNO + real(TNOKi), parameter :: wg(15) = (/ & + 0.007968192496166605615465883474674_TNOKi, 0.018466468311090959142302131912047_TNOKi, & + 0.028784707883323369349719179611292_TNOKi, 0.038799192569627049596801936446348_TNOKi, & + 0.048402672830594052902938140422808_TNOKi, 0.057493156217619066481721689402056_TNOKi, & + 0.065974229882180495128128515115962_TNOKi, 0.073755974737705206268243850022191_TNOKi, & + 0.080755895229420215354694938460530_TNOKi, 0.086899787201082979802387530715126_TNOKi, & + 0.092122522237786128717632707087619_TNOKi, 0.096368737174644259639468626351810_TNOKi, & + 0.099593420586795267062780282103569_TNOKi, 0.101762389748405504596428952168554_TNOKi, & + 0.102852652893558840341285636705415_TNOKi /) contains @@ -62,36 +135,144 @@ function SPL_integrate(Omega,limits,ISSUCTION, & real(ReKi), intent(in ) :: EdgeVelAll(2) !< real(ReKi) :: integrand !< integrand result - real(TNOKi) :: answer !< value returned from qk61, NOTE the typing + real(TNOKi) :: answer !< value returned from qk61_ctx, NOTE the typing + + !> All state needed by the integrands. This is a local variable (not module data), which is what makes this + !! routine safe to call from more than one thread at a time. + type(TNO_ContextType) :: ctx ! local variables that are ignored real(TNOKi) :: abserr,resabs,resasc !< accuracy estimates and residuals. Currently ignored - ! Set module values from input - ISSUCTION_TNO = ISSUCTION - Omega_TNO = real(Omega,TNOKi) + ! Set context values from input + ctx%IsSuction = ISSUCTION + ctx%Omega = real(Omega,TNOKi) ! Mach number of segment - Mach_TNO = real(Mach,TNOKi) + ctx%Mach = real(Mach,TNOKi) ! Atmospheric values - co = real(SpdSound, TNOKi) - rho = real(AirDens, TNOKi) - nu = real(KinVisc, TNOKi) + ctx%co = real(SpdSound, TNOKi) + ctx%rho = real(AirDens, TNOKi) + ctx%nu = real(KinVisc, TNOKi) ! Blade node values - Cf = real(Cfall, TNOKi) - d99 = real(d99all, TNOKi) - edgevel = real(ABS(EdgeVelAll),TNOKi) + ctx%Cf = real(Cfall, TNOKi) + ctx%d99 = real(d99all, TNOKi) + ctx%edgevel = real(ABS(EdgeVelAll),TNOKi) - call slatec_qk61(f_int2,limits(1),limits(2),answer,abserr,resabs,resasc) + call qk61_ctx(f_int2,limits(1),limits(2),ctx,answer,abserr,resabs,resasc) integrand = real( answer, ReKi ) end function SPL_integrate +!================================================================================================================================== +!> Integrate FUNC over (a,b) using the 61-point Gauss-Kronrod rule, passing CTX through to the integrand. +!! +!! This is a direct port of the SLATEC routine qk61 (see modules/nwtc-library/src/NetLib/slatec/qk61.f), changed only +!! so that the integrand takes a context argument. The accumulation order is deliberately identical to the original so +!! that the computed RESULT is bit-for-bit the same as before this port. +!! +!! Two things the original could not offer: +!! 1. an explicit interface for the integrand (qk61 used "external f", an implicit interface), and +!! 2. a way to give the integrand its parameters other than global data. +!! +!! Declared RECURSIVE because the TNO model nests two integrations: the outer integral over k1 calls f_int2, which +!! calls Pressure, which integrates over x2. The original relied on per-file compiler flags (-frecursive / +!! -assume recursion, see modules/nwtc-library/CMakeLists.txt) to make qk61's locals automatic; stating RECURSIVE in +!! the source expresses that requirement in the code instead of in the build system. +recursive subroutine qk61_ctx(func, a, b, ctx, result, abserr, resabs, resasc) + procedure(TNO_Integrand) :: func !< integrand, evaluated as func(x, ctx) + real(TNOKi), intent(in ) :: a !< lower limit of integration + real(TNOKi), intent(in ) :: b !< upper limit of integration + type(TNO_ContextType), intent(inout) :: ctx !< context handed to the integrand + real(TNOKi), intent( out) :: result !< approximation to the integral, from the 61-point Kronrod rule + real(TNOKi), intent( out) :: abserr !< estimate of the modulus of the absolute error + real(TNOKi), intent( out) :: resabs !< approximation to the integral of abs(func) + real(TNOKi), intent( out) :: resasc !< approximation to the integral of abs(func-i/(b-a)) + + ! Local variables (names kept from qk61 to keep this reviewable against the original) + real(TNOKi) :: absc ! abscissa + real(TNOKi) :: centr ! mid point of the interval + real(TNOKi) :: dhlgth ! abs(hlgth) + real(TNOKi) :: epmach ! the largest relative spacing + real(TNOKi) :: fc ! function value at the mid point + real(TNOKi) :: fsum + real(TNOKi) :: fval1, fval2 ! function values + real(TNOKi) :: fv1(30), fv2(30) + real(TNOKi) :: hlgth ! half-length of the interval + real(TNOKi) :: resg ! result of the 30-point Gauss rule + real(TNOKi) :: resk ! result of the 61-point Kronrod rule + real(TNOKi) :: reskh ! approximation to the mean value of func over (a,b), i.e. to i/(b-a) + real(TNOKi) :: uflow ! the smallest positive magnitude + integer :: j, jtw, jtwm1 + + ! r1mach(4)/d1mach(4) is the largest relative spacing and r1mach(1)/d1mach(1) the smallest positive magnitude; + ! EPSILON and TINY are the standard intrinsics for exactly those quantities, so the SLATEC *1mach dependency is + ! not needed here. + epmach = epsilon(1.0_TNOKi) + uflow = tiny(1.0_TNOKi) + + centr = 0.5_TNOKi*(b+a) + hlgth = 0.5_TNOKi*(b-a) + dhlgth = abs(hlgth) + + ! Compute the 61-point Kronrod approximation to the integral, and estimate the absolute error. + resg = 0.0_TNOKi + fc = func(centr, ctx) + resk = wgk(31)*fc + resabs = abs(resk) + + do j = 1,15 + jtw = j*2 + absc = hlgth*xgk(jtw) + fval1 = func(centr-absc, ctx) + fval2 = func(centr+absc, ctx) + fv1(jtw) = fval1 + fv2(jtw) = fval2 + fsum = fval1+fval2 + resg = resg+wg(j)*fsum + resk = resk+wgk(jtw)*fsum + resabs = resabs+wgk(jtw)*(abs(fval1)+abs(fval2)) + end do + + do j = 1,15 + jtwm1 = j*2-1 + absc = hlgth*xgk(jtwm1) + fval1 = func(centr-absc, ctx) + fval2 = func(centr+absc, ctx) + fv1(jtwm1) = fval1 + fv2(jtwm1) = fval2 + fsum = fval1+fval2 + resk = resk+wgk(jtwm1)*fsum + resabs = resabs+wgk(jtwm1)*(abs(fval1)+abs(fval2)) + end do + + reskh = resk*0.5_TNOKi + resasc = wgk(31)*abs(fc-reskh) + do j = 1,30 + resasc = resasc+wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh)) + end do + + result = resk*hlgth + resabs = resabs*dhlgth + resasc = resasc*dhlgth + abserr = abs((resk-resg)*hlgth) + if (resasc /= 0.0_TNOKi .and. abserr /= 0.0_TNOKi) & + abserr = resasc*min(1.0_TNOKi,(200.0_TNOKi*abserr/resasc)**1.5_TNOKi) + if (resabs > uflow/(50.0_TNOKi*epmach)) & + abserr = max((epmach*50.0_TNOKi)*resabs,abserr) + +end subroutine qk61_ctx +!================================================================================================================================== + + +FUNCTION f_int1(x2, ctx) result(f) + REAL(TNOKi), intent(in) :: x2 + type(TNO_ContextType), intent(inout) :: ctx + REAL(TNOKi):: f -FUNCTION f_int1(x2) REAL(TNOKi):: alpha REAL(TNOKi):: alpha_gauss REAL(TNOKi):: Cfin @@ -111,23 +292,21 @@ FUNCTION f_int1(x2) REAL(TNOKi):: Uc REAL(TNOKi):: Uo REAL(TNOKi):: W - REAL(TNOKi), intent(in) :: x2 - REAL(TNOKi):: f_int1 ! changed and being multiplied with edge velocity taken from xfoil output - ! Uo=Mach_TNO*co ISSUCTION_TNO use edgevel(1) + ! Uo=ctx%Mach*ctx%co ctx%IsSuction use ctx%edgevel(1) !constants from xfoil - if (ISSUCTION_TNO) then + if (ctx%IsSuction) then alpha = 0.45 ! = 0.3 pressure, = 0.45 suction - Cfin = Cf(1) - delta = d99(1) - Uo=Mach_TNO*co*edgevel(1) + Cfin = ctx%Cf(1) + delta = ctx%d99(1) + Uo=ctx%Mach*ctx%co*ctx%edgevel(1) else alpha = 0.30 - Cfin = Cf(2) - delta = d99(2) - Uo=Mach_TNO*co*edgevel(2) + Cfin = ctx%Cf(2) + delta = ctx%d99(2) + Uo=ctx%Mach*ctx%co*ctx%edgevel(2) endif ! bjj: Bail out (contributing nothing to the integral) instead of producing NaN/Inf or killing the program. ! - Cf <= 0 used to execute a bare "stop", which aborts without OpenFAST's error handling and without closing @@ -138,39 +317,39 @@ FUNCTION f_int1(x2) ! so U evaluates to 0*(-Inf) = NaN. ! - L is zero at x2 = 0 (the lower integration limit), which makes ke = sqrt(pi)/L infinite. if (Cfin .le. 0. .or. delta .le. 0. .or. x2 .le. 0.) then - f_int1 = 0. + f = 0. RETURN endif u_star = Uo*sqrt(Cfin/2.) if (u_star .le. 0.) then - f_int1 = 0. + f = 0. RETURN endif L = 0.085*delta*tanh(kappa*x2/(0.085*delta)) if (L .le. 0.) then - f_int1 = 0. + f = 0. RETURN endif if (x2 .gt. delta)then U = Uo dudx = 0. - f_int1 = 0. + f = 0. RETURN else W = 1.-cos(pi*x2/delta); - U = u_star*(1./kappa*log(u_star*x2/nu) +Cnuk+ (Uo/u_star-1./kappa*log(u_star*delta/nu)-Cnuk)*0.5*W) - dudx = u_star*(1./(kappa*x2)+(Uo/u_star-1./kappa*log(u_star*delta/nu)-Cnuk)* & + U = u_star*(1./kappa*log(u_star*x2/ctx%nu) +Cnuk+ (Uo/u_star-1./kappa*log(u_star*delta/ctx%nu)-Cnuk)*0.5*W) + dudx = u_star*(1./(kappa*x2)+(Uo/u_star-1./kappa*log(u_star*delta/ctx%nu)-Cnuk)* & 0.5*(pi/delta)*sin(pi*x2/delta)) endif ke=sqrt(pi)/L*0.4213560764 !gamma(5./6.)/gamma(1./3.) - k1_hat = k1/ke - k3_hat = k3/ke + k1_hat = ctx%k1/ke + k3_hat = ctx%k3/ke Nut = (L*kappa)**2.*abs(dudx) kT = sqrt((Nut*dudx)**2./Cmu) @@ -179,54 +358,57 @@ FUNCTION f_int1(x2) Uc = 0.7*U alpha_gauss = 0.05*Uc/L - phim = 1./(alpha_gauss*sqrt(pi))*exp(-((Omega_TNO-Uc*k1)/alpha_gauss)**2.) + phim = 1./(alpha_gauss*sqrt(pi))*exp(-((ctx%Omega-Uc*ctx%k1)/alpha_gauss)**2.) phi22 = 4./9./pi*1/ke**2.*(k1_hat**2.+k3_hat**2.)/(1.+k1_hat**2.+k3_hat**2.)**(7./3.) - f_int1 = L*ums*(dudx)**2*phi22*phim*exp(-2*abs(k)*x2) + f = L*ums*(dudx)**2*phi22*phim*exp(-2*abs(ctx%k)*x2) RETURN END FUNCTION f_int1 -FUNCTION f_int2(k1_in) ! changed name from 'int2' to avoid conflicts with intrinsic of same name +FUNCTION f_int2(k1_in, ctx) result(f) ! changed name from 'int2' to avoid conflicts with intrinsic of same name REAL (TNOKi), intent(in) :: k1_in - REAL (TNOKi) :: f_int2 + type(TNO_ContextType), intent(inout) :: ctx + REAL (TNOKi) :: f - ! bjj: The lower integration limit passed to slatec_qk61 is exactly zero. The 61-point Gauss-Kronrod rule does not + ! bjj: The lower integration limit passed to qk61_ctx is exactly zero. The 61-point Gauss-Kronrod rule does not ! evaluate the integrand at the interval end points, so k1_in is not zero in practice, but guard the 1/k1_in here ! (and the k1**2/(k1**2+k3**2) = 0/0 in Pressure) so this does not depend on the internals of the quadrature rule. if (k1_in .le. 0.0_TNOKi) then - f_int2 = 0.0_TNOKi + f = 0.0_TNOKi RETURN endif - f_int2 = Omega_TNO/co/k1_in*Pressure(k1_in) + f = ctx%Omega/ctx%co/k1_in*Pressure(k1_in, ctx) RETURN END FUNCTION f_int2 -FUNCTION Pressure(k1_in) +FUNCTION Pressure(k1_in, ctx) result(P) ! Variables + REAL(TNOKi), intent(in) :: k1_in + type(TNO_ContextType), intent(inout) :: ctx + real(TNOKi) :: P + REAL(TNOKi) :: a,b,answer REAL(TNOKi) :: abserr,resabs,resasc - REAL(TNOKi) :: k1_in - real(TNOKi) :: Pressure - ! Set variables used in f_int1 - k1 = k1_in + ! Set wavenumbers used in f_int1 + ctx%k1 = k1_in - a = 0.0_TNOKi !1e-4*d99(1) - IF (ISSUCTION_TNO)THEN - b = d99(1) + a = 0.0_TNOKi !1e-4*ctx%d99(1) + IF (ctx%IsSuction)THEN + b = ctx%d99(1) ELSE - b = d99(2) + b = ctx%d99(2) ENDIF - k3 = 0. - k= sqrt(k1**2+k3**2) + ctx%k3 = 0. + ctx%k = sqrt(ctx%k1**2+ctx%k3**2) - CALL slatec_qk61(f_int1,a,b,answer,abserr,resabs,resasc) + CALL qk61_ctx(f_int1,a,b,ctx,answer,abserr,resabs,resasc) - Pressure = 4.0_TNOKi*rho**2 * k1**2 / (k1**2 + k3**2)*answer + P = 4.0_TNOKi*ctx%rho**2 * ctx%k1**2 / (ctx%k1**2 + ctx%k3**2)*answer RETURN END FUNCTION Pressure From ff32e31cb875007a166b42a1abd1c3247496448f Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Fri, 11 Sep 2026 16:41:22 -0600 Subject: [PATCH 12/12] Build: remove unused files, fix VS ADI C Binding - after the previous commit fixing reentrant code in AeroAcoustics, the ScaLAPACK routines are not needed in the NWTC Library build. I removed them from the cmake file and Visual Studio project. - In rebuilding, I found and fixed an error with the Visual Studio AeroDyn_Inflow_C_Binding project sending the wrong argument to run the registry, so it was generating the ADI instead of ADI_c_bindings Types file. - The SubDyn_Types.f90 file had a different comment when I regenerated the file. --- modules/nwtc-library/CMakeLists.txt | 28 +++---- modules/subdyn/src/SubDyn_Types.f90 | 4 +- .../AeroDyn_Inflow_C_Binding.vfproj | 20 ++--- vs-build/modules/NWTC-Library.vfproj | 73 ++----------------- 4 files changed, 31 insertions(+), 94 deletions(-) diff --git a/modules/nwtc-library/CMakeLists.txt b/modules/nwtc-library/CMakeLists.txt index c40457fcad..51a43b70f8 100644 --- a/modules/nwtc-library/CMakeLists.txt +++ b/modules/nwtc-library/CMakeLists.txt @@ -106,20 +106,20 @@ set(NWTCLIBS_SOURCES src/NetLib/scalapack/NWTC_ScaLAPACK.f90 # NetLib SLATEC sources - src/NetLib/slatec/NWTC_SLATEC.f90 - src/NetLib/slatec/dqk61.f - src/NetLib/slatec/qk61.f - src/NetLib/slatec/d1mach.f - src/NetLib/slatec/r1mach.f - src/NetLib/slatec/xercnt.f - src/NetLib/slatec/xerhlt.f - src/NetLib/slatec/xerprn.f - src/NetLib/slatec/xersve.f - src/NetLib/slatec/fdump.f - src/NetLib/slatec/i1mach.f - src/NetLib/slatec/j4save.f - src/NetLib/slatec/xgetua.f - src/NetLib/slatec/xermsg.f + #src/NetLib/slatec/NWTC_SLATEC.f90 + #src/NetLib/slatec/dqk61.f + #src/NetLib/slatec/qk61.f + #src/NetLib/slatec/d1mach.f + #src/NetLib/slatec/r1mach.f + #src/NetLib/slatec/xercnt.f + #src/NetLib/slatec/xerhlt.f + #src/NetLib/slatec/xerprn.f + #src/NetLib/slatec/xersve.f + #src/NetLib/slatec/fdump.f + #src/NetLib/slatec/i1mach.f + #src/NetLib/slatec/j4save.f + #src/NetLib/slatec/xgetua.f + #src/NetLib/slatec/xermsg.f ) get_filename_component(FCNAME ${CMAKE_Fortran_COMPILER} NAME) diff --git a/modules/subdyn/src/SubDyn_Types.f90 b/modules/subdyn/src/SubDyn_Types.f90 index 9fee52b5f2..5bdaa7d288 100644 --- a/modules/subdyn/src/SubDyn_Types.f90 +++ b/modules/subdyn/src/SubDyn_Types.f90 @@ -170,7 +170,7 @@ MODULE SubDyn_Types REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: MemberSpin !< Member spin angle about its axis - for rectangular members [rad] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: MemberDivSize !< Optional maximum element length for each member [m] INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: MemberNDiv !< Resolved number of finite elements per member [-] - INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: MemberElemStart !< First element index for each member in p%Elems [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: MemberElemStart !< First element index for each member in p%Elems [-] CHARACTER(ChanLen) , DIMENSION(:), ALLOCATABLE :: SSOutList !< List of Output Channels [-] LOGICAL :: OutCOSM = .false. !< Output Cos-matrices Flag [-] LOGICAL :: TabDelim = .false. !< Generate a tab-delimited output file in OutJckF-Flag [-] @@ -196,7 +196,7 @@ MODULE SubDyn_Types REAL(R8Ki) , DIMENSION(:,:), ALLOCATABLE :: K !< System stiffness matrix [-] REAL(R8Ki) , DIMENSION(:,:), ALLOCATABLE :: M !< System mass matrix [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: ElemProps !< Element properties(A, L, Ixx, Iyy, Jzz, Shear, Kappa, E, G, Rho, DirCos(1,1), DirCos(2, 1), ....., DirCos(3, 3) ) [-] - INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: MemberNodes !< Member number and list of nodes making up a member (>2 if subdivided) [-] + INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: MemberNodes !< Member number and endpoint node IDs (interior member nodes reconstructed from connectivity) [-] INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: NodesConnN !< Nodes that connect to a common node [-] INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: NodesConnE !< Elements that connect to a common node [-] LOGICAL :: SSSum = .false. !< SubDyn Summary File Flag [-] diff --git a/vs-build/c-bindings/AeroDyn_Inflow_C_Binding.vfproj b/vs-build/c-bindings/AeroDyn_Inflow_C_Binding.vfproj index 126c78c66f..b61ec0074c 100644 --- a/vs-build/c-bindings/AeroDyn_Inflow_C_Binding.vfproj +++ b/vs-build/c-bindings/AeroDyn_Inflow_C_Binding.vfproj @@ -74,23 +74,23 @@ - - + + + + + - + - + - + - - - - - + + diff --git a/vs-build/modules/NWTC-Library.vfproj b/vs-build/modules/NWTC-Library.vfproj index bf5d4164da..db76bb4258 100644 --- a/vs-build/modules/NWTC-Library.vfproj +++ b/vs-build/modules/NWTC-Library.vfproj @@ -115,12 +115,8 @@ - - - - - - + + @@ -147,70 +143,11 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + - - - - - - - - - - - - - - - - - - - - - - - - - + + - - - - - -