From 02d93217803adc82625f681c6929b6452fb004c5 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 08:12:48 -0600 Subject: [PATCH 01/13] AeroDyn/OLAF Registry: add xPts/yPts/zPts and list-file fields to GridOutType --- modules/aerodyn/src/FVW_Registry.txt | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/modules/aerodyn/src/FVW_Registry.txt b/modules/aerodyn/src/FVW_Registry.txt index 45a7aff474..2d9ce90812 100644 --- a/modules/aerodyn/src/FVW_Registry.txt +++ b/modules/aerodyn/src/FVW_Registry.txt @@ -28,6 +28,12 @@ typedef ^ ^ IntKi typedef ^ ^ ReKi uGrid {:}{:}{:}{:} - - "Grid velocity 3 x nz x ny x nx" - typedef ^ ^ ReKi omGrid {:}{:}{:}{:} - - "Grid vorticity 3 x nz x ny x nx" - typedef ^ ^ DbKi tLastOutput - - - "Last output time" - +typedef ^ ^ ReKi xPts {:} - - "Explicit x coordinates (non-equidistant grid, if used)" m +typedef ^ ^ ReKi yPts {:} - - "Explicit y coordinates (non-equidistant grid, if used)" m +typedef ^ ^ ReKi zPts {:} - - "Explicit z coordinates (non-equidistant grid, if used)" m +typedef ^ ^ CHARACTER(1024) xListFile - - - "File with explicit x coordinates (empty if equidistant)" - +typedef ^ ^ CHARACTER(1024) yListFile - - - "File with explicit y coordinates (empty if equidistant)" - +typedef ^ ^ CHARACTER(1024) zListFile - - - "File with explicit z coordinates (empty if equidistant)" - ##################### Segments ############### typedef FVW/FVW T_Sgmt ReKi Points :: - - "Points delimiting the segments" - From 1463748db34cc0c1b737fab0604643847ce48815 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 08:30:55 -0600 Subject: [PATCH 02/13] AeroDyn/OLAF: parse a number or quoted filename in ReadGridOut --- modules/aerodyn/src/FVW_IO.f90 | 43 ++++++++++++++++++++++++++-------- 1 file changed, 33 insertions(+), 10 deletions(-) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 4cde5064d4..014d36f607 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -327,19 +327,42 @@ subroutine ReadGridOut(sLine, GridOut) return endif endif - ! x,y,z + ! x ErrMsg2='Error reading OLAF "x" inputs for grid outputs line: '//trim(sLine) - if (.not. is_numeric(StrArray( 6), GridOut%xStart) ) return - if (.not. is_numeric(StrArray( 7), GridOut%xEnd ) ) return - if (.not. is_integer(StrArray( 8), GridOut%nx ) ) return + GridOut%xListFile = '' + if ( is_numeric(StrArray(6), GridOut%xStart) ) then + if (.not. is_numeric(StrArray(7), GridOut%xEnd) ) return + if (.not. is_integer(StrArray(8), GridOut%nx ) ) return + else + GridOut%xListFile = StrArray(6) + GridOut%xStart = 0.0_ReKi + GridOut%xEnd = 0.0_ReKi + GridOut%nx = -1 + endif + ! y ErrMsg2='Error reading OLAF "y" inputs for grid outputs line: '//trim(sLine) - if (.not. is_numeric(StrArray( 9), GridOut%yStart) ) return - if (.not. is_numeric(StrArray(10), GridOut%yEnd ) ) return - if (.not. is_integer(StrArray(11), GridOut%ny ) ) return + GridOut%yListFile = '' + if ( is_numeric(StrArray(9), GridOut%yStart) ) then + if (.not. is_numeric(StrArray(10), GridOut%yEnd) ) return + if (.not. is_integer(StrArray(11), GridOut%ny ) ) return + else + GridOut%yListFile = StrArray(9) + GridOut%yStart = 0.0_ReKi + GridOut%yEnd = 0.0_ReKi + GridOut%ny = -1 + endif + ! z ErrMsg2='Error reading OLAF "z" inputs for grid outputs line: '//trim(sLine) - if (.not. is_numeric(StrArray(12), GridOut%zStart) ) return - if (.not. is_numeric(StrArray(13), GridOut%zEnd ) ) return - if (.not. is_integer(StrArray(14), GridOut%nz ) ) return + GridOut%zListFile = '' + if ( is_numeric(StrArray(12), GridOut%zStart) ) then + if (.not. is_numeric(StrArray(13), GridOut%zEnd) ) return + if (.not. is_integer(StrArray(14), GridOut%nz ) ) return + else + GridOut%zListFile = StrArray(12) + GridOut%zStart = 0.0_ReKi + GridOut%zEnd = 0.0_ReKi + GridOut%nz = -1 + endif ! Success ErrStat2=ErrID_None ErrMsg2='' From 1ae78909c2d969110bfcca9c2ac94aed21185cd4 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 08:55:38 -0600 Subject: [PATCH 03/13] AeroDyn/OLAF: call ResolveGridAxis in FVW_ReadInputFile and require GridType=1 for list-based grids --- modules/aerodyn/src/FVW_IO.f90 | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 014d36f607..5a16dea07e 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -100,13 +100,28 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) allocate(m%GridOutputs(p%nGridOut), stat=ErrStat2); CALL ReadCom (UnIn,FileName, 'GridOutHeaders', ErrStat2,ErrMsg2); if(Failed()) return CALL ReadCom (UnIn,FileName, 'GridOutUnits', ErrStat2,ErrMsg2); if(Failed()) return - do i =1, p%nGridOut + do i =1, p%nGridOut ErrMsg2='Error reading OLAF grid outputs line '//trim(num2lstr(i)) read(UnIn, fmt='(A)', iostat=ErrStat2) sLine ; if(Failed()) return call ReadGridOut(sLine, m%GridOutputs(i)); if(Failed()) return + ! Resolve each axis to an explicit coordinate array (regardless of whether it is a list file or a range) + call ResolveGridAxis(m%GridOutputs(i)%xStart, m%GridOutputs(i)%xEnd, m%GridOutputs(i)%nx, & + m%GridOutputs(i)%xListFile, PriPath, m%GridOutputs(i)%xPts, ErrStat2, ErrMsg2); if(Failed()) return + call ResolveGridAxis(m%GridOutputs(i)%yStart, m%GridOutputs(i)%yEnd, m%GridOutputs(i)%ny, & + m%GridOutputs(i)%yListFile, PriPath, m%GridOutputs(i)%yPts, ErrStat2, ErrMsg2); if(Failed()) return + call ResolveGridAxis(m%GridOutputs(i)%zStart, m%GridOutputs(i)%zEnd, m%GridOutputs(i)%nz, & + m%GridOutputs(i)%zListFile, PriPath, m%GridOutputs(i)%zPts, ErrStat2, ErrMsg2); if(Failed()) return + ! Error checking if (Check(m%GridOutputs(i)%nx<1, 'Grid output nx needs to be >=1')) return if (Check(m%GridOutputs(i)%ny<1, 'Grid output ny needs to be >=1')) return if (Check(m%GridOutputs(i)%nz<1, 'Grid output nz needs to be >=1')) return + ! Vorticity requires equidistant spacing. GridType must be 1 when using list files. + if (m%GridOutputs(i)%type==idGridVelVorticity) then + if (Check( len_trim(m%GridOutputs(i)%xListFile)>0 .or. & + len_trim(m%GridOutputs(i)%yListFile)>0 .or. & + len_trim(m%GridOutputs(i)%zListFile)>0, & + 'Grid "'//trim(m%GridOutputs(i)%name)//'": vorticity requires equidistant spacing. nGridOut must be 1 when using list files.')) return + endif enddo endif From aa525a00b2e84894a134691fb613076bb9a78aa3 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 09:19:43 -0600 Subject: [PATCH 04/13] AeroDyn/OLAF: New ResolveGridAxis for list-based or range --- modules/aerodyn/src/FVW_IO.f90 | 67 +++++++++++++++++++++++++++++++++- 1 file changed, 66 insertions(+), 1 deletion(-) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 5a16dea07e..756db5686d 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -120,7 +120,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) if (Check( len_trim(m%GridOutputs(i)%xListFile)>0 .or. & len_trim(m%GridOutputs(i)%yListFile)>0 .or. & len_trim(m%GridOutputs(i)%zListFile)>0, & - 'Grid "'//trim(m%GridOutputs(i)%name)//'": vorticity requires equidistant spacing. nGridOut must be 1 when using list files.')) return + 'Grid "'//trim(m%GridOutputs(i)%name)//'": vorticity requires equidistant spacing. GridType must be 1 when using list files.')) return endif enddo endif @@ -383,6 +383,71 @@ subroutine ReadGridOut(sLine, GridOut) ErrMsg2='' end subroutine ReadGridOut + ! Resolve one grid axis to an explicit, strictly increasing coordinate array, + ! either read from ListFile (if given) or built from an equidistant Start/End/n range. + subroutine ResolveGridAxis(AStart, AEnd, n, ListFile, Pts) + real(ReKi), intent(in) :: AStart, AEnd + integer(IntKi), intent(inout) :: n + character(*), intent(in) :: ListFile + real(ReKi), allocatable, intent(out) :: Pts(:) + character(1024) :: FullFile + integer(IntKi) :: UnList, j, IOS + real(ReKi) :: val + + if (len_trim(ListFile) == 0) then + ! Equidistant points, expressed as an explicit array + allocate(Pts(max(n,1))) + do j = 1, n + Pts(j) = AStart + (AEnd - AStart) * real(j-1,ReKi) / real(max(n-1,1),ReKi) + enddo + ErrStat2 = ErrID_None; ErrMsg2 = '' + return + endif + + ! Explicit list from file + FullFile = ListFile + if (PathIsRelative(FullFile)) FullFile = trim(PriPath)//trim(FullFile) + call GetNewUnit(UnList) + call OpenFInpFile(UnList, trim(FullFile), ErrStat2, ErrMsg2) + if (ErrStat2 /= ErrID_None) return + + ! Count valid lines + n = 0 + do + read(UnList, *, iostat=IOS) val + if (IOS /= 0) exit + n = n + 1 + enddo + if (n < 1) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = 'Grid point list file "'//trim(FullFile)//'" contains no valid values.' + close(UnList) + return + endif + + ! Read values + rewind(UnList) + allocate(Pts(n)) + do j = 1, n + read(UnList, *) Pts(j) + enddo + close(UnList) + + ! Validate strictly increasing, no duplicates + do j = 2, n + if (Pts(j) <= Pts(j-1)) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = 'Grid point list file "'//trim(FullFile)//'" must be strictly increasing (no duplicates); '// & + 'value at line '//trim(Num2LStr(j))//' ('//trim(Num2LStr(Pts(j)))//') is not greater than '// & + 'the previous value ('//trim(Num2LStr(Pts(j-1)))//').' + return + endif + enddo + + ErrStat2 = ErrID_None + ErrMsg2 = '' + end subroutine ResolveGridAxis + END SUBROUTINE FVW_ReadInputFile From 11a716246fc28e97a1f3b309c261fb9d669a006f Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 09:34:58 -0600 Subject: [PATCH 05/13] AeroDyn/OLAF: export grid velocity as RectilinearGrid (list-based) or StructuredPoints (range) --- modules/aerodyn/src/FVW_IO.f90 | 33 ++++++++++++++++++++++----------- 1 file changed, 22 insertions(+), 11 deletions(-) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 756db5686d..b05f5007a0 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -600,6 +600,7 @@ subroutine WrVTK_FVW_Grid(p, m, iGrid, FileRootName, VTKcount, Twidth, HubOrient character(255) :: Label character(Twidth) :: Tstr ! string for current VTK write-out step (padded with zeros) real(ReKi), dimension(3) :: dx + logical :: bListBased type(GridOutType), pointer :: g type(VTK_Misc) :: mvtk @@ -613,19 +614,29 @@ subroutine WrVTK_FVW_Grid(p, m, iGrid, FileRootName, VTKcount, Twidth, HubOrient g => m%GridOutputs(iGrid) Label=trim(g%name) Filename = TRIM(FileRootName)//'.'//trim(Label)//'.'//Tstr//'.vtk' + bListBased = len_trim(g%xListFile)>0 .or. len_trim(g%yListFile)>0 .or. len_trim(g%zListFile)>0 if ( vtk_new_ascii_file(trim(filename),Label,mvtk) ) then - dx(1) = (g%xEnd- g%xStart)/max(g%nx-1,1) - dx(2) = (g%yEnd- g%yStart)/max(g%ny-1,1) - dx(3) = (g%zEnd- g%zStart)/max(g%nz-1,1) - call vtk_dataset_structured_points((/g%xStart, g%yStart, g%zStart/),dx,(/g%nx,g%ny,g%nz/),mvtk) - call vtk_point_data_init(mvtk) - call vtk_point_data_vector(g%uGrid(1:3,:,:,:),'Velocity',mvtk) - ! Compute vorticity on the fly - if (g%type==idGridVelVorticity) then - call curl_regular_grid(g%uGrid, g%omgrid, 1,1,1, g%nx,g%ny,g%nz, dx(1),dx(2),dx(3)) - call vtk_point_data_vector(g%omGrid(1:3,:,:,:),'Vorticity',mvtk) + if (bListBased) then + ! List-based grid + ! At least one axis is non-equidistant: RectilinearGrid (explicit coordinates) + ! NOTE: vorticity is blocked for this case at read time. GridType = 1 (velocity only). + call vtk_dataset_rectilinear(g%xPts, g%yPts, g%zPts, mvtk) + call vtk_point_data_init(mvtk) + call vtk_point_data_vector(g%uGrid(1:3,:,:,:),'Velocity',mvtk) + else + ! Equidistant grid: StructuredPoints (implicit coordinates) + dx(1) = (g%xEnd- g%xStart)/max(g%nx-1,1) + dx(2) = (g%yEnd- g%yStart)/max(g%ny-1,1) + dx(3) = (g%zEnd- g%zStart)/max(g%nz-1,1) + call vtk_dataset_structured_points((/g%xStart, g%yStart, g%zStart/),dx,(/g%nx,g%ny,g%nz/),mvtk) + call vtk_point_data_init(mvtk) + call vtk_point_data_vector(g%uGrid(1:3,:,:,:),'Velocity',mvtk) + ! Compute vorticity on the fly + if (g%type==idGridVelVorticity) then + call curl_regular_grid(g%uGrid, g%omgrid, 1,1,1, g%nx,g%ny,g%nz, dx(1),dx(2),dx(3)) + call vtk_point_data_vector(g%omGrid(1:3,:,:,:),'Vorticity',mvtk) + endif endif - ! call vtk_close_file(mvtk) endif From 231c6eab9b94f7bbeff53702b638f931fa024291 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 09:38:54 -0600 Subject: [PATCH 06/13] FVW_IO.f90: trim trailing whitespace --- modules/aerodyn/src/FVW_IO.f90 | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index b05f5007a0..1ea8dad016 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -18,7 +18,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None ! Local variables character(1024) :: PriPath ! the path to the primary input file - character(1024) :: sDummy, sLine ! string to temporarially hold value of read line + character(1024) :: sDummy, sLine ! string to temporarially hold value of read line integer(IntKi) :: UnIn, i integer(IntKi) :: ErrStat2 character(ErrMsgLen) :: ErrMsg2 @@ -26,7 +26,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) ErrMsg = "" Inp%SrcPnlFile = '' ! TODO registry init for empty strings ! Open file - CALL GetNewUnit( UnIn ) + CALL GetNewUnit( UnIn ) CALL OpenFInpfile(UnIn, TRIM(FileName), ErrStat2, ErrMsg2) if (Check( ErrStat2 /= ErrID_None , 'Could not open input file')) return CALL GetPath( FileName, PriPath ) ! Input files will be relative to the path where the primary input file is located. @@ -128,7 +128,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) ! --- Advanced Options ! NOTE: no error handling since this is for debug ! Default options are typically "true" - CALL ReadCom(UnIn,FileName, '=== Separator' ,ErrStat2,ErrMsg2); + CALL ReadCom(UnIn,FileName, '=== Separator' ,ErrStat2,ErrMsg2); CALL ReadCom(UnIn,FileName, '--- Advanced options header' ,ErrStat2,ErrMsg2); if(ErrStat2==ErrID_None) then call WrScr(' - Reading advanced options for OLAF:') @@ -218,7 +218,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) if (Check(Inp%WingRegParam<0 , 'Wing regularization parameter (WakeRegParam) should be positive')) return if (Check(Inp%CoreSpreadEddyVisc<0 , 'Core spreading eddy viscosity (CoreSpreadEddyVisc) should be positive')) return - ! Removing the shed vorticity is a dangerous option if this is done too close to the blades. + ! Removing the shed vorticity is a dangerous option if this is done too close to the blades. ! To be safe, we will no matter what ensure that the last segments of NW are 0 if FWShedVorticity is False (see PackPanelsToSegments) ! Still we force the user to be responsible. if (Check((.not.(Inp%FWShedVorticity)) .and. Inp%nNWPanels<30, '`FWShedVorticity` should be true if `nNWPanels`<30. Alternatively, use a larger number of NWPanels ')) return @@ -231,7 +231,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) CONTAINS logical function Failed() - call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'FVW_ReadInputFile') + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'FVW_ReadInputFile') Failed = ErrStat >= AbortErrLev if (Failed) call CleanUp() end function Failed @@ -304,7 +304,7 @@ subroutine ReadGridOut(sLine, GridOut) ErrStat2=ErrID_Fatal ErrMsg2='Error reading OLAF grid outputs line: '//trim(sLine) ! Name - GridOut%name =StrArray(1) + GridOut%name =StrArray(1) ! Type if (.not. is_integer(StrArray(2), GridOut%type ) ) then ErrMsg2=trim(ErrMsg2)//NewLine//'GridType needs to be an integer.' @@ -315,7 +315,7 @@ subroutine ReadGridOut(sLine, GridOut) if ( index(StrArray(3), "DEFAULT" ) == 1 ) then GridOut%tStart = 0.0_ReKi else - if (.not. is_numeric(StrArray(3), GridOut%tStart) ) then + if (.not. is_numeric(StrArray(3), GridOut%tStart) ) then ErrMsg2=trim(ErrMsg2)//NewLine//'TStart needs to be numeric or "default".' return endif @@ -490,7 +490,7 @@ subroutine WrVTK_FVW(p, x, z, m, FileRootName, VTKcount, Twidth, bladeFrame, Hub Call ProgAbort('Programming error in WrVTK_FVW call: Cannot use the WrVTK_FVW with bladeFrame==TRUE without the optional arguments of HubOrientation and HubPosition') endif endif - + if (DEV_VERSION) then print*,'------------------------------------------------------------------------------' print'(A,L1,A,I0,A,I0,A,I0)','VTK Output - First call ',m%FirstCall, ' nNW:',m%nNW,' nFW:',m%nFW,' i:',VTKCount @@ -502,7 +502,7 @@ subroutine WrVTK_FVW(p, x, z, m, FileRootName, VTKcount, Twidth, bladeFrame, Hub write(Tstr, '(i' // trim(Num2LStr(Twidth)) //'.'// trim(Num2LStr(Twidth)) // ')') VTKcount ! --------------------------------------------------------------------------------} - ! --- Blade + ! --- Blade ! --------------------------------------------------------------------------------{ ! --- Blade Quarter chord points (AC) do iW=1,p%VTKBlades @@ -530,7 +530,7 @@ subroutine WrVTK_FVW(p, x, z, m, FileRootName, VTKcount, Twidth, bladeFrame, Hub ! call WrVTK_Lattice(FileName, mvtk, m%W(iW)%r_LL(1:3,:,:), m%W(iW)%Gamma_LL(:), bladeFrame=bladeFrame) ! enddo ! --------------------------------------------------------------------------------} - ! --- Near wake + ! --- Near wake ! --------------------------------------------------------------------------------{ ! --- Near wake panels do iW=1,p%VTKBlades @@ -548,7 +548,7 @@ subroutine WrVTK_FVW(p, x, z, m, FileRootName, VTKcount, Twidth, bladeFrame, Hub endif enddo ! --------------------------------------------------------------------------------} - ! --- Far wake + ! --- Far wake ! --------------------------------------------------------------------------------{ ! --- Far wake panels do iW=1,p%VTKBlades @@ -572,7 +572,7 @@ subroutine WrVTK_FVW(p, x, z, m, FileRootName, VTKcount, Twidth, bladeFrame, Hub endif if (nSeg>0) then Filename = TRIM(FileRootName)//'.AllSeg.'//Tstr//'.vtk' - CALL WrVTK_Segments(Filename, mvtk, m%Sgmt%Points(:,1:nSegP), m%Sgmt%Connct(:,1:nSeg), m%Sgmt%Gamma(1:nSeg), m%Sgmt%Epsilon(1:nSeg), bladeFrame) + CALL WrVTK_Segments(Filename, mvtk, m%Sgmt%Points(:,1:nSegP), m%Sgmt%Connct(:,1:nSeg), m%Sgmt%Gamma(1:nSeg), m%Sgmt%Epsilon(1:nSeg), bladeFrame) endif if (p%SrcPnl%n>0) then @@ -683,14 +683,14 @@ subroutine WrVTK_Panels(filename, mvtk, p, m, z) endsubroutine WrVTK_Panels -subroutine WrVTK_Segments(filename, mvtk, SegPoints, SegConnct, SegGamma, SegEpsilon, bladeFrame) +subroutine WrVTK_Segments(filename, mvtk, SegPoints, SegConnct, SegGamma, SegEpsilon, bladeFrame) use VTK character(len=*),intent(in) :: filename type(VTK_Misc), intent(inout) :: mvtk !< miscvars for VTK output - real(ReKi), dimension(:,:), intent(in) :: SegPoints !< - integer(IntKi), dimension(:,:), intent(in) :: SegConnct !< - real(ReKi), dimension(:) , intent(in) :: SegGamma !< - real(ReKi), dimension(:) , intent(in) :: SegEpsilon !< + real(ReKi), dimension(:,:), intent(in) :: SegPoints !< + integer(IntKi), dimension(:,:), intent(in) :: SegConnct !< + real(ReKi), dimension(:) , intent(in) :: SegGamma !< + real(ReKi), dimension(:) , intent(in) :: SegEpsilon !< logical, intent(in ) :: bladeFrame !< Output in blade coordinate frame if ( vtk_new_ascii_file(filename,'Sgmt',mvtk) ) then call vtk_dataset_polydata(SegPoints(1:3,:),mvtk,bladeFrame) From b204b29e3de18c2415be0c075b209de4e2264784 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 09:47:56 -0600 Subject: [PATCH 07/13] AeroDyn/OLAF: use resolved xPts/yPts/zPts in SetRequestedWindPoints --- modules/aerodyn/src/FVW_Subs.f90 | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/modules/aerodyn/src/FVW_Subs.f90 b/modules/aerodyn/src/FVW_Subs.f90 index 29b310b274..904fff13cd 100644 --- a/modules/aerodyn/src/FVW_Subs.f90 +++ b/modules/aerodyn/src/FVW_Subs.f90 @@ -572,7 +572,7 @@ subroutine SetRequestedWindPoints(r_wind, x, p, m) type(FVW_MiscVarType), intent(in ), target :: m !< Initial misc/optimization variables integer(IntKi) :: iP_start,iP_end ! Current index of point, start and end of range integer(IntKi) :: iGrid,i,j,k,iW - real(ReKi) :: xP,yP,zP,dx,dy,dz + real(ReKi) :: xP,yP,zP type(GridOutType), pointer :: g ! Using array reshaping to ensure a given near or far wake point is always at the same location in the array. @@ -608,15 +608,12 @@ subroutine SetRequestedWindPoints(r_wind, x, p, m) iP_start=iP_end+1 do iGrid=1,p%nGridOut g => m%GridOutputs(iGrid) - dx = (g%xEnd- g%xStart)/max(g%nx-1,1) - dy = (g%yEnd- g%yStart)/max(g%ny-1,1) - dz = (g%zEnd- g%zStart)/max(g%nz-1,1) do k=1,g%nz - zP = g%zStart + (k-1)*dz + zP = g%zPts(k) do j=1,g%ny - yP = g%yStart + (j-1)*dy + yP = g%yPts(j) do i=1,g%nx - xP = g%xStart + (i-1)*dx + xP = g%xPts(i) r_wind(1:3,iP_start) = (/xP,yP,zP/) iP_start=iP_start+1 enddo From ea53472ffdeab0cfd8da295aa338985e4953d43e Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 09:50:48 -0600 Subject: [PATCH 08/13] AeroDyn/OLAF: use resolved xPts/yPts/zPts in InducedVelocitiesAll_OnGrid --- modules/aerodyn/src/FVW_Subs.f90 | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/modules/aerodyn/src/FVW_Subs.f90 b/modules/aerodyn/src/FVW_Subs.f90 index 904fff13cd..e9b60f81d0 100644 --- a/modules/aerodyn/src/FVW_Subs.f90 +++ b/modules/aerodyn/src/FVW_Subs.f90 @@ -1135,7 +1135,7 @@ subroutine InducedVelocitiesAll_OnGrid(g, p, x, m, ErrStat, ErrMsg) ! Local variables integer(IntKi) :: nCPs, iHeadP integer(IntKi) :: i,j,k - real(ReKi) :: xP,yP,zP,dx,dy,dz + real(ReKi) :: xP,yP,zP ! TODO new options type(T_Tree) :: Tree type(T_Panl) :: Panl @@ -1149,15 +1149,12 @@ subroutine InducedVelocitiesAll_OnGrid(g, p, x, m, ErrStat, ErrMsg) nCPs = g%nx * g%ny * g%nz allocate(CPs(3, nCPs), stat=ErrStat) iHeadP=1 - dx = (g%xEnd- g%xStart)/max(g%nx-1,1) - dy = (g%yEnd- g%yStart)/max(g%ny-1,1) - dz = (g%zEnd- g%zStart)/max(g%nz-1,1) do k=1,g%nz - zP = g%zStart + (k-1)*dz + zP = g%zPts(k) do j=1,g%ny - yP = g%yStart + (j-1)*dy + yP = g%yPts(j) do i=1,g%nx - xP = g%xStart + (i-1)*dx + xP = g%xPts(i) CPs(1:3,iHeadP) = (/xP,yP,zP/) iHeadP=iHeadP+1 enddo From edb2fc0ae3aeb1ffb2f42c336fadae69d80b8e72 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 10:01:17 -0600 Subject: [PATCH 09/13] FVW_Subs.f90: trim trailing whitespace --- modules/aerodyn/src/FVW_Subs.f90 | 120 +++++++++++++++---------------- 1 file changed, 60 insertions(+), 60 deletions(-) diff --git a/modules/aerodyn/src/FVW_Subs.f90 b/modules/aerodyn/src/FVW_Subs.f90 index e9b60f81d0..3c9061e537 100644 --- a/modules/aerodyn/src/FVW_Subs.f90 +++ b/modules/aerodyn/src/FVW_Subs.f90 @@ -504,7 +504,7 @@ subroutine find_nan_1D(array, varname) tot=tot+1 endif enddo - if (found) then + if (found) then print*,'OLAF NAN ',trim(varname),tot,n STOP endif @@ -526,7 +526,7 @@ subroutine find_nan_2D(array, varname) tot=tot+1 endif enddo - if (found) then + if (found) then print*,'OLAF NAN ',trim(varname),tot,n STOP endif @@ -784,8 +784,8 @@ subroutine FVW_InitStates( x, p, ErrStat, ErrMsg ) allocate(x%W(p%nWings)) do iW=1,p%nWings - call AllocAry( x%W(iW)%Gamma_NW, p%W(iW)%nSpan , p%nNWMax , 'NW Panels Circulation', ErrStat2, ErrMsg2 );call SetErrStat ( ErrStat2, ErrMsg2, ErrStat,ErrMsg,'FVW_InitStates' ); - call AllocAry( x%W(iW)%Gamma_FW, FWnSpan , p%nFWMax , 'FW Panels Circulation', ErrStat2, ErrMsg2 );call SetErrStat ( ErrStat2, ErrMsg2, ErrStat,ErrMsg,'FVW_InitStates' ); + call AllocAry( x%W(iW)%Gamma_NW, p%W(iW)%nSpan , p%nNWMax , 'NW Panels Circulation', ErrStat2, ErrMsg2 );call SetErrStat ( ErrStat2, ErrMsg2, ErrStat,ErrMsg,'FVW_InitStates' ); + call AllocAry( x%W(iW)%Gamma_FW, FWnSpan , p%nFWMax , 'FW Panels Circulation', ErrStat2, ErrMsg2 );call SetErrStat ( ErrStat2, ErrMsg2, ErrStat,ErrMsg,'FVW_InitStates' ); call AllocAry( x%W(iW)%Eps_NW , 3, p%W(iW)%nSpan , p%nNWMax , 'NW Panels Reg Param' , ErrStat2, ErrMsg2 );call SetErrStat ( ErrStat2, ErrMsg2, ErrStat,ErrMsg,'FVW_InitStates' ); call AllocAry( x%W(iW)%Eps_FW , 3, FWnSpan , p%nFWMax , 'FW Panels Reg Param' , ErrStat2, ErrMsg2 );call SetErrStat ( ErrStat2, ErrMsg2, ErrStat,ErrMsg,'FVW_InitStates' ); ! set x%W(iW)%r_NW and x%W(iW)%r_FW to (0,0,0) so that InflowWind can shortcut the calculations @@ -794,10 +794,10 @@ subroutine FVW_InitStates( x, p, ErrStat, ErrMsg ) if (ErrStat >= AbortErrLev) return x%W(iW)%r_NW = 0.0_ReKi x%W(iW)%r_FW = 0.0_ReKi - x%W(iW)%Gamma_NW = 0.0_ReKi ! First call of calcoutput, states might not be set + x%W(iW)%Gamma_NW = 0.0_ReKi ! First call of calcoutput, states might not be set x%W(iW)%Gamma_FW = 0.0_ReKi ! NOTE, these values might be mapped from z%W(iW)%Gamma_LL at init - x%W(iW)%Eps_NW = 0.001_ReKi - x%W(iW)%Eps_FW = 0.001_ReKi + x%W(iW)%Eps_NW = 0.001_ReKi + x%W(iW)%Eps_FW = 0.001_ReKi enddo end subroutine FVW_InitStates @@ -836,7 +836,7 @@ subroutine FVW_InitMiscVarsPostParam( p, m, ErrStat, ErrMsg ) bWakeNeedsPart = p%VelocityMethod(1)==idVelocityPart .or. p%VelocityMethod(1)==idVelocityTreePart bLLNeedsPart = p%VelocityMethod(2)==idVelocityPart .or. p%VelocityMethod(2)==idVelocityTreePart if (bLLNeedsPart .or. bWakeNeedsPart) then - nPart = 0 + nPart = 0 if (bWakeNeedsPart) nPart = max(nPart, nSeg * p%PartPerSegment(1)) if (bLLNeedsPart) nPart = max(nPart, nSeg * p%PartPerSegment(2)) call AllocAry( m%Part%P , 3, nPart, 'PartP' , ErrStat2, ErrMsg2 ); if(Failed())return; m%Part%P = -999999_ReKi; @@ -851,7 +851,7 @@ subroutine FVW_InitMiscVarsPostParam( p, m, ErrStat, ErrMsg ) call AllocAry( m%Uind , 3, nCPs, 'Uind' , ErrStat2, ErrMsg2 ); if(Failed())return; m%Uind= -999999_ReKi; contains logical function Failed() - call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'FVW_InitMiscVarsPostParam') + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'FVW_InitMiscVarsPostParam') Failed = ErrStat >= AbortErrLev end function Failed end subroutine FVW_InitMiscVarsPostParam @@ -1205,7 +1205,7 @@ subroutine SegmentsToPartWrap(Sgmt, nSeg, PartPerSegment, RegFunction, Part, all else ! check that we have enough space if (.not. allocated(Part%P)) then - print*,'>>> PartP not allocated'; + print*,'>>> PartP not allocated'; STOP endif if (size(Part%P,2) Initialize / allocated main variables for source panels. +!> Initialize / allocated main variables for source panels. !! If an non empty input file is provided, the panels points and connectivity are read !! Otherwise, Points and IDs should be provided in "p" -!! -!! Acknowledgements: +!! +!! Acknowledgements: !! The original implementation of the source panel method was funded by Accelerate Wind, !! and implemented by E. Branlard. !! For more acknowledgements, visit: https://openfast.readthedocs.io/en/main/source/acknowledgements.html -!! +!! subroutine srcPnl_init(p, m, z, errStat, errMsg, filename) use VTK !, only: ReadVTK_PD_info, ReadVTK_PD_fields type(T_SrcPanlParam), intent(inout) :: p @@ -1894,7 +1894,7 @@ subroutine srcPnl_init(p, m, z, errStat, errMsg, filename) ! --- Compute influence matrix ! For now, panels don't move so we compute this only once, otherwise, put this in FVW_CalcConstrStateResidual - call srcPnl_build_mat(p, m%AI, m%UUI) + call srcPnl_build_mat(p, m%AI, m%UUI) ! --- Factorization call linalg_factor(m%AI, m%IPIV, errStat2, errMsg2); if(Failed()) return @@ -1912,27 +1912,27 @@ subroutine srcPnl_geometry(Panl, errStat, errMsg) type(T_SrcPanlParam), intent(inout) :: Panl integer(IntKi) , intent(out) :: errStat !< Error status of the operation character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None - real(ReKi) :: alpha !< - real(ReKi) :: d1 !< - real(ReKi) :: DLastRingTE !< - real(ReKi) :: eta0 !< - integer(IntKi), dimension(4) :: IDs !< - real(ReKi), dimension(3) :: P1 !< - real(ReKi), dimension(3) :: P1e !< - real(ReKi), dimension(3) :: P1es !< - real(ReKi), dimension(3) :: P1p !< - real(ReKi), dimension(3) :: P2 !< - real(ReKi), dimension(3) :: P2e !< - real(ReKi), dimension(3) :: P2es !< - real(ReKi), dimension(3) :: P2p !< - real(ReKi), dimension(3) :: P3 !< - real(ReKi), dimension(3) :: P3e !< - real(ReKi), dimension(3) :: P3es !< - real(ReKi), dimension(3) :: P3p !< - real(ReKi), dimension(3) :: P4 !< - real(ReKi), dimension(3) :: P4e !< - real(ReKi), dimension(3) :: P4es !< - real(ReKi), dimension(3) :: P4p !< + real(ReKi) :: alpha !< + real(ReKi) :: d1 !< + real(ReKi) :: DLastRingTE !< + real(ReKi) :: eta0 !< + integer(IntKi), dimension(4) :: IDs !< + real(ReKi), dimension(3) :: P1 !< + real(ReKi), dimension(3) :: P1e !< + real(ReKi), dimension(3) :: P1es !< + real(ReKi), dimension(3) :: P1p !< + real(ReKi), dimension(3) :: P2 !< + real(ReKi), dimension(3) :: P2e !< + real(ReKi), dimension(3) :: P2es !< + real(ReKi), dimension(3) :: P2p !< + real(ReKi), dimension(3) :: P3 !< + real(ReKi), dimension(3) :: P3e !< + real(ReKi), dimension(3) :: P3es !< + real(ReKi), dimension(3) :: P3p !< + real(ReKi), dimension(3) :: P4 !< + real(ReKi), dimension(3) :: P4e !< + real(ReKi), dimension(3) :: P4es !< + real(ReKi), dimension(3) :: P4p !< real(ReKi), dimension(3) :: T1 real(ReKi), dimension(3) :: T2 real(ReKi), dimension(3) :: Ptmp !< temp for computing norm of delta points @@ -1941,7 +1941,7 @@ subroutine srcPnl_geometry(Panl, errStat, errMsg) real(ReKi) :: norm_T1 real(ReKi) :: norm_T2 integer :: ip - real(ReKi) :: xi0 !< + real(ReKi) :: xi0 !< integer(IntKi) :: errStat2 !< temporary Error status character(ErrMsgLen) :: errMsg2 !< temporary Error message errStat = ErrID_None @@ -1964,7 +1964,7 @@ subroutine srcPnl_geometry(Panl, errStat, errMsg) call AllocAry(Panl%eta , 4, Panl%n,'eta ' ,errStat2,errMsg2); if(Failed())return call AllocAry(Panl%xi , 4, Panl%n,'xi ' ,errStat2,errMsg2); if(Failed())return call AllocAry(Panl%Area , Panl%n,'Area ' ,errStat2,errMsg2); if(Failed())return - + do ip = 1, Panl%n IDs = Panl%IDs(:,ip) P1 = Panl%P(:,IDs(1)) @@ -1974,10 +1974,10 @@ subroutine srcPnl_geometry(Panl, errStat, errMsg) Panl%Pmid(1:3,ip) = (P1+P2+P3+p4)/4 ! that's hess's barred coordinates T1 = P3-P1 T2 = P4-P2 - ! maximum diagonal + ! maximum diagonal norm_T1 = sqrt(T1(1)**2+ T1(2)**2+ T1(3)**2) norm_T2 = sqrt(T2(1)**2+ T2(2)**2+ T2(3)**2) - ! flat panel coordinate system + ! flat panel coordinate system Ptmp(1) = T2(2) * T1(3) - T2(3) * T1(2) Ptmp(2) = T2(3) * T1(1) - T2(1) * T1(3) Ptmp(3) = T2(1) * T1(2) - T2(2) * T1(1) @@ -1995,22 +1995,22 @@ subroutine srcPnl_geometry(Panl, errStat, errMsg) Panl%R_g2p(2, 1:3, ip) = T2 Panl%R_g2p(3, 1:3, ip) = Panl%Normal(:,ip) Mat = Panl%R_g2p(:,:,ip) - ! Projection of the surface into a flat panel - Hess primed coordinates + ! Projection of the surface into a flat panel - Hess primed coordinates d1 = dot_product(Panl%Normal(:,ip),Panl%Pmid(:,ip)-P1) P1p = P1+Panl%Normal(:,ip)*(-1)**(1-1)*d1 P2p = P2+Panl%Normal(:,ip)*(-1)**(2-1)*d1 P3p = P3+Panl%Normal(:,ip)*(-1)**(3-1)*d1 P4p = P4+Panl%Normal(:,ip)*(-1)**(4-1)*d1 - !Coordinates of flat panel points in panel coordinate system - Hess starred coordinates with greek letters - ! the transformation is such that the zeta coordinate will always be zero + !Coordinates of flat panel points in panel coordinate system - Hess starred coordinates with greek letters + ! the transformation is such that the zeta coordinate will always be zero P1es = matmul(Mat,(P1p-Panl%Pmid(:,ip))) P2es = matmul(Mat,(P2p-Panl%Pmid(:,ip))) P3es = matmul(Mat,(P3p-Panl%Pmid(:,ip))) P4es = matmul(Mat,(P4p-Panl%Pmid(:,ip))) - ! Coordinates of the centroid + ! Coordinates of the centroid xi0 = 1._ReKi/3._ReKi*1.0_ReKi/(P2es(2)-P4es(2)) * (P4es(1)*(P1es(2)-P2es(2))+P2es(1)*(P4es(2)-P1es(2) )) eta0 = -1._ReKi/3._ReKi * P1es(2) - ! Coordinates based on centroid - Hess greek letters coordinates + ! Coordinates based on centroid - Hess greek letters coordinates P1e = P1es-(/ xi0,eta0,0.0_ReKi /) P2e = P2es-(/ xi0,eta0,0.0_ReKi /) P3e = P3es-(/ xi0,eta0,0.0_ReKi /) @@ -2019,7 +2019,7 @@ subroutine srcPnl_geometry(Panl, errStat, errMsg) Panl%eta(:,ip) = (/ P1e(2), P2e(2), P3e(2), P4e(2)/) ! Centroid in reference frame Panl%Pcent(:,ip) = Panl%Pmid(:,ip) + matmul( (/ xi0,eta0,0.0_ReKi /),Mat) - ! Area + ! Area Panl%Area(ip) = 0.5_ReKi*(Panl%xi(3,ip)-Panl%xi(1,ip))*(Panl%eta(2,ip)-Panl%eta(4,ip)) end do ! Loop on panels contains @@ -2035,8 +2035,8 @@ subroutine srcPnl_build_mat(Panl, AI, UUI) real(ReKi), dimension(:,:), intent(out) :: AI !< (nCPs x nPanels) Self Induced Velocities matrix along normal real(ReKi), dimension(:,:,:), intent(out) :: UUI !< (3 x nCPs x nPanels) Unit induced velocity ! Variables - real(ReKi), parameter :: UnitIntensity=1.0_ReKi !< - real(ReKi), dimension(3) :: Uind_tmp !< + real(ReKi), parameter :: UnitIntensity=1.0_ReKi !< + real(ReKi), dimension(3) :: Uind_tmp !< integer :: icp, ip !< loop variables if (OLAF_PROFILING) call tic('SrcPanel build matrix') !$OMP PARALLEL DEFAULT(shared) @@ -2047,12 +2047,12 @@ subroutine srcPnl_build_mat(Panl, AI, UUI) ! ---- loop on all panls do ip = 1, Panl%n call ui_quad_src_11(Panl%Pcent(:,icp), UnitIntensity, Panl%xi(1:4,ip), Panl%eta(1:4,ip), Panl%Pcent(1:3,ip), Panl%R_g2p(1:3,1:3,ip), Uind_tmp) - ! AI= Vi . N + ! AI= Vi . N AI(icp, ip) = dot_product(Uind_tmp, Panl%Normal(1:3,icp)) UUI(:, icp, ip) = Uind_tmp - end do - end do - !$OMP END DO + end do + end do + !$OMP END DO !$OMP END PARALLEL if (OLAF_PROFILING) call toc() end subroutine srcPnl_build_mat @@ -2079,7 +2079,7 @@ subroutine srcPnl_ExtVelocities_OnPanels(u, p, x, m, errStat, errMsg) m%SrcPnl%Uext(1:3,:) = 0.0_ReKi ! Due to side effects of ui_ functions ! Convert Panels to segments, segments to particles, particles to tree call InducedVelocitiesAll_Init(p, x, m, m%Sgmt, m%Part, Tree, Panl, errStat, errMsg, allocPart=.false.) - ! We don't want the influence of panels so we nullify + ! We don't want the influence of panels so we nullify nullify(Panl%p_Src); nullify(Panl%m_Src) call InducedVelocitiesAll_Calc(p%SrcPnl%Pcent(1:3,:), p%SrcPnl%n, m%SrcPnl%Uext, p, m%Sgmt, m%Part, Tree, Panl, errStat, errMsg) call InducedVelocitiesAll_End(p, Tree, m%Part, Panl, errStat, errMsg, deallocPart=.false.) @@ -2139,7 +2139,7 @@ subroutine srcPnl_calcOutput(p, m, z, rho) !, errStat, errMsg ! Static pressure ps = 1/2 rho Utot**2 ! Reference pressure qinf = 1/2 rho Uwnd**2 ! Pressure Coefficient Cp = ps-pinf/qinf (by convention) - ! Pressure force F = (ps-p0) A e_n + ! Pressure force F = (ps-p0) A e_n if (Uwnd_norm2>0) then Cp = max( 1-(norm2(m%Utot(1:3,ip))**2)/Uwnd_norm2, -10._ReKi) ! Bernoulli. else @@ -2167,7 +2167,7 @@ subroutine linalg_factor(AA, IPIV, errStat, errMsg) n = size(AA,1) m = n call LAPACK_GETRF(m, n, AA, IPIV, errStat, errMsg) -endsubroutine +endsubroutine subroutine linalg_solve(AFact, RHS, IPIV, errStat, errMsg) use NWTC_LAPACK, only : LAPACK_GETRS @@ -2179,7 +2179,7 @@ subroutine linalg_solve(AFact, RHS, IPIV, errStat, errMsg) integer :: n n = size(AFact,1) call LAPACK_GETRS('N', n, AFact, IPIV, RHS, errStat, errMsg) -end subroutine +end subroutine !> Solve A x = B subroutine linalg_solveWrap(AA, B, X, errStat, errMsg) From fdab91ba29cc71e8bd22e1fd82f4e246daf59fa7 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 12:40:20 -0600 Subject: [PATCH 10/13] FVW_Types.f90: updated FVW registry outputs --- modules/aerodyn/src/FVW_Types.f90 | 66 +++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/modules/aerodyn/src/FVW_Types.f90 b/modules/aerodyn/src/FVW_Types.f90 index 53ff476c2f..9f0b0b1aa7 100644 --- a/modules/aerodyn/src/FVW_Types.f90 +++ b/modules/aerodyn/src/FVW_Types.f90 @@ -56,6 +56,12 @@ MODULE FVW_Types REAL(ReKi) , DIMENSION(:,:,:,:), ALLOCATABLE :: uGrid !< Grid velocity 3 x nz x ny x nx [-] REAL(ReKi) , DIMENSION(:,:,:,:), ALLOCATABLE :: omGrid !< Grid vorticity 3 x nz x ny x nx [-] REAL(DbKi) :: tLastOutput = 0.0_R8Ki !< Last output time [-] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: xPts !< Explicit x coordinates (non-equidistant grid, if used) [m] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: yPts !< Explicit y coordinates (non-equidistant grid, if used) [m] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: zPts !< Explicit z coordinates (non-equidistant grid, if used) [m] + CHARACTER(1024) :: xListFile !< File with explicit x coordinates (empty if equidistant) [-] + CHARACTER(1024) :: yListFile !< File with explicit y coordinates (empty if equidistant) [-] + CHARACTER(1024) :: zListFile !< File with explicit z coordinates (empty if equidistant) [-] END TYPE GridOutType ! ======================= ! ========= T_Sgmt ======= @@ -466,6 +472,45 @@ subroutine FVW_CopyGridOutType(SrcGridOutTypeData, DstGridOutTypeData, CtrlCode, DstGridOutTypeData%omGrid = SrcGridOutTypeData%omGrid end if DstGridOutTypeData%tLastOutput = SrcGridOutTypeData%tLastOutput + if (allocated(SrcGridOutTypeData%xPts)) then + LB(1:1) = lbound(SrcGridOutTypeData%xPts) + UB(1:1) = ubound(SrcGridOutTypeData%xPts) + if (.not. allocated(DstGridOutTypeData%xPts)) then + allocate(DstGridOutTypeData%xPts(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstGridOutTypeData%xPts.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstGridOutTypeData%xPts = SrcGridOutTypeData%xPts + end if + if (allocated(SrcGridOutTypeData%yPts)) then + LB(1:1) = lbound(SrcGridOutTypeData%yPts) + UB(1:1) = ubound(SrcGridOutTypeData%yPts) + if (.not. allocated(DstGridOutTypeData%yPts)) then + allocate(DstGridOutTypeData%yPts(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstGridOutTypeData%yPts.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstGridOutTypeData%yPts = SrcGridOutTypeData%yPts + end if + if (allocated(SrcGridOutTypeData%zPts)) then + LB(1:1) = lbound(SrcGridOutTypeData%zPts) + UB(1:1) = ubound(SrcGridOutTypeData%zPts) + if (.not. allocated(DstGridOutTypeData%zPts)) then + allocate(DstGridOutTypeData%zPts(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstGridOutTypeData%zPts.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstGridOutTypeData%zPts = SrcGridOutTypeData%zPts + end if + DstGridOutTypeData%xListFile = SrcGridOutTypeData%xListFile + DstGridOutTypeData%yListFile = SrcGridOutTypeData%yListFile + DstGridOutTypeData%zListFile = SrcGridOutTypeData%zListFile end subroutine subroutine FVW_DestroyGridOutType(GridOutTypeData, ErrStat, ErrMsg) @@ -481,6 +526,15 @@ subroutine FVW_DestroyGridOutType(GridOutTypeData, ErrStat, ErrMsg) if (allocated(GridOutTypeData%omGrid)) then deallocate(GridOutTypeData%omGrid) end if + if (allocated(GridOutTypeData%xPts)) then + deallocate(GridOutTypeData%xPts) + end if + if (allocated(GridOutTypeData%yPts)) then + deallocate(GridOutTypeData%yPts) + end if + if (allocated(GridOutTypeData%zPts)) then + deallocate(GridOutTypeData%zPts) + end if end subroutine subroutine FVW_PackGridOutType(RF, Indata) @@ -505,6 +559,12 @@ subroutine FVW_PackGridOutType(RF, Indata) call RegPackAlloc(RF, InData%uGrid) call RegPackAlloc(RF, InData%omGrid) call RegPack(RF, InData%tLastOutput) + call RegPackAlloc(RF, InData%xPts) + call RegPackAlloc(RF, InData%yPts) + call RegPackAlloc(RF, InData%zPts) + call RegPack(RF, InData%xListFile) + call RegPack(RF, InData%yListFile) + call RegPack(RF, InData%zListFile) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -533,6 +593,12 @@ subroutine FVW_UnPackGridOutType(RF, OutData) call RegUnpackAlloc(RF, OutData%uGrid); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%omGrid); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%tLastOutput); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%xPts); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%yPts); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%zPts); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%xListFile); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%yListFile); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%zListFile); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine FVW_CopyT_Sgmt(SrcT_SgmtData, DstT_SgmtData, CtrlCode, ErrStat, ErrMsg) From 242241e824a8a2a7c2bac5f28718140fed30ee72 Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 13:07:20 -0600 Subject: [PATCH 11/13] AeroDyn/OLAF: fix ResolveGridAxis call (5 arguments) --- modules/aerodyn/src/FVW_IO.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 1ea8dad016..cc43223ba7 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -106,11 +106,11 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) call ReadGridOut(sLine, m%GridOutputs(i)); if(Failed()) return ! Resolve each axis to an explicit coordinate array (regardless of whether it is a list file or a range) call ResolveGridAxis(m%GridOutputs(i)%xStart, m%GridOutputs(i)%xEnd, m%GridOutputs(i)%nx, & - m%GridOutputs(i)%xListFile, PriPath, m%GridOutputs(i)%xPts, ErrStat2, ErrMsg2); if(Failed()) return + m%GridOutputs(i)%xListFile, m%GridOutputs(i)%xPts); if(Failed()) return call ResolveGridAxis(m%GridOutputs(i)%yStart, m%GridOutputs(i)%yEnd, m%GridOutputs(i)%ny, & - m%GridOutputs(i)%yListFile, PriPath, m%GridOutputs(i)%yPts, ErrStat2, ErrMsg2); if(Failed()) return + m%GridOutputs(i)%yListFile, m%GridOutputs(i)%yPts); if(Failed()) return call ResolveGridAxis(m%GridOutputs(i)%zStart, m%GridOutputs(i)%zEnd, m%GridOutputs(i)%nz, & - m%GridOutputs(i)%zListFile, PriPath, m%GridOutputs(i)%zPts, ErrStat2, ErrMsg2); if(Failed()) return + m%GridOutputs(i)%zListFile, m%GridOutputs(i)%zPts); if(Failed()) return ! Error checking if (Check(m%GridOutputs(i)%nx<1, 'Grid output nx needs to be >=1')) return if (Check(m%GridOutputs(i)%ny<1, 'Grid output ny needs to be >=1')) return From 7244c06479cb90881ae1bf2b39922f1f6b071a1b Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 14:21:05 -0600 Subject: [PATCH 12/13] AeroDyn/OLAF docs: document non-equidistant grid output (list-file) --- docs/source/user/aerodyn-olaf/InputFiles.rst | 53 ++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/docs/source/user/aerodyn-olaf/InputFiles.rst b/docs/source/user/aerodyn-olaf/InputFiles.rst index 3c6fabd8ad..afc81f7b85 100644 --- a/docs/source/user/aerodyn-olaf/InputFiles.rst +++ b/docs/source/user/aerodyn-olaf/InputFiles.rst @@ -312,6 +312,59 @@ of a box of shape 5x20x30 and dimension 1200x300x295. The grid contains both th The two other grids are vertical and horizontal planes containing only the velocity. +Non-equidistant grid points +~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +By default, grid points along each axis (X, Y, Z) are equidistant, defined by +the ``Start``, ``End``, and ``n`` columns of the grid output table. + +Alternatively, a ``Start`` cell can be given as a quoted filename instead of a +number. In that case, the corresponding axis uses an explicit, user-defined +list of coordinates read from that file, and the ``End``/``n`` columns for +that axis are ignored. + +This choice is made independently for each axis, so a single grid can mix +equidistant and list-defined axes (e.g., equidistant in X and Z, list-defined +in Y), and a filename can be given for ``XStart``, ``YStart``, ``ZStart``, or +any combination of the three. + +The referenced file must: + +- contain exactly one coordinate value per line, +- list values in strictly ascending order (no duplicates), +- contain no comments or blank lines, +- be located in the same directory as the main OLAF input file. + +Example, requesting a non-equidistant Y-axis:: + + GridName GridType TStart TEnd DTOut XStart XEnd nX YStart YEnd nY ZStart ZEnd nZ + (-) (-) (s) (s) (s) (m) (m) (-) (m) (m) (-) (m) (m) (-) + "Yline" 1 default default default 90 90 1 "Ypoints.dat" - - 100 100 1 + +with ``Ypoints.dat`` (located next to the OLAF input file) containing:: + + -125 + -100 + -75 + -50 + -37.5 + -25 + -12.5 + 0 + 12.5 + 25 + 37.5 + 50 + 75 + 100 + 125 + +.. note:: + Vorticity output (``GridType=2``) requires equidistant spacing + and is not available for a grid that uses a list-defined axis. + Use ``GridType=1`` (velocity only) for such grids. + + Advanced Options ~~~~~~~~~~~~~~~~ From 4fcd5a1da046f0f8b85b58a3461b124f290a012f Mon Sep 17 00:00:00 2001 From: Roger Bergua Date: Fri, 31 Jul 2026 21:55:12 -0600 Subject: [PATCH 13/13] AeroDyn/OLAF: make error reporting explicit in ResolveGridAxis and guard all allocations --- modules/aerodyn/src/FVW_IO.f90 | 68 ++++++++++++++++++++++++---------- 1 file changed, 48 insertions(+), 20 deletions(-) diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index cc43223ba7..a7c252bf04 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -106,16 +106,16 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) call ReadGridOut(sLine, m%GridOutputs(i)); if(Failed()) return ! Resolve each axis to an explicit coordinate array (regardless of whether it is a list file or a range) call ResolveGridAxis(m%GridOutputs(i)%xStart, m%GridOutputs(i)%xEnd, m%GridOutputs(i)%nx, & - m%GridOutputs(i)%xListFile, m%GridOutputs(i)%xPts); if(Failed()) return + m%GridOutputs(i)%xListFile, m%GridOutputs(i)%xPts, ErrStat2, ErrMsg2); if(Failed()) return call ResolveGridAxis(m%GridOutputs(i)%yStart, m%GridOutputs(i)%yEnd, m%GridOutputs(i)%ny, & - m%GridOutputs(i)%yListFile, m%GridOutputs(i)%yPts); if(Failed()) return + m%GridOutputs(i)%yListFile, m%GridOutputs(i)%yPts, ErrStat2, ErrMsg2); if(Failed()) return call ResolveGridAxis(m%GridOutputs(i)%zStart, m%GridOutputs(i)%zEnd, m%GridOutputs(i)%nz, & - m%GridOutputs(i)%zListFile, m%GridOutputs(i)%zPts); if(Failed()) return + m%GridOutputs(i)%zListFile, m%GridOutputs(i)%zPts, ErrStat2, ErrMsg2); if(Failed()) return ! Error checking if (Check(m%GridOutputs(i)%nx<1, 'Grid output nx needs to be >=1')) return if (Check(m%GridOutputs(i)%ny<1, 'Grid output ny needs to be >=1')) return if (Check(m%GridOutputs(i)%nz<1, 'Grid output nz needs to be >=1')) return - ! Vorticity requires equidistant spacing. GridType must be 1 when using list files. + ! Vorticity (GridType=2) requires equidistant spacing. GridType must be 1 when using list files. if (m%GridOutputs(i)%type==idGridVelVorticity) then if (Check( len_trim(m%GridOutputs(i)%xListFile)>0 .or. & len_trim(m%GridOutputs(i)%yListFile)>0 .or. & @@ -385,31 +385,53 @@ end subroutine ReadGridOut ! Resolve one grid axis to an explicit, strictly increasing coordinate array, ! either read from ListFile (if given) or built from an equidistant Start/End/n range. - subroutine ResolveGridAxis(AStart, AEnd, n, ListFile, Pts) + subroutine ResolveGridAxis(AStart, AEnd, n, ListFile, Pts, ErrStat, ErrMsg) real(ReKi), intent(in) :: AStart, AEnd integer(IntKi), intent(inout) :: n character(*), intent(in) :: ListFile real(ReKi), allocatable, intent(out) :: Pts(:) - character(1024) :: FullFile - integer(IntKi) :: UnList, j, IOS - real(ReKi) :: val + integer(IntKi), intent(out) :: ErrStat + character(*), intent(out) :: ErrMsg + ! Locals + character(1024) :: FullFile + integer(IntKi) :: UnList, j, IOS + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + real(ReKi) :: val + + ErrStat = ErrID_None + ErrMsg = '' if (len_trim(ListFile) == 0) then ! Equidistant points, expressed as an explicit array - allocate(Pts(max(n,1))) + if (n < 1) then + call SetErrStat(ErrID_Fatal, 'ResolveGridAxis: grid axis definition requires n >= 1 (got '//trim(Num2LStr(n))//').', & + ErrStat, ErrMsg, 'ResolveGridAxis') + return + endif + allocate(Pts(n), stat=IOS) + if (IOS /= 0) then + call SetErrStat(ErrID_Fatal, 'ResolveGridAxis: error allocating grid axis points array (n='//trim(Num2LStr(n))//').', & + ErrStat, ErrMsg, 'ResolveGridAxis') + return + endif do j = 1, n Pts(j) = AStart + (AEnd - AStart) * real(j-1,ReKi) / real(max(n-1,1),ReKi) enddo - ErrStat2 = ErrID_None; ErrMsg2 = '' return endif ! Explicit list from file FullFile = ListFile if (PathIsRelative(FullFile)) FullFile = trim(PriPath)//trim(FullFile) + call GetNewUnit(UnList) call OpenFInpFile(UnList, trim(FullFile), ErrStat2, ErrMsg2) - if (ErrStat2 /= ErrID_None) return + if (ErrStat2 /= ErrID_None) then + call SetErrStat(ErrID_Fatal, 'ResolveGridAxis: could not open grid point list file "'//trim(FullFile)//'": '//trim(ErrMsg2), & + ErrStat, ErrMsg, 'ResolveGridAxis') + return + endif ! Count valid lines n = 0 @@ -418,16 +440,23 @@ subroutine ResolveGridAxis(AStart, AEnd, n, ListFile, Pts) if (IOS /= 0) exit n = n + 1 enddo + if (n < 1) then - ErrStat2 = ErrID_Fatal - ErrMsg2 = 'Grid point list file "'//trim(FullFile)//'" contains no valid values.' + call SetErrStat(ErrID_Fatal, 'ResolveGridAxis: grid point list file "'//trim(FullFile)//'" contains no valid values.', & + ErrStat, ErrMsg, 'ResolveGridAxis') close(UnList) return endif ! Read values rewind(UnList) - allocate(Pts(n)) + allocate(Pts(n), stat=IOS) + if (IOS /= 0) then + call SetErrStat(ErrID_Fatal, 'ResolveGridAxis: error allocating grid axis points array (n='//trim(Num2LStr(n))//').', & + ErrStat, ErrMsg, 'ResolveGridAxis') + close(UnList) + return + endif do j = 1, n read(UnList, *) Pts(j) enddo @@ -436,16 +465,15 @@ subroutine ResolveGridAxis(AStart, AEnd, n, ListFile, Pts) ! Validate strictly increasing, no duplicates do j = 2, n if (Pts(j) <= Pts(j-1)) then - ErrStat2 = ErrID_Fatal - ErrMsg2 = 'Grid point list file "'//trim(FullFile)//'" must be strictly increasing (no duplicates); '// & - 'value at line '//trim(Num2LStr(j))//' ('//trim(Num2LStr(Pts(j)))//') is not greater than '// & - 'the previous value ('//trim(Num2LStr(Pts(j-1)))//').' + call SetErrStat(ErrID_Fatal, & + 'ResolveGridAxis: grid point list file "'//trim(FullFile)//'" must be strictly increasing (no duplicates); '// & + 'value at line '//trim(Num2LStr(j))//' ('//trim(Num2LStr(Pts(j)))//') is not greater than '// & + 'the previous value ('//trim(Num2LStr(Pts(j-1)))//').', & + ErrStat, ErrMsg, 'ResolveGridAxis') return endif enddo - ErrStat2 = ErrID_None - ErrMsg2 = '' end subroutine ResolveGridAxis END SUBROUTINE FVW_ReadInputFile