diff --git a/.gitignore b/.gitignore index 09166cd156..74be88339b 100644 --- a/.gitignore +++ b/.gitignore @@ -64,3 +64,8 @@ varcache # Python cache files openfast_io/dist/ openfast_io/openfast_io/_version.py + +# docker with VScode config +*.code-workspace +docker-compose.yml +.devcontainer/ 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/glue-codes/fast-farm/src/FAST_Farm_IO.f90 b/glue-codes/fast-farm/src/FAST_Farm_IO.f90 index 32fce6d53c..d954a4ea39 100644 --- a/glue-codes/fast-farm/src/FAST_Farm_IO.f90 +++ b/glue-codes/fast-farm/src/FAST_Farm_IO.f90 @@ -1047,6 +1047,7 @@ SUBROUTINE Farm_ValidateInput( p, WD_InitInp, AWAE_InitInp, ErrStat, ErrMsg ) CHARACTER(*), PARAMETER :: RoutineName = 'Farm_ValidateInput' INTEGER(IntKi) :: n_disDT_dt character(60) :: tmpStr + logical :: file_exists ! Flag for file existence check ErrStat = ErrID_None ErrMsg = "" @@ -1066,8 +1067,11 @@ SUBROUTINE Farm_ValidateInput( p, WD_InitInp, AWAE_InitInp, ErrStat, ErrMsg ) ! TODO : Verify that the DLL file exists ! --- SHARED MOORING SYSTEM --- - ! TODO : Verify that p%MD_FileName file exists - if ((p%DT_mooring <= 0.0_ReKi) .or. (p%DT_mooring > p%DT_high)) CALL SetErrStat(ErrID_Fatal,'DT_mooring must be greater than zero and no greater than dt_high.',ErrStat,ErrMsg,RoutineName) + if (p%MooringMod == 3) then + inquire(file=trim(p%MD_FileName), exist=file_exists) + if (.not. file_exists) call SetErrStat(ErrID_Fatal,'Cannot find MoorDyn input file '//trim(p%MD_FileName),ErrStat,ErrMsg,RoutineName) + if ((p%DT_mooring <= 0.0_ReKi) .or. (p%DT_mooring > p%DT_high)) CALL SetErrStat(ErrID_Fatal,'DT_mooring must be greater than zero and no greater than dt_high.',ErrStat,ErrMsg,RoutineName) + end if ! --- AMBIENT WIND: INFLOWWIND MODULE --- [used only for Mod_AmbWind=2 or 3] --- ! FIXME: this really should be checked with the turbine specific size diameter -- maybe relocate this check to AWAE or in FF after initializing all turbines? diff --git a/modules/aerodyn/src/AeroAcoustics.f90 b/modules/aerodyn/src/AeroAcoustics.f90 index 9142a8eea4..708ab8fe41 100644 --- a/modules/aerodyn/src/AeroAcoustics.f90 +++ b/modules/aerodyn/src/AeroAcoustics.f90 @@ -165,6 +165,7 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) character(*), parameter :: RoutineName = 'SetParameters' REAL(ReKi) :: val1,val10,f2,f4, dist1, dist10 REAL(ReKi) :: BladeSpanUsedForNoise + REAL(ReKi) :: LastElemPct ! Initialize variables for this routine ErrStat = ErrID_None @@ -271,7 +272,22 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) p%BlSpn = InitInp%BlSpn p%BlChord = InitInp%BlChord - p%startnode = max(1, p%NumBlNds - 1) + ! Calculate last element size as percentage of blade span + IF (p%NumBlNds > 1) THEN + LastElemPct = 100.0 * (p%BlSpn(p%NumBlNds,1) - p%BlSpn(p%NumBlNds-1,1)) / p%BlSpn(p%NumBlNds,1) + ELSE + LastElemPct = 100.0 ! Single node means element spans entire blade + ENDIF + + IF (InputFileData%AA_Bl_Prcntge .lt. LastElemPct) THEN + CALL SetErrStat(ErrID_Warn, 'BldPrcnt is smaller than the last blade element size. '// & + 'OpenFAST will move on assuming the last blade element, which is the minimum blade span used for noise calculations. '// & + 'Either increase BldPrcnt in your aeroacoustic input file '// & + 'or refine the aerodynamic mesh in your blade AeroDyn input file.', ErrStat2, ErrMsg2, RoutineName ) + if(Failed()) return + endif + + p%startnode = min(p%NumBlNds, 2) BladeSpanUsedForNoise = p%BlSpn(p%NumBlNds,1)*(1.0 - InputFileData%AA_Bl_Prcntge/100.0) do j=p%NumBlNds-1,2,-1 IF ( p%BlSpn(j,1) .lt. BladeSpanUsedForNoise )THEN @@ -279,7 +295,7 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) exit ! exit the loop endif enddo - p%startnode = max(min(p%NumBlNds,2),p%startnode) + p%BlElemSpn = 0; DO I = 1,p%numBlades diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 4cde5064d4..a037c85dc3 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,146 @@ 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, *, iostat=IOS) Pts(j) + if (IOS /= 0) then + call SetErrStat(ErrID_Fatal, 'ResolveGridAxis: error reading grid point #'//trim(Num2LStr(j))//' from grid point list file "'//trim(FullFile)//'" (iostat='//trim(Num2LStr(IOS))//').', & + ErrStat, ErrMsg, 'ResolveGridAxis') + close(UnList) + return + end if + 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 +524,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 +536,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 +564,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 +582,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 +606,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 +634,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 +648,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 +717,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) diff --git a/modules/inflowwind/src/IfW_FlowField.f90 b/modules/inflowwind/src/IfW_FlowField.f90 index 922b831843..43f21947e0 100644 --- a/modules/inflowwind/src/IfW_FlowField.f90 +++ b/modules/inflowwind/src/IfW_FlowField.f90 @@ -78,9 +78,19 @@ subroutine IfW_FlowField_GetVelAcc(FF, IStart, Time, PositionXYZ, VelocityUVW, A ! Determine if acceleration should be calculated and returned OutputAccel = allocated(AccelUVW) - if (OutputAccel .and. .not. FF%AccFieldValid) then - call SetErrStat(ErrID_Fatal, "Accel output requested, but accel field is not valid", & - ErrStat, ErrMsg, RoutineName) + ! Cubic velocity interpolation also requires a valid acceleration field, + ! since its formula uses the derivative data even when acceleration output is not requested. + if ((OutputAccel .or. FF%VelInterpCubic) .and. .not. FF%AccFieldValid) then + if (OutputAccel .and. FF%VelInterpCubic) then + call SetErrStat(ErrID_Fatal, "Acceleration output and cubic velocity interpolation both require a valid acceleration field, but the acceleration field is not valid", & + ErrStat, ErrMsg, RoutineName) + else if (OutputAccel) then + call SetErrStat(ErrID_Fatal, "Acceleration output requested, but the acceleration field is not valid", & + ErrStat, ErrMsg, RoutineName) + else ! FF%VelInterpCubic + call SetErrStat(ErrID_Fatal, "Cubic velocity interpolation requires a valid acceleration field, but the acceleration field is not valid", & + ErrStat, ErrMsg, RoutineName) + end if return end if @@ -224,8 +234,10 @@ subroutine IfW_FlowField_GetVelAcc(FF, IStart, Time, PositionXYZ, VelocityUVW, A ! Calculate grid cells for interpolation, returns velocity and acceleration ! components at corners of grid cell containing time and position. Also - ! returns interpolation values Xi. - call Grid3DField_GetCell(FF%Grid3D, Time, Position(:, i), OutputAccel, GridExceedAllow, & + ! returns interpolation values Xi. AccCell is required whenever cubic + ! velocity interpolation is used (its formula uses the derivative data), + ! not only when acceleration is explicitly requested as an output. + call Grid3DField_GetCell(FF%Grid3D, Time, Position(:, i), OutputAccel .or. FF%VelInterpCubic, GridExceedAllow, & VelCell, AccCell, Xi, Is3D, TmpErrStat, TmpErrMsg) if (TmpErrStat >= AbortErrLev) then call SetErrStat(TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName) @@ -1442,9 +1454,18 @@ subroutine IfW_Grid3DField_CalcAccel(G3D, ErrStat, ErrMsg) call SetErrStat(TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return - ! If number of time grids is 1 or 2, set all accelerations to zero, return - if (G3D%NTGrids < 3) then + ! If number of time steps is 1 or 2, set all accelerations to zero, return + if (G3D%NSteps < 3) then G3D%Acc = 0.0_SiKi + ! Also allocate/zero tower acceleration so GetCellInTower has a valid array to reference + if (G3D%NTGrids > 0) then + call AllocAry(G3D%AccTower, size(G3D%VelTower, dim=1), & + size(G3D%VelTower, dim=2), size(G3D%VelTower, dim=3), & + 'tower wind acceleration data.', TmpErrStat, TmpErrMsg) + call SetErrStat(TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return + G3D%AccTower = 0.0_SiKi + end if return end if @@ -1477,16 +1498,14 @@ subroutine IfW_Grid3DField_CalcAccel(G3D, ErrStat, ErrMsg) call SetErrStat(TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return - ! If number of time grids is 1 or 2, set all accelerations to zero - if (G3D%NTGrids < 3) then - G3D%Acc = 0.0_SiKi - else ! Otherwise, calculate acceleration at each grid point - do iz = 1, G3D%NTGrids - do ic = 1, G3D%NComp - call CalcCubicSplineDeriv(G3D%NSteps, G3D%DTime, G3D%VelTower(ic, iz, :), G3D%AccTower(ic, iz, :)) - end do + ! Calculate acceleration at each tower grid point. NSteps < 3 is already handled by the + ! early return above, and each tower height is interpolated independently in time, so no + ! minimum tower height count is required here. + do iz = 1, G3D%NTGrids + do ic = 1, G3D%NComp + call CalcCubicSplineDeriv(G3D%NSteps, G3D%DTime, G3D%VelTower(ic, iz, :), G3D%AccTower(ic, iz, :)) end do - end if + end do contains diff --git a/modules/inflowwind/tests/inflowwind_utest.F90 b/modules/inflowwind/tests/inflowwind_utest.F90 index 758c25de13..30c2a5e258 100644 --- a/modules/inflowwind/tests/inflowwind_utest.F90 +++ b/modules/inflowwind/tests/inflowwind_utest.F90 @@ -3,6 +3,7 @@ program inflowwind_utest use testdrive, only: run_testsuite, new_testsuite, testsuite_type use test_bladed_wind, only: test_bladed_wind_suite +use test_grid3d_field, only: test_grid3d_field_suite use test_hawc_wind, only: test_hawc_wind_suite use test_outputs, only: test_outputs_suite use test_steady_wind, only: test_steady_wind_suite @@ -21,6 +22,7 @@ program inflowwind_utest testsuites = [ & new_testsuite("Bladed Wind", test_bladed_wind_suite), & + new_testsuite("Grid3D Field", test_grid3d_field_suite), & new_testsuite("HAWC Wind", test_hawc_wind_suite), & new_testsuite("Outputs", test_outputs_suite), & new_testsuite("Steady Wind", test_steady_wind_suite), & diff --git a/modules/inflowwind/tests/test_grid3d_field.F90 b/modules/inflowwind/tests/test_grid3d_field.F90 new file mode 100644 index 0000000000..16a0441e43 --- /dev/null +++ b/modules/inflowwind/tests/test_grid3d_field.F90 @@ -0,0 +1,218 @@ +module test_grid3d_field + +use, intrinsic :: ieee_arithmetic, only: ieee_is_finite +use testdrive, only: new_unittest, unittest_type, error_type, check +use ifw_test_tools +use IfW_FlowField +use IfW_FlowField_Types +use NWTC_Library + +implicit none +private +public :: test_grid3d_field_suite + +integer(IntKi), parameter :: NY = 4, NZ = 4, NT = 10 +real(ReKi), parameter :: DTIME = 0.1_ReKi + +contains + +!> Collect all exported unit tests +subroutine test_grid3d_field_suite(testsuite) + type(unittest_type), allocatable, intent(out) :: testsuite(:) + testsuite = [ & + new_unittest("test_grid3d_cubic_vel_only", test_grid3d_cubic_vel_only), & + new_unittest("test_grid3d_calcaccel_no_tower", test_grid3d_calcaccel_no_tower), & + new_unittest("test_grid3d_calcaccel_few_tower_points", test_grid3d_calcaccel_few_tower_points) & + ] +end subroutine + +!> Reproduces a bug in IfW_FlowField_GetVelAcc: when VelInterpCubic is enabled and the +!! caller does not request acceleration output (AccelUVW left unallocated), the local +!! AccCell array used by the cubic Hermite velocity formula is never populated, so the +!! returned velocity is computed from uninitialized memory. This test builds a minimal +!! Grid3D flow field with a velocity that is a known linear ramp in time (spatially +!! uniform), so the exact correct answer at any query time is known, and checks that +!! querying velocity without requesting acceleration gives the same (finite, correct) +!! answer as querying with acceleration requested. +subroutine test_grid3d_cubic_vel_only(error) + type(error_type), allocatable, intent(out) :: error + + type(FlowFieldType) :: FF + real(ReKi), allocatable :: Position(:, :), VelNoAcc(:, :), VelWithAcc(:, :) + real(ReKi), allocatable :: AccelUVW(:, :) + integer(IntKi) :: it, iy, iz + integer(IntKi) :: TmpErrStat + character(ErrMsgLen) :: TmpErrMsg + real(DbKi) :: QueryTime + real(ReKi) :: Expected + + ! Build minimal Grid3D field: NY x NZ spatial points, NT time steps, no tower grid. + FF%FieldType = Grid3D_FieldType + FF%VelInterpCubic = .true. + FF%RotateWindBox = .false. + + FF%Grid3D%NComp = 3 + FF%Grid3D%NYGrids = NY + FF%Grid3D%NZGrids = NZ + FF%Grid3D%NTGrids = 0 + FF%Grid3D%NSteps = NT + FF%Grid3D%DTime = DTIME + FF%Grid3D%Rate = 1.0_ReKi/DTIME + FF%Grid3D%YHWid = 5.0_ReKi + FF%Grid3D%ZHWid = 5.0_ReKi + FF%Grid3D%GridBase = 5.0_ReKi + FF%Grid3D%InvDY = real(NY - 1, ReKi)/(2.0_ReKi*FF%Grid3D%YHWid) + FF%Grid3D%InvDZ = real(NZ - 1, ReKi)/(2.0_ReKi*FF%Grid3D%ZHWid) + FF%Grid3D%MeanWS = 8.0_ReKi + FF%Grid3D%InvMWS = 1.0_ReKi/FF%Grid3D%MeanWS + FF%Grid3D%InitXPosition = 0.0_ReKi + FF%Grid3D%TotalTime = real(NT - 1, ReKi)*DTIME + FF%Grid3D%Periodic = .false. + FF%Grid3D%InterpTower = .false. + + allocate (FF%Grid3D%Vel(3, NY, NZ, NT)) + FF%Grid3D%Vel = 0.0_SiKi + ! U component is a spatially-uniform linear ramp in time: U(t) = t (seconds -> m/s) + do it = 1, NT + do iz = 1, NZ + do iy = 1, NY + FF%Grid3D%Vel(1, iy, iz, it) = real((it - 1), SiKi)*real(DTIME, SiKi) + end do + end do + end do + + ! Set the true time-derivative (dU/dt = 1, dV/dt = dW/dt = 0) directly rather than via + ! IfW_Grid3DField_CalcAccel so this test isolates the IfW_FlowField_GetVelAcc behavior. + ! (test_grid3d_calcaccel_no_tower below covers IfW_Grid3DField_CalcAccel itself.) + allocate (FF%Grid3D%Acc(3, NY, NZ, NT)) + FF%Grid3D%Acc = 0.0_SiKi + FF%Grid3D%Acc(1, :, :, :) = 1.0_SiKi + FF%AccFieldValid = .true. + + ! Query at a position centered in the grid, at a time between samples + allocate (Position(3, 1)) + Position(:, 1) = [0.0_ReKi, 0.0_ReKi, FF%Grid3D%GridBase + FF%Grid3D%ZHWid] + QueryTime = 0.25_DbKi + Expected = real(QueryTime, ReKi) + + ! Case 1: acceleration NOT requested (AccelUVW left unallocated) - buggy path + allocate (VelNoAcc(3, 1)) + call IfW_FlowField_GetVelAcc(FF, 1, QueryTime, Position, VelNoAcc, AccelUVW, TmpErrStat, TmpErrMsg) + call check(error, TmpErrStat, ErrID_None, message='GetVelAcc (no accel) error: '//trim(TmpErrMsg)); if (allocated(error)) return + + call check(error, ieee_is_finite(VelNoAcc(1, 1)), message='Velocity (no accel requested) is not finite (NaN/Inf)'); if (allocated(error)) return + call check(error, VelNoAcc(1, 1), Expected, thr=1.0e-3_ReKi); if (allocated(error)) return + + ! Case 2: acceleration requested (AccelUVW allocated) - reference path + allocate (VelWithAcc(3, 1)) + allocate (AccelUVW(3, 1)) + call IfW_FlowField_GetVelAcc(FF, 1, QueryTime, Position, VelWithAcc, AccelUVW, TmpErrStat, TmpErrMsg) + call check(error, TmpErrStat, ErrID_None, message='GetVelAcc (with accel) error: '//trim(TmpErrMsg)); if (allocated(error)) return + + call check(error, ieee_is_finite(VelWithAcc(1, 1)), message='Velocity (accel requested) is not finite (NaN/Inf)'); if (allocated(error)) return + call check(error, VelWithAcc(1, 1), Expected, thr=1.0e-3_ReKi); if (allocated(error)) return + + ! The two calls must agree: requesting acceleration must not change the velocity result + call check(error, VelNoAcc(1, 1), VelWithAcc(1, 1), thr=1.0e-6_ReKi); if (allocated(error)) return + +end subroutine + +!> Reproduces a bug in IfW_Grid3DField_CalcAccel: it checks G3D%NTGrids (the tower-grid +!! point count) instead of G3D%NSteps (the number of time samples) to decide whether to +!! compute real cubic-spline time derivatives. For the common case of no tower file +!! (NTGrids=0), this unconditionally zeroes G3D%Acc regardless of how many time steps +!! are actually available, silently defeating cubic-in-time interpolation. This test +!! uses a spatially-uniform linear velocity ramp in time (dU/dt = 1 exactly), so the +!! correct computed derivative is known, and checks that it is recovered when there is +!! no tower grid but plenty of time steps. +subroutine test_grid3d_calcaccel_no_tower(error) + type(error_type), allocatable, intent(out) :: error + + type(Grid3DFieldType) :: G3D + integer(IntKi) :: it, iy, iz + integer(IntKi) :: TmpErrStat + character(ErrMsgLen) :: TmpErrMsg + + G3D%NComp = 3 + G3D%NYGrids = NY + G3D%NZGrids = NZ + G3D%NTGrids = 0 + G3D%NSteps = NT + G3D%DTime = DTIME + G3D%Periodic = .false. + + allocate (G3D%Vel(3, NY, NZ, NT)) + G3D%Vel = 0.0_SiKi + do it = 1, NT + do iz = 1, NZ + do iy = 1, NY + G3D%Vel(1, iy, iz, it) = real((it - 1), SiKi)*real(DTIME, SiKi) + end do + end do + end do + + call IfW_Grid3DField_CalcAccel(G3D, TmpErrStat, TmpErrMsg) + call check(error, TmpErrStat, ErrID_None, message='CalcAccel error: '//trim(TmpErrMsg)); if (allocated(error)) return + + ! Interior time points should recover close to the true derivative (dU/dt = 1); a + ! natural cubic spline has some boundary-driven error, but zero here would indicate + ! the NTGrids-vs-NSteps bug has regressed (computation skipped entirely). + call check(error, real(G3D%Acc(1, 2, 2, 5), ReKi), 1.0_ReKi, thr=0.01_ReKi); if (allocated(error)) return + +end subroutine + +!> Reproduces the mirror-image bug in IfW_Grid3DField_CalcAccel's tower branch: it checked +!! G3D%NTGrids (tower height count) against the same <3 threshold meant for time-step counts, +!! zeroing out tower acceleration whenever there were only 1 or 2 tower grid points, even though +!! each tower height's time derivative is computed independently and only needs G3D%NSteps >= 3 +!! (already guaranteed by the earlier NSteps<3 early return). This test uses NTGrids=2 with a +!! spatially-uniform linear velocity ramp in time (dU/dt = 1 exactly) and checks that the tower +!! acceleration is actually computed rather than silently zeroed. +subroutine test_grid3d_calcaccel_few_tower_points(error) + type(error_type), allocatable, intent(out) :: error + + integer(IntKi), parameter :: NTWR = 2 + type(Grid3DFieldType) :: G3D + integer(IntKi) :: it, iy, iz + integer(IntKi) :: TmpErrStat + character(ErrMsgLen) :: TmpErrMsg + + G3D%NComp = 3 + G3D%NYGrids = NY + G3D%NZGrids = NZ + G3D%NTGrids = NTWR + G3D%NSteps = NT + G3D%DTime = DTIME + G3D%Periodic = .false. + + allocate (G3D%Vel(3, NY, NZ, NT)) + G3D%Vel = 0.0_SiKi + do it = 1, NT + do iz = 1, NZ + do iy = 1, NY + G3D%Vel(1, iy, iz, it) = real((it - 1), SiKi)*real(DTIME, SiKi) + end do + end do + end do + + allocate (G3D%VelTower(3, NTWR, NT)) + G3D%VelTower = 0.0_SiKi + do it = 1, NT + do iz = 1, NTWR + G3D%VelTower(1, iz, it) = real((it - 1), SiKi)*real(DTIME, SiKi) + end do + end do + + call IfW_Grid3DField_CalcAccel(G3D, TmpErrStat, TmpErrMsg) + call check(error, TmpErrStat, ErrID_None, message='CalcAccel error: '//trim(TmpErrMsg)); if (allocated(error)) return + + call check(error, allocated(G3D%AccTower), message='AccTower was not allocated'); if (allocated(error)) return + + ! Interior time point at either tower height should recover close to the true derivative + ! (dU/dt = 1); zero here would indicate the NTGrids<3 tower guard bug has regressed. + call check(error, real(G3D%AccTower(1, 1, 5), ReKi), 1.0_ReKi, thr=0.01_ReKi); if (allocated(error)) return + call check(error, real(G3D%AccTower(1, 2, 5), ReKi), 1.0_ReKi, thr=0.01_ReKi); if (allocated(error)) return + +end subroutine + +end module diff --git a/modules/nwtc-library/src/GridInterp.f90 b/modules/nwtc-library/src/GridInterp.f90 index 769b905fca..22ba5a9de1 100644 --- a/modules/nwtc-library/src/GridInterp.f90 +++ b/modules/nwtc-library/src/GridInterp.f90 @@ -302,7 +302,7 @@ Subroutine GridInterpSetup3D( position, p, m, ErrStat, ErrMsg ) character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None character(*), parameter :: RoutineName = 'GridInterpSetup3D' - integer(IntKi) :: dim,i,j,k + integer(IntKi) :: dim,j,k integer(IntKi) :: support real(ReKi) :: N1D(4,3) real(ReKi) :: isopc ! isoparametric coordinates @@ -320,9 +320,7 @@ Subroutine GridInterpSetup3D( position, p, m, ErrStat, ErrMsg ) do k = 1,4 do j = 1,4 - do i = 1,4 - m%N3D(i,j,k) = N1D(i,1)*N1D(j,2)*N1D(k,3) - end do + m%N3D(:,j,k) = N1D(:,1) * (N1D(j,2)*N1D(k,3)) end do end do @@ -343,9 +341,10 @@ Subroutine GridInterpSetup4D( position, p, m, ErrStat, ErrMsg ) character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None character(*), parameter :: RoutineName = 'GridInterpSetup4D' - integer(IntKi) :: dim,i,j,k,l + integer(IntKi) :: dim,j,k,l integer(IntKi) :: support real(ReKi) :: N1D(4,4) + real(ReKi) :: Nkl real(ReKi) :: isopc ! isoparametric coordinates integer(IntKi) :: ErrStat2 character(ErrMsgLen) :: ErrMsg2 @@ -361,10 +360,9 @@ Subroutine GridInterpSetup4D( position, p, m, ErrStat, ErrMsg ) do l = 1,4 do k = 1,4 + Nkl = N1D(k,3)*N1D(l,4) do j = 1,4 - do i = 1,4 - m%N4D(i,j,k,l) = N1D(i,1)*N1D(j,2)*N1D(k,3)*N1D(l,4) - end do + m%N4D(:,j,k,l) = N1D(:,1) * (N1D(j,2)*Nkl) end do end do end do @@ -386,7 +384,7 @@ Subroutine GridInterpSetupN( position, p, m, ErrStat, ErrMsg ) character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None character(*), parameter :: RoutineName = 'GridInterpSetupN' - integer(IntKi) :: dim,i,j,k + integer(IntKi) :: dim,j,k integer(IntKi) :: support real(ReKi) :: N1D(4,3) real(ReKi) :: N1Ddx(4,2:3) @@ -409,10 +407,8 @@ Subroutine GridInterpSetupN( position, p, m, ErrStat, ErrMsg ) ! Need two sets of weights for d(.)/dx and d(.)/dy. Borrow m%N4D for this. do k = 1,4 do j = 1,4 - do i = 1,4 - m%N4D(i,j,k,1) = N1D(i,1)*N1Ddx(j,2)*N1D (k,3) - m%N4D(i,j,k,2) = N1D(i,1)*N1D (j,2)*N1Ddx(k,3) - end do + m%N4D(:,j,k,1) = N1D(:,1) * (N1Ddx(j,2)*N1D (k,3)) + m%N4D(:,j,k,2) = N1D(:,1) * (N1D (j,2)*N1Ddx(k,3)) end do end do @@ -436,17 +432,21 @@ function GridInterp3DR4( data, m ) character(*), parameter :: RoutineName = 'GridInterp3DR4' real(SiKi) :: GridInterp3DR4 - integer(IntKi) :: i,j,k + real(SiKi) :: acc(4) + integer(IntKi) :: i,j,k,jj,kk ! interpolate - GridInterp3DR4 = 0.0_SiKi + acc = 0.0_SiKi do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - GridInterp3DR4 = GridInterp3DR4 + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + acc(i) = acc(i) + m%N3D(i,j,k) * data( m%Indx(i,1), jj, kk ) end do end do end do + GridInterp3DR4 = acc(1) + acc(2) + acc(3) + acc(4) end function GridInterp3DR4 @@ -456,17 +456,21 @@ function GridInterp3DR8( data, m ) character(*), parameter :: RoutineName = 'GridInterp3DR8' real(DbKi) :: GridInterp3DR8 - integer(IntKi) :: i,j,k + real(DbKi) :: acc(4) + integer(IntKi) :: i,j,k,jj,kk ! interpolate - GridInterp3DR8 = 0.0_DbKi + acc = 0.0_DbKi do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - GridInterp3DR8 = GridInterp3DR8 + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + acc(i) = acc(i) + m%N3D(i,j,k) * data( m%Indx(i,1), jj, kk ) end do end do end do + GridInterp3DR8 = acc(1) + acc(2) + acc(3) + acc(4) end function GridInterp3DR8 @@ -482,18 +486,23 @@ function GridInterp3DVecR4( data, m ) character(*), parameter :: RoutineName = 'GridInterp3DVecR4' integer(IntKi), parameter :: vDim = 3 integer(IntKi) :: i,j,k,vi + integer(IntKi) :: jj,kk + real(SiKi) :: acc(4) real(SiKi) :: GridInterp3DVecR4(vDim) ! interpolate - GridInterp3DVecR4 = 0.0_SiKi - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp3DVecR4(vi) = GridInterp3DVecR4(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + do vi = 1,vDim + acc = 0.0_SiKi + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N3D(i,j,k) * data( m%Indx(i,1), jj, kk, vi ) end do end do end do + GridInterp3DVecR4(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp3DVecR4 @@ -505,18 +514,23 @@ function GridInterp3DVecR8( data, m ) character(*), parameter :: RoutineName = 'GridInterp3DVecR8' integer(IntKi), parameter :: vDim = 3 integer(IntKi) :: i,j,k,vi + integer(IntKi) :: jj,kk + real(DbKi) :: acc(4) real(DbKi) :: GridInterp3DVecR8(vDim) ! interpolate - GridInterp3DVecR8 = 0.0_DbKi - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp3DVecR8(vi) = GridInterp3DVecR8(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + do vi = 1,vDim + acc = 0.0_DbKi + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N3D(i,j,k) * data( m%Indx(i,1), jj, kk, vi ) end do end do end do + GridInterp3DVecR8(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp3DVecR8 @@ -533,18 +547,23 @@ function GridInterp3DVec6R4( data, m ) character(*), parameter :: RoutineName = 'GridInterp3DVec6R4' integer(IntKi), parameter :: vDim = 6 integer(IntKi) :: i,j,k,vi + integer(IntKi) :: jj,kk + real(SiKi) :: acc(4) real(SiKi) :: GridInterp3DVec6R4(vDim) ! interpolate - GridInterp3DVec6R4 = 0.0_SiKi - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp3DVec6R4(vi) = GridInterp3DVec6R4(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + do vi = 1,vDim + acc = 0.0_SiKi + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N3D(i,j,k) * data( m%Indx(i,1), jj, kk, vi ) end do end do end do + GridInterp3DVec6R4(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp3DVec6R4 @@ -556,18 +575,23 @@ function GridInterp3DVec6R8( data, m ) character(*), parameter :: RoutineName = 'GridInterp3DVec6R8' integer(IntKi), parameter :: vDim = 6 integer(IntKi) :: i,j,k,vi + integer(IntKi) :: jj,kk + real(DbKi) :: acc(4) real(DbKi) :: GridInterp3DVec6R8(vDim) ! interpolate - GridInterp3DVec6R8 = 0.0_DbKi - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp3DVec6R8(vi) = GridInterp3DVec6R8(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + do vi = 1,vDim + acc = 0.0_DbKi + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N3D(i,j,k) * data( m%Indx(i,1), jj, kk, vi ) end do end do end do + GridInterp3DVec6R8(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp3DVec6R8 @@ -583,19 +607,24 @@ function GridInterp4DR4( data, m ) character(*), parameter :: RoutineName = 'GridInterp4DR4' real(SiKi) :: GridInterp4DR4 - integer(IntKi) :: i,j,k,l + real(SiKi) :: acc(4) + integer(IntKi) :: i,j,k,l,jj,kk,ll ! interpolate - GridInterp4DR4 = 0.0_SiKi + acc = 0.0_SiKi do l = 1,4 + ll = m%Indx(l,4) do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - GridInterp4DR4 = GridInterp4DR4 + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4) ) + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll ) end do end do end do end do + GridInterp4DR4 = acc(1) + acc(2) + acc(3) + acc(4) end function GridInterp4DR4 @@ -605,19 +634,24 @@ function GridInterp4DR8( data, m ) character(*), parameter :: RoutineName = 'GridInterp4DR8' real(DbKi) :: GridInterp4DR8 - integer(IntKi) :: i,j,k,l + real(DbKi) :: acc(4) + integer(IntKi) :: i,j,k,l,jj,kk,ll ! interpolate - GridInterp4DR8 = 0.0_DbKi + acc = 0.0_DbKi do l = 1,4 + ll = m%Indx(l,4) do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - GridInterp4DR8 = GridInterp4DR8 + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4) ) + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll ) end do end do end do end do + GridInterp4DR8 = acc(1) + acc(2) + acc(3) + acc(4) end function GridInterp4DR8 @@ -633,20 +667,26 @@ function GridInterp4DVecR4( data, m ) character(*), parameter :: RoutineName = 'GridInterp4DVecR4' integer(IntKi), parameter :: vDim = 3 integer(IntKi) :: i,j,k,l,vi + integer(IntKi) :: jj,kk,ll + real(SiKi) :: acc(4) real(SiKi) :: GridInterp4DVecR4(vDim) ! interpolate - GridInterp4DVecR4 = 0.0_SiKi - do l = 1,4 - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp4DVecR4(vi) = GridInterp4DVecR4(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + do vi = 1,vDim + acc = 0.0_SiKi + do l = 1,4 + ll = m%Indx(l,4) + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll, vi ) end do end do end do end do + GridInterp4DVecR4(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp4DVecR4 @@ -658,20 +698,26 @@ function GridInterp4DVecR8( data, m ) character(*), parameter :: RoutineName = 'GridInterp4DVecR8' integer(IntKi), parameter :: vDim = 3 integer(IntKi) :: i,j,k,l,vi + integer(IntKi) :: jj,kk,ll + real(DbKi) :: acc(4) real(DbKi) :: GridInterp4DVecR8(vDim) ! interpolate - GridInterp4DVecR8 = 0.0_DbKi - do l = 1,4 - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp4DVecR8(vi) = GridInterp4DVecR8(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + do vi = 1,vDim + acc = 0.0_DbKi + do l = 1,4 + ll = m%Indx(l,4) + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll, vi ) end do end do end do end do + GridInterp4DVecR8(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp4DVecR8 @@ -688,20 +734,26 @@ function GridInterp4DVec6R4( data, m ) character(*), parameter :: RoutineName = 'GridInterp4DVec6R4' integer(IntKi), parameter :: vDim = 6 integer(IntKi) :: i,j,k,l,vi + integer(IntKi) :: jj,kk,ll + real(SiKi) :: acc(4) real(SiKi) :: GridInterp4DVec6R4(vDim) ! interpolate - GridInterp4DVec6R4 = 0.0_SiKi - do l = 1,4 - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp4DVec6R4(vi) = GridInterp4DVec6R4(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + do vi = 1,vDim + acc = 0.0_SiKi + do l = 1,4 + ll = m%Indx(l,4) + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll, vi ) end do end do end do end do + GridInterp4DVec6R4(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp4DVec6R4 @@ -713,20 +765,26 @@ function GridInterp4DVec6R8( data, m ) character(*), parameter :: RoutineName = 'GridInterp4DVec6R8' integer(IntKi), parameter :: vDim = 6 integer(IntKi) :: i,j,k,l,vi + integer(IntKi) :: jj,kk,ll + real(DbKi) :: acc(4) real(DbKi) :: GridInterp4DVec6R8(vDim) ! interpolate - GridInterp4DVec6R8 = 0.0_DbKi - do l = 1,4 - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp4DVec6R8(vi) = GridInterp4DVec6R8(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + do vi = 1,vDim + acc = 0.0_DbKi + do l = 1,4 + ll = m%Indx(l,4) + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll, vi ) end do end do end do end do + GridInterp4DVec6R8(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp4DVec6R8 @@ -743,20 +801,26 @@ function GridInterp4DVecNR4( vDim, data, m ) character(*), parameter :: RoutineName = 'GridInterp4DVecNR4' integer(IntKi) :: i,j,k,l,vi + integer(IntKi) :: jj,kk,ll + real(SiKi) :: acc(4) real(SiKi) :: GridInterp4DVecNR4(vDim) ! interpolate - GridInterp4DVecNR4 = 0.0_SiKi - do l = 1,4 - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp4DVecNR4(vi) = GridInterp4DVecNR4(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + do vi = 1,vDim + acc = 0.0_SiKi + do l = 1,4 + ll = m%Indx(l,4) + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll, vi ) end do end do end do end do + GridInterp4DVecNR4(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp4DVecNR4 @@ -768,20 +832,26 @@ function GridInterp4DVecNR8( vDim, data, m ) character(*), parameter :: RoutineName = 'GridInterp4DVecNR8' integer(IntKi) :: i,j,k,l,vi + integer(IntKi) :: jj,kk,ll + real(DbKi) :: acc(4) real(DbKi) :: GridInterp4DVecNR8(vDim) ! interpolate - GridInterp4DVecNR8 = 0.0_DbKi - do l = 1,4 - do k = 1,4 - do j = 1,4 - do i = 1,4 - do vi = 1,vDim - GridInterp4DVecNR8(vi) = GridInterp4DVecNR8(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + do vi = 1,vDim + acc = 0.0_DbKi + do l = 1,4 + ll = m%Indx(l,4) + do k = 1,4 + kk = m%Indx(k,3) + do j = 1,4 + jj = m%Indx(j,2) + do i = 1,4 + acc(i) = acc(i) + m%N4D(i,j,k,l) * data( m%Indx(i,1), jj, kk, ll, vi ) end do end do end do end do + GridInterp4DVecNR8(vi) = acc(1) + acc(2) + acc(3) + acc(4) end do end function GridInterp4DVecNR8 @@ -799,21 +869,25 @@ function GridInterpNR4( data, p, m ) character(*), parameter :: RoutineName = 'GridInterpNR4' real(SiKi) :: GridInterpNR4(3) real(SiKi) :: dZetadx, dZetady - integer(IntKi) :: i,j,k + real(SiKi) :: accx(4), accy(4), d + integer(IntKi) :: i,j,k,jj,kk ! interpolate slope - dZetadx = 0.0_SiKi - dZetady = 0.0_SiKi + accx = 0.0_SiKi + accy = 0.0_SiKi do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - dZetadx = dZetadx + m%N4D(i,j,k,1) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - dZetady = dZetady + m%N4D(i,j,k,2) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + d = data( m%Indx(i,1), jj, kk ) + accx(i) = accx(i) + m%N4D(i,j,k,1) * d + accy(i) = accy(i) + m%N4D(i,j,k,2) * d end do end do end do - dZetadx = dZetadx / p%delta(2) - dZetady = dZetady / p%delta(3) + dZetadx = ( accx(1) + accx(2) + accx(3) + accx(4) ) / p%delta(2) + dZetady = ( accy(1) + accy(2) + accy(3) + accy(4) ) / p%delta(3) GridInterpNR4 = [-dZetadx,-dZetady,1.0_SiKi] GridInterpNR4 = GridInterpNR4 / TwoNorm(GridInterpNR4) @@ -828,21 +902,25 @@ function GridInterpNR8( data, p, m ) character(*), parameter :: RoutineName = 'GridInterpNR8' real(DbKi) :: GridInterpNR8(3) real(DbKi) :: dZetadx, dZetady - integer(IntKi) :: i,j,k + real(DbKi) :: accx(4), accy(4), d + integer(IntKi) :: i,j,k,jj,kk ! interpolate slope - dZetadx = 0.0_DbKi - dZetady = 0.0_DbKi + accx = 0.0_DbKi + accy = 0.0_DbKi do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - dZetadx = dZetadx + m%N4D(i,j,k,1) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - dZetady = dZetady + m%N4D(i,j,k,2) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + d = data( m%Indx(i,1), jj, kk ) + accx(i) = accx(i) + m%N4D(i,j,k,1) * d + accy(i) = accy(i) + m%N4D(i,j,k,2) * d end do end do end do - dZetadx = dZetadx / p%delta(2) - dZetady = dZetady / p%delta(3) + dZetadx = ( accx(1) + accx(2) + accx(3) + accx(4) ) / p%delta(2) + dZetady = ( accy(1) + accy(2) + accy(3) + accy(4) ) / p%delta(3) GridInterpNR8 = (/-dZetadx,-dZetady,1.0_DbKi/) GridInterpNR8 = GridInterpNR8 / TwoNorm(GridInterpNR8) @@ -861,19 +939,25 @@ function GridInterpSR4( data, p, m ) character(*), parameter :: RoutineName = 'GridInterpSR4' real(SiKi) :: GridInterpSR4(2) - integer(IntKi) :: i,j,k,dir + real(SiKi) :: accx(4), accy(4), d + integer(IntKi) :: i,j,k,jj,kk ! interpolate slope - GridInterpSR4 = 0.0_SiKi + accx = 0.0_SiKi + accy = 0.0_SiKi do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - do dir = 1,2 - GridInterpSR4(dir) = GridInterpSR4(dir) + m%N4D(i,j,k,dir) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - end do + d = data( m%Indx(i,1), jj, kk ) + accx(i) = accx(i) + m%N4D(i,j,k,1) * d + accy(i) = accy(i) + m%N4D(i,j,k,2) * d end do end do end do + GridInterpSR4(1) = accx(1) + accx(2) + accx(3) + accx(4) + GridInterpSR4(2) = accy(1) + accy(2) + accy(3) + accy(4) GridInterpSR4 = GridInterpSR4 / p%delta(2:3) end function GridInterpSR4 @@ -885,19 +969,25 @@ function GridInterpSR8( data, p, m ) character(*), parameter :: RoutineName = 'GridInterpSR8' real(DbKi) :: GridInterpSR8(2) - integer(IntKi) :: i,j,k,dir + real(DbKi) :: accx(4), accy(4), d + integer(IntKi) :: i,j,k,jj,kk ! interpolate slope - GridInterpSR8 = 0.0_DbKi + accx = 0.0_DbKi + accy = 0.0_DbKi do k = 1,4 + kk = m%Indx(k,3) do j = 1,4 + jj = m%Indx(j,2) do i = 1,4 - do dir = 1,2 - GridInterpSR8(dir) = GridInterpSR8(dir) + m%N4D(i,j,k,dir) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - end do + d = data( m%Indx(i,1), jj, kk ) + accx(i) = accx(i) + m%N4D(i,j,k,1) * d + accy(i) = accy(i) + m%N4D(i,j,k,2) * d end do end do end do + GridInterpSR8(1) = accx(1) + accx(2) + accx(3) + accx(4) + GridInterpSR8(2) = accy(1) + accy(2) + accy(3) + accy(4) GridInterpSR8 = GridInterpSR8 / p%delta(2:3) end function GridInterpSR8 diff --git a/modules/nwtc-library/src/SysFlangLinux.f90 b/modules/nwtc-library/src/SysFlangLinux.f90 index 8805951090..ee97ce4702 100644 --- a/modules/nwtc-library/src/SysFlangLinux.f90 +++ b/modules/nwtc-library/src/SysFlangLinux.f90 @@ -54,9 +54,9 @@ MODULE SysSubs MODULE PROCEDURE NWTC_gammaR8 END INTERFACE - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 6 ! The I/O unit for the console. Unit 6 causes ADAMS to crash. - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .TRUE. ! A flag to tell the program that keyboard input is allowed in the environment. CHARACTER(*), PARAMETER :: NewLine = ACHAR(10) ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}] CHARACTER(*), PARAMETER :: OS_Desc = 'GNU Fortran for Linux' ! Description of the language/OS diff --git a/modules/nwtc-library/src/SysGnuLinux.f90 b/modules/nwtc-library/src/SysGnuLinux.f90 index 6eddb9b410..44d7fe73ea 100644 --- a/modules/nwtc-library/src/SysGnuLinux.f90 +++ b/modules/nwtc-library/src/SysGnuLinux.f90 @@ -55,9 +55,9 @@ MODULE SysSubs MODULE PROCEDURE NWTC_gammaR8 END INTERFACE - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 6 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .TRUE. ! A flag to tell the program that keyboard input is allowed in the environment. CHARACTER(*), PARAMETER :: NewLine = ACHAR(10) ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}] CHARACTER(*), PARAMETER :: OS_Desc = 'GNU Fortran for Linux' ! Description of the language/OS diff --git a/modules/nwtc-library/src/SysGnuWin.f90 b/modules/nwtc-library/src/SysGnuWin.f90 index 1422989cb2..cad5fda7d4 100644 --- a/modules/nwtc-library/src/SysGnuWin.f90 +++ b/modules/nwtc-library/src/SysGnuWin.f90 @@ -55,9 +55,9 @@ MODULE SysSubs MODULE PROCEDURE NWTC_gammaR8 END INTERFACE - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 6 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .TRUE. ! A flag to tell the program that keyboard input is allowed in the environment. CHARACTER(*), PARAMETER :: NewLine = ACHAR(10) ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}] CHARACTER(*), PARAMETER :: OS_Desc = 'GNU Fortran for Windows' ! Description of the language/OS diff --git a/modules/nwtc-library/src/SysIFL.f90 b/modules/nwtc-library/src/SysIFL.f90 index bf427e8fb6..17e7db0e28 100644 --- a/modules/nwtc-library/src/SysIFL.f90 +++ b/modules/nwtc-library/src/SysIFL.f90 @@ -55,9 +55,9 @@ MODULE SysSubs MODULE PROCEDURE NWTC_gammaR8 END INTERFACE - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 6 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .TRUE. ! A flag to tell the program that keyboard input is allowed in the environment. CHARACTER(*), PARAMETER :: NewLine = ACHAR(10) ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}] CHARACTER(*), PARAMETER :: OS_Desc = 'Intel Fortran for Linux' ! Description of the language/OS diff --git a/modules/nwtc-library/src/SysIVF.f90 b/modules/nwtc-library/src/SysIVF.f90 index 41af20386a..0286de427d 100644 --- a/modules/nwtc-library/src/SysIVF.f90 +++ b/modules/nwtc-library/src/SysIVF.f90 @@ -20,19 +20,19 @@ MODULE SysSubs ! This module contains routines with system-specific logic and references, including all references to the console unit, CU. ! It also contains standard (but not system-specific) routines it uses. - ! SysIVF.f90 is specifically for the Intel Visual Fortran for Windows compiler. + ! SysIVF.f90 is specifically for the Intel Visual Fortran for Windows compiler. ! It contains the following routines: ! FUNCTION FileSize( Unit ) ! Returns the size (in bytes) of an open file. ! FUNCTION Is_NaN( DblNum ) ! Please use IEEE_IS_NAN() instead ! FUNCTION NWTC_ERF( x ) - ! FUNCTION NWTC_gamma( x ) ! Returns the gamma value of its argument. + ! FUNCTION NWTC_gamma( x ) ! Returns the gamma value of its argument. ! SUBROUTINE FlushOut ( Unit ) ! SUBROUTINE GET_CWD( DirName, Status ) ! SUBROUTINE MKDIR( new_directory_path ) ! SUBROUTINE OpenCon ! SUBROUTINE OpenUnfInpBEFile ( Un, InFile, RecLen, Error ) ! SUBROUTINE ProgExit ( StatCode ) - ! SUBROUTINE Set_IEEE_Constants( NaN_D, Inf_D, NaN, Inf, NaN_S, Inf_S ) + ! SUBROUTINE Set_IEEE_Constants( NaN_D, Inf_D, NaN, Inf, NaN_S, Inf_S ) ! SUBROUTINE UsrAlarm ! SUBROUTINE WrNR ( Str ) ! SUBROUTINE WrOver ( Str ) @@ -43,11 +43,11 @@ MODULE SysSubs USE NWTC_Base IMPLICIT NONE - + INTERFACE NWTC_ERF ! Returns the ERF value of its argument MODULE PROCEDURE NWTC_ERFR4 MODULE PROCEDURE NWTC_ERFR8 - END INTERFACE + END INTERFACE INTERFACE NWTC_gamma ! Returns the gamma value of its argument ! note: gamma is part of the F08 standard, but may not be implemented everywhere... @@ -55,9 +55,9 @@ MODULE SysSubs MODULE PROCEDURE NWTC_gammaR8 END INTERFACE - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 7 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .TRUE. ! A flag to tell the program that keyboard input is allowed in the environment. CHARACTER(*), PARAMETER :: NewLine = ACHAR(10) ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}]; Note: NewLine change to ACHAR(10) here on Windows to fix issues with C/Fortran interoperability using WrScr CHARACTER(*), PARAMETER :: OS_Desc = 'Intel Visual Fortran for Windows'! Description of the language/OS @@ -115,20 +115,20 @@ FUNCTION Is_NaN( DblNum ) RETURN END FUNCTION Is_NaN ! ( DblNum ) !======================================================================= -!> This function returns the ERF value of its argument. The result has a value equal -!! to the Gauss error function: +!> This function returns the ERF value of its argument. The result has a value equal +!! to the Gauss error function: !! \f{equation}{ !! \mathrm{erf}(x)=\frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}\,dt !! \f} \n !! Use NWTC_ERF (syssubs::nwtc_erf) instead of directly calling a specific routine in the generic interface. FUNCTION NWTC_ERFR4( x ) - ! Returns the ERF value of its argument. The result has a value equal - ! to the error function: 2/pi * integral_from_0_to_x of e^(-t^2) dt. + ! Returns the ERF value of its argument. The result has a value equal + ! to the error function: 2/pi * integral_from_0_to_x of e^(-t^2) dt. - REAL(SiKi), INTENT(IN) :: x !< input + REAL(SiKi), INTENT(IN) :: x !< input REAL(SiKi) :: NWTC_ERFR4 !< \f$\mathrm{erf}(x)\f$ - + NWTC_ERFR4 = ERF( x ) END FUNCTION NWTC_ERFR4 @@ -136,14 +136,14 @@ END FUNCTION NWTC_ERFR4 !> \copydoc syssubs::nwtc_erfr4 FUNCTION NWTC_ERFR8( x ) - REAL(R8Ki), INTENT(IN) :: x ! input + REAL(R8Ki), INTENT(IN) :: x ! input REAL(R8Ki) :: NWTC_ERFR8 ! this function - + NWTC_ERFR8 = ERF( x ) END FUNCTION NWTC_ERFR8 !======================================================================= -!> Returns the gamma value of its argument. The result has a value equal +!> Returns the gamma value of its argument. The result has a value equal !! to a processor-dependent approximation to the gamma function of x: !! \f{equation}{ !! \mathrm{ \Gamma }(x) = \int_0^\infty t^{x-1}e^{-t}\,dt @@ -153,9 +153,9 @@ FUNCTION NWTC_GammaR4( x ) ! \mathrm{ \Gamma }(x) = \int_0^\infinity t^{x-1}e^{-t}\,dt - REAL(SiKi), INTENT(IN) :: x !< input + REAL(SiKi), INTENT(IN) :: x !< input REAL(SiKi) :: NWTC_GammaR4 !< \f$\mathrm{\Gamma}(x)\f$ - + NWTC_GammaR4 = gamma( x ) END FUNCTION NWTC_GammaR4 @@ -163,9 +163,9 @@ END FUNCTION NWTC_GammaR4 !> \copydoc syssubs::nwtc_gammar4 FUNCTION NWTC_GammaR8( x ) - REAL(R8Ki), INTENT(IN) :: x ! input + REAL(R8Ki), INTENT(IN) :: x ! input REAL(R8Ki) :: NWTC_GammaR8 ! result - + NWTC_GammaR8 = gamma( x ) END FUNCTION NWTC_GammaR8 @@ -230,7 +230,7 @@ SUBROUTINE OpenCon() END SUBROUTINE OpenCon !======================================================================= !> This routine opens a binary input file with data stored in Big Endian format (created on a UNIX machine.) -!! Data are stored in RecLen-byte records. (This routine is used for +!! Data are stored in RecLen-byte records. (This routine is used for SUBROUTINE OpenUnfInpBEFile ( Un, InFile, RecLen, Error ) IMPLICIT NONE @@ -265,7 +265,7 @@ SUBROUTINE ProgExit ( StatCode ) INTEGER, INTENT(IN) :: StatCode ! The status code to pass to the OS. CALL EXIT ( StatCode ) - + ! IF ( StatCode == 0 ) THEN ! STOP 0 ! ELSE @@ -276,12 +276,12 @@ SUBROUTINE ProgExit ( StatCode ) ! END IF END SUBROUTINE ProgExit ! ( StatCode ) !======================================================================= -!> This routine sets the values of NaN_D, Inf_D, NaN, Inf (IEEE -!! values for not-a-number and infinity in sindle and double -!! precision) This uses standard F03 intrinsic routines, +!> This routine sets the values of NaN_D, Inf_D, NaN, Inf (IEEE +!! values for not-a-number and infinity in single and double +!! precision) This uses standard F03 intrinsic routines, !! however Gnu has not yet implemented it, so we've placed this !! routine in the system-specific code. -SUBROUTINE Set_IEEE_Constants( NaN_D, Inf_D, NaN, Inf, NaN_S, Inf_S ) +SUBROUTINE Set_IEEE_Constants( NaN_D, Inf_D, NaN, Inf, NaN_S, Inf_S ) USE, INTRINSIC :: ieee_arithmetic ! use this for compilers that have implemented ieee_arithmetic from F03 standard (otherwise see logic in SysGnu*.f90) @@ -290,7 +290,7 @@ SUBROUTINE Set_IEEE_Constants( NaN_D, Inf_D, NaN, Inf, NaN_S, Inf_S ) REAL(ReKi), INTENT(inout) :: Inf !< IEEE value for NaN (not-a-number) REAL(ReKi), INTENT(inout) :: NaN !< IEEE value for Inf (infinity) - + REAL(SiKi), INTENT(inout) :: Inf_S !< IEEE value for NaN (not-a-number) in single precision REAL(SiKi), INTENT(inout) :: NaN_S !< IEEE value for Inf (infinity) in single precision @@ -298,12 +298,12 @@ SUBROUTINE Set_IEEE_Constants( NaN_D, Inf_D, NaN, Inf, NaN_S, Inf_S ) Inf_D = ieee_value(0.0_DbKi, ieee_positive_inf) NaN = ieee_value(0.0_ReKi, ieee_quiet_nan) - Inf = ieee_value(0.0_ReKi, ieee_positive_inf) + Inf = ieee_value(0.0_ReKi, ieee_positive_inf) NaN_S = ieee_value(0.0_SiKi, ieee_quiet_nan) Inf_S = ieee_value(0.0_SiKi, ieee_positive_inf) -END SUBROUTINE Set_IEEE_Constants +END SUBROUTINE Set_IEEE_Constants !======================================================================= !> This routine generates an alarm to warn the user that something went wrong. SUBROUTINE UsrAlarm() @@ -327,11 +327,11 @@ END SUBROUTINE WrNR ! ( Str ) SUBROUTINE WrOver ( Str ) CHARACTER(*), INTENT(IN) :: Str !< The string to write to the screen. - + INTEGER :: MaxLen !< maximum length of string to be written to the screen (cannot exceed ConRecL here) ! When the file is opened using CARRIAGECONTROL='FORTRAN', the "+" character allows writing over the previous line. However, the Fortran carriage control has been deleted from the Fortran standard. - + MaxLen = min(ConRecL-1, len(Str)) if (MaxLen > 0) then WRITE (CU,'("+",A)') Str(1:MaxLen) @@ -404,11 +404,11 @@ SUBROUTINE LoadDynamicLibProc ( DLL, ErrStat, ErrMsg ) ErrStat = ErrID_None ErrMsg = '' - + ! Get the procedure addresses: do i=1,NWTC_MAX_DLL_PROC if ( len_trim( DLL%ProcName(i) ) > 0 ) then - + ProcAddr = GetProcAddress( DLL%FileAddr, TRIM(DLL%ProcName(i))//C_NULL_CHAR ) !the "C_NULL_CHAR" converts the Fortran string to a C-type string (i.e., adds //CHAR(0) to the end) DLL%ProcAddr(i) = TRANSFER(ProcAddr, DLL%ProcAddr(i)) !convert INTEGER(LPVOID) to INTEGER(C_FUNPTR) [used only for compatibility with gfortran] @@ -417,7 +417,7 @@ SUBROUTINE LoadDynamicLibProc ( DLL, ErrStat, ErrMsg ) ErrMsg = 'The procedure '//TRIM(DLL%ProcName(i))//' in file '//TRIM(DLL%FileName)//' could not be loaded.' RETURN END IF - + end if end do @@ -435,12 +435,12 @@ SUBROUTINE FreeDynamicLib ( DLL, ErrStat, ErrMsg ) CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None INTEGER(HANDLE) :: FileAddr ! The address of file FileName. (RETURN value from LoadLibrary in kernel32.f90) INTEGER(BOOL) :: Success ! Whether or not the call to FreeLibrary was successful - + ErrStat = ErrID_None ErrMsg = '' - + IF ( DLL%FileAddr == INT(0,C_INTPTR_T) ) RETURN - + FileAddr = TRANSFER(DLL%FileAddr, FileAddr) !convert INTEGER(C_INTPTR_T) to INTEGER(HANDLE) [used only for compatibility with gfortran] ! Free the DLL: diff --git a/modules/nwtc-library/src/SysIVF_Labview.f90 b/modules/nwtc-library/src/SysIVF_Labview.f90 index 51a57144f5..4b2d433723 100644 --- a/modules/nwtc-library/src/SysIVF_Labview.f90 +++ b/modules/nwtc-library/src/SysIVF_Labview.f90 @@ -72,9 +72,9 @@ MODULE SysSubs !======================================================================= - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 7 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .FALSE. ! A flag to tell the program that keyboard input is allowed in the environment. diff --git a/modules/nwtc-library/src/SysMatlabLinuxGnu.f90 b/modules/nwtc-library/src/SysMatlabLinuxGnu.f90 index 1072ca8e3b..e80a4ae599 100644 --- a/modules/nwtc-library/src/SysMatlabLinuxGnu.f90 +++ b/modules/nwtc-library/src/SysMatlabLinuxGnu.f90 @@ -58,9 +58,9 @@ MODULE SysSubs MODULE PROCEDURE NWTC_gammaR8 END INTERFACE - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 6 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .FALSE. ! A flag to tell the program that keyboard input is allowed in the environment. CHARACTER(*), PARAMETER :: NewLine = ACHAR(10) ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}] CHARACTER(*), PARAMETER :: OS_Desc = 'GNU Fortran for Linux with Matlab' ! Description of the language/OS diff --git a/modules/nwtc-library/src/SysMatlabLinuxIntel.f90 b/modules/nwtc-library/src/SysMatlabLinuxIntel.f90 index 303206a9a4..2d55265fcc 100644 --- a/modules/nwtc-library/src/SysMatlabLinuxIntel.f90 +++ b/modules/nwtc-library/src/SysMatlabLinuxIntel.f90 @@ -58,9 +58,9 @@ MODULE SysSubs MODULE PROCEDURE NWTC_gammaR8 END INTERFACE - INTEGER, PARAMETER :: ConRecL = 120 ! The record length for console output. + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 6 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .FALSE. ! A flag to tell the program that keyboard input is allowed in the environment. CHARACTER(*), PARAMETER :: NewLine = ACHAR(10) ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}] CHARACTER(*), PARAMETER :: OS_Desc = 'Intel Fortran for Linux with Matlab' ! Description of the language/OS diff --git a/modules/nwtc-library/src/SysMatlabWindows.f90 b/modules/nwtc-library/src/SysMatlabWindows.f90 index 417a0e9ef6..c082043b8d 100644 --- a/modules/nwtc-library/src/SysMatlabWindows.f90 +++ b/modules/nwtc-library/src/SysMatlabWindows.f90 @@ -47,8 +47,9 @@ MODULE SysSubs !======================================================================= + INTEGER, PARAMETER :: ConRecL = 180 ! The record length for console output (maximum number of characters that can be written in WrOver(), must be larger than MaxWrScrLen) INTEGER, PUBLIC :: CU = 6 ! The I/O unit for the console (Can be changed with SetConsoleUnit subroutine) - INTEGER, PARAMETER :: MaxWrScrLen = 256 ! The maximum number of characters allowed to be written to a line in WrScr + INTEGER, PARAMETER :: MaxWrScrLen = ConRecL-1 ! The maximum number of characters allowed to be written to a line in WrScr, must be smaller than ConRecL LOGICAL, PARAMETER :: KBInputOK = .FALSE. ! A flag to tell the program that keyboard input is allowed in the environment. diff --git a/openfast_io/openfast_io/FAST_writer.py b/openfast_io/openfast_io/FAST_writer.py index 022bc289d0..d97d1304b0 100644 --- a/openfast_io/openfast_io/FAST_writer.py +++ b/openfast_io/openfast_io/FAST_writer.py @@ -619,7 +619,7 @@ def write_ElastoDynBlade(self, bldInd = None): blade_file = os.path.join(self.FAST_runDirectory,self.fst_vt['ElastoDyn']['BldFile1']) else: EDbld_dict = self.fst_vt['ElastoDynBlade'][bldInd] - blade_file = os.path.join(self.FAST_runDirectory,self.fst_vt['ElastoDyn']['BldFile'+(bldInd+1)]) + blade_file = os.path.join(self.FAST_runDirectory,self.fst_vt['ElastoDyn']['BldFile'+ str(bldInd+1)]) f = open(blade_file, 'w') diff --git a/reg_tests/r-test b/reg_tests/r-test index a7372cb071..174a8ef776 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit a7372cb071ec8ecbaf7f367d269e54822430b882 +Subproject commit 174a8ef7764a845782932122ab6ec5992c0f72e9 diff --git a/unit_tests/CMakeLists.txt b/unit_tests/CMakeLists.txt index 00b002f863..c85dae1018 100644 --- a/unit_tests/CMakeLists.txt +++ b/unit_tests/CMakeLists.txt @@ -51,6 +51,7 @@ add_test(NAME beamdyn_utest COMMAND beamdyn_utest) add_executable(inflowwind_utest ${PROJECT_SOURCE_DIR}/modules/inflowwind/tests/inflowwind_utest.F90 ${PROJECT_SOURCE_DIR}/modules/inflowwind/tests/test_bladed_wind.F90 + ${PROJECT_SOURCE_DIR}/modules/inflowwind/tests/test_grid3d_field.F90 ${PROJECT_SOURCE_DIR}/modules/inflowwind/tests/test_hawc_wind.F90 ${PROJECT_SOURCE_DIR}/modules/inflowwind/tests/test_outputs.F90 ${PROJECT_SOURCE_DIR}/modules/inflowwind/tests/test_steady_wind.F90