From d2b4c995a795e851ce370bdddeb6bef1b710073c Mon Sep 17 00:00:00 2001 From: Emmanuel Branlard Date: Sun, 23 Feb 2025 14:14:14 -0500 Subject: [PATCH] SP: adding source panels --- .../ExampleFiles/ExampleFile--OLAF.dat | 6 +- docs/source/user/aerodyn-olaf/InputFiles.rst | 112 +++ modules/aerodyn/src/AeroDyn.f90 | 35 +- modules/aerodyn/src/AeroDyn_Driver_Subs.f90 | 27 +- modules/aerodyn/src/AeroDyn_Inflow.f90 | 47 +- modules/aerodyn/src/FVW.f90 | 126 +++- modules/aerodyn/src/FVW_BiotSavart.f90 | 176 ++++- modules/aerodyn/src/FVW_IO.f90 | 81 ++- modules/aerodyn/src/FVW_Registry.txt | 73 +- modules/aerodyn/src/FVW_Subs.f90 | 576 ++++++++++++++-- modules/aerodyn/src/FVW_Tests.f90 | 646 ++++++++++-------- modules/aerodyn/src/FVW_Types.f90 | 646 +++++++++++++++++- modules/aerodyn/src/FVW_VortexTools.f90 | 51 +- modules/aerodyn/tests/test_FVW_testsuite.F90 | 42 ++ modules/nwtc-library/src/VTK.f90 | 321 ++++++++- 15 files changed, 2513 insertions(+), 452 deletions(-) create mode 100644 modules/aerodyn/tests/test_FVW_testsuite.F90 diff --git a/docs/source/user/aerodyn-olaf/ExampleFiles/ExampleFile--OLAF.dat b/docs/source/user/aerodyn-olaf/ExampleFiles/ExampleFile--OLAF.dat index bd9d6a46d3..4fdfb0ad97 100644 --- a/docs/source/user/aerodyn-olaf/ExampleFiles/ExampleFile--OLAF.dat +++ b/docs/source/user/aerodyn-olaf/ExampleFiles/ExampleFile--OLAF.dat @@ -43,4 +43,8 @@ default VTK_fps - Frame rate for VTK output (frames per second) {"all 0 nGridOut - Number of grid outputs GridName GridType TStart TEnd DTGrid XStart XEnd nX YStart YEnd nY ZStart ZEnd nZ (-) (-) (s) (s) (s) (m) (m) (-) (m) (m) (-) (m) (m) (-) ------------------------------------------------------------------------------------------------- +=============================================================================================== +--------------------------- ADVANCED OPTIONS -------------------------------------------------- +=============================================================================================== +! Advanced options may be placed here in arbitrary order, with the regular format: +! Value Key - Comment diff --git a/docs/source/user/aerodyn-olaf/InputFiles.rst b/docs/source/user/aerodyn-olaf/InputFiles.rst index a841fbf814..e67d3d3592 100644 --- a/docs/source/user/aerodyn-olaf/InputFiles.rst +++ b/docs/source/user/aerodyn-olaf/InputFiles.rst @@ -312,6 +312,118 @@ 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. +Advanced Options +~~~~~~~~~~~~~~~~ + + +Advanced options (typically used for developers or beta features) can be placed at the end of the OLAF input file: + + +- These options use the regular format: `Value Key - Comment`. +- They can be placed in arbitrary order. +- They can be commented out with a `!` character at the beginning of the line. +- Blank lines or unsupported options are ignored (display to screen). +- If not provided their default value is being used, but the keyword "default" is not available. +- We recall that these are **beta features** and should mostly be used by developers. + + +The end of the OLAF input file would look as follows: + +.. code:: + + [...] + GridName GridType TStart TEnd DTGrid XStart XEnd nX YStart YEnd nY ZStart ZEnd nZ + (-) (-) (s) (s) (s) (m) (m) (-) (m) (m) (-) (m) (m) (-) + =============================================================================================== + --------------------------- ADVANCED OPTIONS -------------------------------------------------- + =============================================================================================== + ! Advanced options may be placed here in arbitrary order, with the regular format: + ! Value1 Key1 - Comment1 + ! Value2 Key2 - Comment2 + ! Lines starting with `!` are ignored, empty lines are ignored + [...] etc + +Currently, the advanced options supported as as follows: + +.. code:: + + =============================================================================================== + --------------------------- ADVANCED OPTIONS -------------------------------------------------- + =============================================================================================== + "Panels.vtk" SrcPnlFile - Name of VTK file containing source panels {default: ""} + 1 nSrcPnlUpdate - How often do src panel updates (in time steps of OLAF), {default: 1} + True Induction - Compute induction, {default: True} + True InductionAtCP - Compute induced velocities at nodes or CP, {default: True} + True WakeAtTE - Start the wake at the trailing edge, or at the LL, {default: True} + False DStallOnWake - Dynamic stall has influence on wake, {default: False} + 0.75 kFrozenNWStart - Fraction of wake induced velocity at start of frozen wake, {default: 0.75} + 0.5 kFrozenNWEnd - Fraction of wake induced velocity at end of frozen wake, {default: 0.5} + 0.0 zGround - Ground height {default: 0.0} + 0.1 zGroundPush - Ground push back {default: 0.1} + + + +**SrcPanlFile** [string] specifies the name of the VTK file used for the source panel method. +The VTK file should be a legacy ASCII VTK file, with DATASET POLYDATA POINTS and POLYGONS. +A sample VTK file is provided below for two panels forming a regular grid in the XY plane. +The connectivity of the polygons needs to be such that the normal points away from the body and +into the fluid. In the example below, the normals are in the z direction which would be typical +for a bottom wall with fluid above. +When use **WrVTK**, OLAF will write separate VTK files containing various information related to the source panels, such as the pressure coefficient, area, force per area, normals. Looking at the orientation of the normals is extremely important (In Paraview, select 3D Glyph, Arrows, and select the `Normals` output, scaled by the `Normals`). +The BodyID CELL_DATA can be used to separate different patches, which can help the post processing. This is optional, OLAF doesn't use it currently, but it is written back in the output files of OLAF. + +Curious readers may look at the unittest `Test_SrcPnl_Sphere` of OLAF that test the pressure coefficient +about a sphere. + +.. code:: + + # vtk DataFile Version 2.0 + Comment + ASCII + DATASET POLYDATA + POINTS 6 double + 0.0 0.0 0.0 + 0.0 25.0 0.0 + 0.0 50.0 0.0 + 10.0 0.0 0.0 + 10.0 25.0 0.0 + 10.0 50.0 0.0 + + POLYGONS 2 10 + 4 0 1 4 3 + 4 1 2 5 4 + + CELL_DATA 2 + SCALARS BodyID int + LOOKUP_TABLE default + 0 + 0 + + + +**nSrcPnlUpdate** [int] Define how often the source panel updates (in time steps of OLAF), the default is `1`, at each time step. + +**Induction** [switch] Compute induced velocities, otherwise all induced velocity are 0 (no wake, no panels, etc) ! Default is `True`. + +**InductionAtCP** [switch] When performing the lifting line calculations, compute induced velocities at nodes (False) or at control points (CP, True). Default is `True`. + +**WakeAtTE** [switch] Start the wake at the trailing edge (True), or directly at the LL (False, no chordwise panel). Default is `True`. + +**DStallOnWake** [switch] Include the influence of the dynamic stall on the wake (True), i.e. use the dynamic Cl to update the circulation of the lifting line and wake panel. Default is `False`. + +**kFrozenNWStart** [float] Fraction of wake induced velocity at start of frozen wake, default is `0.75`. See OLAF theory related to Frozen NW, :numref:`sec:vortconvfrozen`. + +**kFrozenNWEnd** [float] Fraction of wake induced velocity at end of frozen wake, default is `0.5`. See OLAF theory related to Frozen NW, :numref:`sec:vortconvfrozen`. + +**zGround** [float] Height below which vortex points are not allowed to be and if any are present they will be pushed back above the ground at the height defined by **zGroundPush**. Default is `0.0`. +For MHK, the sea bed location, based on the water depth is added to **zGround**. + +**zGroundPush** [float] Ground push back, see **zGround**. Default is `0.1`. + + + + + AeroDyn Input File -------------------- Input file modifications diff --git a/modules/aerodyn/src/AeroDyn.f90 b/modules/aerodyn/src/AeroDyn.f90 index c1324965ef..70bf290949 100644 --- a/modules/aerodyn/src/AeroDyn.f90 +++ b/modules/aerodyn/src/AeroDyn.f90 @@ -385,7 +385,12 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut ! set the rest of the parameters - p%Skew_Mod = InputFileData%Skew_Mod + ! NOTE: some parameters need to be set even if no rotors + p%Skew_Mod = InputFileData%Skew_Mod + p%UA_Flag = InputFileData%UA_Init%UAMod > UA_None + p%CompAeroMaps = InitInp%CompAeroMaps + p%DT = InputFileData%DTAero + p%Wake_Mod = InputFileData%Wake_Mod do iR = 1, nRotors p%rotors(iR)%AeroProjMod = AeroProjMod(iR) call WrScr(' AeroDyn: projMod: '//trim(num2lstr(p%rotors(iR)%AeroProjMod))) @@ -1444,9 +1449,8 @@ subroutine SetParameters( InitInp, InputFileData, RotData, p, p_AD, ErrStat, Err p%MHK = InitInp%MHK - p_AD%DT = InputFileData%DTAero - p_AD%Wake_Mod = InputFileData%Wake_Mod p%DBEMT_Mod = InputFileData%DBEMT_Mod + p%TwrPotent = InputFileData%TwrPotent p%TwrShadow = InputFileData%TwrShadow p%TwrAero = InputFileData%TwrAero @@ -1736,9 +1740,11 @@ subroutine AD_End( u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg ) if (p%Wake_Mod == WakeMod_FVW ) then if ( p%UA_Flag ) then - do iW=1,p%FVW%nWings - call UA_End(m%FVW%W(iW)%p_UA) - enddo + if (allocated(m%FVW%W)) then + do iW=1,p%FVW%nWings + call UA_End(m%FVW%W(iW)%p_UA) + enddo + endif end if call FVW_End( m%FVW_u, p%FVW, x%FVW, xd%FVW, z%FVW, OtherState%FVW, m%FVW_y, m%FVW, ErrStat, ErrMsg ) @@ -5096,11 +5102,25 @@ SUBROUTINE Init_OLAF( InputFileData, u_AD, u, p, x, xd, z, OtherState, m, ErrSta allocate(InitInp%W(nWings) , STAT = ErrStat2); ErrMsg2='Allocate W'; if(Failed()) return allocate(InitInp%WingsMesh(nWings), STAT = ErrStat2); ErrMsg2='Allocate Wings Mesh'; if(Failed()) return + ! --- Default inputs (not per retor) + ! UA and inputs that are not per rotor + InitInp%UA_Flag = p%UA_Flag + ! Important inputs if no rotors are present! + if (size(p%rotors)==0) then + InitInp%numBladeNodes = 0 + InitInp%AirDens = InputFileData%AirDens + InitInp%KinVisc = InputFileData%KinVisc + InitInp%MHK = 0 ! TODO this is an initinp of AeroDyn + InitInp%WtrDpth = 0.0_ReKi ! TODO this is an initinp of AeroDyn + InitInp%RootName = p%RootName(1:len_trim(p%RootName)-2) ! Removing "AD" + endif + ! --- Inputs per wings/blades iW_incr=0 do iR=1, size(p%rotors) InitInp%numBladeNodes = p%rotors(iR)%numBlNds ! TODO TODO TODO per wing + InitInp%AirDens = p%rotors(iR)%AirDens InitInp%KinVisc = p%rotors(iR)%KinVisc InitInp%MHK = p%rotors(iR)%MHK InitInp%WtrDpth = p%rotors(iR)%WtrDpth @@ -5163,9 +5183,6 @@ SUBROUTINE Init_OLAF( InputFileData, u_AD, u, p, x, xd, z, OtherState, m, ErrSta if(Failed()) return enddo ! iB, blades - - ! Unsteady Aero Data - InitInp%UA_Flag = p%UA_Flag call UA_CopyInitInput(InputFileData%UA_Init, InitInp%UA_Init, MESH_NEWCOPY, ErrStat2, ErrMsg2) iW_incr = iW_incr+p%rotors(iR)%numBlades diff --git a/modules/aerodyn/src/AeroDyn_Driver_Subs.f90 b/modules/aerodyn/src/AeroDyn_Driver_Subs.f90 index 844e6c8d89..b93f0760f6 100644 --- a/modules/aerodyn/src/AeroDyn_Driver_Subs.f90 +++ b/modules/aerodyn/src/AeroDyn_Driver_Subs.f90 @@ -540,7 +540,7 @@ subroutine Init_ADI_ForDriver(iCase, ADI, dvr, FED, dt, needInitIW, errStat, err InitInp%AD%defPvap = dvr%Pvap InitInp%AD%WtrDpth = dvr%WtrDpth InitInp%AD%MSL2SWL = dvr%MSL2SWL - InitInp%AD%CompSeaSt = dvr%SS_InitInp%CompSeaSt + InitInp%AD%CompSeaSt = dvr%SS_InitInp%CompSeaSt /=0 ! Init data per rotor allocate(InitInp%AD%rotors(dvr%numTurbines), stat=errStat) if (errStat/=0) then @@ -1037,7 +1037,7 @@ subroutine Dvr_ReadInputFile(fileName, dvr, errStat, errMsg ) call ParseCom(FileInfo_In, CurLine, Line, errStat2, errMsg2, unEc); if (Failed()) return call ParseVar(FileInfo_In, CurLine, "compInflow", dvr%IW_InitInp%compInflow , errStat2, errMsg2, unEc); if (Failed()) return call ParseVar(FileInfo_In, CurLine, "InflowFile", dvr%IW_InitInp%InputFile, errStat2, errMsg2, unEc, IsPath=.true.); if (Failed()) return - if (dvr%IW_InitInp%compInflow==0) then + if (dvr%IW_InitInp%compInflow/=1) then call ParseVar(FileInfo_In, CurLine, "HWindSpeed", dvr%IW_InitInp%HWindSpeed , errStat2, errMsg2, unEc); if (Failed()) return call ParseVar(FileInfo_In, CurLine, "RefHt" , dvr%IW_InitInp%RefHt , errStat2, errMsg2, unEc); if (Failed()) return call ParseVar(FileInfo_In, CurLine, "PLExp" , dvr%IW_InitInp%PLExp , errStat2, errMsg2, unEc); if (Failed()) return @@ -1543,6 +1543,9 @@ subroutine Dvr_InitializeOutputs(nWT, out, numSteps, errStat, errMsg) end do ! i write (out%unOutFile(iWT),'()') enddo + if (nWT==0) then + ! Special case, no turbine, TODO + endif endif ! --- Binary @@ -1744,7 +1747,11 @@ subroutine Dvr_WriteOutputs(nt, t, dvr, out, yADI, SeaSt, errStat, errMsg) END IF ! Packing all outputs excpet time into one array - nAD = size(yADI%AD%rotors(1)%WriteOutput) + if (size(yADI%AD%rotors)>=1) then + nAD = size(yADI%AD%rotors(1)%WriteOutput) + else + nAD = 0 + endif nIW = size(yADI%IW_WriteOutput) nSS = size(SeaSt%y%WriteOutput) nDV = out%nDvrOutputs @@ -1872,11 +1879,15 @@ subroutine setVTKParameters(DVR_Outs, dvr, ADI, errStat, errMsg, dirname) ! Get radius for ground (blade length + hub radius): GroundRad = MaxBladeLength + MaxTwrLength+ DVR_Outs%VTKHubRad ! write the ground or seabed reference polygon: - RefPoint(1:2) = dvr%WT(1)%originInit(1:2) - do iWT=2,dvr%numTurbines - RefPoint(1:2) = RefPoint(1:2) + dvr%WT(iWT)%originInit(1:2) - end do - RefPoint(1:2) = RefPoint(1:2) / dvr%numTurbines + if (dvr%numTurbines>0) then + RefPoint(1:2) = dvr%WT(1)%originInit(1:2) + do iWT=2,dvr%numTurbines + RefPoint(1:2) = RefPoint(1:2) + dvr%WT(iWT)%originInit(1:2) + end do + RefPoint(1:2) = RefPoint(1:2) / dvr%numTurbines + else + RefPoint(1:3) = 0.0_ReKi + endif RefPoint(3) = 0.0_ReKi RefLengths = GroundRad + sqrt((WorldBoxMax(1)-WorldBoxMin(1))**2 + (WorldBoxMax(2)-WorldBoxMin(2))**2) diff --git a/modules/aerodyn/src/AeroDyn_Inflow.f90 b/modules/aerodyn/src/AeroDyn_Inflow.f90 index 90a5bbd9c1..62bead89f4 100644 --- a/modules/aerodyn/src/AeroDyn_Inflow.f90 +++ b/modules/aerodyn/src/AeroDyn_Inflow.f90 @@ -86,12 +86,18 @@ subroutine ADI_Init(InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut InitOut%Ver = InitOut_AD%ver ! Add writeoutput units and headers to driver, same for all cases and rotors! !TODO: this header is too short if we add more rotors. Should also add a rotor identifier - call concatOutputHeaders(InitOut%WriteOutputHdr, InitOut%WriteOutputUnt, InitOut_AD%rotors(1)%WriteOutputHdr, InitOut_AD%rotors(1)%WriteOutputUnt, errStat2, errMsg2); if(Failed()) return + if (size(InitOut_AD%rotors)>=1) then + call concatOutputHeaders(InitOut%WriteOutputHdr, InitOut%WriteOutputUnt, InitOut_AD%rotors(1)%WriteOutputHdr, InitOut_AD%rotors(1)%WriteOutputUnt, errStat2, errMsg2); if(Failed()) return + endif ! --- Initialize grouped outputs - !TODO: assumes one rotor - p%NumOuts = p%AD%rotors(1)%NumOuts + p%AD%rotors(1)%BldNd_TotNumOuts + m%IW%p%NumOuts - call AllocAry(y%WriteOutput, p%NumOuts, 'WriteOutput', errStat2, errMsg2); if (Failed()) return + if (size(InitOut_AD%rotors)>=1) then + !TODO: assumes one rotor + p%NumOuts = p%AD%rotors(1)%NumOuts + p%AD%rotors(1)%BldNd_TotNumOuts + m%IW%p%NumOuts + call AllocAry(y%WriteOutput, p%NumOuts, 'WriteOutput', errStat2, errMsg2); if (Failed()) return + else + p%NumOuts = m%IW%p%NumOuts + endif ! --- Initialize outputs call AllocAry(y%IW_WriteOutput, size(m%IW%y%WriteOutput),'IW_WriteOutput', errStat2, errMsg2); if(Failed()) return @@ -298,21 +304,26 @@ subroutine ADI_CalcOutput(t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg) y%PLExp = m%IW%PLExp ! --- Set outputs - !TODO: this assumes one rotor!!! - associate(AD_NumOuts => p%AD%rotors(1)%NumOuts + p%AD%rotors(1)%BldNd_TotNumOuts, & - IW_NumOuts => m%IW%p%NumOuts) - y%WriteOutput(1:IW_NumOuts) = y%IW_WriteOutput(1:IW_NumOuts) - y%WriteOutput(IW_NumOuts+1:p%NumOuts) = y%AD%rotors(1)%WriteOutput(1:AD_NumOuts) - end associate - - !---------------------------------------------------------------------------- - ! Store hub height velocity calculated in CalcOutput - !---------------------------------------------------------------------------- + if (size(p%AD%rotors)>=1) then + !TODO: this assumes one rotor!!! + associate(AD_NumOuts => p%AD%rotors(1)%NumOuts + p%AD%rotors(1)%BldNd_TotNumOuts, & + IW_NumOuts => m%IW%p%NumOuts) + y%WriteOutput(1:IW_NumOuts) = y%IW_WriteOutput(1:IW_NumOuts) + y%WriteOutput(IW_NumOuts+1:p%NumOuts) = y%AD%rotors(1)%WriteOutput(1:AD_NumOuts) + end associate + + !---------------------------------------------------------------------------- + ! Store hub height velocity calculated in CalcOutput + !---------------------------------------------------------------------------- + + if (p%storeHHVel) then + do iWT = 1, size(u%AD%rotors) + y%HHVel(:,iWT) = m%AD%Inflow(1)%RotInflow(iWT)%InflowOnHub(:,1) + end do + endif - if (p%storeHHVel) then - do iWT = 1, size(u%AD%rotors) - y%HHVel(:,iWT) = m%AD%Inflow(1)%RotInflow(iWT)%InflowOnHub(:,1) - end do + else + y%WriteOutput(1:p%NumOuts) = y%IW_WriteOutput(1:m%IW%p%NumOuts) endif contains diff --git a/modules/aerodyn/src/FVW.f90 b/modules/aerodyn/src/FVW.f90 index 2c55f4dc27..cadbde92f3 100644 --- a/modules/aerodyn/src/FVW.f90 +++ b/modules/aerodyn/src/FVW.f90 @@ -106,6 +106,12 @@ subroutine FVW_Init(AFInfo, InitInp, u, p, x, xd, z, OtherState, y, m, Interval, p%iNWStart=1 endif + ! Source panels (parameters, misc and constraints) + ! NOTE: needs to be before FVW_InitMiscVars otherwise number of panels not accounted for + if (len_trim(InputFileData%SrcPnlFile)>0) then + call srcPnl_init(p%SrcPnl, m%SrcPnl, z%SrcPnl, ErrStat2, ErrMsg2, filename=InputFileData%SrcPnlFile); if(Failed()) return + endif + ! Initialize Misc Vars (may depend on input file) CALL FVW_InitMiscVars( p, m, ErrStat2, ErrMsg2 ); if(Failed()) return @@ -148,8 +154,19 @@ subroutine FVW_Init(AFInfo, InitInp, u, p, x, xd, z, OtherState, y, m, Interval, interval = InitInp%DTAero ! important, gluecode and UA, needs proper interval call UA_Init_Wrapper(AFInfo, InitInp, interval, p, x, xd, OtherState, m, ErrStat2, ErrMsg2); if (Failed()) return + ! --- Other States + OtherState%Initialized = .true. + OtherState%ShedScale = 1.0_ReKi ! Will be overriden + ! TODO: should be otherstate + !m%FirstCall = .True. + !m%nNW = p%iNWStart-1 ! Number of active nearwake panels + !m%nFW = 0 ! Number of active farwake panels + !m%iStep = 0 ! Current step number + !m%VTKstep = -1 ! Counter of VTK outputs + !m%VTKlastTime = -HUGE(1.0_DbKi) + !m%OldWakeTime = -HUGE(1.0_DbKi) + ! Framework types unused - OtherState%Dummy = 0 xd%Dummy = 0 InitOut%Dummy = 0 CONTAINS @@ -267,6 +284,8 @@ subroutine FVW_InitMiscVars( p, m, ErrStat, ErrMsg ) endif m%GridOutputs(iGrid)%tLastOutput = -HUGE(1.0_DbKi) enddo + ! Panels + nMax = nMax + p%SrcPnl%n call AllocAry( m%r_wind, 3, nMax, 'Requested wind points', ErrStat2, ErrMsg2 );call SetErrStat ( ErrStat2, ErrMsg2, ErrStat,ErrMsg,RoutineName ) m%r_wind = 0.0_ReKi ! set to zero so InflowWind can shortcut calculations @@ -362,6 +381,7 @@ subroutine FVW_SetParametersFromInputs( InitInp, p, ErrStat, ErrMsg ) ! p%nWings = size(InitInp%WingsMesh) p%DTaero = InitInp%DTaero ! AeroDyn Time step + p%AirDens = InitInp%AirDens ! Air Density p%KinVisc = InitInp%KinVisc ! Kinematic air viscosity p%MHK = InitInp%MHK ! MHK flag p%WtrDpth = InitInp%WtrDpth ! Water depth @@ -421,6 +441,8 @@ SUBROUTINE FVW_SetParametersFromInputFile( InputFileData, p, ErrStat, ErrMsg ) character(ErrMsgLen) :: ErrMsg2 ErrStat = ErrID_None ErrMsg = "" + errStat2 = ErrID_None + errMsg2 = "" ! Set parameters from input file p%IntMethod = InputFileData%IntMethod @@ -451,26 +473,32 @@ SUBROUTINE FVW_SetParametersFromInputFile( InputFileData, p, ErrStat, ErrMsg ) do iW=1,p%nWings if (allocated(p%W(iW)%PrescribedCirculation)) deallocate(p%W(iW)%PrescribedCirculation) if (InputFileData%CircSolvMethod==idCircPrescribed) then - call AllocAry(p%W(iW)%PrescribedCirculation, p%W(iW)%nSpan, 'Prescribed Circulation', ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat,ErrMsg,'FVW_SetParameters'); + call AllocAry(p%W(iW)%PrescribedCirculation, p%W(iW)%nSpan, 'Prescribed Circulation', ErrStat2, ErrMsg2); if(Failed()) return p%W(iW)%PrescribedCirculation = -999999_ReKi; if (.not. allocated(p%W(iW)%s_CP)) then - ErrMsg = 'Spanwise coordinate not allocated.' - ErrStat = ErrID_Fatal - return + ErrMsg2 = 'Spanwise coordinate not allocated.' + ErrStat2 = ErrID_Fatal + if (Failed()) return endif - call ReadAndInterpGamma(trim(InputFileData%CirculationFile), p%W(iW)%s_CP(1:p%W(iW)%nSpan), p%W(iW)%s_LL(p%W(iW)%nSpan+1), p%W(iW)%PrescribedCirculation, ErrStat2, ErrMsg2) - call SetErrStat ( ErrStat2, ErrMsg2, ErrStat,ErrMsg,'FVW_SetParameters' ); + call ReadAndInterpGamma(trim(InputFileData%CirculationFile), p%W(iW)%s_CP(1:p%W(iW)%nSpan), p%W(iW)%s_LL(p%W(iW)%nSpan+1), p%W(iW)%PrescribedCirculation, ErrStat2, ErrMsg2); if(Failed()) return endif enddo - + if(Failed()) return +contains + logical function Failed() + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'FVW_SetParametersFromInputFile') + Failed = ErrStat >= AbortErrLev + end function Failed end subroutine FVW_SetParametersFromInputFile !---------------------------------------------------------------------------------------------------------------------------------- -!> This routine is called at the end of the simulation. -subroutine FVW_FinalWrite(u, p, x, z, m, ErrStat, ErrMsg) +!> This routine is called at the end of the simulation. +!! NOTE: we don't want to call this if OLAF is not fully initialized as some variables might be unallocated +subroutine FVW_FinalWrite(u, p, x, z, OtherState, m, ErrStat, ErrMsg) type(FVW_InputType), intent(in ) :: u !< System inputs type(FVW_ParameterType), intent(in ) :: p !< Parameters type(FVW_ContinuousStateType), intent(in ) :: x !< Continuous states type(FVW_ConstraintStateType), intent(in ) :: z !< Constraint states + type(FVW_OtherStateType), intent(inout) :: OtherState !< Input: Other states at t; Output: at t+DTaero type(FVW_MiscVarType), intent(inout) :: m !< Misc/optimization variables integer(IntKi), intent( out) :: ErrStat !< Error status of the operation character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None @@ -479,7 +507,7 @@ subroutine FVW_FinalWrite(u, p, x, z, m, ErrStat, ErrMsg) ErrStat = ErrID_None ErrMsg = "" ! Place any last minute operations or calculations here: - if (p%WrVTK>0 .and. m%VTKstep0 .and. m%VTKstep This is a tight coupling routine for solving for the residual of the constraint state functions. +!> Solves for the residual of the constraint state functions. +!! Constraints are: +!! - Lifting line circulation +!! - Intensity of the source panels +!! In theory, we should solve for both at the same time, or iterate between the two +!! Here, we do it sequentially, and only once. subroutine FVW_CalcConstrStateResidual( t, u, p, x, xd, z_guess, OtherState, m, z_out, AFInfo, ErrStat, ErrMsg, iLabel) real(DbKi), intent(in ) :: t !< Current simulation time in seconds type(FVW_InputType), intent(in ) :: u !< Inputs at t type(FVW_ParameterType), intent(in ) :: p !< Parameters - type(FVW_ContinuousStateType), intent(in ) :: x !< Continuous states at t + type(FVW_ContinuousStateType), intent(inout) :: x !< Continuous states at t. NOTE: inout because of LL mapping for panels type(FVW_DiscreteStateType), intent(in ) :: xd !< Discrete states at t type(FVW_ConstraintStateType), intent(in ) :: z_guess !< Constraint states at t (possibly a guess) type(FVW_OtherStateType), intent(in ) :: OtherState !< Other states at t @@ -1336,7 +1367,9 @@ subroutine FVW_CalcConstrStateResidual( t, u, p, x, xd, z_guess, OtherState, m, integer(IntKi), intent(in ) :: iLabel integer(IntKi), intent( OUT) :: ErrStat !< Error status of the operation character(*), intent( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None - integer :: iW + integer :: iW + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 ! Initialize ErrStat ErrStat = ErrID_None @@ -1350,16 +1383,45 @@ subroutine FVW_CalcConstrStateResidual( t, u, p, x, xd, z_guess, OtherState, m, ! Solve for the residual of the constraint state functions here: !z%residual = 0.0_ReKi !z%W(iW)%Gamma_LL = 0.0_ReKi + !if (.not.allocated(z_out%W)) allocate(z_out%W(p%nWings)) allocate(z_out%W(p%nWings)) do iW= 1,p%nWings - call AllocAry( z_out%W(iW)%Gamma_LL, p%W(iW)%nSpan, 'Lifting line Circulation', ErrStat, ErrMsg ); + call AllocAry( z_out%W(iW)%Gamma_LL, p%W(iW)%nSpan, 'Lifting line Circulation', ErrStat2, ErrMsg2); z_out%W(iW)%Gamma_LL = -999999_ReKi; enddo - CALL Wings_ComputeCirculation(t, z_out, z_guess, p, x, m, AFInfo, ErrStat, ErrMsg, iLabel) + CALL Wings_ComputeCirculation(t, z_out, z_guess, p, x, m, AFInfo, ErrStat2, ErrMsg2, iLabel); if(Failed()) return + + ! ---- Source panels + if (p%SrcPnl%n>0) then + if (mod(m%iStep, p%nSrcPnlUpdate)==0 .or. m%iStep<3) then + ! NOTE: panels don't move for now, so no need to recompute influence matrix + ! call srcPnl_build_mat(p%SrcPnl, m%SrcPnl%AI) + + ! Map LL circulation to states + call Map_LL_NW(p, m, z_out, x, OtherState%ShedScale, errStat2, errMsg2); if(Failed()) return + + ! Compute Wind (Uwnd) and induced velocity (Uext) from elements different than panels + ! NOTE: Both are summed and stored into Uext + call srcPnl_ExtVelocities_OnPanels(u, p, x, m, errStat2, errMsg2); if(Failed()) return + !m%SrcPnl%Uext = 0.0_ReKi ! HACK + !m%SrcPnl%Uwnd = (/1.0, 0., 0./) + + ! Compute the value of the source panel (sigma) ! assumes that Uext and AI were computed before + if (.not.allocated(z_out%SrcPnl%Sigma)) allocate(z_out%SrcPnl%Sigma(p%SrcPnl%n)) + z_out%SrcPnl%Sigma = 0.0_ReKi + call srcPnl_solve(p%SrcPnl, m%SrcPnl, z_out%SrcPnl, errStat2, errMsg2); if(Failed()) return + else + z_out%SrcPnl%Sigma = m%SrcPnl%RHS ! Use previous intensities + endif + endif - if(.false.) print*,OtherState%Dummy !unused var if(.false.) print*,xd%Dummy +contains + logical function Failed() + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'FVW_CalcConstrStateResidual') + Failed = ErrStat >= AbortErrLev + end function Failed end subroutine FVW_CalcConstrStateResidual @@ -1445,6 +1507,11 @@ subroutine FVW_CalcOutput(t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg) ! Compute induced velocity at AD nodes call CalcOutputForAD(u,p,x,y,m, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + + ! --- Source Panels: Compute velocity, pressure, loads on the source panels + if (p%SrcPnl%n>0) then + call srcPnl_calcOutput(p%SrcPnl, m%SrcPnl, z%SrcPnl, p%AirDens) + endif ! Write some info to screen when major milestone achieved if (m%iStep == p%nNWFree .and. p%nNWFree0) then @@ -1553,6 +1620,7 @@ subroutine WriteVTKOutputs(t, force, VTKstep, u, p, x, z, m, ErrStat, ErrMsg) endif enddo endif + if (OLAF_PROFILING) call toc() end subroutine WriteVTKOutputs !---------------------------------------------------------------------------------------------------------------------------------- diff --git a/modules/aerodyn/src/FVW_BiotSavart.f90 b/modules/aerodyn/src/FVW_BiotSavart.f90 index 99dc7d0aca..36c73d1dd9 100644 --- a/modules/aerodyn/src/FVW_BiotSavart.f90 +++ b/modules/aerodyn/src/FVW_BiotSavart.f90 @@ -2,12 +2,13 @@ !! NOTE: these functions should be independent of the framework types module FVW_BiotSavart - use NWTC_Library, only: ReKi, IntKi + use NWTC_Library, only: ReKi, IntKi, Pi, EqualRealNos use OMP_LIB implicit none real(ReKi),parameter :: PRECISION_UI = epsilon(1.0_ReKi)/100 !< NOTE assuming problem of size 1 + real(ReKi),parameter :: PRECISION_EPS = epsilon(1.0_ReKi) !< Machine Precision For the given ReKi for problems of scale 1! real(ReKi),parameter :: MIN_EXP_VALUE=-10.0_ReKi real(ReKi),parameter :: MINDENOM=0.0_ReKi ! real(ReKi),parameter :: MINDENOM=1e-15_ReKi @@ -24,6 +25,7 @@ module FVW_BiotSavart integer(IntKi), parameter, dimension(3) :: idRegPartVALID = (/idRegNone,idRegExp,idRegCompact/) real(ReKi),parameter :: fourpi_inv = 0.25_ReKi / ACOS(-1.0_Reki ) + real(ReKi),parameter :: fourpi = 4.00_ReKi * ACOS(-1.0_Reki ) contains @@ -429,4 +431,176 @@ subroutine ui_quad_n1(CPs, nCPs, P1, P2, P3, P4, Gamm, RegFunction, RegParam, Ui !OMP END PARALLEL end subroutine ui_quad_n1 + +subroutine ui_quad_src_11(CP, Sigma, xi, eta, RefPoint, R_g2p, UI) + real(ReKi), intent(in) :: Sigma !< Source panel intensity + real(ReKi), dimension(3), intent(in) :: CP !< Control Point + real(ReKi), dimension(3), intent(out) :: UI !< Induced velocity + real(ReKi), dimension(3), intent(in) :: RefPoint !< Coordinate of panel origin in ref coordinates + real(ReKi), dimension(4), intent(in) :: xi !< Panel points coordinates + real(ReKi), dimension(4), intent(in) :: eta !< Panel points coordinates + real(ReKi), dimension(3,3), intent(in) :: R_g2p !< 3 x 3, global 2 panel + real(ReKi),parameter :: eps_quadsource=1e-6_ReKi !!!!!!!!!!!!!!!!!! !< Used if z coordinate close to zero + real(ReKi), dimension(3,3) :: tA !< + real(ReKi) :: d12, d23, d34, d41 !< + real(ReKi) :: m12, m23, m34, m41 !< + real(ReKi) :: xi1, xi2, xi3, xi4 !< + real(ReKi) :: eta1, eta2, eta3, eta4 !< + real(ReKi) :: e1, e2, e3, e4 !< + real(ReKi) :: h1, h2, h3, h4 !< + real(ReKi) :: r1, r2, r3, r4 !< + real(ReKi) :: RJ12, RJ23, RJ34, RJ41 !< + real(ReKi) :: TAN12, TAN23, TAN34, TAN41 !< + real(ReKi), dimension(3) :: Vp !< + real(ReKi), dimension(3) :: DP !< + real(ReKi), dimension(3) :: DPp !< + xi1=xi(1) + xi2=xi(2) + xi3=xi(3) + xi4=xi(4) + eta1=eta(1) + eta2=eta(2) + eta3=eta(3) + eta4=eta(4) + !param that are constant for each panel, distances and slopes - The slopes can be if divided by zero NaN => security required + d12 = sqrt((xi2-xi1)**2+(eta2-eta1)**2) + d23 = sqrt((xi3-xi2)**2+(eta3-eta2)**2) + d34 = sqrt((xi4-xi3)**2+(eta4-eta3)**2) + d41 = sqrt((xi1-xi4)**2+(eta1-eta4)**2) + + ! transform control points in panel coordinate system using matrix + DP(1:3) = CP(1:3)-RefPoint(1:3) + DPp = matmul(R_g2p, DP) ! transfo in element coordinate system, noted x,y,z, but in fact xi eta zeta + ! scalars + r1 = sqrt((DPp(1)-xi1)**2 + (DPp(2)-eta1)**2 + DPp(3)**2) + r2 = sqrt((DPp(1)-xi2)**2 + (DPp(2)-eta2)**2 + DPp(3)**2) + r3 = sqrt((DPp(1)-xi3)**2 + (DPp(2)-eta3)**2 + DPp(3)**2) + r4 = sqrt((DPp(1)-xi4)**2 + (DPp(2)-eta4)**2 + DPp(3)**2) + ! + e1 = DPp(3)**2 + (DPp(1)-xi1)**2 + e2 = DPp(3)**2 + (DPp(1)-xi2)**2 + e3 = DPp(3)**2 + (DPp(1)-xi3)**2 + e4 = DPp(3)**2 + (DPp(1)-xi4)**2 + ! + h1 = (DPp(2)-eta1)*(DPp(1)-xi1) + h2 = (DPp(2)-eta2)*(DPp(1)-xi2) + h3 = (DPp(2)-eta3)*(DPp(1)-xi3) + h4 = (DPp(2)-eta4)*(DPp(1)-xi4) + ! Velocities in element frame + ! --- Log term + ! Security - Katz Plotkin page 608 appendix D Code 11 + if ( r1+r2-d12<=0.0_ReKi .or. d12 <=0.0_ReKi ) then + RJ12=0._ReKi + else + RJ12=1/d12 * log((r1+r2-d12)/(r1+r2+d12)) + endif + + if ( r2+r3-d23<=0.0_ReKi .or. d23 <=0.0_ReKi ) then + RJ23=0._ReKi + else + RJ23=1/d23 * log((r2+r3-d23)/(r2+r3+d23)) + endif + + if ( r3+r4-d34<=0.0_ReKi .or. d34 <=0.0_ReKi ) then + RJ34=0._ReKi + else + RJ34=1/d34 * log((r3+r4-d34)/(r3+r4+d34)) + endif + + if ( r4+r1-d41<=0.0_ReKi .or. d41 <=0.0_ReKi ) then + RJ41=0._ReKi + else + RJ41=1/d41 * log((r4+r1-d41)/(r4+r1+d41)) + endif + ! --- Tan term + ! 12 + if (EqualRealNos(xi2,xi1)) then ! Security - Hess 1962 - page 47 - bottom + TAN12=0._ReKi + else + m12=(eta2-eta1)/(xi2-xi1) + if( abs(DPp(3)) Induced velocity by several flat quadrilateral source panels on multiple control points (CPs) +subroutine ui_quad_src_nn(CPs, Sigmas, xi, eta, RefPoint, R_g2p, UI, nCPs, nPanels) + integer, intent(in) :: nCPs + integer, intent(in) :: nPanels + real(ReKi), dimension(3,nCPs), intent(in) :: CPs !< 3 x nCPs, coordinates of the control points + real(ReKi), dimension(3,nCPs), intent(inout) :: UI !< 3 x nCPs, induced velocity on control points (side effects) + real(ReKi), dimension(nPanels), intent(in) :: Sigmas !< np, intensities of Panels + real(ReKi), dimension(3,nPanels), intent(in) :: RefPoint !< 3 x np , coordinate of panel origin in ref coordinates + real(ReKi), dimension(4,nPanels), intent(in) :: xi !< 4 x np + real(ReKi), dimension(4,nPanels), intent(in) :: eta !< 4 x np + real(ReKi), dimension(3,3,nPanels), intent(in) :: R_g2p !< 3 x 3 x np, transformation matrix global to panel + real(ReKi) :: Uind_tmp(3) !< + real(ReKi) :: Uind_cum(3) !< + integer :: ip, icp !< loop index + !$OMP PARALLEL DEFAULT(SHARED) + !$OMP DO PRIVATE(icp, Uind_cum, Uind_tmp, ip) schedule(runtime) + do icp=1,nCPs ! loop on Control Points + Uind_cum = 0.0_ReKi + do ip=1,nPanels !loop on panels + call ui_quad_src_11(CPs(:,icp), Sigmas(ip), xi(:,ip), eta(:,ip), RefPoint(:,ip), R_g2p(:,:,ip), Uind_tmp) + Uind_cum = Uind_cum + Uind_tmp + enddo + UI(1:3,icp) = UI(1:3,icp) + Uind_cum + end do ! control points + !$OMP END DO + !$OMP END PARALLEL +end subroutine ui_quad_src_nn + +elemental real(ReKi) function signit(ref, val) + real(ReKi),intent(in) ::ref + real(ReKi),intent(in) ::val + if ( abs(val)>PRECISION_EPS ) then + signit = sign(ref, val) + else + signit = 1.0_ReKi + endif +endfunction + end module FVW_BiotSavart diff --git a/modules/aerodyn/src/FVW_IO.f90 b/modules/aerodyn/src/FVW_IO.f90 index 48580e578e..4cde5064d4 100644 --- a/modules/aerodyn/src/FVW_IO.f90 +++ b/modules/aerodyn/src/FVW_IO.f90 @@ -24,6 +24,7 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) character(ErrMsgLen) :: ErrMsg2 ErrStat = ErrID_None ErrMsg = "" + Inp%SrcPnlFile = '' ! TODO registry init for empty strings ! Open file CALL GetNewUnit( UnIn ) CALL OpenFInpfile(UnIn, TRIM(FileName), ErrStat2, ErrMsg2) @@ -117,37 +118,50 @@ SUBROUTINE FVW_ReadInputFile( FileName, p, m, Inp, ErrStat, ErrMsg ) if(ErrStat2==ErrID_None) then call WrScr(' - Reading advanced options for OLAF:') do while(ErrStat2==ErrID_None) - read(UnIn, '(A)', iostat=ErrStat2) sDummy + read(UnIn, '(A)', iostat=ErrStat2) sLine if (ErrStat2/=ErrID_None) exit + sDummy = sLine call Conv2UC(sDummy) ! to uppercase if (index(sDummy, '!') == 1 .or. index(sDummy, '=') == 1 .or. index(sDummy, '#') == 1) then ! pass comment lines elseif (index(sDummy, 'INDUCTIONATCP')>1) then - read(sDummy, '(L1)') p%InductionAtCP + read(sLine, '(L1)') p%InductionAtCP print*,' >>> InductionAtCP ',p%InductionAtCP elseif (index(sDummy, 'WAKEATTE')>1) then - read(sDummy, '(L1)') p%WakeAtTE + read(sLine, '(L1)') p%WakeAtTE print*,' >>> WakeAtTE ',p%WakeAtTE elseif (index(sDummy, 'DSTALLONWAKE')>1) then - read(sDummy, '(L1)') p%DStallOnWake + read(sLine, '(L1)') p%DStallOnWake print*,' >>> DStallOnWake ',p%DStallOnWake elseif (index(sDummy, 'INDUCTION')>1) then - read(sDummy, '(L1)') p%Induction + read(sLine, '(L1)') p%Induction print*,' >>> Induction ',p%Induction elseif (index(sDummy, 'KFROZENNWEND')>1) then - read(sDummy, *) p%kFrozenNWEnd + read(sLine, *) p%kFrozenNWEnd print*,' >>> kFrozenNWEnd ',p%kFrozenNWEnd elseif (index(sDummy, 'KFROZENNWSTART')>1) then - read(sDummy, *) p%kFrozenNWStart + read(sLine, *) p%kFrozenNWStart print*,' >>> kFrozenNWStart ',p%kFrozenNWStart elseif (index(sDummy, 'VELOCITYMETHODLL')>1) then - read(sDummy, *) Inp%VelocityMethod(2) + read(sLine, *) Inp%VelocityMethod(2) print*,' >>> VelocityMethod ',Inp%VelocityMethod elseif (index(sDummy, 'TREEBRANCHFACTORLL')>1) then - read(sDummy, *) Inp%TreeBranchFactor(2) + read(sLine, *) Inp%TreeBranchFactor(2) print*,' >>> TreeBranchFactor ',Inp%TreeBranchFactor + elseif (index(sDummy, 'SRCPNLFILE')>1) then + read(sLine, *) Inp%SrcPnlFile + print*,' >>> SrcPnlFile ',trim(Inp%SrcPnlFile) + elseif (index(sDummy, 'ZGROUNDPUSH')>1) then + read(sLine, *) p%zGroundPush + print*,' >>> zGroundPush ',p%zGroundPush + elseif (index(sDummy, 'ZGROUND')>1) then + read(sLine, *) p%zGround + print*,' >>> zGround ',p%zGround + elseif (index(sDummy, 'NSRCPNLUPDATE')>1) then + read(sLine, *) p%nSrcPnlUpdate + print*,' >>> nSrcPnlUpdate ',p%nSrcPnlUpdate else - print*,'[WARN] Line ignored: '//trim(sDummy) + print*,'[WARN] Line ignored: '//trim(sLine) endif enddo endif @@ -453,8 +467,15 @@ subroutine WrVTK_FVW(p, x, z, m, FileRootName, VTKcount, Twidth, bladeFrame, Hub nSeg = 2*nSeg nSegP = 2*nSegP endif - 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) + 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) + endif + + if (p%SrcPnl%n>0) then + Filename = TRIM(FileRootName)//'.SrcPnl.'//Tstr//'.vtk' + CALL WrVTK_Panels(Filename, mvtk, p%SrcPnl, m%SrcPnl, z%SrcPnl) + endif if(.false.) print*,z%W(1)%Gamma_LL(1) ! unused var for now end subroutine WrVTK_FVW @@ -511,6 +532,42 @@ subroutine WrVTK_FVW_Grid(p, m, iGrid, FileRootName, VTKcount, Twidth, HubOrient end subroutine WrVTK_FVW_Grid +subroutine WrVTK_Panels(filename, mvtk, p, m, z) + use VTK + character(len=*), intent(in) :: filename + type(VTK_Misc), intent(inout) :: mvtk !< miscvars for VTK output + type(T_SrcPanlParam), intent(in) :: p + type(T_SrcPanlMisc) , intent(in) :: m + type(T_SrcPanlVar) , intent(in) :: z + if ( vtk_new_ascii_file(trim(filename), 'SourcePanels '//trim(p%Comment), mvtk) ) then + call vtk_dataset_polydata(p%P, mvtk, bladeFrame=.false.) + call vtk_quad(p%IDs-1, mvtk) + call vtk_cell_data_init(mvtk) + if (.not. allocated(z%Sigma)) then ! First time step + call vtk_cell_data_scalar(m%RHS ,'Intensities', mvtk) + else + call vtk_cell_data_scalar(z%Sigma ,'Intensities', mvtk) + endif + call vtk_cell_data_scalar(p%BodyIDs,'BodyID' , mvtk) + call vtk_cell_data_scalar(m%Cp ,'Cp' , mvtk) + call vtk_cell_data_scalar(m%p ,'p' , mvtk) + call vtk_cell_data_scalar(p%Area ,'Area' , mvtk) + call vtk_cell_data_vector(m%F ,'F' , mvtk) + call vtk_cell_data_vector(m%FpA ,'F/A' , mvtk) + call vtk_cell_data_vector(p%Pmid ,'Pmid' , mvtk) + call vtk_cell_data_vector(p%Pcent ,'Pcent' , mvtk) + call vtk_cell_data_vector(m%Uwnd ,'Uwnd' , mvtk) + call vtk_cell_data_vector(m%Uext ,'Uext' , mvtk) + call vtk_cell_data_vector(m%Uind ,'Uind' , mvtk) + call vtk_cell_data_vector(m%Utot ,'Utot' , mvtk) + call vtk_cell_data_vector(p%Normal ,'Normal' , mvtk) + !call vtk_cell_data_vector(p%Ts ,'Ts',mvtk) + !call vtk_cell_data_vector(p%Tc ,'Tc',mvtk) + !call vtk_cell_data_vector(p%To ,'To',mvtk) + call vtk_close_file(mvtk) + endif +endsubroutine WrVTK_Panels + subroutine WrVTK_Segments(filename, mvtk, SegPoints, SegConnct, SegGamma, SegEpsilon, bladeFrame) use VTK diff --git a/modules/aerodyn/src/FVW_Registry.txt b/modules/aerodyn/src/FVW_Registry.txt index bdfef5796e..45a7aff474 100644 --- a/modules/aerodyn/src/FVW_Registry.txt +++ b/modules/aerodyn/src/FVW_Registry.txt @@ -22,9 +22,9 @@ typedef ^ ^ ReKi typedef ^ ^ ReKi xEnd - - - "xEnd" - typedef ^ ^ ReKi yEnd - - - "yEnd" - typedef ^ ^ ReKi zEnd - - - "zEnd" - -typedef ^ ^ IntKi nx - - - "nx" - -typedef ^ ^ IntKi ny - - - "ny" - -typedef ^ ^ IntKi nz - - - "nz" - +typedef ^ ^ IntKi nx - 0 - "nx" - +typedef ^ ^ IntKi ny - 0 - "ny" - +typedef ^ ^ IntKi nz - 0 - "nz" - 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" - @@ -35,15 +35,43 @@ typedef ^ ^ IntKi typedef ^ ^ ReKi Gamma : - - "Segment circulations" - typedef ^ ^ ReKi Epsilon : - - "Segment regularization parameter" - typedef ^ ^ IntKi RegFunction - - - "Type of regularizaion function (LambOseen, Vatistas, see FVW_BiotSavart)" - -typedef ^ ^ IntKi nAct - - - "Number of active segments" - -typedef ^ ^ IntKi nActP - - - "Number of active segment points" - +typedef ^ ^ IntKi nAct - 0 - "Number of active segments" - +typedef ^ ^ IntKi nActP - 0 - "Number of active segment points" - # TODO add tree ##################### Particles ############### typedef FVW/FVW T_Part ReKi P :: - - "Particle Points" - typedef ^ ^ ReKi Alpha :: - - "Particle intensity 3 x nP" - typedef ^ ^ ReKi RegParam : - - "Particle regularization parameter" - typedef ^ ^ IntKi RegFunction - - - "Type of regularizaion function (FVW_BiotSavart)" - -typedef ^ ^ IntKi nAct - - - "Number of active particles <=nP" - +typedef ^ ^ IntKi nAct - 0 - "Number of active particles <=nP" - +##################### Panels ############### +typedef FVW/FVW T_SrcPanlParam IntKi n - 0 - "Number of panels" - +typedef ^ ^ CHARACTER(1024) Comment - "" - "Area" - +typedef ^ ^ ReKi Area : - - "Area" - +typedef ^ ^ ReKi P :: - - "Points database" - +typedef ^ ^ IntKi IDs :: - - "Points connectivity ID for each panel 4xnPn" - +typedef ^ ^ IntKi BodyIDs : - - "Body ID for each panel" - +typedef ^ ^ ReKi Pcent :: - - "Centroid point" - +typedef ^ ^ ReKi Pmid :: - - "Mid point" - +typedef ^ ^ ReKi Normal :: - - "Normal to panel" - +typedef ^ ^ ReKi xi :: - - "Xi coordinate for each panel point" - +typedef ^ ^ ReKi eta :: - - "eta coordinate for each panel point" - +typedef ^ ^ ReKi R_g2p ::: - - "Transformation matrix global to Panel" - + +typedef FVW/FVW T_SrcPanlVar ReKi Sigma : - - "Intensities" - + +typedef FVW/FVW T_SrcPanlMisc ReKi AI :: - - "Influence matrix to solve for intensity" - +typedef ^ ^ ReKi UUI ::: - - "Unit induced velocity" - +typedef ^ ^ ReKi Uext :: - - "External inflow velocity at panel control points" - +typedef ^ ^ ReKi Uwnd :: - - "Inflow velocity vector at panel control points" - +typedef ^ ^ ReKi Uind :: - - "Induced velocity vector at panel control points" - +typedef ^ ^ ReKi Utot :: - - "Full velocity (include induced) vector at panel control points" - +typedef ^ ^ ReKi Cp : - - "Pressure coefficient" - +typedef ^ ^ ReKi p : - - "Pressure (1/2 rho Utot^2)" - +typedef ^ ^ ReKi F :: - - "Force from pressure" - +typedef ^ ^ ReKi FpA :: - - "Force per area" - +typedef ^ ^ ReKi RHS : - - "Right hand side" - +typedef ^ ^ IntKi IPIV : - - "Pivot for factorization" - ##################### Registry for FVW ############### # ..... PARAMETERS ............. @@ -51,24 +79,24 @@ typedef FVW/FVW Wng_ParameterType ReKi typedef ^ ^ ReKi chord_CP : - - "Chord on LL cp " m typedef ^ ^ ReKi s_LL : - - "Spanwise coordinate of LL elements" m typedef ^ ^ ReKi s_CP : - - "Spanwise coordinate of LL CP" m -typedef ^ ^ IntKi iRotor - - - "Index of rotor the wing belong to" - +typedef ^ ^ IntKi iRotor - 0 - "Index of rotor the wing belong to" - typedef ^ ^ IntKi AFindx :: - - "Index to the airfoils from AD15 [BladeNode,BladeIndex=1]" - -typedef ^ ^ IntKi nSpan - - - "TODO, should be defined per wing. Number of spanwise element" - +typedef ^ ^ IntKi nSpan - 0 - "TODO, should be defined per wing. Number of spanwise element" - typedef ^ ^ ReKi PrescribedCirculation : - - "Prescribed circulation on all lifting lines" "m/s" #FVW_ParameterType -typedef FVW/FVW ParameterType IntKi nRotors - - - "Number of Wings" - -typedef ^ ^ IntKi nWings - - - "Number of Wings" - +typedef FVW/FVW ParameterType IntKi nRotors - 0 - "Number of Rotors" - +typedef ^ ^ IntKi nWings - 0 - "Number of Wings" - #typedef ^ ^ IntKi Rot2Wings : - - "Index mapping from wings to rotors" - typedef ^ ^ Wng_ParameterType W : - - "Wings parameters" - typedef ^ ^ IntKi Bld2Wings :: - - "Index mapping from blades to wings" - -typedef ^ ^ IntKi iNWStart - - - "Index where NW start in r_NW. (iNWStart=2, the first panel contains the lifting line panel, otherwise, start at 1)" - -typedef ^ ^ IntKi nNWMax - - - "Maximum number of nw panels, per wing" - -typedef ^ ^ IntKi nNWFree - - - "Number of nw panels that are free, per wing" - -typedef ^ ^ IntKi nFWMax - - - "Maximum number of fw panels, per wing" - -typedef ^ ^ IntKi nFWFree - - - "Number of fw panels that are free, per wing" - +typedef ^ ^ IntKi iNWStart - 0 - "Index where NW start in r_NW. (iNWStart=2, the first panel contains the lifting line panel, otherwise, start at 1)" - +typedef ^ ^ IntKi nNWMax - 0 - "Maximum number of nw panels, per wing" - +typedef ^ ^ IntKi nNWFree - 0 - "Number of nw panels that are free, per wing" - +typedef ^ ^ IntKi nFWMax - 0 - "Maximum number of fw panels, per wing" - +typedef ^ ^ IntKi nFWFree - 0 - "Number of fw panels that are free, per wing" - typedef ^ ^ Logical FWShedVorticity - - - "Include shed vorticity in the far wake" - typedef ^ ^ IntKi IntMethod - - - "Integration Method (1=RK4, 2=AB4, 3=ABM4, 5=Euler1)" - typedef ^ ^ ReKi FreeWakeStart - - - "Time when wake starts convecting (rolling up)" s @@ -92,6 +120,7 @@ typedef ^ ^ ReKi typedef ^ ^ IntKi PartPerSegment 2 - - "Number of particles per segment, e.g. for tree method, for full wake and lifting line" typedef ^ ^ DbKi DTaero - - - "Time interval for calls calculations" s typedef ^ ^ DbKi DTfvw - - - "Time interval for calculating wake induced velocities" s +typedef ^ ^ ReKi AirDens - - - "Air density" kg/m^3 typedef ^ ^ ReKi KinVisc - - - "Kinematic air viscosity" m^2/s typedef ^ ^ IntKi MHK - - - "MHK flag" - typedef ^ ^ ReKi WtrDpth - - - "Water depth" m @@ -103,7 +132,7 @@ typedef ^ ^ IntKi typedef ^ ^ CHARACTER(1024) RootName - - - "RootName for writing output files" - typedef ^ ^ CHARACTER(1024) VTK_OutFileRoot - - - "Rootdirectory for writing VTK files" - typedef ^ ^ CHARACTER(1024) VTK_OutFileBase - - - "Basename for writing VTK files" - -typedef ^ ^ IntKi nGridOut - - - "Number of VTK grid to output" - +typedef ^ ^ IntKi nGridOut - 0 - "Number of VTK grid to output" - # Parameters advanced options typedef ^ ^ Logical InductionAtCP - .true. - "Compute induced velocities at nodes or CP" typedef ^ ^ Logical WakeAtTE - .true. - "Start the wake at the trailing edge, or at the LL" @@ -111,6 +140,11 @@ typedef ^ ^ Logical typedef ^ ^ Logical Induction - .true. - "Compute induction" typedef ^ ^ ReKi kFrozenNWStart - 0.75 - "Fraction of wake induced velocity at start of frozen wake. 1 seems too strong." typedef ^ ^ ReKi kFrozenNWEnd - 0.5 - "Fraction of wake induced velocity at end of frozen wake" +typedef ^ ^ ReKi zGround - 0.0 - "Ground height" +typedef ^ ^ ReKi zGroundPush - 0.1 - "Distance above ground where vortices are pushed back" +typedef ^ ^ IntKi nSrcPnlUpdate - 1 - "How often do src panel updates (in time steps of OLAF)" +# Parameters panels +typedef ^ ^ T_SrcPanlParam SrcPnl - - - "Source panel parameters" - #.......... ContinuousStateType ...... typedef FVW/FVW Wng_ContinuousStateType ReKi Gamma_NW :: - - "Circulation of the near wake panels ( nSpan x nNW )" - @@ -204,6 +238,7 @@ typedef ^ ^ LOGICAL # Element storage (buffers) typedef ^ ^ T_Sgmt Sgmt - - - "Segments storage" - typedef ^ ^ T_Part Part - - - "Particle storage" - +typedef ^ ^ T_SrcPanlMisc SrcPnl - - - "Source panels storage" - # Wake rollup storage (buffer) typedef ^ ^ ReKi CPs :: - - "Control points used for wake rollup computation" - typedef ^ ^ ReKi Uind :: - - "Induced velocities obtained at control points" - @@ -237,11 +272,13 @@ typedef FVW/FVW Wng_ConstraintStateType Reki # FVW_ConstraintStateType typedef FVW/FVW ConstraintStateType Wng_ConstraintStateType W : - - "rotors constr. states" - typedef ^ ^ Reki residual - - - "Residual" - +typedef ^ ^ T_SrcPanlVar SrcPnl - - - "Source panel constraints" - # ....... OtherStateType ............ # FVW_OtherStateType -typedef FVW/FVW OtherStateType IntKi Dummy - - - "Empty to satisfy framework" +typedef FVW/FVW OtherStateType ReKi ShedScale - - - "Scale shed vorticity at begining of simulation" typedef ^ ^ UA_OtherStateType UA : - - "other states for UnsteadyAero for each wing" - +typedef ^ ^ LOGICAL Initialized - .false. - "True if OLAF is initialized" - #.......... InitInputType ...... @@ -260,6 +297,7 @@ typedef ^ ^ Wng_InitInputType typedef ^ ^ MeshType WingsMesh : - - "Input Mesh defining position and orientation of wings (nSpan+1) " - typedef ^ ^ IntKi numBladeNodes - - - "Number of nodes on each blade" - typedef ^ ^ DbKi DTaero - - - "Time interval for calls (from AD15)" s +typedef ^ ^ ReKi AirDens - - - "Air density" kg/m^3 typedef ^ ^ ReKi KinVisc - - - "Kinematic air viscosity" m^2/s typedef ^ ^ IntKi MHK - - - "MHK flag" - typedef ^ ^ ReKi WtrDpth - - - "Water depth" m @@ -302,6 +340,7 @@ typedef ^ ^ IntKi typedef ^ ^ IntKi VTKBlades - - - "Outputs VTk for each blade 0=no blade, 1=Bld 1" - typedef ^ ^ DbKi DTvtk - - - "Requested timestep between VTK outputs (calculated from the VTK_fps read in)" s typedef ^ ^ IntKi VTKCoord - - - "Switch for VTK outputs coordinate system" - +typedef ^ ^ CHARACTER(1024) SrcPnlFile - - '' "Input file for source panels" - #.......... InitOutputType ...... # FVW_InitOutputType diff --git a/modules/aerodyn/src/FVW_Subs.f90 b/modules/aerodyn/src/FVW_Subs.f90 index c874f07c8c..4bafacfcbd 100644 --- a/modules/aerodyn/src/FVW_Subs.f90 +++ b/modules/aerodyn/src/FVW_Subs.f90 @@ -6,6 +6,11 @@ module FVW_SUBS use FVW_BiotSavart implicit none + type :: T_Panl + type(T_SrcPanlParam), pointer :: p_Src =>null() + type(T_SrcPanlMisc) , pointer :: m_Src =>null() + type(T_SrcPanlVar) , pointer :: z_Src =>null() + endtype ! --- Module parameters ! Circulation solving methods @@ -105,7 +110,7 @@ subroutine Output_Gamma(CP, Gamma_LL, iW, iStep, iLabel, iIter) norm=sqrt(CP(1,i)**2+CP(2,i)**2+CP(3,i)**2) write(iUnit,'(E14.7,A,E14.7,A,E14.7,A,E14.7,A,E14.7)') norm,',', CP(1,i),',',CP(2,i),',',CP(3,i),',', Gamma_LL(i) enddo - close(iUnit) + if (iUnit>0) close(iUnit) endsubroutine Output_Gamma ! ===================================================================================== !> Read a delimited file containing a circulation and interpolate it on the requested Control Points @@ -500,7 +505,7 @@ subroutine find_nan_1D(array, varname) endif enddo if (found) then - print*,'>>>>>>>>>>>>> NAN ',trim(varname),tot,n + print*,'OLAF NAN ',trim(varname),tot,n STOP endif end subroutine @@ -522,7 +527,7 @@ subroutine find_nan_2D(array, varname) endif enddo if (found) then - print*,'>>>>>>>>>>>>> NAN ',trim(varname),tot,n + print*,'OLAF NAN ',trim(varname),tot,n STOP endif end subroutine @@ -547,7 +552,7 @@ subroutine find_nan_3D(array, varname) enddo enddo if (found) then - print*,'>>>>>>>>>>>>> NAN ',trim(varname), tot,n*m + print*,'OLAF: NAN ',trim(varname), tot,n*m STOP endif end subroutine @@ -561,7 +566,7 @@ subroutine find_nan_3D(array, varname) !! ensure that we do not violate requirements in the framework later for changing the size !! of input and output arrays. subroutine SetRequestedWindPoints(r_wind, x, p, m) - real(ReKi), dimension(:,:), allocatable, intent(inout) :: r_wind !< Position where wind is requested + real(ReKi), dimension(:,:), intent(inout) :: r_wind !< Position where wind is requested type(FVW_ContinuousStateType), intent(inout) :: x !< States type(FVW_ParameterType), intent(in ) :: p !< Parameters type(FVW_MiscVarType), intent(in ), target :: m !< Initial misc/optimization variables @@ -574,17 +579,7 @@ subroutine SetRequestedWindPoints(r_wind, x, p, m) ! NOTE: Maximum number of points are passed, whether they "exist" or not. ! NOTE: InflowWind ignores points at (0,0,0) !if (DEV_VERSION) then - ! ! Removing points that don't exist - ! !call print_x_NW_FW(p,m,x,'wind befr') - ! if (m%nNW<=p%nNWMax) then - ! x%W(iW)%r_NW(1:3, 1:p%W(iW)%nSpan+1, m%nNW+2:p%nNWMax+1, 1:p%nWings) = 0.0_ReKi - ! endif - ! if ( ((p%nNWMax<=1) .and. (m%nFW==0)) .or. ((m%nFW>0) .and. (m%nFW<=p%nFWMax))) then - ! x%W(iW)%r_FW(1:3, 1:FWnSpan+1, m%nFW+2:p%nFWMax+1, 1:p%nWings) = 0.0_ReKi - ! else - ! x%W(iW)%r_FW(1:3, 1:FWnSpan+1, m%nFW+1:p%nFWMax+1, 1:p%nWings) = 0.0_ReKi - ! endif - ! !call print_x_NW_FW(p,m,x,'wind after') + ! r_wind = -99999_ReKi !endif iP_end=0 @@ -629,22 +624,18 @@ subroutine SetRequestedWindPoints(r_wind, x, p, m) enddo ! Loop on z enddo ! Loop on grids + ! --- Panels + if (p%SrcPnl%n>0) then + r_wind(1:3,iP_start:iP_start-1+p%SrcPnl%n) = p%SrcPnl%Pcent(1:3,1:p%SrcPnl%n) + endif + !if (DEV_VERSION) then ! ! Additional checks - ! if (any(r_wind(3,:)<=-99999_ReKi)) then - ! call print_x_NW_FW(p,m,x,'wind after') - ! print*,'Error in wind' - ! STOP - ! endif - ! ! Removing points that don't exist - ! if (m%nNW<=p%nNWMax) then - ! x%W(iW)%r_NW(1:3, 1:p%W(iW)%nSpan+1, m%nNW+2:p%nNWMax+1, 1:p%nWings) = -999999.0_ReKi - ! endif - ! if ( ((p%nNWMax<=1) .and. (m%nFW==0)) .or. ((m%nFW>0) .and. (m%nFW<=p%nFWMax))) then - ! x%W(iW)%r_FW(1:3, 1:FWnSpan+1, m%nFW+2:p%nFWMax+1, 1:p%nWings) =-999999.0_ReKi - ! else - ! x%W(iW)%r_FW(1:3, 1:FWnSpan+1, m%nFW+1:p%nFWMax+1, 1:p%nWings) =-999999.0_ReKi - ! endif + ! if (any(r_wind(3,:)<=-99999_ReKi)) then + ! call print_x_NW_FW(p,m,x,'wind after') + ! print*,'Error in wind' + ! STOP + ! endif !endif end subroutine SetRequestedWindPoints @@ -749,6 +740,35 @@ subroutine DistributeRequestedWind_Grid(V_wind, p, m) enddo ! Loop on grids end subroutine DistributeRequestedWind_Grid +subroutine DistributeRequestedWind_Panels(V_wind, p, m) + real(ReKi), dimension(:,:), intent(in ) :: V_wind !< Requested wind, packed + type(FVW_ParameterType), intent(in ) :: p !< Parameters + type(FVW_MiscVarType), target, intent(inout) :: 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 + type(GridOutType), pointer :: g + iP_end=0 + ! --- LL/NW/FW + do iW=1,p%nWings; iP_end = iP_end + p%W(iW)%nSpan; enddo + do iW=1,p%nWings; iP_end = iP_end+(p%W(iW)%nSpan+1)*(p%nNWMax+1) ; enddo + if (p%nFWMax>0) then + do iW=1,p%nWings; iP_end = iP_end+(FWnSpan+1)*(p%nFWMax+1); enddo + endif + ! --- VTK points + iP_end=iP_end+1 + do iGrid=1,p%nGridOut + g => m%GridOutputs(iGrid) + do k=1,g%nz + do j=1,g%ny + do i=1,g%nx + iP_end=iP_end+1 + enddo + enddo + enddo ! Loop on x + enddo ! Loop on grids + iP_start=iP_end + m%SrcPnl%Uwnd(:,1:p%SrcPnl%n) = V_wind(:,iP_start:iP_start+p%SrcPnl%n-1) +end subroutine DistributeRequestedWind_Panels !> Init States @@ -808,31 +828,35 @@ subroutine FVW_InitMiscVarsPostParam( p, m, ErrStat, ErrMsg ) nSeg = nSeg*2 nSegP = nSegP*2 endif - call AllocAry( m%Sgmt%Connct, 4, nSeg , 'SegConnct' , ErrStat2, ErrMsg2 );call SetErrStat(ErrStat2, ErrMsg2, ErrStat,ErrMsg,RoutineName); m%Sgmt%Connct = -999; - call AllocAry( m%Sgmt%Points, 3, nSegP, 'SegPoints' , ErrStat2, ErrMsg2 );call SetErrStat(ErrStat2, ErrMsg2, ErrStat,ErrMsg,RoutineName); m%Sgmt%Points = -999999_ReKi; - call AllocAry( m%Sgmt%Gamma , nSeg, 'SegGamma' , ErrStat2, ErrMsg2 );call SetErrStat(ErrStat2, ErrMsg2, ErrStat,ErrMsg,RoutineName); m%Sgmt%Gamma = -999999_ReKi; - call AllocAry( m%Sgmt%Epsilon, nSeg, 'SegEpsilon', ErrStat2, ErrMsg2 );call SetErrStat(ErrStat2, ErrMsg2, ErrStat,ErrMsg,RoutineName); m%Sgmt%Epsilon= -999999_ReKi; + call AllocAry( m%Sgmt%Connct, 4, nSeg , 'SegConnct' , ErrStat2, ErrMsg2 ); if(Failed())return; m%Sgmt%Connct = -999; + call AllocAry( m%Sgmt%Points, 3, nSegP, 'SegPoints' , ErrStat2, ErrMsg2 ); if(Failed())return; m%Sgmt%Points = -999999_ReKi; + call AllocAry( m%Sgmt%Gamma , nSeg, 'SegGamma' , ErrStat2, ErrMsg2 ); if(Failed())return; m%Sgmt%Gamma = -999999_ReKi; + call AllocAry( m%Sgmt%Epsilon, nSeg, 'SegEpsilon', ErrStat2, ErrMsg2 ); if(Failed())return; m%Sgmt%Epsilon= -999999_ReKi; m%Sgmt%nAct = -1 ! Active segments m%Sgmt%nActP = -1 m%Sgmt%RegFunction = p%RegFunction - bWakeNeedsPart = p%VelocityMethod(1)==idVelocityPart .or.p%VelocityMethod(1)==idVelocityTreePart + 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 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 );call SetErrStat(ErrStat2, ErrMsg2, ErrStat,ErrMsg,RoutineName); m%Part%P = -999999_ReKi; - call AllocAry( m%Part%Alpha , 3, nPart, 'PartAlpha' , ErrStat2, ErrMsg2 );call SetErrStat(ErrStat2, ErrMsg2, ErrStat,ErrMsg,RoutineName); m%Part%Alpha = -999999_ReKi; - call AllocAry( m%Part%RegParam, nPart, 'PartEpsilon', ErrStat2, ErrMsg2 );call SetErrStat(ErrStat2, ErrMsg2, ErrStat,ErrMsg,RoutineName); m%Part%RegParam= -999999_ReKi; + call AllocAry( m%Part%P , 3, nPart, 'PartP' , ErrStat2, ErrMsg2 ); if(Failed())return; m%Part%P = -999999_ReKi; + call AllocAry( m%Part%Alpha , 3, nPart, 'PartAlpha' , ErrStat2, ErrMsg2 ); if(Failed())return; m%Part%Alpha = -999999_ReKi; + call AllocAry( m%Part%RegParam, nPart, 'PartEpsilon', ErrStat2, ErrMsg2 ); if(Failed())return; m%Part%RegParam= -999999_ReKi; m%Part%nAct = -1 ! Active particles m%Part%RegFunction = p%RegFunction endif ! TODO Figure out Uind, CPs needed for grid - call AllocAry( m%CPs , 3, nCPs, 'CPs' , ErrStat2, ErrMsg2 );call SetErrStat(ErrStat2, ErrMsg2, ErrStat,ErrMsg,RoutineName); m%CPs= -999999_ReKi; - call AllocAry( m%Uind , 3, nCPs, 'Uind' , ErrStat2, ErrMsg2 );call SetErrStat(ErrStat2, ErrMsg2, ErrStat,ErrMsg,RoutineName); m%Uind= -999999_ReKi; - + call AllocAry( m%CPs , 3, nCPs, 'CPs' , ErrStat2, ErrMsg2 ); if(Failed())return; m%CPs= -999999_ReKi; + 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') + Failed = ErrStat >= AbortErrLev + end function Failed end subroutine FVW_InitMiscVarsPostParam !> Count how many segments are needed to represent the Near wake and far wakes, starting at a given depth @@ -1117,10 +1141,12 @@ subroutine InducedVelocitiesAll_OnGrid(g, p, x, m, ErrStat, ErrMsg) real(ReKi) :: xP,yP,zP,dx,dy,dz ! TODO new options type(T_Tree) :: Tree + type(T_Panl) :: Panl real(ReKi), dimension(:,:), allocatable :: CPs ! TODO get rid of me with dedicated functions real(ReKi), dimension(:,:), allocatable :: Uind ! TODO get rid of me with dedicated functions ErrStat= ErrID_None ErrMsg ='' + if (OLAF_PROFILING) call tic('InducedVelocitiesAll_OnGrid') ! --- Packing control points nCPs = g%nx * g%ny * g%nz @@ -1148,9 +1174,9 @@ subroutine InducedVelocitiesAll_OnGrid(g, p, x, m, ErrStat, ErrMsg) ! --- Compute induced velocity ! Convert Panels to segments, segments to particles, particles to tree - call InducedVelocitiesAll_Init(p, x, m, m%Sgmt, m%Part, Tree, ErrStat, ErrMsg, allocPart=.false.) - call InducedVelocitiesAll_Calc(CPs, nCPs, Uind, p, m%Sgmt, m%Part, Tree, ErrStat, ErrMsg) - call InducedVelocitiesAll_End(p, Tree, m%Part, ErrStat, ErrMsg, deallocPart=.false.) + call InducedVelocitiesAll_Init(p, x, m, m%Sgmt, m%Part, Tree, Panl, ErrStat, ErrMsg, allocPart=.false.) + call InducedVelocitiesAll_Calc(CPs, nCPs, Uind, p, m%Sgmt, m%Part, Tree, Panl, ErrStat, ErrMsg) + call InducedVelocitiesAll_End(p, Tree, m%Part, Panl, ErrStat, ErrMsg, deallocPart=.false.) ! --- Unpacking induced velocity points iHeadP=1 @@ -1159,6 +1185,8 @@ subroutine InducedVelocitiesAll_OnGrid(g, p, x, m, ErrStat, ErrMsg) if(allocated(CPs )) deallocate(CPs , stat=ErrStat) if(allocated(Uind)) deallocate(Uind, stat=ErrStat) + if (OLAF_PROFILING) call toc() + end subroutine InducedVelocitiesAll_OnGrid !> Wrapper to setup part from set of segments @@ -1213,13 +1241,14 @@ end subroutine SegmentsToPartWrap !> Perform initialization steps before requesting induced velocities from All vortex elements !! In : x%W(iW)%r_NW, x%W(iW)%r_FW, x%W(iW)%Gamma_NW, x%W(iW)%Gamma_FW !! Out: Tree, Part, m -subroutine InducedVelocitiesAll_Init(p, x, m, Sgmt, Part, Tree, ErrStat, ErrMsg, allocPart) - type(FVW_ParameterType), intent(in ) :: p !< Parameters +subroutine InducedVelocitiesAll_Init(p, x, m, Sgmt, Part, Tree, Panl, ErrStat, ErrMsg, allocPart) + type(FVW_ParameterType), intent(in ), target :: p !< Parameters type(FVW_ContinuousStateType), intent(in ) :: x !< States - type(FVW_MiscVarType), intent(in ) :: m !< Misc + type(FVW_MiscVarType), intent(in ), target :: m !< Misc type(T_Sgmt), intent(inout) :: Sgmt !< Segments type(T_Part), intent(inout) :: Part !< Particle storage if needed type(T_Tree), intent(out) :: Tree !< Tree of particles if needed + type(T_Panl), intent(out ) :: Panl !< Panel storage if needed integer(IntKi), intent( out) :: ErrStat !< Error status of the operation character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None logical, intent(in ) :: allocPart !< allocate particles @@ -1251,16 +1280,22 @@ subroutine InducedVelocitiesAll_Init(p, x, m, Sgmt, Part, Tree, ErrStat, ErrMsg call grow_tree_segment(Tree, nSeg, Sgmt%Points, Sgmt%Connct(:,1:nSeg), Sgmt%Gamma(1:nSeg), p%RegFunction, Sgmt%Epsilon(1:nSeg), 0) endif + ! --- Src + Panl%p_Src => p%SrcPnl + Panl%m_Src => m%SrcPnl + !Panl%z_Src => z%SrcPanel !we will use RHS for sigmas.. + end subroutine InducedVelocitiesAll_Init !> Compute induced velocity on flat CPs -subroutine InducedVelocitiesAll_Calc(CPs, nCPs, Uind, p, Sgmt, Part, Tree, ErrStat, ErrMsg) +subroutine InducedVelocitiesAll_Calc(CPs, nCPs, Uind, p, Sgmt, Part, Tree, Panl, ErrStat, ErrMsg) real(ReKi), dimension(:,:), intent(in) :: CPs !< Control points (3 x nCPs++) integer(IntKi) , intent(in) :: nCPs !< Number of control points on which to compute (nCPs <= size(CPs,2)) real(ReKi), dimension(:,: ) , intent(inout) :: Uind !< Induced velocity vector - Side effects!!! (3 x nCPs++) type(FVW_ParameterType), intent(in ) :: p !< Parameters type(T_Sgmt), intent(in ) :: Sgmt !< Segments type(T_Part), intent(in ) :: Part !< Particle storage if needed + type(T_Panl), intent(in ) :: Panl !< Panel storage if needed type(T_Tree), intent(inout) :: Tree !< Tree of particles if needed integer(IntKi), intent( out) :: ErrStat !< Error status of the operation character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None @@ -1283,15 +1318,24 @@ subroutine InducedVelocitiesAll_Calc(CPs, nCPs, Uind, p, Sgmt, Part, Tree, ErrSt elseif (p%VelocityMethod(iVel)==idVelocityTreeSeg) then call ui_tree_segment(Tree, CPs, nCPs, p%TreeBranchFactor(iVel), Tree%DistanceDirect, Uind, ErrStat, ErrMsg) endif + + ! --- Src Panels + if (associated(Panl%p_Src)) then + if (Panl%p_Src%n>0) then + call ui_quad_src_nn(CPs, Panl%m_Src%RHS, Panl%p_Src%xi, Panl%p_Src%eta, Panl%p_Src%Pcent, Panl%p_Src%R_g2p, Uind, nCPs, Panl%p_Src%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) + endif + endif end subroutine InducedVelocitiesAll_Calc !> Perform termination steps after velocity was requested from all vortex elements !! InOut: Tree, Part, m -subroutine InducedVelocitiesAll_End(p, Tree, Part, ErrStat, ErrMsg, deallocPart) +subroutine InducedVelocitiesAll_End(p, Tree, Part, Panl, ErrStat, ErrMsg, deallocPart) type(FVW_ParameterType), intent(in ) :: p !< Parameters type(T_Tree), intent(inout) :: Tree !< Tree of particles if needed type(T_Part), intent(inout) :: Part !< Particle storage if needed + type(T_Panl), intent(inout) :: Panl !< Panel storage if needed integer(IntKi), intent( out) :: ErrStat !< Error status of the operation character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None logical, intent(in ) :: deallocPart @@ -1314,6 +1358,10 @@ subroutine InducedVelocitiesAll_End(p, Tree, Part, ErrStat, ErrMsg, deallocPart) call cut_tree(Tree) ! We do not deallocate segment endif + ! Src Panels (we nullify only) + nullify(Panl%p_Src) + nullify(Panl%m_Src) + nullify(Panl%z_Src) end subroutine InducedVelocitiesAll_End @@ -1333,6 +1381,7 @@ subroutine WakeInducedVelocities(p, x, m, ErrStat, ErrMsg) integer(IntKi) :: nFWEff ! Number of farwake panels that are free at current time step integer(IntKi) :: nNWEff ! Number of nearwake panels that are free at current time step type(T_Tree) :: Tree + type(T_Panl) :: Panl if (OLAF_PROFILING) call tic('WakeInduced Calc') ErrStat= ErrID_None ErrMsg ='' @@ -1347,9 +1396,9 @@ subroutine WakeInducedVelocities(p, x, m, ErrStat, ErrMsg) ! Convert Panels to segments, segments to particles, particles to tree m%Uind=0.0_ReKi ! very important due to side effects of ui_* methods m%Uind(:,nCPs+1:)=1000.0_ReKi ! TODO For debugging only - call InducedVelocitiesAll_Init(p, x, m, m%Sgmt, m%Part, Tree, ErrStat, ErrMsg, allocPart=.false.) - call InducedVelocitiesAll_Calc(m%CPs, nCPs, m%Uind, p, m%Sgmt, m%Part, Tree, ErrStat, ErrMsg) - call InducedVelocitiesAll_End(p, Tree, m%Part, ErrStat, ErrMsg, deallocPart=.false.) + call InducedVelocitiesAll_Init(p, x, m, m%Sgmt, m%Part, Tree, Panl, ErrStat, ErrMsg, allocPart=.false.) + call InducedVelocitiesAll_Calc(m%CPs, nCPs, m%Uind, p, m%Sgmt, m%Part, Tree, Panl, ErrStat, ErrMsg) + call InducedVelocitiesAll_End(p, Tree, m%Part, Panl, ErrStat, ErrMsg, deallocPart=.false.) call UnPackInducedVelocity() if (DEV_VERSION) then @@ -1524,6 +1573,10 @@ subroutine LiftingLineInducedVelocities(p, x, InductionAtCP, iDepthStart, m, Err !deallocate(Part%P, Part%Alpha, Part%RegParam) call cut_tree(Tree) endif + ! --- Src Panel contribution + if (p%SrcPnl%n>0) then + call ui_quad_src_nn(CPs, m%SrcPnl%RHS, p%SrcPnl%xi, p%SrcPnl%eta, p%SrcPnl%Pcent, p%SrcPnl%R_g2p, Uind, nCPs, p%SrcPnl%n) + endif ! --- Unpack call UnPackLiftingLineVelocities() @@ -1608,13 +1661,11 @@ subroutine FakeGroundEffect(p, x, m, ErrStat, ErrMsg) ErrStat = ErrID_None ErrMsg = "" + GROUND = 1.e-4_ReKi + p%zGround if ( p%MHK /= MHK_None ) then - GROUND = 1.e-4_ReKi - p%WtrDpth - ABOVE_GROUND = 0.1_ReKi - p%WtrDpth - else - GROUND = 1.e-4_ReKi - ABOVE_GROUND = 0.1_ReKi + GROUND = GROUND - p%WtrDpth endif + ABOVE_GROUND = GROUND + p%zGroundPush ! Typically 0.1m above ground nBelow=0 nBelowFW=0 @@ -1742,4 +1793,415 @@ subroutine AlphaVrel_Generic(M_ag, Vstr_g, Vind_g, Vwnd_g, KinVisc, Chord, Vrel end subroutine AlphaVrel_Generic +!> 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" +subroutine srcPnl_init(p, m, z, errStat, errMsg, filename) + use VTK !, only: ReadVTK_PD_info, ReadVTK_PD_fields + type(T_SrcPanlParam), intent(inout) :: p + type(T_SrcPanlMisc) , intent(inout) :: m + type(T_SrcPanlVar) , intent(inout) :: z + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None + character(len=*), intent(in), optional :: filename + integer(IntKi) :: fid !< Unit number for file reading + character(ErrMsgLen) :: descr !< Description in input file + integer(IntKi) :: errStat2 !< temporary Error status + character(ErrMsgLen) :: errMsg2 !< temporary Error message + type(vtk_field), pointer :: fields + type(vtk_field), pointer :: bodyID + integer :: nPanels + nullify(fields) + nullify(bodyID) + fid = 0 + + errStat = ErrID_None + errMsg = "" + fid = 0 + ! --- Read points and connectivity + if (present(filename)) then + if (len_trim(filename)>0) then + !call EllipsoidPanels(17, 15, 1., 1., 1., p%P, p%IDs) + ! --- Read Source Panels + call WrScr(' - Source panel file: '//trim(filename)) + call ReadVTK_PD_info(filename, descr, p%P, p%IDs, fid, errStat2, errMsg2); if(Failed()) return + call ReadVTK_CD_fields(filename, fid, size(p%IDs, 2), fields, errStat2, errMsg2); if(Failed()) return + if (errStat2==ErrID_Warn) then + call WrScr('[WARN] Fields not found in VTK file, see '//trim(errMsg)) + errMsg='' + endif + if (fid>0) close(fid) + + ! --- Extract Body ID + !call VTK_printFields(fields) + call VTK_getField(fields, 'BODYID', bodyID) ! If not found, we'll replace by dummy ID + nPanels = size(p%IDs, 2) + call AllocAry(p%BodyIDs, nPanels, 'BodyIDs', errStat2, errMsg2); if(Failed()) return + if (associated(bodyID)) then + p%BodyIDs(1:nPanels) = bodyID%values(1,1:nPanels) + else + call WrScr('[WARN] Body IDs not found in source panel input file') + p%BodyIDs(:) = 1 + endif + call VTK_destroyFields(fields) + print'(A,I0,A,I0)',' - Number of panels: ', nPanels, ' , Points per panels: ',size(p%IDs,1) + p%IDs = p%IDs+1 ! VTK index starts at 0 + + else + errStat2 = ErrID_Fatal + errMsg2 = 'Filename should not be of length zero if provided.' + if (Failed()) return + endif + endif + if (.not. allocated(p%P) .or. .not.allocated(p%IDs)) then + errStat2 = ErrID_Fatal + errMsg2 = 'Points and connectivity (`IDs`) should be allocated (when an input file for source panels is not provided).' + if (Failed()) return + endif + + ! --- geometrical parameters (Area, normal vector, transformation matrices, etc) + p%Comment = '' + p%n = 0 + call srcPnl_geometry(p, errStat2, errMsg2); if(Failed()) return + + if (p%n>0) then + ! --- Allocations + call AllocAry(m%AI , p%n, p%n, 'AI ', errStat2, errMsg2); if(Failed())return + call AllocAry(m%UUI , 3, p%n, p%n, 'UUI ', errStat2, errMsg2); if(Failed())return + call AllocAry(m%Uwnd , 3 , p%n, 'Uwnd', errStat2, errMsg2); if(Failed())return + call AllocAry(m%Uext , 3 , p%n, 'Uext', errStat2, errMsg2); if(Failed())return + call AllocAry(m%Uind , 3 , p%n, 'Uind', errStat2, errMsg2); if(Failed())return + call AllocAry(m%Utot , 3 , p%n, 'Utot', errStat2, errMsg2); if(Failed())return + call AllocAry(m%F , 3 ,p%n, 'Fs ', errStat2, errMsg2); if(Failed())return + call AllocAry(m%FpA , 3 ,p%n, 'FpA ', errStat2, errMsg2); if(Failed())return + call AllocAry(m%Cp ,p%n, 'Cps ', errStat2, errMsg2); if(Failed())return + call AllocAry(m%p ,p%n, 'ps ', errStat2, errMsg2); if(Failed())return + call AllocAry(m%RHS ,p%n, 'RHS ', errStat2, errMsg2); if(Failed())return + call AllocAry(m%IPIV ,p%n, 'IPIV', errStat2, errMsg2); if(Failed())return + call AllocAry(z%Sigma ,p%n, 'Sigm', errStat2, errMsg2); if(Failed())return + m%AI = -999999_ReKi + m%UUI = -999999_ReKi + m%Uwnd = 0_ReKi + m%Uext = 0_ReKi + m%Uind = -999999_ReKi + m%Utot = -999999_ReKi + m%F = 0_ReKi + m%Cp = 0_ReKi + m%p = 0_ReKi + m%RHS = 0_ReKi + m%IPIV = -999999_ReKi + z%Sigma= 0_ReKi + + ! --- 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) + + ! --- Factorization + call linalg_factor(m%AI, m%IPIV, errStat2, errMsg2); if(Failed()) return + endif +contains + logical function Failed() + call SeterrStat(errStat2, errMsg2, errStat, errMsg, 'srcPnl_init') + Failed = errStat >= AbortErrLev + if (Failed .and. fid>0) close(fid) + end function +end subroutine srcPnl_init + +!> Geometrical parameters (Area, normal vector, transformation matrices, etc) +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), dimension(3) :: T1 + real(ReKi), dimension(3) :: T2 + real(ReKi), dimension(3) :: Ptmp !< temp for computing norm of delta points + real(ReKi), dimension(3,3) :: Mat + real(ReKi) :: norm_1 + real(ReKi) :: norm_T1 + real(ReKi) :: norm_T2 + integer :: ip + real(ReKi) :: xi0 !< + integer(IntKi) :: errStat2 !< temporary Error status + character(ErrMsgLen) :: errMsg2 !< temporary Error message + errStat = ErrID_None + errMsg = "" + Panl%n = size(Panl%IDs, 2) ! Number of panels + + ! --- Allocate storage for geometrical variables + if (allocated(Panl%R_g2p)) deallocate(Panl%R_g2p) + if (allocated(Panl%Pmid)) deallocate(Panl%Pmid) + if (allocated(Panl%Pcent)) deallocate(Panl%Pcent) + if (allocated(Panl%Normal))deallocate(Panl%Normal) + if (allocated(Panl%eta)) deallocate(Panl%eta) + if (allocated(Panl%xi)) deallocate(Panl%xi) + if (allocated(Panl%Area)) deallocate(Panl%Area) + + call AllocAry(Panl%R_g2p, 3, 3 ,Panl%n,'Normal',errStat2,errMsg2); if(Failed())return + call AllocAry(Panl%Pmid , 3, Panl%n,'Pmid ' ,errStat2,errMsg2); if(Failed())return + call AllocAry(Panl%Pcent , 3, Panl%n,'Pcent' ,errStat2,errMsg2); if(Failed())return + call AllocAry(Panl%Normal, 3, Panl%n,'Normal',errStat2,errMsg2); if(Failed())return + 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)) + P2 = Panl%P(:,IDs(2)) + P3 = Panl%P(:,IDs(3)) + P4 = Panl%P(:,IDs(4)) + Panl%Pmid(1:3,ip) = (P1+P2+P3+p4)/4 ! that's hess's barred coordinates + T1 = P3-P1 + T2 = P4-P2 + ! 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 + 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) + norm_1=sqrt(Ptmp(1)**2+ Ptmp(2)**2+ Ptmp(3)**2) + Panl%Normal(1:3,ip) = Ptmp/norm_1 !z or zeta axis + ! CP will be centroid ref (at the end, the middle) for Source panels + T1 = T1/norm_T1 ! x or xi axis + ! T2 = cross(N,T1) + T2(1) = Panl%Normal(2,ip) * T1(3) - Panl%Normal(3,ip) * T1(2) + T2(2) = Panl%Normal(3,ip) * T1(1) - Panl%Normal(1,ip) * T1(3) + T2(3) = Panl%Normal(1,ip) * T1(2) - Panl%Normal(2,ip) * T1(1) + norm_T2 = sqrt(T2(1)**2+ T2(2)**2+ T2(3)**2) + T2 = T2/norm_T2 ! y or eta axis + Panl%R_g2p(1, 1:3, ip) = T1 + 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 + 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 + 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 + 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 + P1e = P1es-(/ xi0,eta0,0.0_ReKi /) + P2e = P2es-(/ xi0,eta0,0.0_ReKi /) + P3e = P3es-(/ xi0,eta0,0.0_ReKi /) + P4e = P4es-(/ xi0,eta0,0.0_ReKi /) + Panl%xi(:,ip) = (/ P1e(1), P2e(1), P3e(1), P4e(1)/) + 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 + 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 + logical function Failed() + call SeterrStat(errStat2, errMsg2, errStat, errMsg, 'srcPnl_geometry') + Failed = errStat >= AbortErrLev + end function +end subroutine srcPnl_geometry + +subroutine srcPnl_build_mat(Panl, AI, UUI) + ! Arguments + type(T_SrcPanlParam), intent(in) :: Panl + 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 !< + integer :: icp, ip !< loop variables + if (OLAF_PROFILING) call tic('SrcPanel build matrix') + !$OMP PARALLEL DEFAULT(shared) + !$OMP DO PRIVATE(icp, ip, Uind_tmp) schedule(runtime) + ! Loop on control points + do icp = 1, Panl%n + ! Control Point (for a given panel) + ! ---- 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(icp, ip) = dot_product(Uind_tmp, Panl%Normal(1:3,icp)) + UUI(:, icp, ip) = Uind_tmp + end do + end do + !$OMP END DO + !$OMP END PARALLEL + if (OLAF_PROFILING) call toc() +end subroutine srcPnl_build_mat + + +!> Compute Wind + induced velocities from all vortex elements except panels on panel points +!! In : Control points p%SrcPnl%Pcent +!! Out: induced velocity m%SrcPnl%Uext +subroutine srcPnl_ExtVelocities_OnPanels(u, p, x, m, errStat, errMsg) + type(FVW_InputType), intent(in ) :: u !< Inputs + type(FVW_ParameterType), intent(in ) :: p !< Parameters + type(FVW_ContinuousStateType), intent(in ) :: x !< States + type(FVW_MiscVarType), intent(inout) :: m !< Initial misc/optimization variables + integer(IntKi), intent( out) :: errStat !< Error status of the operation + character(*), intent( out) :: errMsg !< Error message if ErrStat /= ErrID_None + integer :: icpp + ! TODO new options + type(T_Tree) :: Tree + type(T_Panl) :: Panl + errStat= ErrID_None + errMsg ='' + if (OLAF_PROFILING) call tic('SrcPanel ExtVelocities') + ! --- Induced velocities from all but panels + 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 + 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.) + + ! Compute external inflow m%SrcPnl%Uwnd + CALL DistributeRequestedWind_Panels(u%V_wind, p, m) + ! Combine Free stream and external velocity + m%SrcPnl%Uext = m%SrcPnl%Uext + m%SrcPnl%Uwnd + + if (OLAF_PROFILING) call toc() +end subroutine srcPnl_ExtVelocities_OnPanels + +!> Compute the value of the source panel (sigma) +!! assumes that Uext and AI were computed before +subroutine srcPnl_solve(p, m, z, errStat, errMsg) + type(T_SrcPanlParam), intent(in ) :: p + type(T_SrcPanlMisc) , intent(inout) :: m + type(T_SrcPanlVar) , intent(inout) :: z + integer(IntKi) , intent(out) :: errStat ! < Error status of the operation + character(errMsgLen), intent(out) :: errMsg ! < Error message if errStat / = ErrID_None + integer :: icpp + if (OLAF_PROFILING) call tic('SrcPanel Solve') + ! Compute RHS + do icpp = 1, p%n + m%RHS(icpp) = -dot_product(m%Uext(1:3,icpp), p%Normal(1:3,icpp)) + enddo + ! Solve for source intensities Sigma + call linalg_solve(m%AI, m%RHS, m%IPIV, errStat, errMsg) + z%Sigma = m%RHS + if (OLAF_PROFILING) call toc() +end subroutine srcPnl_solve + + +!> Compute velocity, pressure, loads on the source panels +subroutine srcPnl_calcOutput(p, m, z, rho) !, errStat, errMsg + type(T_SrcPanlParam), intent(in ) :: p + type(T_SrcPanlMisc) , intent(inout) :: m + type(T_SrcPanlVar) , intent(in ) :: z + real(ReKi) , intent(in ) :: rho +! integer(IntKi) , intent(out) :: errStat ! < Error status of the operation +! character(errMsgLen), intent(out) :: errMsg ! < Error message if errStat / = ErrID_None + integer :: ip + real(ReKi) :: Uwnd_norm2 + real(ReKi) :: Cp, qinf + if (OLAF_PROFILING) call tic('SrcPanel CalcOutput') + + ! --- Compute induced and total velocity + m%Uind(1,:) = matmul(m%UUI(1,:,:), m%RHS) + m%Uind(2,:) = matmul(m%UUI(2,:,:), m%RHS) + m%Uind(3,:) = matmul(m%UUI(3,:,:), m%RHS) + m%Utot = m%Uind + m%Uext + + ! --- Compute loads + do ip = 1, p%n + Uwnd_norm2 = norm2(m%Uwnd(1:3,ip))**2 + qinf = (0.5_ReKi * rho * Uwnd_norm2) ! TODO use Vinfty... + ! 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 + if (Uwnd_norm2>0) then + Cp = max( 1-(norm2(m%Utot(1:3,ip))**2)/Uwnd_norm2, -10._ReKi) ! Bernoulli. + else + Cp = -10._ReKi !!! + endif + !m%p (ip) = (1.0_ReKi - Cp) * qinf + m%p (ip) = Cp * qinf ! Delta p + m%F (1:3,ip) = - Cp * p%Area(ip) * qinf * p%Normal(1:3,ip) + m%FpA(1:3,ip)= m%F (1:3,ip)/p%Area(ip) + m%Cp(ip) = Cp + enddo + if (OLAF_PROFILING) call toc() +endsubroutine srcPnl_calcOutput + +! --- Linear Solve module +subroutine linalg_factor(AA, IPIV, errStat, errMsg) + use NWTC_LAPACK, only : LAPACK_GETRF + real(ReKi), dimension(:, :), intent(inout) :: AA + integer(IntKi),dimension(:), intent(inout) :: IPIV + integer(IntKi), intent(out) :: errStat !< Error status of the operation + character(*), intent(out) :: errMsg !< Error message if ErrStat /= ErrID_None + integer :: m, n + errStat = ErrID_None + errMsg = '' + n = size(AA,1) + m = n + call LAPACK_GETRF(m, n, AA, IPIV, errStat, errMsg) +endsubroutine + +subroutine linalg_solve(AFact, RHS, IPIV, errStat, errMsg) + use NWTC_LAPACK, only : LAPACK_GETRS + real(ReKi), dimension(:, :), intent(in) :: AFact + real(ReKi), dimension(:), intent(inout) :: RHS + integer(IntKi),dimension(:), intent(inout) :: IPIV + integer(IntKi), intent(out) :: errStat !< Error status of the operation + character(*), intent(out) :: errMsg !< Error message if ErrStat /= ErrID_None + integer :: n + n = size(AFact,1) + call LAPACK_GETRS('N', n, AFact, IPIV, RHS, errStat, errMsg) +end subroutine + +!> Solve A x = B +subroutine linalg_solveWrap(AA, B, X, errStat, errMsg) + real(ReKi), dimension(:, :), intent(in) :: AA !< + real(ReKi), dimension(:), intent(in) :: B !< + real(ReKi), dimension(:), intent(out) :: X !< + integer(IntKi), intent(out) :: errStat !< Error status of the operation + character(*), intent(out) :: errMsg !< Error message if ErrStat /= ErrID_None + real(ReKi), dimension(:, :), allocatable :: AA2 + real(ReKi), dimension(:), allocatable :: B2 + integer(IntKi), dimension(:), allocatable :: IPIV + allocate(AA2(size(AA,1), size(AA,2))) + allocate(B2 (size(B,1))) + allocate(IPIV(size(B,1))) + AA2 = AA + B2 = B + call linalg_factor(AA2, IPIV, errStat, errMsg) + call linalg_solve (AA2, B2, IPIV, errStat, errMsg) + X(:) = B2(:) + deallocate(AA2) + deallocate(B2) + deallocate(IPIV) +endsubroutine + end module FVW_Subs diff --git a/modules/aerodyn/src/FVW_Tests.f90 b/modules/aerodyn/src/FVW_Tests.f90 index e2299e97e1..1edc14e935 100644 --- a/modules/aerodyn/src/FVW_Tests.f90 +++ b/modules/aerodyn/src/FVW_Tests.f90 @@ -1,7 +1,5 @@ module FVW_Tests - use NWTC_Library - use FVW_Types use FVW_Subs use FVW_VortexTools @@ -13,216 +11,184 @@ module FVW_Tests implicit none public :: FVW_RunTests + public :: Test_LinSolve + public :: Test_SrcPnl_Sphere + public :: Test_BiotSavart_SrcPnl + public :: Test_BiotSavart_Sgmt + public :: Test_BiotSavart_Part + public :: Test_BiotSavart_PartTree + public :: Test_SegmentsToPart + public :: FVW_Test_WakeInducedVelocities + + private + character(len=*), parameter :: testname = 'FVW' interface test_equal; module procedure & test_equal_i1, & test_equal_i0 end interface - interface test_almost_equal; module procedure & - test_almost_equal_0, & - test_almost_equal_1, & - test_almost_equal_2 - end interface + interface test_almost_equal; module procedure & + test_almost_equal_0, & + test_almost_equal_1, & + test_almost_equal_2 + end interface contains ! -------------------------------------------------------------------------------- ! --- Helper functions (should be part of NWTC library) ! -------------------------------------------------------------------------------- - subroutine test_success(testname,info,bPrint_in) - character(len=*), intent(in) :: testname - character(len=*), intent(in) :: info - logical, intent(in), optional :: bPrint_in - if(present(bPrint_in)) then - if(bPrint_in) then - write(*,'(A)')'[ OK ] '//trim(testname)//': '//trim(Info) - endif - else - write(*,'(A)')'[ OK ] '//trim(testname)//': '//trim(Info) - endif - end subroutine - - subroutine test_fail(testname,info,bPrint_in,bStop_in) - character(len=*), intent(in) :: testname - character(len=*), intent(in) :: info - logical, intent(in), optional :: bPrint_in - logical, intent(in), optional :: bStop_in - if(present(bPrint_in)) then - if(bPrint_in) then - write(*,'(A)')'[FAIL] '//trim(testname)//': '//trim(Info) - endif - else - write(*,'(A)')'[FAIL] '//trim(testname)//': '//trim(Info) - endif - if(present(bStop_in)) then - if(bStop_in) then - STOP - endif - else - STOP - endif - end subroutine - - subroutine test_equal_i0(testname,Var,iTry,iRef) - ! Arguments - character(len=*), intent(in) :: testname - character(len=*), intent(in) :: Var - integer, intent(in) :: iTry !< - integer, intent(in) :: iRef !< - ! Variables - character(len=255) :: InfoAbs - if(iRef/=iTry) then - write(InfoAbs,'(A,I0,A,I0)') trim(Var),iRef,'/',iTry - call test_fail(testname,InfoAbs) - STOP - else - write(InfoAbs,'(A,A,I0)') trim(Var),' ok ',iRef - call test_success(testname,InfoAbs) - endif - end subroutine - - subroutine test_equal_i1(testname,Var,VecTry,VecRef,bTest,bPrintOnly,bPassed) - ! Arguments - character(len=*), intent(in) :: testname - character(len=*), intent(in) :: Var - integer, dimension(:), intent(in) :: VecTry !< - integer, dimension(:), intent(in) :: VecRef !< - logical, intent(in) :: bTest - logical, intent(in) :: bPrintOnly - logical, intent(out),optional :: bPassed - ! Variables - character(len=255) :: InfoAbs - integer :: i,cpt - ! - cpt=0 - do i=1,size(VecRef) - if(VecRef(i)/=VecTry(i)) then - cpt=cpt+1 - endif - enddo - if(cpt>0) then - write(InfoAbs,'(A,I0)') trim(Var)//' Elements different: ',cpt - if(present(bPassed)) then - bPassed=.false. - endif - else - write(InfoAbs,'(A)') trim(Var)//' reproduced to identity' - if(present(bPassed)) then - bPassed=.true. - endif - endif - if(bPrintOnly) then - print'(A)',trim(InfoAbs) - endif - if(bTest) then - if(cpt>0) then - call test_fail(testname,InfoAbs) - STOP - else - call test_success(testname,InfoAbs) - endif - endif - end subroutine - - subroutine test_almost_equal_0(testname,Var,Ref,Try,MINNORM,bStop,bPrint,bPassed) - ! Arguments - character(len=*), intent(in) :: testname - character(len=*), intent(in) :: Var - real(ReKi), intent(in) :: Ref !< - real(ReKi), intent(in) :: Try !< - real(ReKi), intent(in) :: MINNORM - logical, intent(in) :: bStop - logical, intent(in) :: bPrint - logical, intent(out),optional :: bPassed - ! Variables - character(len=255) :: InfoAbs - real(ReKi) :: delta - integer :: cpt - ! - cpt=0 - delta=abs(Ref-Try) - if(delta>MINNORM) then - write(InfoAbs,'(A,ES8.1E2,A,ES8.1E2,A,I0)') trim(Var)//' tol: ',MINNORM,', mean: ',delta,' - Failed:',cpt - call test_fail(testname,InfoAbs,bPrint,bStop) - else - write(InfoAbs,'(A,ES8.1E2,A,ES8.1E2)') trim(Var)//' tol: ',MINNORM,', mean: ',delta - call test_success(testname,InfoAbs,bPrint) - endif - if(present(bPassed)) then - bPassed=delta>MINNORM - endif - end subroutine - subroutine test_almost_equal_1(testname,Var,VecRef,VecTry,MINNORM,bStop,bPrint,bPassed) - ! Arguments - character(len=*), intent(in) :: testname - character(len=*), intent(in) :: Var - real(ReKi), dimension(:), intent(in) :: VecRef !< - real(ReKi), dimension(:), intent(in) :: VecTry !< - real(ReKi), intent(in) :: MINNORM - logical, intent(in) :: bStop - logical, intent(in) :: bPrint - logical, intent(out),optional :: bPassed - ! Variables - character(len=255) :: InfoAbs - integer :: i,cpt - real(ReKi) :: delta - real(ReKi) :: delta_cum - ! - cpt=0 - delta_cum=0.0_ReKi - do i=1,size(VecRef,1) - delta=abs(VecRef(i)-VecTry(i)) - delta_cum=delta_cum+delta - if(delta>MINNORM) then - cpt=cpt+1 - endif - enddo - delta_cum=delta_cum/size(VecRef) - - if(cpt>0) then - write(InfoAbs,'(A,ES8.1E2,A,ES8.1E2,A,I0)') trim(Var)//' tol: ',MINNORM,', mean: ',delta_cum,' - Failed:',cpt - call test_fail(testname,InfoAbs,bPrint,bStop) - else - write(InfoAbs,'(A,ES8.1E2,A,ES8.1E2)') trim(Var)//' tol: ',MINNORM,', mean: ',delta_cum - call test_success(testname,InfoAbs,bPrint) - endif - if(present(bPassed)) then - bPassed=(cpt==0) - endif - end subroutine - subroutine test_almost_equal_2(testname,Var,VecRef,VecTry,MINNORM,bStop,bPrint,bPassed) - ! Arguments - character(len=*), intent(in) :: testname - character(len=*), intent(in) :: Var - real(ReKi), dimension(:,:), intent(in) :: VecRef !< - real(ReKi), dimension(:,:), intent(in) :: VecTry !< - real(ReKi), intent(in) :: MINNORM - logical, intent(in) :: bStop - logical, intent(in) :: bPrint - logical, intent(out),optional :: bPassed - ! Variables - real(ReKi), dimension(:),allocatable :: VecRef2 !< - real(ReKi), dimension(:),allocatable :: VecTry2 !< - integer :: p, i,j,n1,n2,nCPs - ! - n1 = size(VecRef,1); n2 = size(VecRef,2); nCPs=n1*n2 - allocate ( VecRef2 (n1*n2) ) ; allocate ( VecTry2 (n1*n2) ) - p=0 - do j=1,n2; do i=1,n1 - p=p+1 - VecRef2(p)=VecRef(i,j) - VecTry2(p)=VecTry(i,j) - enddo; enddo; - call test_almost_equal(testname,Var,VecRef2,VecTry2,MINNORM,bStop,bPrint,bPassed) - end subroutine + subroutine test_success(info, errStat, errMsg) + character(len=*), intent(in) :: info + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None + write(*,'(A)')'[ OK ] '//trim(testname)//': '//trim(Info) + errStat = ErrID_None + errMsg = '' + end subroutine + + subroutine test_fail(info, errStat, errMsg) + character(len=*), intent(in) :: info + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None + write(*,'(A)')'[FAIL] '//trim(testname)//': '//trim(Info) + errStat = ErrID_Fatal + errMsg = '[FAIL] '//trim(testname)//': '//trim(Info) + end subroutine + + subroutine test_equal_i0(Var, iTry, iRef, errStat, errMsg) + ! Arguments + character(len=*), intent(in) :: Var + integer, intent(in) :: iTry !< + integer, intent(in) :: iRef !< + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None + ! Variables + character(len=255) :: InfoAbs + if(iRef/=iTry) then + write(InfoAbs,'(A,I0,A,I0)') trim(Var),iRef,'/',iTry + call test_fail(InfoAbs, errStat, errMsg) + else + write(InfoAbs,'(A,A,I0)') trim(Var),' ok ',iRef + call test_success(InfoAbs, errStat, errMsg) + endif + end subroutine + + subroutine test_equal_i1(Var, VecTry, VecRef, errStat, errMsg) + ! Arguments + character(len=*), intent(in) :: Var + integer, dimension(:), intent(in) :: VecTry !< + integer, dimension(:), intent(in) :: VecRef !< + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None + ! Variables + character(len=255) :: InfoAbs + integer :: i,cpt + ! + cpt=0 + do i=1,size(VecRef) + if(VecRef(i)/=VecTry(i)) then + cpt=cpt+1 + endif + enddo + if(cpt>0) then + write(InfoAbs,'(A,I0)') trim(Var)//' Elements different: ',cpt + else + write(InfoAbs,'(A)') trim(Var)//' reproduced to identity' + endif + if(cpt>0) then + call test_fail(InfoAbs, errStat, errMsg) + else + call test_success(InfoAbs, errStat, errMsg) + endif + end subroutine + + subroutine test_almost_equal_0(Var, Ref, Try, MINNORM, errStat, errMsg) + ! Arguments + character(len=*), intent(in) :: Var + real(ReKi), intent(in) :: Ref !< + real(ReKi), intent(in) :: Try !< + real(ReKi), intent(in) :: MINNORM + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None + ! Variables + character(len=255) :: InfoAbs + real(ReKi) :: delta + integer :: cpt + ! + cpt=0 + delta=abs(Ref-Try) + if(delta>MINNORM) then + write(InfoAbs,'(A,ES8.1E2,A,ES8.1E2,A,I0)') trim(Var)//' tol: ',MINNORM,', mean: ',delta,' - Failed:',cpt + call test_fail(InfoAbs, errStat, errMsg) + else + write(InfoAbs,'(A,ES8.1E2,A,ES8.1E2)') trim(Var)//' tol: ',MINNORM,', mean: ',delta + call test_success(InfoAbs, errStat, errMsg) + endif + end subroutine + subroutine test_almost_equal_1(Var, VecRef, VecTry, MINNORM, errStat, errMsg) + ! Arguments + character(len=*), intent(in) :: Var + real(ReKi), dimension(:), intent(in) :: VecRef !< + real(ReKi), dimension(:), intent(in) :: VecTry !< + real(ReKi), intent(in) :: MINNORM + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None + ! Variables + character(len=255) :: InfoAbs + integer :: i,cpt + real(ReKi) :: delta + real(ReKi) :: delta_cum + ! + cpt=0 + delta_cum=0.0_ReKi + do i=1,size(VecRef,1) + delta=abs(VecRef(i)-VecTry(i)) + delta_cum=delta_cum+delta + if(delta>MINNORM) then + cpt=cpt+1 + endif + enddo + delta_cum=delta_cum/size(VecRef) + + if(cpt>0) then + write(InfoAbs,'(A,ES8.1E2,A,ES8.1E2,A,I0)') trim(Var)//' tol: ',MINNORM,', mean: ',delta_cum,' - Failed:',cpt + call test_fail(InfoAbs, errStat, errMsg) + else + write(InfoAbs,'(A,ES8.1E2,A,ES8.1E2)') trim(Var)//' tol: ',MINNORM,', mean: ',delta_cum + call test_success(InfoAbs, errStat, errMsg) + endif + end subroutine + subroutine test_almost_equal_2(Var, VecRef, VecTry, MINNORM, errStat, errMsg) + ! Arguments + character(len=*), intent(in) :: Var + real(ReKi), dimension(:,:), intent(in) :: VecRef !< + real(ReKi), dimension(:,:), intent(in) :: VecTry !< + real(ReKi), intent(in) :: MINNORM + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None + ! Variables + real(ReKi), dimension(:),allocatable :: VecRef2 !< + real(ReKi), dimension(:),allocatable :: VecTry2 !< + integer :: p, i,j,n1,n2,nCPs + ! + n1 = size(VecRef,1); n2 = size(VecRef,2); nCPs=n1*n2 + allocate ( VecRef2 (n1*n2) ) ; allocate ( VecTry2 (n1*n2) ) + p=0 + do j=1,n2; do i=1,n1 + p=p+1 + VecRef2(p)=VecRef(i,j) + VecTry2(p)=VecTry(i,j) + enddo; enddo; + call test_almost_equal(Var,VecRef2,VecTry2,MINNORM, errStat, errMsg) + end subroutine ! --------------------------------------------------------------------------------} ! --- Specific FVW tests ! --------------------------------------------------------------------------------{ !> - subroutine Test_BiotSavart_Sgmt(testname, ErrStat, ErrMsg) - character(len=*), intent(in) :: testname - integer(IntKi) , intent(out) :: ErrStat !< Error status of the operation - character(ErrMsgLen), intent(out) :: ErrMsg !< Error message if ErrStat /= ErrID_None + subroutine Test_BiotSavart_Sgmt(errStat, errMsg) + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None real(ReKi), dimension(3) :: P1,P2,P3,CP real(ReKi), dimension(3) :: U1 real(ReKi) :: SegGamma1 !< Circulation [m^2/s] @@ -239,9 +205,9 @@ subroutine Test_BiotSavart_Sgmt(testname, ErrStat, ErrMsg) real(ReKi), dimension(nSegTot) :: RegParam !< Regularization parameter real(ReKi), dimension(3,nCPsTot) :: Uind_out !< Induced velocity vector - Side effects!!! real(ReKi), dimension(3,4) :: CPs_test !< - ! Initialize ErrStat - ErrStat = ErrID_None - ErrMsg = "" + ! Initialize errStat + errStat = ErrID_None + errMsg = "" ! --- Test that the two functions return the same values P1=(/0. ,0.,-1./) P2=(/0. ,0., 1./) @@ -274,7 +240,7 @@ subroutine Test_BiotSavart_Sgmt(testname, ErrStat, ErrMsg) !print*,'Reg function', RegFunction, 'CP',CP !print*,'Uind_out',Uind_out !print*,'U1 ',U1 - call test_almost_equal(testname,'Uind method1/2', U1, Uind_out(:,1), 1e-4_ReKi, .true.,.true.) + call test_almost_equal('Uind method1/2', U1, Uind_out(:,1), 1e-4_ReKi, errStat, errMsg) !call test_almost_equal('Uind method1/2', U1, Uind_out(:,1), 1e-4, .false.,.true.) enddo enddo @@ -313,16 +279,15 @@ subroutine Test_BiotSavart_Sgmt(testname, ErrStat, ErrMsg) !print*,'Reg function', RegFunction, 'CP',CP !print*,'Uind_out',Uind_out !print*,'U1 ',U1 - call test_almost_equal(testname,'Uind 1seg/2seg', U1, Uind_out(:,1), 1e-4_ReKi, .true.,.true.) + call test_almost_equal('Uind 1seg/2seg', U1, Uind_out(:,1), 1e-4_ReKi, errStat, errMsg) enddo enddo end subroutine !> - subroutine Test_BiotSavart_Part(testname, ErrStat, ErrMsg) - character(len=*), intent(in) :: testname - integer(IntKi) , intent(out) :: ErrStat !< Error status of the operation - character(ErrMsgLen), intent(out) :: ErrMsg !< Error message if ErrStat /= ErrID_None + subroutine Test_BiotSavart_Part(errStat, errMsg) + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None real(ReKi), dimension(3) :: P1,CP real(ReKi), dimension(3) :: U1 real(ReKi), dimension(3) :: PartAlpha1 !< Particle intensity alpha=om.dV [m^3/s] @@ -337,9 +302,9 @@ subroutine Test_BiotSavart_Part(testname, ErrStat, ErrMsg) real(ReKi), dimension(nPart) :: RegParam !< Regularization parameter real(ReKi), dimension(3,nCPs) :: Uind_out !< Induced velocity vector - Side effects!!! real(ReKi), dimension(3,4) :: CPs_test !< - ! Initialize ErrStat - ErrStat = ErrID_None - ErrMsg = "" + ! Initialize errStat + errStat = ErrID_None + errMsg = "" ! --- Test that the two functions return the same values P1=(/0.0, 0.0, -1.0 /) CPs_test(:,1) = (/ 0.0, 0., 0.0 /) ! Middle @@ -368,16 +333,15 @@ subroutine Test_BiotSavart_Part(testname, ErrStat, ErrMsg) !print*,'Reg function', RegFunction, 'CP',CP !print*,'Uind_out',Uind_out !print*,'U1 ',U1 - call test_almost_equal(testname,'Uind part method1/2', U1, Uind_out(:,1), 1e-4_ReKi, .true.,.true.) + call test_almost_equal('Uind part method1/2', U1, Uind_out(:,1), 1e-4_ReKi, errStat, errMsg) enddo enddo end subroutine Test_BiotSavart_Part !> This test compares calls using the tree algorithm and the direct N^2 evaluation - subroutine Test_BiotSavart_PartTree(testname, ErrStat, ErrMsg) - character(len=*), intent(in) :: testname - integer(IntKi) , intent(out) :: ErrStat !< Error status of the operation - character(ErrMsgLen), intent(out) :: ErrMsg !< Error message if ErrStat /= ErrID_None + subroutine Test_BiotSavart_PartTree(errStat, errMsg) + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None type(T_Tree) :: Tree real(ReKi), dimension(3) :: U_ref integer(IntKi) :: i1,i2,i3,k, iCP @@ -392,9 +356,9 @@ subroutine Test_BiotSavart_PartTree(testname, ErrStat, ErrMsg) real(ReKi), dimension(:,:), allocatable :: Uind2 !< Induced velocity vector - Side effects!!! real(ReKi) :: BranchFactor, BranchSmall real(ReKi), dimension(3,5) :: CPs_test !< - ! Initialize ErrStat - ErrStat = ErrID_None - ErrMsg = "" + ! Initialize errStat + errStat = ErrID_None + errMsg = "" BranchFactor = 2.0_ReKi !< Should be above1 BranchSmall = 0.0_ReKi RegFunction = 1 @@ -408,10 +372,10 @@ subroutine Test_BiotSavart_PartTree(testname, ErrStat, ErrMsg) U_ref =0.0_ReKi call grow_tree_part(Tree, nPart, PartPoints, PartAlpha, RegFunction, RegParam, 0) !call print_tree(Tree) - call ui_tree_part(Tree, nCPs, CPs, BranchFactor, BranchSmall, Uind2, ErrStat, ErrMsg) + call ui_tree_part(Tree, nCPs, CPs, BranchFactor, BranchSmall, Uind2, errStat, errMsg) call ui_part_nograd(nCPS, CPs, nPart, PartPoints, PartAlpha, RegFunction, RegParam, Uind1) ! Test - call test_almost_equal(testname,'Uind tree 0 part', U_ref, Uind2(:,1), 1e-4_ReKi, .true.,.true.) + call test_almost_equal('Uind tree 0 part', U_ref, Uind2(:,1), 1e-4_ReKi, errStat, errMsg) call cut_tree(Tree) call dealloc() @@ -424,10 +388,10 @@ subroutine Test_BiotSavart_PartTree(testname, ErrStat, ErrMsg) U_ref =0.0_ReKi call grow_tree_part(Tree, nPart, PartPoints, PartAlpha, RegFunction, RegParam, 0) !call print_tree(Tree) - call ui_tree_part(Tree, nCPs, CPs, BranchFactor, BranchSmall, Uind2, ErrStat, ErrMsg) + call ui_tree_part(Tree, nCPs, CPs, BranchFactor, BranchSmall, Uind2, errStat, errMsg) call ui_part_nograd(nCPS, CPs, nPart, PartPoints, PartAlpha, RegFunction, RegParam, Uind1) ! Test - call test_almost_equal(testname,'Uind tree 1 part', Uind1, Uind2, 1e-4_ReKi, .true.,.true.) + call test_almost_equal('Uind tree 1 part', Uind1, Uind2, 1e-4_ReKi, errStat, errMsg) call cut_tree(Tree) !call print_tree(Tree) call dealloc() @@ -457,16 +421,16 @@ subroutine Test_BiotSavart_PartTree(testname, ErrStat, ErrMsg) do iCP=1,4 CPs(:,1) = CPs_test(:,icp) Uind2=0.0_ReKi; Uind1=0.0_ReKi - call ui_tree_part(Tree, nCPs, CPs, BranchFactor, BranchSmall, Uind2, ErrStat, ErrMsg) + call ui_tree_part(Tree, nCPs, CPs, BranchFactor, BranchSmall, Uind2, errStat, errMsg) call ui_part_nograd(nCPs, CPs, nPart, PartPoints, PartAlpha, RegFunction, RegParam, Uind1) !print*,'Uind',Uind1, Uind2 ! Test - call test_almost_equal(testname,'Uind tree 81 part', Uind1, Uind2, 1e-2_ReKi, .true.,.true.) + call test_almost_equal('Uind tree 81 part', Uind1, Uind2, 1e-2_ReKi, errStat, errMsg) enddo call cut_tree(Tree) ! --- Test that tree ui cannot be called after tree has been cut - call ui_tree_part(Tree, nCPs, CPs, BranchFactor, BranchSmall, Uind2, ErrStat, ErrMsg) - call test_equal(testname,'Err. stat tree cut',ErrStat,ErrID_Fatal) + call ui_tree_part(Tree, nCPs, CPs, BranchFactor, BranchSmall, Uind2, errStat, errMsg) + call test_equal('Err. stat tree cut', errStat, ErrID_Fatal, errStat, errMsg) call dealloc() contains @@ -487,11 +451,49 @@ subroutine dealloc() end subroutine end subroutine Test_BiotSavart_PartTree + subroutine Test_BiotSavart_SrcPnl(errStat, errMsg) + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None + integer, parameter :: ncp=5 + integer, parameter :: np=1 + real(ReKi), dimension(3,ncp) :: CPs + real(ReKi), dimension(3,ncp) :: UI + real(ReKi), dimension(3,ncp) :: UI_ref + real(ReKi), dimension(3,ncp) :: Grad + real(ReKi), dimension(np) :: Sigmas + real(ReKi), dimension(3,np) :: RefPoint + real(ReKi), dimension(4,np) :: xi + real(ReKi), dimension(4,np) :: eta + real(ReKi), dimension(3,3,np) :: TransfoMat + errStat = ErrID_None + errMsg = "" + Sigmas(1) = 2._ReKi + xi (1:4,1)= [-1,0,1,0] + eta(1:4,1)= [0 ,1,0,-1] + TransfoMat = 0.0_ReKi + TransfoMat(1,1,1) = 1._ReKi + TransfoMat(2,2,1) = 1._ReKi + TransfoMat(3,3,1) = 1._ReKi + RefPoint(1:3,1) = [0,0,0] + UI = 0.0_ReKi + CPs(1:3,1)= [ 0.0_ReKi, 0.0_ReKi, 0.0_ReKi] + CPs(1:3,2)= [-1.0_ReKi, 1.0_ReKi, 1.0_ReKi] + CPs(1:3,3)= [ 1.0_ReKi,-1.0_ReKi, 1.0_ReKi] + CPs(1:3,4)= [ 2.0_ReKi, 2.0_ReKi, 2.0_ReKi] + CPs(1:3,5)= [ 0.0_ReKi, 0.0_ReKi, 1.0_ReKi] + call ui_quad_src_nn(CPs, Sigmas, xi, eta, RefPoint, TransfoMat, UI, ncp, np) + + UI_ref(1,:)= [0.0000e+00,-5.676174595996415e-02, 5.676174595996415e-02, 1.508218771767720e-02, 0.000000000000000e+00] + UI_ref(2,:)= [0.0000e+00, 5.676174595996415e-02 ,-5.676174595996415e-02, 1.508218771767721e-02, 0.000000000000000e+00] + UI_ref(3,:)= [1.0000e+00, 6.672740817112868e-02 , 6.672740817112868e-02, 1.571922940762083e-02, 2.163468959387855e-01] + call test_almost_equal('Uind src', UI, UI_ref, 1e-6_ReKi, errStat, errMsg) + end subroutine Test_BiotSavart_SrcPnl + + !> Compares the velocity field obtained from a segment and its convert to particle version - subroutine Test_SegmentsToPart(testname, ErrStat, ErrMsg) - character(len=*), intent(in) :: testname - integer(IntKi) , intent(out) :: ErrStat !< Error status of the operation - character(ErrMsgLen), intent(out) :: ErrMsg !< Error message if ErrStat /= ErrID_None + subroutine Test_SegmentsToPart(errStat, errMsg) + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None real(ReKi), dimension(:,:), allocatable :: PartPoints !< Particle points real(ReKi), dimension(:,:), allocatable :: PartAlpha !< Particle circulation real(ReKi), dimension(:) , allocatable :: PartEpsilon !< Regularization parameter @@ -508,8 +510,10 @@ subroutine Test_SegmentsToPart(testname, ErrStat, ErrMsg) real(ReKi) :: RegParam1 !< integer(IntKi) :: nPartPerSeg, nPart, iHeadP integer(IntKi) :: RegFunctionPart, RegFunctionSeg - ErrStat = ErrID_None - ErrMsg = "" + integer(IntKi) :: errStat2 !< Error status of the operation + character(errMsgLen) :: errMsg2 !< Error message if errStat /= ErrID_None + errStat = ErrID_None + errMsg = "" RegParam1=1.0 ! Creating two aligned segments SegConnct(:,1)=(/1,2/) @@ -544,7 +548,7 @@ subroutine Test_SegmentsToPart(testname, ErrStat, ErrMsg) Uind1 =0.0_ReKi; Uind2 =0.0_ReKi; call ui_seg(1, nCPsTot, CPs, 1, nSegTot, SegPoints, SegConnct, SegGamma, RegFunctionSeg, SegEpsilon, Uind1) call ui_part_nograd(nCPSTot, CPs, nPart, PartPoints, PartAlpha, RegFunctionPart, PartEpsilon, Uind2) - call test_almost_equal(testname,'Uind 10 part/sgmt no reg', Uind1, Uind2, 1e-3_ReKi, .true.,.true.) + call test_almost_equal('Uind 10 part/sgmt no reg', Uind1, Uind2, 1e-3_ReKi, errStat2, errMsg2); if(Failed())return call dealloc() ! --- Test 1 - 2 particles, no regularization @@ -560,7 +564,7 @@ subroutine Test_SegmentsToPart(testname, ErrStat, ErrMsg) Uind1 =0.0_ReKi; Uind2 =0.0_ReKi; call ui_seg(1, nCPsTot, CPs, 1, nSegTot, SegPoints, SegConnct, SegGamma, RegFunctionSeg, SegEpsilon, Uind1) call ui_part_nograd(nCPsTot, CPs, nPart, PartPoints, PartAlpha, RegFunctionPart, PartEpsilon, Uind2) - call test_almost_equal(testname,'Uind 2 part/sgmt noreg', Uind1, Uind2, 3e-1_ReKi, .true.,.true.) + call test_almost_equal('Uind 2 part/sgmt noreg', Uind1, Uind2, 3e-1_ReKi, errStat2, errMsg2); if(Failed())return call dealloc() @@ -584,7 +588,7 @@ subroutine Test_SegmentsToPart(testname, ErrStat, ErrMsg) !print'(A,10F7.3)','Uind2',Uind2(2,:) !print'(A,10F7.3)','Uind1',Uind1(3,:) !print'(A,10F7.3)','Uind2',Uind2(3,:) - call test_almost_equal(testname,'Uind 10 part/sgmt w.reg', Uind1, Uind2, 5e-2_ReKi, .true.,.true.) + call test_almost_equal('Uind 10 part/sgmt w.reg', Uind1, Uind2, 5e-2_ReKi, errStat, errMsg) call dealloc() contains @@ -598,6 +602,10 @@ subroutine alloc(n) subroutine dealloc() deallocate(PartPoints, PartAlpha, PartEpsilon) end subroutine + logical function Failed() + call SeterrStat(errStat2, errMsg2, errStat, errMsg, 'Test_SegmentsToPart') + Failed = errStat >= AbortErrLev + end function end subroutine Test_SegmentsToPart !> @@ -748,25 +756,25 @@ subroutine MeshMe(M,offset) !> Test Wake Induced velocity calcualtion when using nNWMax or nNWFree !! A dummy helical wake is created. The induced velocity is computed on !! either the full wake, or just the "free" wake (which should be way faster) - subroutine FVW_Test_WakeInducedVelocities(ErrStat, ErrMsg) - integer(IntKi) , intent(out) :: ErrStat !< Error status of the operation - character(ErrMsgLen), intent(out) :: ErrMsg !< Error message if ErrStat /= ErrID_None + subroutine FVW_Test_WakeInducedVelocities(errStat, errMsg) + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None type(FVW_ParameterType) :: p !< Parameters type(FVW_ContinuousStateType) :: x !< States type(FVW_MiscVarType) :: m !< Initial misc/optimization variables !type(FVW_VTK_Misc) :: mvtk integer :: iW, j, k, nSpan - integer(IntKi) :: ErrStat2 - character(ErrMsgLen) :: ErrMsg2 - character(*), parameter :: RoutineName = 'FVW_Test_CPUTime' + integer(IntKi) :: errStat2 + character(errMsgLen) :: errMsg2 + character(*), parameter :: RoutineName = 'Test_CPUTime' integer(ReKi), parameter :: nR = 20 real(ReKi), parameter :: R = 100 real(ReKi), parameter :: G = 100 real(ReKi), allocatable, dimension(:,:) :: V1 real(ReKi), allocatable, dimension(:,:) :: V2 real(ReKi) :: t1,t2 - ErrStat = ErrID_None - ErrMsg = "" + errStat = ErrID_None + errMsg = "" ! --- Create a helical wake TODO, put me into FVW_* p%nWings = 3 @@ -784,7 +792,7 @@ subroutine FVW_Test_WakeInducedVelocities(ErrStat, ErrMsg) p%PartPerSegment = 1 allocate(p%W(p%nWings)) p%W(:)%nSpan = nSpan - call FVW_InitStates( x, p, ErrStat, ErrMsg ) + call FVW_InitStates( x, p, errStat, errMsg ) do iW=1,size(x%W); do j=1,size(x%W(iW)%r_NW,2); do k=1,size(x%W(iW)%r_NW,3); @@ -802,61 +810,149 @@ subroutine FVW_Test_WakeInducedVelocities(ErrStat, ErrMsg) enddo allocate(m%W(p%nWings)) do iW = 1,p%nWings - call AllocAry( m%W(iW)%Vind_NW , 3 , nSpan+1 ,p%nNWMax+1, 'Vind on NW ', ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName); m%W(iW)%Vind_NW= -999_ReKi; - call AllocAry( m%W(iW)%Vind_FW , 3 , FWnSpan+1,p%nFWMax+1, 'Vind on FW ', ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName); m%W(iW)%Vind_FW= -999_ReKi; + call AllocAry( m%W(iW)%Vind_NW , 3 , nSpan+1 ,p%nNWMax+1, 'Vind on NW ', errStat2, errMsg2); call SeterrStat(errStat2, errMsg2, errStat, errMsg, RoutineName); m%W(iW)%Vind_NW= -999_ReKi; + call AllocAry( m%W(iW)%Vind_FW , 3 , FWnSpan+1,p%nFWMax+1, 'Vind on FW ', errStat2, errMsg2); call SeterrStat(errStat2, errMsg2, errStat, errMsg, RoutineName); m%W(iW)%Vind_FW= -999_ReKi; enddo - call FVW_InitMiscVarsPostParam( p, m, ErrStat2, ErrMsg2) ! Alloc Sgmt, CPs, Uind + call FVW_InitMiscVarsPostParam( p, m, errStat2, errMsg2) ! Alloc Sgmt, CPs, Uind ! --- Compute induced velocity on full wake allocate(V1(3,nSpan+1)) p%nNWFree=p%nNWMax call cpu_time(t1) - call WakeInducedVelocities(p, x, m, ErrStat2, ErrMsg2); + call WakeInducedVelocities(p, x, m, errStat2, errMsg2); call cpu_time(t2) - print*,'Ellapsed time',t2-t1 + !print*,'Ellapsed time',t2-t1 V1 = m%W(1)%Vind_NW(:,:,1) ! --- Compute induced velocity on free wake only allocate(V2(3,nSpan+1)) p%nNWFree=int(p%nNWMax/5) call cpu_time(t1) - call WakeInducedVelocities(p, x, m, ErrStat2, ErrMsg2); + call WakeInducedVelocities(p, x, m, errStat2, errMsg2); call cpu_time(t2) - print*,'Ellapsed time',t2-t1 + !print*,'Ellapsed time',t2-t1 V2 = m%W(1)%Vind_NW(:,:,1) - !print*,'>>>Vx mean ',sum(abs((V1(1,:))))/nSpan - !print*,'>>>Vx mean ',sum(abs((V2(1,:))))/nSpan - !print*,'>>> Vx error',sum(abs(V1(1,:)-V2(1,:)))/nSpan - !print*,'>>> Vy error',sum(abs(V1(2,:)-V2(2,:)))/nSpan - !print*,'>>> Vz error',sum(abs(V1(3,:)-V2(3,:)))/nSpan - !call WrVTK_Segments('_TEST.vtk', mvtk, m%Sgmt%Points(:,:), m%Sgmt%Connct(:,:), m%Sgmt%Gamma(:), m%Sgmt%Epsilon(:), .false.) - call test_almost_equal(RoutineName,'Uind nNW/nNWFree', V1, V2, 1e-6_ReKi, .true.,.true.) + call test_almost_equal('Uind nNW/nNWFree', V1, V2, 1e-6_ReKi, errStat, errMsg) deallocate(V1) deallocate(V2) - call FVW_DestroyParam(p, ErrStat2, ErrMsg2) - call FVW_DestroyContState(x, ErrStat2, ErrMsg2) - call FVW_DestroyMisc(m, ErrStat2, ErrMsg2) + call FVW_DestroyParam(p, errStat2, errMsg2) + call FVW_DestroyContState(x, errStat2, errMsg2) + call FVW_DestroyMisc(m, errStat2, errMsg2) end subroutine FVW_Test_WakeInducedVelocities + + !> Test the resolution of a system A x = b + subroutine Test_LinSolve(errStat, errMsg) + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None + real(ReKi) :: AA(3,3) + real(ReKi) :: B(3) + real(ReKi) :: B2(3) + real(ReKi) :: X(3) + errStat = ErrID_None + errMsg = "" + ! --- Test for non singular system + AA(1,:) = (/ 0., 1., 1./) + AA(2,:) = (/ 1., 2., 0./) + AA(3,:) = (/ 0., 0., 3./) + B(1:3) = (/ 1., 1., 1./) + call linalg_solveWrap(AA, B, X, errStat, errMsg) + B2 = matmul(AA, X) + call test_almost_equal('LinSolve 3x3 ', B, B2, 1e-8, errStat, errMsg) + + end subroutine Test_LinSolve + + !> Test that the pressure coefficient on a sphere (made of source panel) match potential flow theory + subroutine Test_SrcPnl_Sphere(errStat, errMsg) + use VTK, only: vtk_misc_init + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None + ! + type(T_SrcPanlParam) :: p_SrcPnl + type(T_SrcPanlVar) :: z_SrcPnl + type(T_SrcPanlMisc) :: m_SrcPnl + integer(IntKi) :: errStat2 + character(errMsgLen) :: errMsg2 + integer, parameter :: nMain = 17 !< + integer, parameter :: nScnd = 15 !< + !integer, parameter :: nMain =67 !< + !integer, parameter :: nScnd =65 !< + integer :: ncp + integer :: np , icpp, ip + real(ReKi), dimension(:,:), allocatable :: CPs + real(ReKi), dimension(:,:), allocatable :: UI + real(ReKi), dimension(3) :: Uwnd + real(ReKi) :: Uwnd_norm2 + real(ReKi) :: rho, Cp + character(32) :: label + type(VTK_Misc) :: mvtk + errStat = ErrID_None + errMsg = "" + + call EllipsoidPanels(nMain, nScnd, 1., 1., 1., p_SrcPnl%P, p_SrcPnl%IDs) + call srcPnl_init(p_SrcPnl, m_SrcPnl, z_SrcPnl, errStat2, errMsg2); if(Failed()) return + label = 'n1='//trim(num2lstr(nMain))//' - n2='//trim(num2lstr(nScnd)) + p_SrcPnl%Comment = ' - Sphere - '//trim(label) + + ! --- Compute Uext + do icpp = 1, p_SrcPnl%n + m_SrcPnl%Uwnd(:, icpp) = (/1.0, 0., 0./) + enddo + m_SrcPnl%Uext = m_SrcPnl%Uwnd ! + Uind_other + + ! Compute the value of the source panel (sigma) ! assumes that Uext and AI were computed before + call srcPnl_solve(p_SrcPnl, m_SrcPnl, z_SrcPnl, errStat2, errMsg2); if(Failed()) return + + ! Compute velocity, pressure, loads on the source panels + call srcPnl_calcOutput(p_SrcPnl, m_SrcPnl, z_SrcPnl, 1.225_ReKi) !, errStat, errMsg + + call vtk_misc_init(mvtk) + call WrVTK_Panels('../_Sphere_'//trim(label)//'.vtk', mvtk, p_SrcPnl, m_SrcPnl, z_SrcPnl) + + !print*,'>>>> CP Min Max',(/minval(m_SrcPnl%Cp), maxval(m_SrcPnl%Cp)/) + !print*,'>>>> P Min Max',(/minval(m_SrcPnl%p), maxval(m_SrcPnl%p)/) + if (nMain==17) then + call test_almost_equal('Cp sphere', (/minval(m_SrcPnl%Cp), maxval(m_SrcPnl%Cp)/), (/-1.2987,1.0000/), 1e-3_ReKi, errStat2, errMsg2);if(Failed())return + !call test_almost_equal('p sphere' , (/minval(m_SrcPnl%p), maxval(m_SrcPnl%p)/) , (/0.000,1.4079/), 1e-3_ReKi, errStat2, errMsg2);if(Failed())return + call test_almost_equal('p sphere' , (/minval(m_SrcPnl%p), maxval(m_SrcPnl%p)/) , (/-0.7955,0.6125/), 1e-3_ReKi, errStat2, errMsg2);if(Failed())return + else if (nMain==67) then + call test_almost_equal('Cp sphere', (/minval(m_SrcPnl%Cp), maxval(m_SrcPnl%Cp)/), (/-1.26739,1.0000/), 1e-3_ReKi, errStat2, errMsg2);if(Failed())return + !call test_almost_equal('p sphere' , (/minval(m_SrcPnl%p), maxval(m_SrcPnl%p)/) , (/0.000,1.38878/), 1e-3_ReKi, errStat2, errMsg2);if(Failed())return + call test_almost_equal('p sphere' , (/minval(m_SrcPnl%p), maxval(m_SrcPnl%p)/) , (/-0.7955,0.6125/), 1e-3_ReKi, errStat2, errMsg2);if(Failed())return + endif + contains + logical function Failed() + call SeterrStat(errStat2, errMsg2, errStat, errMsg, 'Test_SrcPnl_Sphere') + Failed = errStat >= AbortErrLev + end function + end subroutine Test_SrcPnl_Sphere + !> Main test function - subroutine FVW_RunTests(ErrStat,ErrMsg) - integer(IntKi) , intent(out) :: ErrStat !< Error status of the operation - character(ErrMsgLen), intent(out) :: ErrMsg !< Error message if ErrStat /= ErrID_None - integer(IntKi) :: ErrStat2 - character(ErrMsgLen) :: ErrMsg2 - character(len=255) :: testname - ! Initialize ErrStat - ErrStat = ErrID_None - ErrMsg = "" - testname='FVW' - call Test_BiotSavart_Sgmt(testname, ErrStat2, ErrMsg2) - call Test_BiotSavart_Part(testname, ErrStat2, ErrMsg2) - call Test_BiotSavart_PartTree(testname, ErrStat2, ErrMsg2) - call Test_SegmentsToPart(testname, ErrStat2, ErrMsg2) - call FVW_Test_WakeInducedVelocities(ErrStat2, ErrMsg2) + !! NOTE: Potentially edit ../tests/test_FVW_testsuite.f90 as well + subroutine FVW_RunTests(errStat,errMsg) + integer(IntKi) , intent(out) :: errStat !< Error status of the operation + character(errMsgLen), intent(out) :: errMsg !< Error message if errStat /= ErrID_None + integer(IntKi) :: errStat2 + character(errMsgLen) :: errMsg2 + ! Initialize errStat + errStat = ErrID_None + errMsg = "" + call Test_LinSolve (errStat2,errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, 'FVW_RunTests') + call Test_SrcPnl_Sphere (errStat2,errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, 'FVW_RunTests') + call Test_BiotSavart_SrcPnl (errStat2,errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, 'FVW_RunTests') + call Test_BiotSavart_Sgmt (errStat2,errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, 'FVW_RunTests') + call Test_BiotSavart_Part (errStat2,errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, 'FVW_RunTests') + call Test_BiotSavart_PartTree (errStat2,errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, 'FVW_RunTests') + call Test_SegmentsToPart (errStat2,errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, 'FVW_RunTests') + call FVW_Test_WakeInducedVelocities(errStat2,errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, 'FVW_RunTests') + contains + logical function Failed() + call SetErrStat(errStat2, errMsg2, errStat, errMsg, 'FVW_RunTests') + Failed = errStat >= AbortErrLev + end function end subroutine FVW_RunTests end module FVW_Tests diff --git a/modules/aerodyn/src/FVW_Types.f90 b/modules/aerodyn/src/FVW_Types.f90 index 2ef8bc6c5c..53ff476c2f 100644 --- a/modules/aerodyn/src/FVW_Types.f90 +++ b/modules/aerodyn/src/FVW_Types.f90 @@ -50,9 +50,9 @@ MODULE FVW_Types REAL(ReKi) :: xEnd = 0.0_ReKi !< xEnd [-] REAL(ReKi) :: yEnd = 0.0_ReKi !< yEnd [-] REAL(ReKi) :: zEnd = 0.0_ReKi !< zEnd [-] - INTEGER(IntKi) :: nx = 0_IntKi !< nx [-] - INTEGER(IntKi) :: ny = 0_IntKi !< ny [-] - INTEGER(IntKi) :: nz = 0_IntKi !< nz [-] + INTEGER(IntKi) :: nx = 0 !< nx [-] + INTEGER(IntKi) :: ny = 0 !< ny [-] + INTEGER(IntKi) :: nz = 0 !< nz [-] 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 [-] @@ -65,8 +65,8 @@ MODULE FVW_Types REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: Gamma !< Segment circulations [-] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: Epsilon !< Segment regularization parameter [-] INTEGER(IntKi) :: RegFunction = 0_IntKi !< Type of regularizaion function (LambOseen, Vatistas, see FVW_BiotSavart) [-] - INTEGER(IntKi) :: nAct = 0_IntKi !< Number of active segments [-] - INTEGER(IntKi) :: nActP = 0_IntKi !< Number of active segment points [-] + INTEGER(IntKi) :: nAct = 0 !< Number of active segments [-] + INTEGER(IntKi) :: nActP = 0 !< Number of active segment points [-] END TYPE T_Sgmt ! ======================= ! ========= T_Part ======= @@ -75,32 +75,69 @@ MODULE FVW_Types REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Alpha !< Particle intensity 3 x nP [-] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: RegParam !< Particle regularization parameter [-] INTEGER(IntKi) :: RegFunction = 0_IntKi !< Type of regularizaion function (FVW_BiotSavart) [-] - INTEGER(IntKi) :: nAct = 0_IntKi !< Number of active particles <=nP [-] + INTEGER(IntKi) :: nAct = 0 !< Number of active particles <=nP [-] END TYPE T_Part ! ======================= +! ========= T_SrcPanlParam ======= + TYPE, PUBLIC :: T_SrcPanlParam + INTEGER(IntKi) :: n = 0 !< Number of panels [-] + CHARACTER(1024) :: Comment !< Area [-] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: Area !< Area [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: P !< Points database [-] + INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: IDs !< Points connectivity ID for each panel 4xnPn [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: BodyIDs !< Body ID for each panel [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Pcent !< Centroid point [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Pmid !< Mid point [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Normal !< Normal to panel [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: xi !< Xi coordinate for each panel point [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: eta !< eta coordinate for each panel point [-] + REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: R_g2p !< Transformation matrix global to Panel [-] + END TYPE T_SrcPanlParam +! ======================= +! ========= T_SrcPanlVar ======= + TYPE, PUBLIC :: T_SrcPanlVar + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: Sigma !< Intensities [-] + END TYPE T_SrcPanlVar +! ======================= +! ========= T_SrcPanlMisc ======= + TYPE, PUBLIC :: T_SrcPanlMisc + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: AI !< Influence matrix to solve for intensity [-] + REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: UUI !< Unit induced velocity [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Uext !< External inflow velocity at panel control points [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Uwnd !< Inflow velocity vector at panel control points [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Uind !< Induced velocity vector at panel control points [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Utot !< Full velocity (include induced) vector at panel control points [-] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: Cp !< Pressure coefficient [-] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: p !< Pressure (1/2 rho Utot^2) [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: F !< Force from pressure [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: FpA !< Force per area [-] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: RHS !< Right hand side [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: IPIV !< Pivot for factorization [-] + END TYPE T_SrcPanlMisc +! ======================= ! ========= Wng_ParameterType ======= TYPE, PUBLIC :: Wng_ParameterType REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: chord_LL !< Chord of each blade element from input file [idx1=BladeNode, idx2=Blade number] [-] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: chord_CP !< Chord on LL cp [m] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: s_LL !< Spanwise coordinate of LL elements [m] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: s_CP !< Spanwise coordinate of LL CP [m] - INTEGER(IntKi) :: iRotor = 0_IntKi !< Index of rotor the wing belong to [-] + INTEGER(IntKi) :: iRotor = 0 !< Index of rotor the wing belong to [-] INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: AFindx !< Index to the airfoils from AD15 [BladeNode,BladeIndex=1] [-] - INTEGER(IntKi) :: nSpan = 0_IntKi !< TODO, should be defined per wing. Number of spanwise element [-] + INTEGER(IntKi) :: nSpan = 0 !< TODO, should be defined per wing. Number of spanwise element [-] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: PrescribedCirculation !< Prescribed circulation on all lifting lines [m/s] END TYPE Wng_ParameterType ! ======================= ! ========= FVW_ParameterType ======= TYPE, PUBLIC :: FVW_ParameterType - INTEGER(IntKi) :: nRotors = 0_IntKi !< Number of Wings [-] - INTEGER(IntKi) :: nWings = 0_IntKi !< Number of Wings [-] + INTEGER(IntKi) :: nRotors = 0 !< Number of Rotors [-] + INTEGER(IntKi) :: nWings = 0 !< Number of Wings [-] TYPE(Wng_ParameterType) , DIMENSION(:), ALLOCATABLE :: W !< Wings parameters [-] INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: Bld2Wings !< Index mapping from blades to wings [-] - INTEGER(IntKi) :: iNWStart = 0_IntKi !< Index where NW start in r_NW. (iNWStart=2, the first panel contains the lifting line panel, otherwise, start at 1) [-] - INTEGER(IntKi) :: nNWMax = 0_IntKi !< Maximum number of nw panels, per wing [-] - INTEGER(IntKi) :: nNWFree = 0_IntKi !< Number of nw panels that are free, per wing [-] - INTEGER(IntKi) :: nFWMax = 0_IntKi !< Maximum number of fw panels, per wing [-] - INTEGER(IntKi) :: nFWFree = 0_IntKi !< Number of fw panels that are free, per wing [-] + INTEGER(IntKi) :: iNWStart = 0 !< Index where NW start in r_NW. (iNWStart=2, the first panel contains the lifting line panel, otherwise, start at 1) [-] + INTEGER(IntKi) :: nNWMax = 0 !< Maximum number of nw panels, per wing [-] + INTEGER(IntKi) :: nNWFree = 0 !< Number of nw panels that are free, per wing [-] + INTEGER(IntKi) :: nFWMax = 0 !< Maximum number of fw panels, per wing [-] + INTEGER(IntKi) :: nFWFree = 0 !< Number of fw panels that are free, per wing [-] LOGICAL :: FWShedVorticity = .false. !< Include shed vorticity in the far wake [-] INTEGER(IntKi) :: IntMethod = 0_IntKi !< Integration Method (1=RK4, 2=AB4, 3=ABM4, 5=Euler1) [-] REAL(ReKi) :: FreeWakeStart = 0.0_ReKi !< Time when wake starts convecting (rolling up) [s] @@ -124,6 +161,7 @@ MODULE FVW_Types INTEGER(IntKi) , DIMENSION(1:2) :: PartPerSegment = 0_IntKi !< Number of particles per segment, e.g. for tree method, for full wake and lifting line [-] REAL(DbKi) :: DTaero = 0.0_R8Ki !< Time interval for calls calculations [s] REAL(DbKi) :: DTfvw = 0.0_R8Ki !< Time interval for calculating wake induced velocities [s] + REAL(ReKi) :: AirDens = 0.0_ReKi !< Air density [kg/m^3] REAL(ReKi) :: KinVisc = 0.0_ReKi !< Kinematic air viscosity [m^2/s] INTEGER(IntKi) :: MHK = 0_IntKi !< MHK flag [-] REAL(ReKi) :: WtrDpth = 0.0_ReKi !< Water depth [m] @@ -134,13 +172,17 @@ MODULE FVW_Types CHARACTER(1024) :: RootName !< RootName for writing output files [-] CHARACTER(1024) :: VTK_OutFileRoot !< Rootdirectory for writing VTK files [-] CHARACTER(1024) :: VTK_OutFileBase !< Basename for writing VTK files [-] - INTEGER(IntKi) :: nGridOut = 0_IntKi !< Number of VTK grid to output [-] + INTEGER(IntKi) :: nGridOut = 0 !< Number of VTK grid to output [-] LOGICAL :: InductionAtCP = .true. !< Compute induced velocities at nodes or CP [-] LOGICAL :: WakeAtTE = .true. !< Start the wake at the trailing edge, or at the LL [-] LOGICAL :: DStallOnWake = .false. !< Dynamic stall has influence on wake [-] LOGICAL :: Induction = .true. !< Compute induction [-] REAL(ReKi) :: kFrozenNWStart = 0.75 !< Fraction of wake induced velocity at start of frozen wake. 1 seems too strong. [-] REAL(ReKi) :: kFrozenNWEnd = 0.5 !< Fraction of wake induced velocity at end of frozen wake [-] + REAL(ReKi) :: zGround = 0.0 !< Ground height [-] + REAL(ReKi) :: zGroundPush = 0.1 !< Distance above ground where vortices are pushed back [-] + INTEGER(IntKi) :: nSrcPnlUpdate = 1 !< How often do src panel updates (in time steps of OLAF) [-] + TYPE(T_SrcPanlParam) :: SrcPnl !< Source panel parameters [-] END TYPE FVW_ParameterType ! ======================= ! ========= Wng_ContinuousStateType ======= @@ -237,6 +279,7 @@ MODULE FVW_Types LOGICAL :: UA_Flag = .false. !< logical flag indicating whether to use UnsteadyAero [-] TYPE(T_Sgmt) :: Sgmt !< Segments storage [-] TYPE(T_Part) :: Part !< Particle storage [-] + TYPE(T_SrcPanlMisc) :: SrcPnl !< Source panels storage [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: CPs !< Control points used for wake rollup computation [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Uind !< Induced velocities obtained at control points [-] TYPE(GridOutType) , DIMENSION(:), ALLOCATABLE :: GridOutputs !< Number of VTK grid to output [-] @@ -278,12 +321,14 @@ MODULE FVW_Types TYPE, PUBLIC :: FVW_ConstraintStateType TYPE(Wng_ConstraintStateType) , DIMENSION(:), ALLOCATABLE :: W !< rotors constr. states [-] REAL(ReKi) :: residual = 0.0_ReKi !< Residual [-] + TYPE(T_SrcPanlVar) :: SrcPnl !< Source panel constraints [-] END TYPE FVW_ConstraintStateType ! ======================= ! ========= FVW_OtherStateType ======= TYPE, PUBLIC :: FVW_OtherStateType - INTEGER(IntKi) :: Dummy = 0_IntKi !< Empty to satisfy framework [-] + REAL(ReKi) :: ShedScale = 0.0_ReKi !< Scale shed vorticity at begining of simulation [-] TYPE(UA_OtherStateType) , DIMENSION(:), ALLOCATABLE :: UA !< other states for UnsteadyAero for each wing [-] + LOGICAL :: Initialized = .false. !< True if OLAF is initialized [-] END TYPE FVW_OtherStateType ! ======================= ! ========= Wng_InitInputType ======= @@ -304,6 +349,7 @@ MODULE FVW_Types TYPE(MeshType) , DIMENSION(:), ALLOCATABLE :: WingsMesh !< Input Mesh defining position and orientation of wings (nSpan+1) [-] INTEGER(IntKi) :: numBladeNodes = 0_IntKi !< Number of nodes on each blade [-] REAL(DbKi) :: DTaero = 0.0_R8Ki !< Time interval for calls (from AD15) [s] + REAL(ReKi) :: AirDens = 0.0_ReKi !< Air density [kg/m^3] REAL(ReKi) :: KinVisc = 0.0_ReKi !< Kinematic air viscosity [m^2/s] INTEGER(IntKi) :: MHK = 0_IntKi !< MHK flag [-] REAL(ReKi) :: WtrDpth = 0.0_ReKi !< Water depth [m] @@ -345,6 +391,7 @@ MODULE FVW_Types INTEGER(IntKi) :: VTKBlades = 0_IntKi !< Outputs VTk for each blade 0=no blade, 1=Bld 1 [-] REAL(DbKi) :: DTvtk = 0.0_R8Ki !< Requested timestep between VTK outputs (calculated from the VTK_fps read in) [s] INTEGER(IntKi) :: VTKCoord = 0_IntKi !< Switch for VTK outputs coordinate system [-] + CHARACTER(1024) :: SrcPnlFile !< Input file for source panels [-] END TYPE FVW_InputFile ! ======================= ! ========= FVW_InitOutputType ======= @@ -702,6 +749,523 @@ subroutine FVW_UnPackT_Part(RF, OutData) call RegUnpack(RF, OutData%nAct); if (RegCheckErr(RF, RoutineName)) return end subroutine +subroutine FVW_CopyT_SrcPanlParam(SrcT_SrcPanlParamData, DstT_SrcPanlParamData, CtrlCode, ErrStat, ErrMsg) + type(T_SrcPanlParam), intent(in) :: SrcT_SrcPanlParamData + type(T_SrcPanlParam), intent(inout) :: DstT_SrcPanlParamData + integer(IntKi), intent(in ) :: CtrlCode + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + integer(B4Ki) :: LB(3), UB(3) + integer(IntKi) :: ErrStat2 + character(*), parameter :: RoutineName = 'FVW_CopyT_SrcPanlParam' + ErrStat = ErrID_None + ErrMsg = '' + DstT_SrcPanlParamData%n = SrcT_SrcPanlParamData%n + DstT_SrcPanlParamData%Comment = SrcT_SrcPanlParamData%Comment + if (allocated(SrcT_SrcPanlParamData%Area)) then + LB(1:1) = lbound(SrcT_SrcPanlParamData%Area) + UB(1:1) = ubound(SrcT_SrcPanlParamData%Area) + if (.not. allocated(DstT_SrcPanlParamData%Area)) then + allocate(DstT_SrcPanlParamData%Area(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlParamData%Area.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlParamData%Area = SrcT_SrcPanlParamData%Area + end if + if (allocated(SrcT_SrcPanlParamData%P)) then + LB(1:2) = lbound(SrcT_SrcPanlParamData%P) + UB(1:2) = ubound(SrcT_SrcPanlParamData%P) + if (.not. allocated(DstT_SrcPanlParamData%P)) then + allocate(DstT_SrcPanlParamData%P(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlParamData%P.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlParamData%P = SrcT_SrcPanlParamData%P + end if + if (allocated(SrcT_SrcPanlParamData%IDs)) then + LB(1:2) = lbound(SrcT_SrcPanlParamData%IDs) + UB(1:2) = ubound(SrcT_SrcPanlParamData%IDs) + if (.not. allocated(DstT_SrcPanlParamData%IDs)) then + allocate(DstT_SrcPanlParamData%IDs(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlParamData%IDs.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlParamData%IDs = SrcT_SrcPanlParamData%IDs + end if + if (allocated(SrcT_SrcPanlParamData%BodyIDs)) then + LB(1:1) = lbound(SrcT_SrcPanlParamData%BodyIDs) + UB(1:1) = ubound(SrcT_SrcPanlParamData%BodyIDs) + if (.not. allocated(DstT_SrcPanlParamData%BodyIDs)) then + allocate(DstT_SrcPanlParamData%BodyIDs(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlParamData%BodyIDs.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlParamData%BodyIDs = SrcT_SrcPanlParamData%BodyIDs + end if + if (allocated(SrcT_SrcPanlParamData%Pcent)) then + LB(1:2) = lbound(SrcT_SrcPanlParamData%Pcent) + UB(1:2) = ubound(SrcT_SrcPanlParamData%Pcent) + if (.not. allocated(DstT_SrcPanlParamData%Pcent)) then + allocate(DstT_SrcPanlParamData%Pcent(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlParamData%Pcent.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlParamData%Pcent = SrcT_SrcPanlParamData%Pcent + end if + if (allocated(SrcT_SrcPanlParamData%Pmid)) then + LB(1:2) = lbound(SrcT_SrcPanlParamData%Pmid) + UB(1:2) = ubound(SrcT_SrcPanlParamData%Pmid) + if (.not. allocated(DstT_SrcPanlParamData%Pmid)) then + allocate(DstT_SrcPanlParamData%Pmid(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlParamData%Pmid.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlParamData%Pmid = SrcT_SrcPanlParamData%Pmid + end if + if (allocated(SrcT_SrcPanlParamData%Normal)) then + LB(1:2) = lbound(SrcT_SrcPanlParamData%Normal) + UB(1:2) = ubound(SrcT_SrcPanlParamData%Normal) + if (.not. allocated(DstT_SrcPanlParamData%Normal)) then + allocate(DstT_SrcPanlParamData%Normal(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlParamData%Normal.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlParamData%Normal = SrcT_SrcPanlParamData%Normal + end if + if (allocated(SrcT_SrcPanlParamData%xi)) then + LB(1:2) = lbound(SrcT_SrcPanlParamData%xi) + UB(1:2) = ubound(SrcT_SrcPanlParamData%xi) + if (.not. allocated(DstT_SrcPanlParamData%xi)) then + allocate(DstT_SrcPanlParamData%xi(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlParamData%xi.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlParamData%xi = SrcT_SrcPanlParamData%xi + end if + if (allocated(SrcT_SrcPanlParamData%eta)) then + LB(1:2) = lbound(SrcT_SrcPanlParamData%eta) + UB(1:2) = ubound(SrcT_SrcPanlParamData%eta) + if (.not. allocated(DstT_SrcPanlParamData%eta)) then + allocate(DstT_SrcPanlParamData%eta(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlParamData%eta.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlParamData%eta = SrcT_SrcPanlParamData%eta + end if + if (allocated(SrcT_SrcPanlParamData%R_g2p)) then + LB(1:3) = lbound(SrcT_SrcPanlParamData%R_g2p) + UB(1:3) = ubound(SrcT_SrcPanlParamData%R_g2p) + if (.not. allocated(DstT_SrcPanlParamData%R_g2p)) then + allocate(DstT_SrcPanlParamData%R_g2p(LB(1):UB(1),LB(2):UB(2),LB(3):UB(3)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlParamData%R_g2p.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlParamData%R_g2p = SrcT_SrcPanlParamData%R_g2p + end if +end subroutine + +subroutine FVW_DestroyT_SrcPanlParam(T_SrcPanlParamData, ErrStat, ErrMsg) + type(T_SrcPanlParam), intent(inout) :: T_SrcPanlParamData + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + character(*), parameter :: RoutineName = 'FVW_DestroyT_SrcPanlParam' + ErrStat = ErrID_None + ErrMsg = '' + if (allocated(T_SrcPanlParamData%Area)) then + deallocate(T_SrcPanlParamData%Area) + end if + if (allocated(T_SrcPanlParamData%P)) then + deallocate(T_SrcPanlParamData%P) + end if + if (allocated(T_SrcPanlParamData%IDs)) then + deallocate(T_SrcPanlParamData%IDs) + end if + if (allocated(T_SrcPanlParamData%BodyIDs)) then + deallocate(T_SrcPanlParamData%BodyIDs) + end if + if (allocated(T_SrcPanlParamData%Pcent)) then + deallocate(T_SrcPanlParamData%Pcent) + end if + if (allocated(T_SrcPanlParamData%Pmid)) then + deallocate(T_SrcPanlParamData%Pmid) + end if + if (allocated(T_SrcPanlParamData%Normal)) then + deallocate(T_SrcPanlParamData%Normal) + end if + if (allocated(T_SrcPanlParamData%xi)) then + deallocate(T_SrcPanlParamData%xi) + end if + if (allocated(T_SrcPanlParamData%eta)) then + deallocate(T_SrcPanlParamData%eta) + end if + if (allocated(T_SrcPanlParamData%R_g2p)) then + deallocate(T_SrcPanlParamData%R_g2p) + end if +end subroutine + +subroutine FVW_PackT_SrcPanlParam(RF, Indata) + type(RegFile), intent(inout) :: RF + type(T_SrcPanlParam), intent(in) :: InData + character(*), parameter :: RoutineName = 'FVW_PackT_SrcPanlParam' + if (RF%ErrStat >= AbortErrLev) return + call RegPack(RF, InData%n) + call RegPack(RF, InData%Comment) + call RegPackAlloc(RF, InData%Area) + call RegPackAlloc(RF, InData%P) + call RegPackAlloc(RF, InData%IDs) + call RegPackAlloc(RF, InData%BodyIDs) + call RegPackAlloc(RF, InData%Pcent) + call RegPackAlloc(RF, InData%Pmid) + call RegPackAlloc(RF, InData%Normal) + call RegPackAlloc(RF, InData%xi) + call RegPackAlloc(RF, InData%eta) + call RegPackAlloc(RF, InData%R_g2p) + if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine FVW_UnPackT_SrcPanlParam(RF, OutData) + type(RegFile), intent(inout) :: RF + type(T_SrcPanlParam), intent(inout) :: OutData + character(*), parameter :: RoutineName = 'FVW_UnPackT_SrcPanlParam' + integer(B4Ki) :: LB(3), UB(3) + integer(IntKi) :: stat + logical :: IsAllocAssoc + if (RF%ErrStat /= ErrID_None) return + call RegUnpack(RF, OutData%n); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Comment); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%Area); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%P); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%IDs); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%BodyIDs); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%Pcent); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%Pmid); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%Normal); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%xi); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%eta); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%R_g2p); if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine FVW_CopyT_SrcPanlVar(SrcT_SrcPanlVarData, DstT_SrcPanlVarData, CtrlCode, ErrStat, ErrMsg) + type(T_SrcPanlVar), intent(in) :: SrcT_SrcPanlVarData + type(T_SrcPanlVar), intent(inout) :: DstT_SrcPanlVarData + integer(IntKi), intent(in ) :: CtrlCode + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + integer(B4Ki) :: LB(1), UB(1) + integer(IntKi) :: ErrStat2 + character(*), parameter :: RoutineName = 'FVW_CopyT_SrcPanlVar' + ErrStat = ErrID_None + ErrMsg = '' + if (allocated(SrcT_SrcPanlVarData%Sigma)) then + LB(1:1) = lbound(SrcT_SrcPanlVarData%Sigma) + UB(1:1) = ubound(SrcT_SrcPanlVarData%Sigma) + if (.not. allocated(DstT_SrcPanlVarData%Sigma)) then + allocate(DstT_SrcPanlVarData%Sigma(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlVarData%Sigma.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlVarData%Sigma = SrcT_SrcPanlVarData%Sigma + end if +end subroutine + +subroutine FVW_DestroyT_SrcPanlVar(T_SrcPanlVarData, ErrStat, ErrMsg) + type(T_SrcPanlVar), intent(inout) :: T_SrcPanlVarData + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + character(*), parameter :: RoutineName = 'FVW_DestroyT_SrcPanlVar' + ErrStat = ErrID_None + ErrMsg = '' + if (allocated(T_SrcPanlVarData%Sigma)) then + deallocate(T_SrcPanlVarData%Sigma) + end if +end subroutine + +subroutine FVW_PackT_SrcPanlVar(RF, Indata) + type(RegFile), intent(inout) :: RF + type(T_SrcPanlVar), intent(in) :: InData + character(*), parameter :: RoutineName = 'FVW_PackT_SrcPanlVar' + if (RF%ErrStat >= AbortErrLev) return + call RegPackAlloc(RF, InData%Sigma) + if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine FVW_UnPackT_SrcPanlVar(RF, OutData) + type(RegFile), intent(inout) :: RF + type(T_SrcPanlVar), intent(inout) :: OutData + character(*), parameter :: RoutineName = 'FVW_UnPackT_SrcPanlVar' + integer(B4Ki) :: LB(1), UB(1) + integer(IntKi) :: stat + logical :: IsAllocAssoc + if (RF%ErrStat /= ErrID_None) return + call RegUnpackAlloc(RF, OutData%Sigma); if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine FVW_CopyT_SrcPanlMisc(SrcT_SrcPanlMiscData, DstT_SrcPanlMiscData, CtrlCode, ErrStat, ErrMsg) + type(T_SrcPanlMisc), intent(in) :: SrcT_SrcPanlMiscData + type(T_SrcPanlMisc), intent(inout) :: DstT_SrcPanlMiscData + integer(IntKi), intent(in ) :: CtrlCode + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + integer(B4Ki) :: LB(3), UB(3) + integer(IntKi) :: ErrStat2 + character(*), parameter :: RoutineName = 'FVW_CopyT_SrcPanlMisc' + ErrStat = ErrID_None + ErrMsg = '' + if (allocated(SrcT_SrcPanlMiscData%AI)) then + LB(1:2) = lbound(SrcT_SrcPanlMiscData%AI) + UB(1:2) = ubound(SrcT_SrcPanlMiscData%AI) + if (.not. allocated(DstT_SrcPanlMiscData%AI)) then + allocate(DstT_SrcPanlMiscData%AI(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlMiscData%AI.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlMiscData%AI = SrcT_SrcPanlMiscData%AI + end if + if (allocated(SrcT_SrcPanlMiscData%UUI)) then + LB(1:3) = lbound(SrcT_SrcPanlMiscData%UUI) + UB(1:3) = ubound(SrcT_SrcPanlMiscData%UUI) + if (.not. allocated(DstT_SrcPanlMiscData%UUI)) then + allocate(DstT_SrcPanlMiscData%UUI(LB(1):UB(1),LB(2):UB(2),LB(3):UB(3)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlMiscData%UUI.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlMiscData%UUI = SrcT_SrcPanlMiscData%UUI + end if + if (allocated(SrcT_SrcPanlMiscData%Uext)) then + LB(1:2) = lbound(SrcT_SrcPanlMiscData%Uext) + UB(1:2) = ubound(SrcT_SrcPanlMiscData%Uext) + if (.not. allocated(DstT_SrcPanlMiscData%Uext)) then + allocate(DstT_SrcPanlMiscData%Uext(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlMiscData%Uext.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlMiscData%Uext = SrcT_SrcPanlMiscData%Uext + end if + if (allocated(SrcT_SrcPanlMiscData%Uwnd)) then + LB(1:2) = lbound(SrcT_SrcPanlMiscData%Uwnd) + UB(1:2) = ubound(SrcT_SrcPanlMiscData%Uwnd) + if (.not. allocated(DstT_SrcPanlMiscData%Uwnd)) then + allocate(DstT_SrcPanlMiscData%Uwnd(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlMiscData%Uwnd.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlMiscData%Uwnd = SrcT_SrcPanlMiscData%Uwnd + end if + if (allocated(SrcT_SrcPanlMiscData%Uind)) then + LB(1:2) = lbound(SrcT_SrcPanlMiscData%Uind) + UB(1:2) = ubound(SrcT_SrcPanlMiscData%Uind) + if (.not. allocated(DstT_SrcPanlMiscData%Uind)) then + allocate(DstT_SrcPanlMiscData%Uind(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlMiscData%Uind.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlMiscData%Uind = SrcT_SrcPanlMiscData%Uind + end if + if (allocated(SrcT_SrcPanlMiscData%Utot)) then + LB(1:2) = lbound(SrcT_SrcPanlMiscData%Utot) + UB(1:2) = ubound(SrcT_SrcPanlMiscData%Utot) + if (.not. allocated(DstT_SrcPanlMiscData%Utot)) then + allocate(DstT_SrcPanlMiscData%Utot(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlMiscData%Utot.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlMiscData%Utot = SrcT_SrcPanlMiscData%Utot + end if + if (allocated(SrcT_SrcPanlMiscData%Cp)) then + LB(1:1) = lbound(SrcT_SrcPanlMiscData%Cp) + UB(1:1) = ubound(SrcT_SrcPanlMiscData%Cp) + if (.not. allocated(DstT_SrcPanlMiscData%Cp)) then + allocate(DstT_SrcPanlMiscData%Cp(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlMiscData%Cp.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlMiscData%Cp = SrcT_SrcPanlMiscData%Cp + end if + if (allocated(SrcT_SrcPanlMiscData%p)) then + LB(1:1) = lbound(SrcT_SrcPanlMiscData%p) + UB(1:1) = ubound(SrcT_SrcPanlMiscData%p) + if (.not. allocated(DstT_SrcPanlMiscData%p)) then + allocate(DstT_SrcPanlMiscData%p(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlMiscData%p.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlMiscData%p = SrcT_SrcPanlMiscData%p + end if + if (allocated(SrcT_SrcPanlMiscData%F)) then + LB(1:2) = lbound(SrcT_SrcPanlMiscData%F) + UB(1:2) = ubound(SrcT_SrcPanlMiscData%F) + if (.not. allocated(DstT_SrcPanlMiscData%F)) then + allocate(DstT_SrcPanlMiscData%F(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlMiscData%F.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlMiscData%F = SrcT_SrcPanlMiscData%F + end if + if (allocated(SrcT_SrcPanlMiscData%FpA)) then + LB(1:2) = lbound(SrcT_SrcPanlMiscData%FpA) + UB(1:2) = ubound(SrcT_SrcPanlMiscData%FpA) + if (.not. allocated(DstT_SrcPanlMiscData%FpA)) then + allocate(DstT_SrcPanlMiscData%FpA(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlMiscData%FpA.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlMiscData%FpA = SrcT_SrcPanlMiscData%FpA + end if + if (allocated(SrcT_SrcPanlMiscData%RHS)) then + LB(1:1) = lbound(SrcT_SrcPanlMiscData%RHS) + UB(1:1) = ubound(SrcT_SrcPanlMiscData%RHS) + if (.not. allocated(DstT_SrcPanlMiscData%RHS)) then + allocate(DstT_SrcPanlMiscData%RHS(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlMiscData%RHS.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlMiscData%RHS = SrcT_SrcPanlMiscData%RHS + end if + if (allocated(SrcT_SrcPanlMiscData%IPIV)) then + LB(1:1) = lbound(SrcT_SrcPanlMiscData%IPIV) + UB(1:1) = ubound(SrcT_SrcPanlMiscData%IPIV) + if (.not. allocated(DstT_SrcPanlMiscData%IPIV)) then + allocate(DstT_SrcPanlMiscData%IPIV(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstT_SrcPanlMiscData%IPIV.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstT_SrcPanlMiscData%IPIV = SrcT_SrcPanlMiscData%IPIV + end if +end subroutine + +subroutine FVW_DestroyT_SrcPanlMisc(T_SrcPanlMiscData, ErrStat, ErrMsg) + type(T_SrcPanlMisc), intent(inout) :: T_SrcPanlMiscData + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + character(*), parameter :: RoutineName = 'FVW_DestroyT_SrcPanlMisc' + ErrStat = ErrID_None + ErrMsg = '' + if (allocated(T_SrcPanlMiscData%AI)) then + deallocate(T_SrcPanlMiscData%AI) + end if + if (allocated(T_SrcPanlMiscData%UUI)) then + deallocate(T_SrcPanlMiscData%UUI) + end if + if (allocated(T_SrcPanlMiscData%Uext)) then + deallocate(T_SrcPanlMiscData%Uext) + end if + if (allocated(T_SrcPanlMiscData%Uwnd)) then + deallocate(T_SrcPanlMiscData%Uwnd) + end if + if (allocated(T_SrcPanlMiscData%Uind)) then + deallocate(T_SrcPanlMiscData%Uind) + end if + if (allocated(T_SrcPanlMiscData%Utot)) then + deallocate(T_SrcPanlMiscData%Utot) + end if + if (allocated(T_SrcPanlMiscData%Cp)) then + deallocate(T_SrcPanlMiscData%Cp) + end if + if (allocated(T_SrcPanlMiscData%p)) then + deallocate(T_SrcPanlMiscData%p) + end if + if (allocated(T_SrcPanlMiscData%F)) then + deallocate(T_SrcPanlMiscData%F) + end if + if (allocated(T_SrcPanlMiscData%FpA)) then + deallocate(T_SrcPanlMiscData%FpA) + end if + if (allocated(T_SrcPanlMiscData%RHS)) then + deallocate(T_SrcPanlMiscData%RHS) + end if + if (allocated(T_SrcPanlMiscData%IPIV)) then + deallocate(T_SrcPanlMiscData%IPIV) + end if +end subroutine + +subroutine FVW_PackT_SrcPanlMisc(RF, Indata) + type(RegFile), intent(inout) :: RF + type(T_SrcPanlMisc), intent(in) :: InData + character(*), parameter :: RoutineName = 'FVW_PackT_SrcPanlMisc' + if (RF%ErrStat >= AbortErrLev) return + call RegPackAlloc(RF, InData%AI) + call RegPackAlloc(RF, InData%UUI) + call RegPackAlloc(RF, InData%Uext) + call RegPackAlloc(RF, InData%Uwnd) + call RegPackAlloc(RF, InData%Uind) + call RegPackAlloc(RF, InData%Utot) + call RegPackAlloc(RF, InData%Cp) + call RegPackAlloc(RF, InData%p) + call RegPackAlloc(RF, InData%F) + call RegPackAlloc(RF, InData%FpA) + call RegPackAlloc(RF, InData%RHS) + call RegPackAlloc(RF, InData%IPIV) + if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine FVW_UnPackT_SrcPanlMisc(RF, OutData) + type(RegFile), intent(inout) :: RF + type(T_SrcPanlMisc), intent(inout) :: OutData + character(*), parameter :: RoutineName = 'FVW_UnPackT_SrcPanlMisc' + integer(B4Ki) :: LB(3), UB(3) + integer(IntKi) :: stat + logical :: IsAllocAssoc + if (RF%ErrStat /= ErrID_None) return + call RegUnpackAlloc(RF, OutData%AI); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%UUI); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%Uext); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%Uwnd); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%Uind); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%Utot); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%Cp); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%p); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%F); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%FpA); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%RHS); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%IPIV); if (RegCheckErr(RF, RoutineName)) return +end subroutine + subroutine FVW_CopyWng_ParameterType(SrcWng_ParameterTypeData, DstWng_ParameterTypeData, CtrlCode, ErrStat, ErrMsg) type(Wng_ParameterType), intent(in) :: SrcWng_ParameterTypeData type(Wng_ParameterType), intent(inout) :: DstWng_ParameterTypeData @@ -921,6 +1485,7 @@ subroutine FVW_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%PartPerSegment = SrcParamData%PartPerSegment DstParamData%DTaero = SrcParamData%DTaero DstParamData%DTfvw = SrcParamData%DTfvw + DstParamData%AirDens = SrcParamData%AirDens DstParamData%KinVisc = SrcParamData%KinVisc DstParamData%MHK = SrcParamData%MHK DstParamData%WtrDpth = SrcParamData%WtrDpth @@ -938,6 +1503,12 @@ subroutine FVW_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%Induction = SrcParamData%Induction DstParamData%kFrozenNWStart = SrcParamData%kFrozenNWStart DstParamData%kFrozenNWEnd = SrcParamData%kFrozenNWEnd + DstParamData%zGround = SrcParamData%zGround + DstParamData%zGroundPush = SrcParamData%zGroundPush + DstParamData%nSrcPnlUpdate = SrcParamData%nSrcPnlUpdate + call FVW_CopyT_SrcPanlParam(SrcParamData%SrcPnl, DstParamData%SrcPnl, CtrlCode, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return end subroutine subroutine FVW_DestroyParam(ParamData, ErrStat, ErrMsg) @@ -963,6 +1534,8 @@ subroutine FVW_DestroyParam(ParamData, ErrStat, ErrMsg) if (allocated(ParamData%Bld2Wings)) then deallocate(ParamData%Bld2Wings) end if + call FVW_DestroyT_SrcPanlParam(ParamData%SrcPnl, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end subroutine subroutine FVW_PackParam(RF, Indata) @@ -1012,6 +1585,7 @@ subroutine FVW_PackParam(RF, Indata) call RegPack(RF, InData%PartPerSegment) call RegPack(RF, InData%DTaero) call RegPack(RF, InData%DTfvw) + call RegPack(RF, InData%AirDens) call RegPack(RF, InData%KinVisc) call RegPack(RF, InData%MHK) call RegPack(RF, InData%WtrDpth) @@ -1029,6 +1603,10 @@ subroutine FVW_PackParam(RF, Indata) call RegPack(RF, InData%Induction) call RegPack(RF, InData%kFrozenNWStart) call RegPack(RF, InData%kFrozenNWEnd) + call RegPack(RF, InData%zGround) + call RegPack(RF, InData%zGroundPush) + call RegPack(RF, InData%nSrcPnlUpdate) + call FVW_PackT_SrcPanlParam(RF, InData%SrcPnl) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -1085,6 +1663,7 @@ subroutine FVW_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%PartPerSegment); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%DTaero); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%DTfvw); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%AirDens); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%KinVisc); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%MHK); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WtrDpth); if (RegCheckErr(RF, RoutineName)) return @@ -1102,6 +1681,10 @@ subroutine FVW_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%Induction); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%kFrozenNWStart); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%kFrozenNWEnd); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%zGround); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%zGroundPush); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%nSrcPnlUpdate); if (RegCheckErr(RF, RoutineName)) return + call FVW_UnpackT_SrcPanlParam(RF, OutData%SrcPnl) ! SrcPnl end subroutine subroutine FVW_CopyWng_ContinuousStateType(SrcWng_ContinuousStateTypeData, DstWng_ContinuousStateTypeData, CtrlCode, ErrStat, ErrMsg) @@ -2386,6 +2969,9 @@ subroutine FVW_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) call FVW_CopyT_Part(SrcMiscData%Part, DstMiscData%Part, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return + call FVW_CopyT_SrcPanlMisc(SrcMiscData%SrcPnl, DstMiscData%SrcPnl, CtrlCode, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return if (allocated(SrcMiscData%CPs)) then LB(1:2) = lbound(SrcMiscData%CPs) UB(1:2) = ubound(SrcMiscData%CPs) @@ -2462,6 +3048,8 @@ subroutine FVW_DestroyMisc(MiscData, ErrStat, ErrMsg) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) call FVW_DestroyT_Part(MiscData%Part, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + call FVW_DestroyT_SrcPanlMisc(MiscData%SrcPnl, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (allocated(MiscData%CPs)) then deallocate(MiscData%CPs) end if @@ -2512,6 +3100,7 @@ subroutine FVW_PackMisc(RF, Indata) call RegPack(RF, InData%UA_Flag) call FVW_PackT_Sgmt(RF, InData%Sgmt) call FVW_PackT_Part(RF, InData%Part) + call FVW_PackT_SrcPanlMisc(RF, InData%SrcPnl) call RegPackAlloc(RF, InData%CPs) call RegPackAlloc(RF, InData%Uind) call RegPack(RF, allocated(InData%GridOutputs)) @@ -2566,6 +3155,7 @@ subroutine FVW_UnPackMisc(RF, OutData) call RegUnpack(RF, OutData%UA_Flag); if (RegCheckErr(RF, RoutineName)) return call FVW_UnpackT_Sgmt(RF, OutData%Sgmt) ! Sgmt call FVW_UnpackT_Part(RF, OutData%Part) ! Part + call FVW_UnpackT_SrcPanlMisc(RF, OutData%SrcPnl) ! SrcPnl call RegUnpackAlloc(RF, OutData%CPs); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%Uind); if (RegCheckErr(RF, RoutineName)) return if (allocated(OutData%GridOutputs)) deallocate(OutData%GridOutputs) @@ -3092,6 +3682,9 @@ subroutine FVW_CopyConstrState(SrcConstrStateData, DstConstrStateData, CtrlCode, end do end if DstConstrStateData%residual = SrcConstrStateData%residual + call FVW_CopyT_SrcPanlVar(SrcConstrStateData%SrcPnl, DstConstrStateData%SrcPnl, CtrlCode, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return end subroutine subroutine FVW_DestroyConstrState(ConstrStateData, ErrStat, ErrMsg) @@ -3114,6 +3707,8 @@ subroutine FVW_DestroyConstrState(ConstrStateData, ErrStat, ErrMsg) end do deallocate(ConstrStateData%W) end if + call FVW_DestroyT_SrcPanlVar(ConstrStateData%SrcPnl, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end subroutine subroutine FVW_PackConstrState(RF, Indata) @@ -3133,6 +3728,7 @@ subroutine FVW_PackConstrState(RF, Indata) end do end if call RegPack(RF, InData%residual) + call FVW_PackT_SrcPanlVar(RF, InData%SrcPnl) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -3159,6 +3755,7 @@ subroutine FVW_UnPackConstrState(RF, OutData) end do end if call RegUnpack(RF, OutData%residual); if (RegCheckErr(RF, RoutineName)) return + call FVW_UnpackT_SrcPanlVar(RF, OutData%SrcPnl) ! SrcPnl end subroutine subroutine FVW_CopyOtherState(SrcOtherStateData, DstOtherStateData, CtrlCode, ErrStat, ErrMsg) @@ -3174,7 +3771,7 @@ subroutine FVW_CopyOtherState(SrcOtherStateData, DstOtherStateData, CtrlCode, Er character(*), parameter :: RoutineName = 'FVW_CopyOtherState' ErrStat = ErrID_None ErrMsg = '' - DstOtherStateData%Dummy = SrcOtherStateData%Dummy + DstOtherStateData%ShedScale = SrcOtherStateData%ShedScale if (allocated(SrcOtherStateData%UA)) then LB(1:1) = lbound(SrcOtherStateData%UA) UB(1:1) = ubound(SrcOtherStateData%UA) @@ -3191,6 +3788,7 @@ subroutine FVW_CopyOtherState(SrcOtherStateData, DstOtherStateData, CtrlCode, Er if (ErrStat >= AbortErrLev) return end do end if + DstOtherStateData%Initialized = SrcOtherStateData%Initialized end subroutine subroutine FVW_DestroyOtherState(OtherStateData, ErrStat, ErrMsg) @@ -3222,7 +3820,7 @@ subroutine FVW_PackOtherState(RF, Indata) integer(B4Ki) :: i1 integer(B4Ki) :: LB(1), UB(1) if (RF%ErrStat >= AbortErrLev) return - call RegPack(RF, InData%Dummy) + call RegPack(RF, InData%ShedScale) call RegPack(RF, allocated(InData%UA)) if (allocated(InData%UA)) then call RegPackBounds(RF, 1, lbound(InData%UA), ubound(InData%UA)) @@ -3232,6 +3830,7 @@ subroutine FVW_PackOtherState(RF, Indata) call UA_PackOtherState(RF, InData%UA(i1)) end do end if + call RegPack(RF, InData%Initialized) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -3244,7 +3843,7 @@ subroutine FVW_UnPackOtherState(RF, OutData) integer(IntKi) :: stat logical :: IsAllocAssoc if (RF%ErrStat /= ErrID_None) return - call RegUnpack(RF, OutData%Dummy); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%ShedScale); if (RegCheckErr(RF, RoutineName)) return if (allocated(OutData%UA)) deallocate(OutData%UA) call RegUnpack(RF, IsAllocAssoc); if (RegCheckErr(RF, RoutineName)) return if (IsAllocAssoc) then @@ -3258,6 +3857,7 @@ subroutine FVW_UnPackOtherState(RF, OutData) call UA_UnpackOtherState(RF, OutData%UA(i1)) ! UA end do end if + call RegUnpack(RF, OutData%Initialized); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine FVW_CopyWng_InitInputType(SrcWng_InitInputTypeData, DstWng_InitInputTypeData, CtrlCode, ErrStat, ErrMsg) @@ -3409,6 +4009,7 @@ subroutine FVW_CopyInitInput(SrcInitInputData, DstInitInputData, CtrlCode, ErrSt end if DstInitInputData%numBladeNodes = SrcInitInputData%numBladeNodes DstInitInputData%DTaero = SrcInitInputData%DTaero + DstInitInputData%AirDens = SrcInitInputData%AirDens DstInitInputData%KinVisc = SrcInitInputData%KinVisc DstInitInputData%MHK = SrcInitInputData%MHK DstInitInputData%WtrDpth = SrcInitInputData%WtrDpth @@ -3480,6 +4081,7 @@ subroutine FVW_PackInitInput(RF, Indata) end if call RegPack(RF, InData%numBladeNodes) call RegPack(RF, InData%DTaero) + call RegPack(RF, InData%AirDens) call RegPack(RF, InData%KinVisc) call RegPack(RF, InData%MHK) call RegPack(RF, InData%WtrDpth) @@ -3527,6 +4129,7 @@ subroutine FVW_UnPackInitInput(RF, OutData) end if call RegUnpack(RF, OutData%numBladeNodes); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%DTaero); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%AirDens); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%KinVisc); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%MHK); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WtrDpth); if (RegCheckErr(RF, RoutineName)) return @@ -3575,6 +4178,7 @@ subroutine FVW_CopyInputFile(SrcInputFileData, DstInputFileData, CtrlCode, ErrSt DstInputFileData%VTKBlades = SrcInputFileData%VTKBlades DstInputFileData%DTvtk = SrcInputFileData%DTvtk DstInputFileData%VTKCoord = SrcInputFileData%VTKCoord + DstInputFileData%SrcPnlFile = SrcInputFileData%SrcPnlFile end subroutine subroutine FVW_DestroyInputFile(InputFileData, ErrStat, ErrMsg) @@ -3623,6 +4227,7 @@ subroutine FVW_PackInputFile(RF, Indata) call RegPack(RF, InData%VTKBlades) call RegPack(RF, InData%DTvtk) call RegPack(RF, InData%VTKCoord) + call RegPack(RF, InData%SrcPnlFile) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -3663,6 +4268,7 @@ subroutine FVW_UnPackInputFile(RF, OutData) call RegUnpack(RF, OutData%VTKBlades); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%DTvtk); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%VTKCoord); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%SrcPnlFile); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine FVW_CopyInitOutput(SrcInitOutputData, DstInitOutputData, CtrlCode, ErrStat, ErrMsg) diff --git a/modules/aerodyn/src/FVW_VortexTools.f90 b/modules/aerodyn/src/FVW_VortexTools.f90 index 7537b8e814..870d19dae4 100644 --- a/modules/aerodyn/src/FVW_VortexTools.f90 +++ b/modules/aerodyn/src/FVW_VortexTools.f90 @@ -8,6 +8,7 @@ module FVW_VortexTools use NWTC_LIBRARY, only: ReKi, IntKi, num2lstr, ErrID_Fatal, ErrID_None, EqualRealNos !< Not desired use NWTC_LIBRARY, only: InterpArray + use NWTC_LIBRARY, only: Pi use FVW_BiotSavart, only: fourpi_inv implicit none @@ -2147,7 +2148,55 @@ subroutine curl_regular_grid(f,rotf,ix,iy,iz,nx,ny,nz,dx,dy,dz) !endif end subroutine curl_regular_grid - + !> Return an ellipsoid made of panels + subroutine EllipsoidPanels(n1, n2, L1, L2, L3, Points, IDs) + integer, intent(in) :: n1, n2 !< Number of points in main and secondary direction (spherical coordinates) + real(ReKi), intent(in) :: L1, L2, L3 !< Main axes of ellipsoid in 3 directions + real(ReKi), dimension(:,:), allocatable, intent(out) :: Points !< 3 x nP + integer, dimension(:,:), allocatable, intent(out) :: IDs !< 4 x nPanels, connectivity + ! Variables + integer :: ip, it, node,e + real(ReKi) :: phi !< + real(ReKi) :: theta !< + real(ReKi) :: dphi !< + real(ReKi) :: dtheta !< + integer :: nphi !< + integer :: ntheta !< + real(ReKi), dimension(3) :: P + ntheta = n2+1 + nphi = n1 + dtheta = pi/(ntheta-1) ! -1 + dphi = 2*pi/nphi + allocate(Points(3, nphi*ntheta )) + allocate(IDs (4, nphi*(ntheta-1))) + IDs = 0.0_ReKi + Points = 0.0_ReKi + + ! Building Points + do it=1,ntheta + theta=-pi/2 + real((it-1),ReKi)*dtheta ! + do ip=1,nphi + phi = (ip-1)*dphi + node = (it-1)*nphi +ip + Points(1, node) = L1*cos(theta)*cos(phi) + Points(2, node) = L2*cos(theta)*sin(phi) + Points(3, node) = L3*sin(theta) + end do + end do + ! Ugly for loops for Panl + do it=1,ntheta-1 + do ip=1,nphi-1 + node = (it-1)*nphi +ip + e = (it-1)*(nphi)+ip + IDs(:, e)=(/node,node+nphi,node+nphi+1,node+1 /) + end do + ! special case, ip=nphi + ip = nphi ! + node = (it-1)*nphi +ip + e = (it-1)*(nphi)+ip + IDs(:, e) = (/node,node+nphi,node+nphi-ip+1,node-ip+1/) + end do + end subroutine EllipsoidPanels ! --- TIC TOC MODULE !> Simpler version of matlab tic diff --git a/modules/aerodyn/tests/test_FVW_testsuite.F90 b/modules/aerodyn/tests/test_FVW_testsuite.F90 new file mode 100644 index 0000000000..17d40598f6 --- /dev/null +++ b/modules/aerodyn/tests/test_FVW_testsuite.F90 @@ -0,0 +1,42 @@ +@test +subroutine test_AD_FVW() + ! test branches + ! - known valid checks for various FVW routines (contained in own module) + ! - known invalid rotation matrix: halve the angle of the diagonal elements + + use pFUnit_mod + use NWTC_Num + use FVW_Tests + + implicit none + + integer(IntKi) :: ErrStat + character(ErrMsgLen) :: ErrMsg + character(1024) :: testname + + ! initialize NWTC_Num constants + call SetConstants() + + ! --- Run all tests at once + call FVW_RunTests (errStat,errMsg); + @assertEqual(0, errStat, 'All FVW tests ') + + ! --- Run individual tests + !call Test_LinSolve (errStat,errMsg); + !@assertEqual(0, errStat, 'Test_LinSolve ') + !call Test_SrcPnl_Sphere (errStat,errMsg); + !@assertEqual(0, errStat, 'Test_SrcPnl_Sphere ') + !call Test_BiotSavart_SrcPnl (errStat,errMsg); + !@assertEqual(0, errStat, 'Test_BiotSavart_SrcPnl ') + !call Test_BiotSavart_Sgmt (errStat,errMsg); + !@assertEqual(0, errStat, 'Test_BiotSavart_Sgmt ') + !call Test_BiotSavart_Part (errStat,errMsg); + !@assertEqual(0, errStat, 'Test_BiotSavart_Part ') + !call Test_BiotSavart_PartTree (errStat,errMsg); + !@assertEqual(0, errStat, 'Test_BiotSavart_PartTree ') + !call Test_SegmentsToPart (errStat,errMsg); + !@assertEqual(0, errStat, 'Test_SegmentsToPart ') + !call FVW_Test_WakeInducedVelocities(errStat,errMsg); + !@assertEqual(0, errStat, 'FVW_Test_WakeInducedVelocities') + +end subroutine test_AD_FVW diff --git a/modules/nwtc-library/src/VTK.f90 b/modules/nwtc-library/src/VTK.f90 index e4ffc614d0..1af0e2c6ed 100644 --- a/modules/nwtc-library/src/VTK.f90 +++ b/modules/nwtc-library/src/VTK.f90 @@ -5,9 +5,10 @@ module VTK use Precision, only: IntKi, SiKi, ReKi - use NWTC_Base, only: ErrID_None, ErrID_Fatal, AbortErrLev, ErrMsgLen, SetErrStat + use NWTC_Base, only: ErrID_None, ErrID_Fatal, ErrID_Warn, AbortErrLev, ErrMsgLen use NWTC_IO, only: GetNewUnit, NewLine, WrScr, ReadStr, OpenFOutFile - use NWTC_IO, only: OpenFinpFile, ReadCom, Conv2UC + use NWTC_IO, only: OpenFinpFile, ReadCom, Conv2UC, num2lstr, AllocAry + use NWTC_IO, only: SetErrStat implicit none @@ -30,6 +31,14 @@ module VTK real(ReKi),dimension(3) :: PO_g END TYPE VTK_Misc + type, public :: vtk_field + character(len=255) :: fieldname ='' !< fieldname + real(ReKi), dimension(:,:), allocatable :: values ! < + type(vtk_field), pointer :: next => null() + end type vtk_field + + + interface vtk_dataset_structured_grid; module procedure & vtk_dataset_structured_grid_flat, & vtk_dataset_structured_grid_grid @@ -48,6 +57,7 @@ module VTK end interface interface vtk_cell_data_scalar; module procedure & vtk_cell_data_scalar_1d,& + vtk_cell_data_scalar_1di,& vtk_cell_data_scalar_2d end interface @@ -73,8 +83,14 @@ module VTK public :: vtk_cell_data_vector public :: WrVTK_SP_header public :: WrVTK_SP_vectors3D - public :: ReadVTK_SP_info + public :: ReadVTK_SP_info !< Structured points public :: ReadVTK_SP_vectors + public :: ReadVTK_PD_info !< Polydata + public :: ReadVTK_CD_fields !< read fields + ! --- VTK fields + public :: VTK_printFields + public :: VTK_getField + public :: VTK_destroyFields ! --- VTK XML routines public :: WrVTK_header public :: WrVTK_footer @@ -356,7 +372,285 @@ SUBROUTINE ReadVTK_SP_vectors( FileName, Un, dims, gridVals, ErrStat, ErrMsg ) end if END SUBROUTINE ReadVTK_SP_vectors - + +!======================================================================= +!> Read the header/points/polygons for a vtk, ascii, polydata (PD) file. +!! NOTE: the file is not closed upon exiting the routine, unless an error occurs. +!! Example: +!! # vtk DataFile Version 2.0 +!! SourcePanels - Ellipsoid - n1=67 - n2=65 +!! ASCII +!! +!! DATASET POLYDATA +!! POINTS 3 double +!! 0.0 0.0 -1.0 +!! 0.0 0.0 0.0 +!! 0.0 0.0 1.0 +!! +!! POLYGONS 1 5 +!! 4 0 1 2 0 +!! +!! CELL_DATA 1 +!! SCALARS Pressure double +!! LOOKUP_TABLE default +!! 1.0 +!! VECTORS Velocity double +!! 0.0 0.0 1.0 +!! + subroutine ReadVTK_PD_info(filename, descr, points, polygons, fid, errStat, errMsg) + character(*) , intent(in ) :: filename !< name of output file + character(1024), intent( out) :: descr !< line describing the contents of the file + real(ReKi) , intent( out), dimension(:,:), allocatable :: points !< Points 3 x nPoints + integer(IntKi) , intent( out), dimension(:,:), allocatable :: polygons !< Polygons connectivity, typically 4 x nPolygons + integer(IntKi) , intent( out) :: fid !< unit number of opened file + integer(IntKi) , intent( out) :: errStat !< error level/status of openfoutfile operation + character(*) , intent( out) :: errMsg !< message when error occurs + integer(IntKi) :: errStat2 ! local error level/status of openfoutfile operation + character(ErrMsgLen) :: errMsg2 ! local message when error occurs + character(1024) :: sdummy1, sdummy2, line + integer(IntKi) :: iline, i + integer(IntKi) :: nPoints, nPoly, nTot, nPerLine + integer(IntKi) :: iSection !< keep track of section we are currently reading 0=header, 2=points, 3=polygons + errStat = ErrID_None + errMsg = '' + call GetNewUnit(fid, errStat2, errMsg2) + call OpenFInpFile(fid, trim(filename), errStat2, errMsg2); if(Failed()) return + ! + descr='' + iline = 0 + iSection = 0 + do + iline = iline+1 + call ReadStr(fid, Filename, line, 'line', 'line', errStat2, errMsg2); if(Failed()) return + call Conv2UC(line) + if (index(adjustl(line), 'DATASET POLYDATA') == 1) then + iSection = 1 + + else if (index(adjustl(line), 'POINTS') == 1) then + if (iSection/=1) then + errStat2 = ErrID_Fatal + errMsg2 = 'Cannot read section `Points` if section `Dataset polydata` not read' + endif + ! --- Reading points + iSection = 2 + read(line, *, iostat=errStat2) sdummy1, nPoints, sdummy2 + errMsg2='Problem reading number of points'; if(Failed()) return + call AllocAry(points, 3, nPoints, 'points', errStat2, errMsg2); if(Failed()) return + do i =1,nPoints + iline = iline + i + read(fid, *, iostat=errStat2) points(:,i) + errMsg2='Problem reading point '//trim(num2lstr(i)); if(Failed()) return + enddo + + else if (index(adjustl(line), 'POLYGONS') == 1) then + if (iSection/=2) then + errStat2 = ErrID_Fatal + errMsg2 = 'Cannot read section `Polygons` if section `Points` not read' + endif + ! --- Reading polygons + iSection = 3 + read(line, *, iostat=errStat2) sdummy1, nPoly, nTot + errMsg2='Problem reading number of polygons'; if(Failed()) return + call AllocAry(polygons, 4, nPoly, 'points', errStat2, errMsg2); if(Failed()) return + do i =1,nPoly + iline = iline + i + read(fid, *, iostat=errStat2) nPerLine, polygons(:,i) + errMsg2='Problem reading polygon '//trim(num2lstr(i)); if(Failed()) return + enddo + exit + else + if (iSection==0) then + descr = trim(descr)//'|'//trim(line) + endif + end if + end do + if (iSection/=3) then + errStat2 = ErrID_Fatal + errMsg2 = 'Not all sections were read properly for VTK polydata file.' + endif + contains + logical function Failed() + if (errStat2/=ErrID_None) then + errStat2 = ErrID_Fatal ! All errors are fatal + errMsg2 = trim(errMsg2)//', on line '//trim(num2lstr(iline))//char(10)//' File: '//trim(filename) + close(fid) + endif + call SeterrStat(errStat2, errMsg2, errStat, errMsg, 'ReadVTK_PD_info') + Failed = errStat >= AbortErrLev + end function + end subroutine ReadVTK_PD_info + + !> Read Cell Data in file, return fields, a list of scalar or vector fields + subroutine ReadVTK_CD_fields(filename, fid, nCells, fields, errStat, errMsg) + use, intrinsic :: iso_fortran_env, only : iostat_end + character(*) , intent(in ) :: filename !< name of output file + integer(IntKi) , intent(in ) :: fid !< unit number of opened file + integer(IntKi) , intent(in ) :: nCells !< Number of cells + type(vtk_field), pointer :: fields !< Fields to be read + integer(IntKi) , intent(out) :: errStat !< error level/status of openfoutfile operation + character(*) , intent(out) :: errMsg !< message when error occurs + integer(IntKi) :: errStat2 ! local error level/status of openfoutfile operation + character(ErrMsgLen) :: errMsg2 ! local message when error occurs + character(1024) :: line + character(128) :: sType, fieldname, datatype + type(vtk_field), pointer :: head + type(vtk_field), pointer :: new_field + errStat = ErrID_None + errMsg = '' + ! Read untill cell_data + call ReadVTK_until(filename, fid, 'CELL_DATA', errStat2, errMsg2); if(Failed()) return + if (errStat2==ErrID_Warn) return ! We return if not found + ! + do + read(fid, '(A)', iostat = errStat2) line + if (errStat2 == iostat_end) then + close(fid) + return + endif + + if (index(adjustl(line), 'SCALARS') > 0) then + call prepare_new_field(1); if(Failed()) return + call read_field(); if(Failed()) return + call prepend_field() + + else if (index(adjustl(line), 'VECTORS') > 0) then + call prepare_new_field(3); if(Failed()) return + call read_field(); if(Failed()) return + call prepend_field() + else + ! pass + endif + end do + ! Should never be reached + contains + logical function Failed() + call SeterrStat(errStat2, errMsg2, errStat, errMsg, 'ReadVTK_CD_fields') + Failed = errStat >= AbortErrLev + if (Failed .and. fid>0) close(fid) + end function + + subroutine prepare_new_field(nDim) + integer, intent(in) :: nDim + read(line, *) sType, fieldname, datatype + call Conv2UC(fieldname) + !print*,'VTK: Cell Data: Type: ',trim(sType),', Name: ',trim(fieldname), ', Type:',trim(datatype) + errMsg2 = 'Error allocating new field' + allocate(new_field, stat = errStat2) + allocate(new_field%values(nDim, nCells), stat = errStat2) + if (errStat2/=0) then + errStat2 = ErrID_Fatal + return + endif + new_field%fieldname = fieldname + ! If the line starts with LOOKUP_TABLE DEFAULT, we ignore the line + errMsg2 = 'Error reading lookup table' + read(fid, '(A)', iostat = errStat2) line + if (errStat2/=0) then + errStat2 = ErrID_Fatal + return + endif + call Conv2UC(line) + if (index(adjustl(line), 'LOOKUP_TABLE DEFAULT') <= 0) then + backspace(fid) + endif + end subroutine + subroutine read_field() + integer :: i + do i = 1, nCells + read(fid, *, iostat=errStat2) new_field%values(:,i) + if (errStat2/=0) then + errStat2= ErrID_Fatal + errMsg2 ='Error reading line '//trim(num2lstr(i))//' of table '//trim(new_field%fieldname)//' in VTK file: '//trim(filename) + return + endif + end do + end subroutine + + subroutine prepend_field() + if (associated(fields)) new_field%next => fields + fields => new_field + nullify(new_field) + end subroutine + + end subroutine ReadVTK_CD_fields + + !> Loop and print the list + subroutine VTK_printFields(list) + type(vtk_field), pointer, intent(in) :: list + type(vtk_field), pointer :: field + type(vtk_field), pointer :: head + integer :: n,m + head => list + if (.not.associated(head)) then + call WrScr('VTK list of fields empty.') + endif + do while (associated(head)) + if (allocated(head%values)) then + n = size(head%values,1) + m = size(head%values,2) + else + n=-1 + m=-1 + endif + print *, 'VTK: Fieldname: ', head%fieldname(1:10), ' alloc: ', allocated(head%values), ' shape:',n, m + head => head%next + end do + end subroutine VTK_printFields + + ! Loop through list of fields and find the field with the requested fieldname + subroutine VTK_getField(list, fieldname, field) + type(vtk_field), pointer, intent(in) :: list + character(*), intent(in) :: fieldname + type(vtk_field), pointer, intent(out) :: field + character(len(fieldname)) :: fieldname_UP + type(vtk_field), pointer :: head + fieldname_UP = adjustl(fieldname) + call Conv2UC(fieldname_UP) + nullify(field) + head => list + do while (associated(head)) + if (trim(adjustl(head%fieldname)) == trim(fieldname_UP)) then + field => head + exit + else + head => head%next + endif + end do + end subroutine + + !> Deallocate the list of fields + subroutine VTK_destroyFields(list) + type(vtk_field), pointer :: list, next + do while (associated(list)) + next => list%next + if (allocated(list%values)) deallocate(list%values) + deallocate(list) + list => next + end do + end subroutine VTK_destroyFields + + + !> Read until an UpperCase Keyword is found + subroutine ReadVTK_until(filename, fid, KeyUpper, errStat, errMsg) + character(*) , intent(in ) :: filename !< name of output file + integer(IntKi) , intent(in ) :: fid !< unit number of opened file + character(*) , intent(in ) :: KeyUpper !< Key in upperCase + integer(IntKi) , intent(out) :: errStat !< error level/status of openfoutfile operation + character(*) , intent(out) :: errMsg !< message when error occurs + integer(IntKi) :: errStat2 ! local error level/status of openfoutfile operation + character(1024) :: line + errStat = ErrID_None + errMsg = '' + line='' + do + read(fid, *, iostat=errStat2) line ! Read each line until "CELL_DATA" is encountered + if (errStat2/=0) exit ! We return from this for loop and fail + if (index(adjustl(line), KeyUpper) > 0) return ! Key was found + end do + errStat = ErrID_Warn + errMsg = 'Key '//trim(KeyUpper)//' not found in VTK file:'//trim(filename) + end subroutine ReadVTK_until + !======================================================================= !> This routine writes out the heading for an vtk, ascii, structured_points dataset file . !! It tries to open a text file for writing and returns the Unit number of the opened file. @@ -901,6 +1195,25 @@ subroutine vtk_cell_data_scalar_1d(D,sname,mvtk) endif end subroutine + subroutine vtk_cell_data_scalar_1di(D,sname,mvtk) + integer(IntKi), dimension(:),intent(in)::D + character(len=*),intent(in) ::sname + type(VTK_Misc),intent(inout) :: mvtk + + if ( mvtk%bFileOpen ) then + if (mvtk%bBinary) then + write(mvtk%vtk_unit)'SCALARS '//trim(sname)//' int 1'//NewLine + write(mvtk%vtk_unit)'LOOKUP_TABLE default'//NewLine + write(mvtk%vtk_unit)D + write(mvtk%vtk_unit)NewLine + else + write(mvtk%vtk_unit,fmt='(A,A,A)') 'SCALARS ', sname, ' int' + write(mvtk%vtk_unit,'(A)') 'LOOKUP_TABLE default' + write(mvtk%vtk_unit,'(1I0)')D + endif + endif + end subroutine + subroutine vtk_cell_data_scalar_2d(D,sname,mvtk) real(ReKi), dimension(:,:),intent(in)::D character(len=*),intent(in) ::sname