Merge rc-5.0.1 to dev (maintenance)

.github/scripts/check_tabs.sh

@@ -0,0 +1,104 @@
+#!/usr/bin/env bash
+#
+# check_tabs.sh -- verify that no Fortran source file contains tab characters.
+#
+# OpenFAST style requires 3-space indentation; tab characters are not allowed
+# because they render inconsistently across editors and break the alignment of
+# continuation lines.
+#
+# Scans every git-tracked *.f90 / *.F90 file, with no exclusions.  When a tab is
+# found in a file generated by the OpenFAST Registry, the report says so and
+# points at the registry input, since the fix belongs in the .txt rather than in
+# the generated Fortran.
+#
+# Exit status:
+#    0  no tab characters found
+#    1  one or more tab characters found
+#    2  the check could not run (not a git checkout, or a file failed to scan)
+#
+set -uo pipefail
+
+TAB=$(printf '\t')
+MARK='--->'
+
+repo_root=$(git rev-parse --show-toplevel) || {
+   echo "ERROR: not inside a git checkout; cannot determine repository root."
+   exit 2
+}
+cd "$repo_root" || exit 2
+
+echo "OpenFAST Fortran source style check"
+echo "==================================="
+echo
+echo "OpenFAST style requires 3 space indentation in all Fortran source files."
+echo "Tab characters are not allowed."
+echo
+
+mapfile -t FILES < <(git ls-files -- '*.f90' '*.F90' | sort)
+
+if [ "${#FILES[@]}" -eq 0 ]; then
+   echo "ERROR: no Fortran source files found -- is this a git checkout?"
+   exit 2
+fi
+
+n_bad_files=0
+n_bad_lines=0
+
+for file in "${FILES[@]}"; do
+   # grep exits 0 on a match, 1 on no match, and >1 on an actual error.
+   # Only "no match" may be skipped; a scan failure must not look like a pass.
+   matches=$(grep -n -- "$TAB" "$file")
+   status=$?
+   if [ "$status" -eq 1 ]; then
+      continue
+   elif [ "$status" -ne 0 ]; then
+      echo "ERROR: failed to scan ${file} (grep exit ${status})."
+      exit 2
+   fi
+
+   n_bad_files=$((n_bad_files + 1))
+   echo "$file"
+
+   while IFS= read -r match; do
+      lineno=${match%%:*}
+      content=${match#*:}
+
+      # Compute the column with parameter expansion rather than awk: "awk -v"
+      # decodes backslash escapes in the assigned value, which shifts the index
+      # on any line containing a backslash (Doxygen markup, format strings).
+      prefix=${content%%$TAB*}
+      column=$(( ${#prefix} + 1 ))
+      n_bad_lines=$((n_bad_lines + 1))
+
+      printf '   line %s, column %s: %s\n' "$lineno" "$column" "${content//$TAB/$MARK}"
+
+      # Inline annotation on the pull request diff.
+      if [ "${GITHUB_ACTIONS:-}" = "true" ]; then
+         echo "::error file=${file},line=${lineno},col=${column}::Tab character found. Use 3 space indentation; tabs are not allowed."
+      fi
+   done <<< "$matches"
+
+   # A tab in a Registry-generated file cannot be fixed in place: it was copied
+   # from a description string in the registry input, so point the developer
+   # there instead.
+   if head -n 1 "$file" 2>/dev/null | grep -q 'STARTOFREGISTRYGENERATEDFILE'; then
+      echo "   NOTE: this file is generated by the OpenFAST Registry. Fix the tab in"
+      echo "         the corresponding registry input (.txt) in the same directory,"
+      echo "         then regenerate this file."
+   fi
+
+   echo
+done
+
+echo "Checked ${#FILES[@]} Fortran source file(s)."
+
+if [ "$n_bad_files" -gt 0 ]; then
+   echo
+   echo "ERROR: found ${n_bad_lines} line(s) containing tab characters in ${n_bad_files} file(s)."
+   echo "       Replace each tab with spaces (3 space indentation) and commit the fix."
+   echo "       Tabs above are shown as '${MARK}'."
+   exit 1
+fi
+
+echo "No tab characters found."
+exit 0

.github/workflows/check-tabs.yml

@@ -0,0 +1,36 @@
+
+name: 'Fortran Source Style'
+
+#-------------------------------------------------------------------------------
+# Notes
+# - OpenFAST requires 3 space indentation in Fortran source; tab characters are
+#   not allowed.  This workflow enforces that rule on all git-tracked *.f90 and
+#   *.F90 files.
+# - No files are excluded from the scan.  A tab in a Registry-generated file is
+#   reported with a note pointing at the registry .txt input, since that is
+#   where it has to be fixed.
+#-------------------------------------------------------------------------------
+
+on:
+  push:
+    paths:
+      - '**.f90'
+      - '**.F90'
+      - '.github/workflows/check-tabs.yml'
+      - '.github/scripts/check_tabs.sh'
+
+  pull_request:
+    types: [opened, synchronize, edited, reopened]
+
+jobs:
+  check-tabs:
+    name: Check for tab characters
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v4
+        with:
+          submodules: false
+
+      - name: Check Fortran source for tab characters
+        run: .github/scripts/check_tabs.sh

cmake/OpenfastFortranOptions.cmake

@@ -133,10 +133,10 @@ macro(set_fast_gfortran)
     set( CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -fcheck=all,no-array-temps -pedantic -fbacktrace -finit-real=inf -finit-integer=9999." )
   endif()
 
-  if(CYGWIN)
-    # increase the default 2MB stack size to 16 MB
+  if(CYGWIN OR MINGW)
+    # increase the default 1-2 MB stack size to 16 MB
     MATH(EXPR stack_size "16 * 1024 * 1024")
-    set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS},--stack,${stack_size}")
+    set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -Wl,--stack,${stack_size}")
   endif()
 
   # Profiling

modules/aerodyn/src/AeroDyn.f90

@@ -6119,8 +6119,9 @@ SUBROUTINE TwrInflArray( p, u, RotInflow, m, Positions, Inflow, ErrStat, ErrMsg
    ! these models are valid for only small tower deflections; check for potential division-by-zero errors:   
    call CheckTwrInfl( u, ErrStat2, ErrMsg2 ); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ); if (ErrStat >= AbortErrLev) return
 
+   ! FirstWarn_TowerStrike is firstprivate so each thread starts from the .false. set above (avoids reading uninitialized memory); ErrStat2/ErrMsg2 are private to avoid a data race
    !$OMP PARALLEL default(shared)
-   !$OMP do private(i,Pos,theta_tower_trans,W_tower,xbar,ybar,zbar,TwrCd,TwrTI,TwrClrnc,FirstWarn_TowerStrike,DisturbInflow,v) schedule(runtime)
+   !$OMP do private(i,Pos,theta_tower_trans,W_tower,xbar,ybar,zbar,TwrCd,TwrTI,TwrClrnc,DisturbInflow,v,ErrStat2,ErrMsg2) firstprivate(FirstWarn_TowerStrike) schedule(runtime)
    do i = 1, size(Positions,2)
       Pos=Positions(1:3,i)
 
@@ -7737,13 +7738,25 @@ subroutine PerturbFlowField(Var, BaseFF, PerturbSign, PerturbFF)
       PerturbFF%Uniform%VelH = BaseFF%Uniform%VelH
       PerturbFF%Uniform%ShrV = BaseFF%Uniform%ShrV
       PerturbFF%PropagationDir = BaseFF%PropagationDir
+      PerturbFF%RotToWind = BaseFF%RotToWind
+      PerturbFF%RotFromWind = BaseFF%RotFromWind
+      PerturbFF%RotateWindBox = BaseFF%RotateWindBox
       select case (Var%DL%Num)
       case (AD_u_HWindSpeed) 
          PerturbFF%Uniform%VelH = BaseFF%Uniform%VelH + Var%Perturb*PerturbSign
       case (AD_u_PLexp) 
          PerturbFF%Uniform%ShrV = BaseFF%Uniform%ShrV + Var%Perturb*PerturbSign
       case (AD_u_PropagationDir) 
          PerturbFF%PropagationDir = BaseFF%PropagationDir + Var%Perturb*PerturbSign
+         PerturbFF%RotToWind(1,:) = [ &
+            cos(-PerturbFF%VFlowAngle)*cos(-PerturbFF%PropagationDir), &
+            cos(-PerturbFF%VFlowAngle)*sin(-PerturbFF%PropagationDir), -sin(-PerturbFF%VFlowAngle)]
+         PerturbFF%RotToWind(2,:) = [-sin(-PerturbFF%PropagationDir), cos(-PerturbFF%PropagationDir), 0.0_ReKi]
+         PerturbFF%RotToWind(3,:) = [ &
+            sin(-PerturbFF%VFlowAngle)*cos(-PerturbFF%PropagationDir), &
+            sin(-PerturbFF%VFlowAngle)*sin(-PerturbFF%PropagationDir), cos(-PerturbFF%VFlowAngle)]
+         PerturbFF%RotFromWind = transpose(PerturbFF%RotToWind)
+         PerturbFF%RotateWindBox = .true.
       end select
    end subroutine
 

modules/aerodyn/src/FVW.f90

@@ -14,7 +14,7 @@ module FVW
    use FVW_IO
    use FVW_Wings
    use FVW_BiotSavart
-   use FVW_VortexTools, only: tic, toc
+   use FVW_VortexTools, only: tic, toc, T_Tree
    use FVW_Tests
    use AirFoilInfo
 
@@ -30,9 +30,6 @@ module FVW
    public   :: FVW_CalcOutput
    public   :: FVW_UpdateStates
 
-   ! parameter for deciding if enough time has elapsed (Wake calculation, and vtk output)
-   real(DbKi), parameter      :: OneMinusEpsilon = 1 - 10000*EPSILON(1.0_DbKi)
-
 contains
 
 !----------------------------------------------------------------------------------------------------------------------------------
@@ -603,7 +600,7 @@ subroutine FVW_UpdateStates( t, n, u, utimes, p, x, xd, z, OtherState, AFInfo, m
       bReevaluation=.True.
    endif
    ! Compute Induced wake effects only if time since last compute is > DTfvw
-   if ( (( t - m%OldWakeTime ) >= p%DTfvw*OneMinusEpsilon) )  then
+   if ( (( t - m%OldWakeTime ) >= p%DTfvw - 0.25_DbKi*p%DTaero) )  then
       m%OldWakeTime = t
       m%ComputeWakeInduced = .TRUE.    ! It's time to update the induced velocities from wake
    else
@@ -1570,6 +1567,9 @@ subroutine WriteVTKOutputs(t, force, VTKstep, u, p, x, z, m, ErrStat, ErrMsg)
    character(*), parameter :: RoutineName = 'WriteVTKOutputs'
    integer(IntKi) :: iW, iGrid
    integer(IntKi) :: nSeg, nSegP
+   logical, allocatable :: bDoGrid(:) ! Flag per output grid: .true. if it needs to be written now
+   type(T_Tree)   :: Tree
+   type(T_Panl)   :: Panl
    ErrStat = ErrID_None
    ErrMsg  = ''
    if (OLAF_PROFILING) call tic('WriteVTKOutputs')
@@ -1584,7 +1584,7 @@ subroutine WriteVTKOutputs(t, force, VTKstep, u, p, x, z, m, ErrStat, ErrMsg)
       do iW=1,p%nWings
          m%W(iW)%Vtot_CP = m%W(iW)%Vind_CP + m%W(iW)%Vwnd_CP - m%W(iW)%Vstr_CP
       enddo
-      if ( force .or. (( t - m%VTKlastTime ) >= p%DTvtk*OneMinusEpsilon ))  then
+      if ( force .or. (( t - m%VTKlastTime ) >= p%DTvtk - 0.25_DbKi*p%DTaero ))  then
          m%VTKlastTime = t
          if ((p%VTKCoord==2).or.(p%VTKCoord==3)) then
             ! Hub reference coordinates, for export only, ALL VTK Will be exported in this coordinate system!
@@ -1612,18 +1612,31 @@ subroutine WriteVTKOutputs(t, force, VTKstep, u, p, x, z, m, ErrStat, ErrMsg)
       ! Distribute the Wind we requested to Inflow wind to storage Misc arrays
       ! TODO ANDY: replace with direct call to inflow wind at Grid points
       CALL DistributeRequestedWind_Grid(u%V_wind, p, m)
+      ! Compute once which grids are due for output now
+      allocate(bDoGrid(p%nGridOut))
       do iGrid=1,p%nGridOut
          bWithinTime   = t>=m%GridOutputs(iGrid)%tStart-p%DTaero/2. .and. t<= m%GridOutputs(iGrid)%tEnd+p%DTaero/2.
-         bTimeToOutput = ( t - m%GridOutputs(iGrid)%tLastOutput) >= m%GridOutputs(iGrid)%DTout * OneMinusEpsilon
-         if (force .or. (bWithinTime .and. bTimeToOutput) )  then
-            ! Compute induced velocity on grid, TODO use the same Tree for all CalcOutput
-            call InducedVelocitiesAll_OnGrid(m%GridOutputs(iGrid), p, x, m, ErrStat2, ErrMsg2);
-            call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
-            m%GridOutputs(iGrid)%tLastOutput = t
-            call WrVTK_FVW_Grid(p, m, iGrid, trim(p%VTK_OutFileBase)//'FVW_Grid', VTKstep, 9)
-            m%VTKstep=VTKstep ! We save the step at which writing occurred
-         endif
+         bTimeToOutput = ( t - m%GridOutputs(iGrid)%tLastOutput) >= m%GridOutputs(iGrid)%DTout - 0.25_DbKi*p%DTaero
+         bDoGrid(iGrid) = force .or. (bWithinTime .and. bTimeToOutput)
       enddo
+      if (any(bDoGrid)) then
+         ! Build the wake segments/tree once and reuse it for all grids below
+         call InducedVelocitiesAll_Init(p, x, m, m%Sgmt, m%Part, Tree, Panl, ErrStat2, ErrMsg2, allocPart=.false.)
+         call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
+         do iGrid=1,p%nGridOut
+            if (bDoGrid(iGrid)) then
+               ! Compute induced velocity on grid, reusing the wake tree built once
+               call InducedVelocitiesAll_OnGrid_Calc(m%GridOutputs(iGrid), p, m%Sgmt, m%Part, Tree, Panl, ErrStat2, ErrMsg2)
+               call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
+               m%GridOutputs(iGrid)%tLastOutput = t
+               call WrVTK_FVW_Grid(p, m, iGrid, trim(p%VTK_OutFileBase)//'FVW_Grid', VTKstep, 9)
+               m%VTKstep=VTKstep ! We save the step at which writing occurred
+            endif
+         enddo
+         call InducedVelocitiesAll_End(p, Tree, m%Part, Panl, ErrStat2, ErrMsg2, deallocPart=.false.)
+         call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
+      endif
+      deallocate(bDoGrid)
    endif
    if (OLAF_PROFILING) call toc()
 end subroutine WriteVTKOutputs

modules/aerodyn/src/FVW_BiotSavart.f90

@@ -10,6 +10,8 @@ module FVW_BiotSavart
    real(ReKi),parameter :: PRECISION_UI  = epsilon(1.0_ReKi)/100 !< NOTE assuming problem of size 1
    real(ReKi),parameter :: PRECISION_EPS =  epsilon(1.0_ReKi) !< Machine Precision For the given ReKi for problems of scale 1!
    real(ReKi),parameter :: MIN_EXP_VALUE=-10.0_ReKi
+   real(ReKi),parameter :: PART_REG_NRAD =  2.0_ReKi             !< Particle exp mollifier treated as 1 beyond this many core radii (matches 2*rc far-field multipole floor)
+   real(ReKi),parameter :: PART_REG_CUT3 =  PART_REG_NRAD**3     !< Corresponding (r/rc)^3 cutoff
    real(ReKi),parameter :: MINDENOM=0.0_ReKi
 !    real(ReKi),parameter :: MINDENOM=1e-15_ReKi
    real(ReKi),parameter :: MINNORM=1e-4
@@ -374,26 +376,35 @@ subroutine ui_part_nograd_11(DeltaP, Alpha, RegFunction, RegParam, Ui)
    real(ReKi),dimension(3) :: C          !< Cross product of Alpha and r
    real(ReKi)              :: E          !< Exponential poart for the mollifider
    real(ReKi)              :: r3_inv     !< 
+   real(ReKi)              :: r2, r3, rc3!< |r|^2, |r|^3, RegParam^3 (reused to avoid recomputing ** intrinsics)
    real(ReKi)              :: rDeltaP    !< norm , distance between point and particle
    real(ReKi)              :: ScalarPart !< the part containing the inverse of the distance, but not 4pi, Mollifier
-   rDeltaP=sqrt(DeltaP(1)**2+ DeltaP(2)**2+ DeltaP(3)**2)! norm
+   r2      = DeltaP(1)**2+ DeltaP(2)**2+ DeltaP(3)**2
+   rDeltaP = sqrt(r2)! norm
    if (rDeltaP<MINNORM) then !--- Exactly on the Singularity 
       Ui(1:3)  = 0.0_ReKi
       return
    else !--- Normal Procedure 
+      r3 = r2*rDeltaP ! |r|^3, reused below
       C(1) = Alpha(2) * DeltaP(3) - Alpha(3) * DeltaP(2)
       C(2) = Alpha(3) * DeltaP(1) - Alpha(1) * DeltaP(3)
       C(3) = Alpha(1) * DeltaP(2) - Alpha(2) * DeltaP(1)
       select case (RegFunction) !
       case (idRegNone) ! No mollification
-         r3_inv     = 1._ReKi/(rDeltaP**3)
+         r3_inv     = 1._ReKi/r3
          ScalarPart = r3_inv*fourpi_inv
       case (idRegExp) ! Exponential mollifier
-         r3_inv     = 1._ReKi/(rDeltaP**3)
-         E          = exp(-rDeltaP**3/RegParam**3)
-         ScalarPart = (1._ReKi-E)*r3_inv*fourpi_inv
+         rc3        = RegParam*RegParam*RegParam
+         r3_inv     = 1._ReKi/r3
+         if (r3 > PART_REG_CUT3*rc3) then ! r > 2*rc: mollifier -> 1 (skip exp), consistent with far-field multipole floor
+            ScalarPart = r3_inv*fourpi_inv
+         else
+            E          = exp(-r3/rc3)
+            ScalarPart = (1._ReKi-E)*r3_inv*fourpi_inv
+         endif
       case (idRegCompact) ! Compact support
-         r3_inv     = 1._ReKi/sqrt(RegParam**6+rDeltaP**6)
+         rc3        = RegParam*RegParam*RegParam
+         r3_inv     = 1._ReKi/sqrt(rc3*rc3+r3*r3)
          ScalarPart = r3_inv*fourpi_inv
       case default 
          print*,'[ERROR] Wrong regularization function for particles',RegFunction

modules/aerodyn/src/FVW_IO.f90

@@ -184,6 +184,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg )
 
    ! --- Validation of inputs
    if (PathIsRelative(Inp%CirculationFile)) Inp%CirculationFile = TRIM(PriPath)//TRIM(Inp%CirculationFile)
+   if (len_trim(Inp%SrcPnlFile)>0 .and. PathIsRelative(Inp%SrcPnlFile)) Inp%SrcPnlFile = TRIM(PriPath)//TRIM(Inp%SrcPnlFile)
 
    if (Check(.not.(ANY(idCircVALID ==Inp%CircSolvMethod)), 'Circulation method (CircSolvMethod) not implemented: '//trim(Num2LStr(Inp%CircSolvMethod)))) return
    if (Check(.not.(ANY(idIntMethodVALID==Inp%IntMethod    )) , 'Time integration method (IntMethod) not yet implemented. Use Euler 1st order method for now.')) return