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 ~~~~~~~~~~~~~~~~ diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 4cde5064d4..a7c252bf04 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. @@ -100,20 +100,35 @@ 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, 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, 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, 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 (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. & + len_trim(m%GridOutputs(i)%zListFile)>0, & + 'Grid "'//trim(m%GridOutputs(i)%name)//'": vorticity requires equidistant spacing. GridType must be 1 when using list files.')) return + endif enddo endif ! --- 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:') @@ -203,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 @@ -216,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 @@ -289,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.' @@ -300,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 @@ -327,24 +342,140 @@ 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='' 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, ErrStat, ErrMsg) + real(ReKi), intent(in) :: AStart, AEnd + integer(IntKi), intent(inout) :: n + character(*), intent(in) :: ListFile + real(ReKi), allocatable, intent(out) :: Pts(:) + 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 + 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 + 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) 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 + do + read(UnList, *, iostat=IOS) val + if (IOS /= 0) exit + n = n + 1 + enddo + + if (n < 1) then + 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), 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 + close(UnList) + + ! Validate strictly increasing, no duplicates + do j = 2, n + if (Pts(j) <= Pts(j-1)) then + 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 + + end subroutine ResolveGridAxis + END SUBROUTINE FVW_ReadInputFile @@ -387,7 +518,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 @@ -399,7 +530,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 @@ -427,7 +558,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 @@ -445,7 +576,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 @@ -469,7 +600,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 @@ -497,6 +628,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 @@ -510,19 +642,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 @@ -569,14 +711,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) 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" - diff --git a/modules/aerodyn/src/FVW_Subs.f90 b/modules/aerodyn/src/FVW_Subs.f90 index 29b310b274..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 @@ -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 @@ -787,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 @@ -797,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 @@ -839,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; @@ -854,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 @@ -1138,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 @@ -1152,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 @@ -1211,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 @@ -1900,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 @@ -1918,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 @@ -1947,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 @@ -1970,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)) @@ -1980,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) @@ -2001,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 /) @@ -2025,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 @@ -2041,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) @@ -2053,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 @@ -2085,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.) @@ -2145,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 @@ -2173,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 @@ -2185,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) 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)