diff --git a/docs/source/user/aerodyn/input.rst b/docs/source/user/aerodyn/input.rst
index 93fc266d3e..66e6b59c6e 100644
--- a/docs/source/user/aerodyn/input.rst
+++ b/docs/source/user/aerodyn/input.rst
@@ -90,7 +90,7 @@ When ``Wake_Mod`` is set to 3, the free vortex wake model is used, also referred
:numref:`OLAF`). ``Wake_Mod`` cannot be set to 3 during linearization analyses.
.. note::
- Link to old inputs: The previous input `WakeMod` is removed, `WakeMod=2` used to mean DBEMT, but this now controled using `DBEMT_Mod`.
+ Link to old inputs: The previous input `WakeMod` is removed, `WakeMod=2` used to mean DBEMT, but this now controlled using `DBEMT_Mod`.
`WakeMod=2` is a placeholder for future induction calculation method.
@@ -328,7 +328,7 @@ Most ``UA_Mod`` will require `AoA34` to be set to true. But when using quasi-ste
This feature is currently not implemented due to a lag between the `dev` and `dev-unstable` branch.
.. note::
- Link to previous inputs: `AFAeroMod=1` implies `AoA34=False`. But to have a fair comparison between steady and unsteady aerodynamics model, it would be best to set `AoA34=True` when when doing quasi-steady aero.
+ Link to previous inputs: `AFAeroMod=1` implies `AoA34=False`. But to have a fair comparison between steady and unsteady aerodynamics model, it would be best to set `AoA34=True` when doing quasi-steady aero.
diff --git a/modules/aerodyn/python-lib/aerodyn_inflow_library.py b/modules/aerodyn/python-lib/aerodyn_inflow_library.py
index 9539ddff68..add1bcf5ac 100644
--- a/modules/aerodyn/python-lib/aerodyn_inflow_library.py
+++ b/modules/aerodyn/python-lib/aerodyn_inflow_library.py
@@ -33,7 +33,7 @@
c_byte,
c_int,
c_double,
- c_float,
+ c_float,
c_char,
c_char_p,
c_wchar,
@@ -98,7 +98,7 @@ def __init__(self, library_path):
self.WtrDpth = 0.0 # Water depth (m)
self.MSL2SWL = 0.0 # Offset between still-water level and mean sea level (m) [positive upward]
- # flags
+ # flags
self.storeHHVel = 1 # 0=false, 1=true
self.transposeDCM= 1 # 0=false, 1=true
self.debuglevel = 0 # 0-4 levels
@@ -149,6 +149,7 @@ def __init__(self, library_path):
self.numMeshPts = 1
self.initMeshPos = np.zeros(shape=(self.numMeshPts,3),dtype=c_float ) # Nx3 array [x,y,z]
self.initMeshOrient = np.zeros(shape=(self.numMeshPts,9),dtype=c_double) # Nx9 array [r11,r12,r13,r21,r22,r23,r31,r32,r33]
+ self.meshPtToBladeNum = np.zeros(shape=(self.numMeshPts),dtype=c_int) # Nx1 array [blade number]
# OutRootName
# If HD writes a file (echo, summary, or other), use this for the
@@ -161,11 +162,11 @@ def _initialize_routines(self):
self.ADI_C_PreInit.argtypes = [
POINTER(c_int), # numTurbines
POINTER(c_int), # transposeDCM
- POINTER(c_int), # debuglevel
+ POINTER(c_int), # debuglevel
POINTER(c_int), # ErrStat_C
POINTER(c_char) # ErrMsg_C
]
- self.ADI_C_PreInit.restype = c_int
+ self.ADI_C_PreInit.restype = c_int
# setup one rotor
self.ADI_C_SetupRotor.argtypes = [
@@ -182,6 +183,7 @@ def _initialize_routines(self):
POINTER(c_int), # numMeshPts
POINTER(c_float), # initMeshPos_flat
POINTER(c_double), # initMeshOrient_flat
+ POINTER(c_int), # meshPtToBladeNum
POINTER(c_int), # ErrStat_C
POINTER(c_char) # ErrMsg_C
]
@@ -194,7 +196,7 @@ def _initialize_routines(self):
POINTER(c_int), # IfW input file passed as string
POINTER(c_char_p), # IfW input file as string
POINTER(c_int), # IfW input file string length
- POINTER(c_char), # OutRootName
+ POINTER(c_char), # OutRootName
POINTER(c_float), # gravity
POINTER(c_float), # defFldDens
POINTER(c_float), # defKinVisc
@@ -203,9 +205,9 @@ def _initialize_routines(self):
POINTER(c_float), # defPvap
POINTER(c_float), # WtrDpth
POINTER(c_float), # MSL2SWL
- POINTER(c_int), # InterpOrder
+ POINTER(c_int), # InterpOrder
POINTER(c_double), # dt
- POINTER(c_double), # tmax
+ POINTER(c_double), # tmax
POINTER(c_int), # storeHHVel
POINTER(c_int), # WrVTK
POINTER(c_int), # WrVTK_Type
@@ -219,16 +221,16 @@ def _initialize_routines(self):
POINTER(c_int), # ErrStat_C
POINTER(c_char) # ErrMsg_C
]
- self.ADI_C_Init.restype = c_int
+ self.ADI_C_Init.restype = c_int
#self.ADI_C_ReInit.argtypes = [
- # POINTER(c_double), # t_initial
+ # POINTER(c_double), # t_initial
# POINTER(c_double), # dt
- # POINTER(c_double), # tmax
+ # POINTER(c_double), # tmax
# POINTER(c_int), # ErrStat_C
# POINTER(c_char) # ErrMsg_C
#]
- #self.ADI_C_ReInit.restype = c_int
+ #self.ADI_C_ReInit.restype = c_int
self.ADI_C_SetRotorMotion.argtypes = [
POINTER(c_int), # iturb
@@ -322,6 +324,7 @@ def adi_setuprotor(self,iturb,isHAWT,turbRefPos):
initRootOrient_flat_c = self.flatOrientArr(self._initNumBlades, self.numBlades,self.initRootOrient, 'RootOrient')
initMeshPos_flat_c = self.flatPosArr( self._initNumMeshPts,self.numMeshPts,self.initMeshPos, 'MeshPos')
initMeshOrient_flat_c = self.flatOrientArr(self._initNumMeshPts,self.numMeshPts,self.initMeshOrient,'MeshOrient')
+ initMeshPtToBladeNum_flat_c = (c_int * len(self.meshPtToBladeNum))(*self.meshPtToBladeNum)
self.ADI_C_SetupRotor(
c_int(iturb), # IN: iturb -- current turbine number
@@ -337,6 +340,7 @@ def adi_setuprotor(self,iturb,isHAWT,turbRefPos):
byref(c_int(self.numMeshPts)), # IN: number of mesh points expected
initMeshPos_flat_c, # IN: initMeshPos -- initial node positions in flat array of 3*numMeshPts
initMeshOrient_flat_c, # IN: initMeshOrient -- initial node orientation DCMs in flat array of 9*numMeshPts
+ initMeshPtToBladeNum_flat_c, # IN: initMeshPtToBladeNum -- initial mesh point to blade number mapping
byref(self.error_status_c), # OUT: ErrStat_C
self.error_message_c # OUT: ErrMsg_C
)
@@ -402,7 +406,7 @@ def adi_init(self, AD_input_string_array, IfW_input_string_array):
)
self.check_error()
-
+
# Initialize output channels
self.numChannels = self._numChannels_c.value
@@ -497,7 +501,7 @@ def adi_getrotorloads(self, iturb, meshFrcMom):
)
self.check_error()
-
+
## Reshape Force/Moment into [N,6]
count = 0
for j in range(0,self.numMeshPts):
@@ -518,7 +522,7 @@ def adi_calcOutput(self, time, outputChannelValues):
# Run ADI_C_CalcOutput
self.ADI_C_CalcOutput(
- byref(c_double(time)), # IN: time at which to calculate output forces
+ byref(c_double(time)), # IN: time at which to calculate output forces
outputChannelValues_c, # OUT: output channel values as described in input file
byref(self.error_status_c), # OUT: ErrStat_C
self.error_message_c # OUT: ErrMsg_C
@@ -535,7 +539,7 @@ def adi_updateStates(self, time, timeNext):
# Run AeroDyn_Inflow_UpdateStates_c
self.ADI_C_UpdateStates(
- byref(c_double(time)), # IN: time at which to calculate output forces
+ byref(c_double(time)), # IN: time at which to calculate output forces
byref(c_double(timeNext)), # IN: time T+dt we are stepping to
byref(self.error_status_c), # OUT: ErrStat_C
self.error_message_c # OUT: ErrMsg_C
@@ -662,7 +666,7 @@ def check_init_hubroot(self):
print("Expecting a 9 element array for initNacelleOrient DCM [r11,r12,r13,r21,r22,r23,r31,r32,r33]")
self.adi_end()
raise Exception("\nAeroDyn terminated prematurely.")
-
+
def check_init_mesh(self):
#print("shape of initMeshPos ", self.initMeshPos.shape)
diff --git a/modules/aerodyn/src/AeroDyn.f90 b/modules/aerodyn/src/AeroDyn.f90
index b2bd2bcd80..a33ce9c153 100644
--- a/modules/aerodyn/src/AeroDyn.f90
+++ b/modules/aerodyn/src/AeroDyn.f90
@@ -360,7 +360,6 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut
endif
enddo
- print*,'AD_Init: AeroProjMod B:',AeroProjMod
! -----------------------------------------------------------------
! Read the AeroDyn blade files, or copy from passed input
@@ -4495,6 +4494,9 @@ SUBROUTINE Init_BEMTmodule( InputFileData, RotInputFileData, u_AD, u, p, p_AD, x
elseif (p%AeroProjMod==APM_LiftingLine) then
Label='Projection: Lifting Line'
else
+ ! Normally we wouldn't want to do a print or STOP, but we
+ ! should never get here unless a programmer made a mistake.
+ ! I'll leave this as is for now. - ADP
print*,'Invalid projection method'
STOP
endif
@@ -4745,7 +4747,7 @@ SUBROUTINE TFin_CalcOutput(p, p_AD, u, m, y, ErrStat, ErrMsg )
V_ind = 0.0_ReKi
elseif(p%TFin%TFinIndMod==TFinIndMod_rotavg) then
! TODO TODO
- print*,'TODO TailFin: compute rotor average induced velocity'
+ call WrScr('TODO TailFin: compute rotor average induced velocity')
V_ind = 0.0_ReKi
else
STOP ! Will never happen
diff --git a/modules/aerodyn/src/AeroDyn_IO.f90 b/modules/aerodyn/src/AeroDyn_IO.f90
index 3cceff5290..4744dee847 100644
--- a/modules/aerodyn/src/AeroDyn_IO.f90
+++ b/modules/aerodyn/src/AeroDyn_IO.f90
@@ -659,6 +659,7 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade
integer(IntKi) :: CurLine !< current entry in FileInfo_In%Lines array
real(ReKi) :: TmpRe5(5) !< temporary 8 number array for reading values in
character(1024) :: sDummy !< temporary string
+ character(1024) :: tmpOutStr !< temporary string for writing to screen
logical :: wakeModProvided, frozenWakeProvided, skewModProvided, AFAeroModProvided, UAModProvided, isLegalComment, firstWarn !< Temporary for legacy purposes
logical :: AoA34_Missing
integer :: UAMod_Old
@@ -1177,19 +1178,20 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade
! pass comment lines
elseif (index(sDummy, 'SECTAVG')>1) then
read(sDummy, '(L1)') InputFileData%SectAvg
- print*,' >>> SectAvg ',InputFileData%SectAvg
+ write(tmpOutStr,*) ' >>> SectAvg ',InputFileData%SectAvg
elseif (index(sDummy, 'SA_PSIBWD')>1) then
read(sDummy, *) InputFileData%SA_PsiBwd
- print*,' >>> SA_PsiBwd ',InputFileData%SA_PsiBwd
+ write(tmpOutStr,*) ' >>> SA_PsiBwd ',InputFileData%SA_PsiBwd
elseif (index(sDummy, 'SA_PSIFWD')>1) then
read(sDummy, *) InputFileData%SA_PsiFwd
- print*,' >>> SA_PsiFwd ',InputFileData%SA_PsiFwd
+ write(tmpOutStr,*) ' >>> SA_PsiFwd ',InputFileData%SA_PsiFwd
elseif (index(sDummy, 'SA_NPERSEC')>1) then
read(sDummy, *) InputFileData%SA_nPerSec
- print*,' >>> SA_nPerSec ',InputFileData%SA_nPerSec
+ write(tmpOutStr,*) ' >>> SA_nPerSec ',InputFileData%SA_nPerSec
else
- print*,'[WARN] AeroDyn Line ignored: '//trim(sDummy)
+ write(tmpOutStr,*) '[WARN] AeroDyn Line ignored: '//trim(sDummy)
endif
+ call WrScr(trim(tmpOutStr))
enddo
@@ -1271,6 +1273,7 @@ end function newInputMissing
!-------------------------------------------------------------------------------------------------
subroutine printNewOldInputs()
+ character(1024) :: tmpStr
! Temporary HACK, for WakeMod=10, 11 or 12 use AeroProjMod 2 (will trigger PolarBEM)
if (InputFileData%Wake_Mod==10) then
call WrScr('[WARN] Wake_Mod=10 is a temporary hack. Setting BEM_Mod to 0')
@@ -1285,23 +1288,23 @@ subroutine printNewOldInputs()
!====== Summary of new AeroDyn options ===============================================================
! NOTE: remove me in future release
call WrScr('-------------- New AeroDyn inputs (with new meaning):')
- write (*,'(A20,I0)') 'Wake_Mod: ' , InputFileData%Wake_Mod
- write (*,'(A20,I0)') 'BEM_Mod: ' , InputFileData%BEM_Mod
- write (*,'(A20,L0)') 'SectAvg: ' , InputFileData%SectAvg
- write (*,'(A20,I0)') 'SectAvgWeighting: ', InputFileData%SA_Weighting
- write (*,'(A20,I0)') 'SectAvgNPoints: ', InputFileData%SA_nPerSec
- write (*,'(A20,I0)') 'DBEMT_Mod:' , InputFileData%DBEMT_Mod
- write (*,'(A20,I0)') 'Skew_Mod: ' , InputFileData%Skew_Mod
- write (*,'(A20,L0)') 'SkewMomCorr:' , InputFileData%SkewMomCorr
- write (*,'(A20,I0)') 'SkewRedistr_Mod:' , InputFileData%SkewRedistr_Mod
- write (*,'(A20,L0)') 'AoA34: ' , InputFileData%AoA34
- write (*,'(A20,I0)') 'UA_Mod: ' , InputFileData%UAMod
+ write (tmpStr,'(A20,I0)') 'Wake_Mod: ' , InputFileData%Wake_Mod; call WrScr(trim(tmpStr))
+ write (tmpStr,'(A20,I0)') 'BEM_Mod: ' , InputFileData%BEM_Mod; call WrScr(trim(tmpStr))
+ write (tmpStr,'(A20,L1)') 'SectAvg: ' , InputFileData%SectAvg; call WrScr(trim(tmpStr))
+ write (tmpStr,'(A20,I0)') 'SectAvgWeighting: ', InputFileData%SA_Weighting; call WrScr(trim(tmpStr))
+ write (tmpStr,'(A20,I0)') 'SectAvgNPoints: ', InputFileData%SA_nPerSec; call WrScr(trim(tmpStr))
+ write (tmpStr,'(A20,I0)') 'DBEMT_Mod:' , InputFileData%DBEMT_Mod; call WrScr(trim(tmpStr))
+ write (tmpStr,'(A20,I0)') 'Skew_Mod: ' , InputFileData%Skew_Mod; call WrScr(trim(tmpStr))
+ write (tmpStr,'(A20,L1)') 'SkewMomCorr:' , InputFileData%SkewMomCorr; call WrScr(trim(tmpStr))
+ write (tmpStr,'(A20,I0)') 'SkewRedistr_Mod:' , InputFileData%SkewRedistr_Mod; call WrScr(trim(tmpStr))
+ write (tmpStr,'(A20,L1)') 'AoA34: ' , InputFileData%AoA34; call WrScr(trim(tmpStr))
+ write (tmpStr,'(A20,I0)') 'UA_Mod: ' , InputFileData%UAMod; call WrScr(trim(tmpStr))
call WrScr('-------------- Old AeroDyn inputs:')
- write (*,'(A20,I0)') 'WakeMod: ', InputFileData%WakeMod
- write (*,'(A20,I0)') 'SkewMod: ', InputFileData%SkewMod
- write (*,'(A20,I0)') 'AFAeroMod:', InputFileData%AFAeroMod
- write (*,'(A20,L0)') 'FrozenWake:', InputFileData%FrozenWake
- write (*,'(A20,I0)') 'UAMod: ', UAMod_Old
+ write (tmpStr,'(A20,I0)') 'WakeMod: ', InputFileData%WakeMod; call WrScr(trim(tmpStr))
+ write (tmpStr,'(A20,I0)') 'SkewMod: ', InputFileData%SkewMod; call WrScr(trim(tmpStr))
+ write (tmpStr,'(A20,I0)') 'AFAeroMod:', InputFileData%AFAeroMod; call WrScr(trim(tmpStr))
+ write (tmpStr,'(A20,L1)') 'FrozenWake:', InputFileData%FrozenWake; call WrScr(trim(tmpStr))
+ write (tmpStr,'(A20,I0)') 'UAMod: ', UAMod_Old; call WrScr(trim(tmpStr))
call WrScr('------------------------------------------------------')
end subroutine printNewOldInputs
@@ -2408,7 +2411,7 @@ subroutine calcCantAngle(f, xi,stencilSize,n,cantAngle)
call differ_stencil ( xi(i), 1, 2, xiIn, cx, info )
!call differ_stencil ( xi(i), 1, 2, fIn, cf, info )
if (info /= 0) then
- print*,'Cant Calc failed at i=',i
+ call WrScr('Cant Calc failed at i='//trim(Num2LStr(i)))
else
cPrime(i) = 0.0
fPrime(i) = 0.0
diff --git a/modules/aerodyn/src/AeroDyn_Inflow_C_Binding.f90 b/modules/aerodyn/src/AeroDyn_Inflow_C_Binding.f90
index 0920c4ff5e..34fb3c046a 100644
--- a/modules/aerodyn/src/AeroDyn_Inflow_C_Binding.f90
+++ b/modules/aerodyn/src/AeroDyn_Inflow_C_Binding.f90
@@ -27,7 +27,6 @@ MODULE AeroDyn_Inflow_C_BINDING
USE NWTC_Library
USE VersionInfo
-
IMPLICIT NONE
SAVE
@@ -36,7 +35,7 @@ MODULE AeroDyn_Inflow_C_BINDING
PUBLIC :: ADI_C_CalcOutput
PUBLIC :: ADI_C_UpdateStates
PUBLIC :: ADI_C_End
- PUBLIC :: ADI_C_PreInit ! Initial call to setup number of turbines
+ PUBLIC :: ADI_C_PreInit ! Initial call to setup number of turbines
PUBLIC :: ADI_C_SetupRotor ! Initial node positions etc for a rotor
PUBLIC :: ADI_C_SetRotorMotion ! Set motions for a given rotor
PUBLIC :: ADI_C_GetRotorLoads ! Retrieve loads for a given rotor
@@ -45,7 +44,6 @@ MODULE AeroDyn_Inflow_C_BINDING
! Version info for display
type(ProgDesc), parameter :: version = ProgDesc( 'AeroDyn-Inflow library', '', '' )
-
!------------------------------------------------------------------------------------
! Debugging: debugverbose -- passed at PreInit
! 0 - none
@@ -98,7 +96,7 @@ MODULE AeroDyn_Inflow_C_BINDING
type(ADI_InputType) :: ADI_u !< ADI inputs -- set by AD_SetInputMotion. Copied as needed (necessary for correction steps)
!------------------------------
! Simulation data
- type(Dvr_SimData) :: Sim !< data about the simulation
+ type(Dvr_SimData) :: Sim !< Data about the simulation
!------------------------------
! Outputs
type(Dvr_Outputs) :: WrOutputsData !< Data for writing outputs to file
@@ -122,7 +120,7 @@ MODULE AeroDyn_Inflow_C_BINDING
! interface and therefore must store it here analogously to how it is handled
! in the OpenFAST glue code.
integer(IntKi) :: n_Global ! global timestep
- integer(IntKi) :: n_VTK ! VTK timestep
+ integer(IntKi) :: n_VTK ! VTK timestep
real(DbKi) :: InputTimePrev ! input time of last UpdateStates call
! Note that we are including the previous state info here (not done in OF this way)
integer(IntKi), parameter :: STATE_LAST = 0 ! Index for previous state (not needed in OF, but necessary here)
@@ -140,6 +138,37 @@ MODULE AeroDyn_Inflow_C_BINDING
! positions passed into this module and what is used inside AD. This is done
! through a pair of meshes for the motion and loads corresponding to the node
! positions passed in.
+
+ ! ========= BladeNodeToMeshPointMapType =======
+ TYPE, PUBLIC :: BladeNodeToMeshPointMapType
+ INTEGER(IntKi), ALLOCATABLE :: BladeNodeToMeshPoint(:) !< Blade node -> structural mesh point mapping (sized by the number of nodes on the blade)
+ END TYPE BladeNodeToMeshPointMapType
+ ! =======================
+ ! ========= BladePtMeshCoordsType =======
+ TYPE, PUBLIC :: BladePtMeshCoordsType
+ REAL(ReKi), DIMENSION(:,:), ALLOCATABLE :: Position !< Position of all blade points (sized by 3 x number of mesh points on the blade [x,y,z])
+ REAL(ReKi), DIMENSION(:,:,:), ALLOCATABLE :: Orient !< Orientation of all blade points (sized by 3 x 3 x number of mesh points on the blade [r11,r12,r13,r21,r22,r23,r31,r32,r33])
+ REAL(ReKi), DIMENSION(:,:), ALLOCATABLE :: Velocity !< Velocity of all blade points (sized by 6 x number of mesh points on the blade [u,v,w,p,q,r])
+ REAL(ReKi), DIMENSION(:,:), ALLOCATABLE :: Accln !< Acceleration of all blade points (sized by 6 x number of mesh points on the blade [udot,vdot,wdot,pdot,qdot,rdot])
+ REAL(ReKi), DIMENSION(:,:), ALLOCATABLE :: Force !< Force of all blade points (sized by 6 x number of mesh points on the blade [Fx,Fy,Fz,Mx,My,Mz])
+ END TYPE BladePtMeshCoordsType
+ ! =======================
+ ! ========= StrucPtsToBladeMapType =======
+ TYPE, PUBLIC :: StrucPtsToBladeMapType
+ INTEGER(IntKi) :: NumBlades ! Number of blades on this rotor
+ INTEGER(IntKi), ALLOCATABLE :: NumMeshPtsPerBlade(:) ! Number of structural mesh points on each blade (sized by the number of blades)
+ INTEGER(IntKi), ALLOCATABLE :: MeshPt_2_BladeNum(:) ! Structural mesh point -> which blade on the rotor it is on (sized by the number of mesh points on the rotor)
+ TYPE(BladeNodeToMeshPointMapType),ALLOCATABLE:: BladeNode_2_MeshPt(:) ! Blade node on blade -> structural mesh point (sized by the number of mesh points on the blade)
+ TYPE(BladePtMeshCoordsType), ALLOCATABLE :: BladePtMeshCoords(:) ! Mesh point coordinates for each blade (sized by the number of blades)
+ END TYPE StrucPtsToBladeMapType
+ ! =======================
+ ! ========= MeshByBladeType =======
+ TYPE, PUBLIC :: MeshByBladeType
+ ! TODO: Sometime we should rename Mesh to BldMesh
+ TYPE(MeshType), ALLOCATABLE :: Mesh(:) ! Mesh for motions/loads of external nodes at each blade (sized by number of blades on the rotor)
+ END TYPE MeshByBladeType
+ ! =======================
+
!------------------------------
! Meshes for external nodes
! These point meshes are merely used to simplify the mapping of motions/loads
@@ -147,31 +176,25 @@ MODULE AeroDyn_Inflow_C_BINDING
! one or multiple points.
! - 1 point -- rigid floating body assumption
! - N points -- flexible structure (either floating or fixed bottom)
+ ! TODO: for clarity, sometime it might be worth renaming BldPt* here to RtrPt* instead
logical :: TransposeDCM !< Transpose DCMs as passed in -- test the vtk outputs to see if needed
integer(IntKi), allocatable :: NumMeshPts(:) ! Number of mesh points we are interfacing motions/loads to/from AD for each rotor
- type(MeshType), allocatable :: BldPtMotionMesh(:) ! mesh for motions of external nodes
- type(MeshType), allocatable :: BldPtLoadMesh(:) ! mesh for loads for external nodes
- type(MeshType), allocatable :: BldPtLoadMesh_tmp(:) ! mesh for loads for external nodes -- temporary
-! type(MeshType), allocatable :: NacMotionMesh(:) ! mesh for motion of nacelle -- TODO: add this mesh for nacelle load transfers
-! type(MeshType), allocatable :: NacLoadMesh(:) ! mesh for loads for nacelle loads -- TODO: add this mesh for nacelle load transfers
+ type(MeshByBladeType), allocatable :: BldPtMotionMesh(:) ! Mesh for motions of external nodes (sized by number of rotors)
+ type(MeshByBladeType), allocatable :: BldPtLoadMesh(:) ! Mesh for loads for external nodes (sized by number of rotors)
+ type(MeshByBladeType), allocatable :: BldPtLoadMesh_tmp(:) ! Mesh for loads for external nodes -- temporary storage for loads (sized by number of rotors)
+ ! type(MeshType), allocatable :: NacMotionMesh(:) ! mesh for motion of nacelle -- TODO: add this mesh for nacelle load transfers
+ ! type(MeshType), allocatable :: NacLoadMesh(:) ! mesh for loads for nacelle loads -- TODO: add this mesh for nacelle load transfers
!------------------------------
! Mesh mapping: motions
! The mapping of motions from the nodes passed in to the corresponding AD meshes
- type(MeshMapType), allocatable :: Map_BldPtMotion_2_AD_Blade(:,:) ! Mesh mapping between input motion mesh for blade
- type(MeshMapType), allocatable :: Map_AD_Nac_2_NacPtLoad(:) ! Mesh mapping between input motion mesh for nacelle
+ ! TODO: sometime restructure the Map_BldPtMotion_2_AD_Blade and Map_AD_BldLoad_P_2_BldPtLoad to 1D and place inside a rotor structure
+ type(MeshMapType), allocatable :: Map_BldPtMotion_2_AD_Blade(:,:) ! Mesh mapping between input motion mesh for blade (sized by the number of blades and number of rotors)
+ type(MeshMapType), allocatable :: Map_AD_Nac_2_NacPtLoad(:) ! Mesh mapping between input motion mesh for nacelle
!------------------------------
! Mesh mapping: loads
! The mapping of loads from the AD meshes to the corresponding external nodes
- type(MeshMapType), allocatable :: Map_AD_BldLoad_P_2_BldPtLoad(:,:) ! Mesh mapping between AD output blade line2 load to BldPtLoad for return
-! type(MeshMapType) :: Map_NacPtMotion_2_AD_Nac(:) ! Mesh mapping between AD output nacelle pt load to NacLoad for return
- ! Motions input (so we don't have to reallocate all the time
- real(ReKi), allocatable :: tmpBldPtMeshPos(:,:) ! temp array. Probably don't need this, but makes conversion from C clearer.
- real(ReKi), allocatable :: tmpBldPtMeshOri(:,:,:) ! temp array. Probably don't need this, but makes conversion from C clearer.
- real(ReKi), allocatable :: tmpBldPtMeshVel(:,:) ! temp array. Probably don't need this, but makes conversion from C clearer.
- real(ReKi), allocatable :: tmpBldPtMeshAcc(:,:) ! temp array. Probably don't need this, but makes conversion from C clearer.
- real(ReKi), allocatable :: tmpBldPtMeshFrc(:,:) ! temp array. Probably don't need this, but makes conversion to C clearer.
- !------------------------------------------------------------------------------------
-
+ type(StrucPtsToBladeMapType), allocatable :: StrucPts_2_Bld_Map(:) ! Array mapping info for structural mesh points to blades, and back (sized by the number of rotors/turbines)
+ type(MeshMapType), allocatable :: Map_AD_BldLoad_P_2_BldPtLoad(:,:) ! Mesh mapping between AD output blade line2 load to BldPtLoad for return (sized by the number of blades and number of rotors)
! NOTE on turbine origin
! The turbine origin is set by TurbOrigin_C during the ADI_C_SetupRotor routine. This is the tower base location. All
@@ -180,7 +203,6 @@ MODULE AeroDyn_Inflow_C_BINDING
! TurbOrigin_C (stored as Sim%WT(iWT)%OriginInit). When the mesh and other points are passed in, they are relative to
! their respective rotor origin.
-
CONTAINS
!> This routine sets the error status in C_CHAR for export to calling code.
@@ -225,7 +247,7 @@ subroutine ADI_C_PreInit(NumTurbines_C,TransposeDCM_in,debuglevel,ErrStat_C,ErrM
character(ErrMsgLen) :: ErrMsg !< aggregated error message
integer :: ErrStat2 !< temporary error status from a call
character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call
- character(*), parameter :: RoutineName = 'ADI_C_PreInit' !< for error handling
+ character(*), parameter :: RoutineName = 'ADI_C_PreInit' !< for error handling
! Initialize error handling
ErrStat = ErrID_None
@@ -255,7 +277,6 @@ subroutine ADI_C_PreInit(NumTurbines_C,TransposeDCM_in,debuglevel,ErrStat_C,ErrM
if (Failed()) return;
endif
-
! Set number of turbines
Sim%NumTurbines = int(NumTurbines_C,IntKi)
@@ -279,27 +300,31 @@ subroutine ADI_C_PreInit(NumTurbines_C,TransposeDCM_in,debuglevel,ErrStat_C,ErrM
Sim%WT(iWT)%NumBlades = -999
enddo
- ! storage for numnber of meshpoints
+ ! Storage for number of meshpoints
if (allocated(NumMeshPts)) deallocate(NumMeshPts)
allocate(NumMeshPts(Sim%NumTurbines),stat=errStat2); if (Failed0('NumMeshPts')) return
NumMeshPts = -999
- ! allocate meshes and meshmappings
+ ! Allocate meshes and mesh mappings
if (allocated(BldPtMotionMesh )) deallocate(BldPtMotionMesh )
if (allocated(BldPtLoadMesh )) deallocate(BldPtLoadMesh )
if (allocated(BldPtLoadMesh_tmp)) deallocate(BldPtLoadMesh_tmp)
- !if (allocated(NacMotionMesh )) deallocate(NacMotionMesh )
- !if (allocated(NacLoadMesh )) deallocate(NacLoadMesh )
- if (allocated(Map_BldPtMotion_2_AD_Blade )) deallocate(Map_BldPtMotion_2_AD_Blade )
- if (allocated(Map_AD_BldLoad_P_2_BldPtLoad)) deallocate(Map_AD_BldLoad_P_2_BldPtLoad)
- !if (allocated(Map_NacPtMotion_2_AD_Nac )) deallocate(Map_NacPtMotion_2_AD_Nac )
+ ! if (allocated(NacMotionMesh )) deallocate(NacMotionMesh )
+ ! if (allocated(NacLoadMesh )) deallocate(NacLoadMesh )
allocate(BldPtMotionMesh( Sim%NumTurbines), STAT=ErrStat2); if (Failed0('BldPtMotionMesh' )) return
allocate(BldPtLoadMesh( Sim%NumTurbines), STAT=ErrStat2); if (Failed0('BldPtLoadMesh' )) return
allocate(BldPtLoadMesh_tmp(Sim%NumTurbines), STAT=ErrStat2); if (Failed0('BldPtLoadMesh_tmp')) return
- !allocate(NacMotionMesh( Sim%NumTurbines), STAT=ErrStat2); if (Failed0('NacMotionMesh' )) return
- !allocate(NacLoadMesh( Sim%NumTurbines), STAT=ErrStat2); if (Failed0('NacLoadMesh' )) return
- !allocate(Map_NacPtMotion_2_AD_Nac(Sim%NumTurbines),STAT=ErrStat2); if (Failed0('Map_AD_BldLoad_P_2_BldPtLoad')) return
+ ! allocate(NacMotionMesh( Sim%NumTurbines), STAT=ErrStat2); if (Failed0('NacMotionMesh' )) return
+ ! allocate(NacLoadMesh( Sim%NumTurbines), STAT=ErrStat2); if (Failed0('NacLoadMesh' )) return
+ if (allocated(Map_BldPtMotion_2_AD_Blade )) deallocate(Map_BldPtMotion_2_AD_Blade )
+ if (allocated(Map_AD_BldLoad_P_2_BldPtLoad )) deallocate(Map_AD_BldLoad_P_2_BldPtLoad)
+ ! if (allocated(Map_NacPtMotion_2_AD_Nac )) deallocate(Map_NacPtMotion_2_AD_Nac )
+ ! allocate(Map_NacPtMotion_2_AD_Nac(Sim%NumTurbines),STAT=ErrStat2); if (Failed0('Map_AD_BldLoad_P_2_BldPtLoad')) returns
+
+ ! Allocate the StrucPtsToBladeMapType array used for mapping structural points to blades of the rotor
+ if (allocated(StrucPts_2_Bld_Map)) deallocate(StrucPts_2_Bld_Map)
+ allocate(StrucPts_2_Bld_Map(Sim%NumTurbines), STAT=ErrStat2); if (Failed0('StrucPts_2_Bld_Map' )) return
call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C)
@@ -341,7 +366,7 @@ subroutine ShowPassedData()
call WrScr("-----------------------------------------------------------")
end subroutine ShowPassedData
-end subroutine ADI_C_PreInit
+end subroutine ADI_C_PreInit
!===============================================================================================================
!--------------------------------------------- AeroDyn Init----------------------------------------------------
@@ -391,7 +416,7 @@ SUBROUTINE ADI_C_Init( ADinputFilePassed, ADinputFileString_C, ADinputFileString
real(c_float), intent(in ) :: VTKNacDim_in(6) !< Nacelle dimension passed in for VTK surface rendering [0,y0,z0,Lx,Ly,Lz] (m)
real(c_float), intent(in ) :: VTKHubrad_in !< Hub radius for VTK surface rendering
integer(c_int), intent(in ) :: wrOuts_C !< Write ADI output file
- real(c_double), intent(in ) :: DT_Outs_C !< Timestep to write output file from ADI
+ real(c_double), intent(in ) :: DT_Outs_C !< Timestep to write output file from ADI
! Output
integer(c_int), intent( out) :: NumChannels_C !< Number of output channels requested from the input file
character(kind=c_char), intent( out) :: OutputChannelNames_C(ChanLen*MaxADIOutputs+1) !< NOTE: if MaxADIOutputs is sufficiently large, we may overrun the buffer on the Python side.
@@ -399,7 +424,7 @@ SUBROUTINE ADI_C_Init( ADinputFilePassed, ADinputFileString_C, ADinputFileString
integer(c_int), intent( out) :: ErrStat_C !< Error status
character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) !< Error message (C_NULL_CHAR terminated)
- ! Local Variable4
+ ! Local variables
character(IntfStrLen) :: OutRootName !< Root name to use for echo files and other
character(IntfStrLen) :: TmpFileName !< Temporary file name if not passing AD or IfW input file contents directly
character(kind=C_char, len=ADinputFileStringLength_C), pointer :: ADinputFileString !< Input file as a single string with NULL chracter separating lines
@@ -409,7 +434,7 @@ SUBROUTINE ADI_C_Init( ADinputFilePassed, ADinputFileString_C, ADinputFileString
character(ErrMsgLen) :: ErrMsg !< aggregated error message
integer(IntKi) :: ErrStat2 !< temporary error status from a call
character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call
- integer(IntKi) :: i,j,k !< generic counters
+ integer(IntKi) :: i,j,k !< generic index variables
integer(IntKi) :: iWT !< current turbine number (iterate through during setup for ADI_Init call)
integer(IntKi) :: AeroProjMod !< for checking that all turbines use the same AeroProjMod
character(*), parameter :: RoutineName = 'ADI_C_Init' !< for error handling
@@ -787,9 +812,9 @@ subroutine SetupFileOutputs()
call concatOutputHeaders(WrOutputsData%WriteOutputHdr, WrOutputsData%WriteOutputUnt, InitOutData%WriteOutputHdr, InitOutData%WriteOutputUnt, errStat2, errMsg2); if(Failed()) return
! allocate output file handling and set formats
- WrOutputsData%outFmt = "ES15.8E2"
+ WrOutputsData%outFmt = "ES15.8E2"
WrOutputsData%delim = TAB
- WrOutputsData%AD_ver = InitOutData%Ver
+ WrOutputsData%AD_ver = InitOutData%Ver
allocate(WrOutputsData%unOutFile(Sim%numTurbines), STAT=ErrStat2); if(Failed0("unOutFile")) return;
WrOutputsData%unOutFile = -1
!FIXME: number of timesteps is incorrect!
@@ -847,82 +872,61 @@ end subroutine ShowPassedData
!! meshes
subroutine SetupMotionLoadsInterfaceMeshes()
integer(IntKi) :: iWT !< current rotor/turbine
- integer(IntKi) :: iNode
- integer(IntKi) :: maxBlades !< find out maximum number of blades on all turbine rotors
+ integer(IntKi) :: iBlade !< current blade
+ integer(IntKi) :: maxBlades !< maximum number of blades on all turbine rotors
+ ! Find out maximum number of blades on all turbine rotors
maxBlades = 0_IntKi
do iWT=1,Sim%NumTurbines
maxBlades = max(maxBlades,Sim%WT(iWT)%NumBlades)
enddo
- allocate(Map_BldPtMotion_2_AD_Blade( maxBlades,Sim%NumTurbines),STAT=ErrStat2); if (Failed0('Map_BldPtMotion_2_AD_Blade' )) return
- allocate(Map_AD_BldLoad_P_2_BldPtLoad(maxBlades,Sim%NumTurbines),STAT=ErrStat2); if (Failed0('Map_AD_BldLoad_P_2_BldPtLoad')) return
+ ! NOTE: storing mappings in 2D this way may increase memory usage slightly if one turbine has many more blades than another. However
+ ! the speed an memory penalties are negligible, so I don't see much reason to change that at this point.
+ allocate(Map_BldPtMotion_2_AD_Blade( maxBlades, Sim%NumTurbines), STAT=ErrStat2); if (Failed0('Map_BldPtMotion_2_AD_Blade' )) return
+ allocate(Map_AD_BldLoad_P_2_BldPtLoad(maxBlades, Sim%NumTurbines), STAT=ErrStat2); if (Failed0('Map_AD_BldLoad_P_2_BldPtLoad')) return
- ! step through all turbine rotors
+ ! Step through all turbine rotors
do iWT=1,Sim%NumTurbines
-
!-------------------------------------------------------------
! Load mesh for blades
- CALL MeshCopy( SrcMesh = BldPtMotionMesh(iWT) ,&
- DestMesh = BldPtLoadMesh(iWT) ,&
- CtrlCode = MESH_SIBLING ,&
- IOS = COMPONENT_OUTPUT ,&
- ErrStat = ErrStat2 ,&
- ErrMess = ErrMsg2 ,&
- Force = .TRUE. ,&
- Moment = .TRUE. )
- if(Failed()) return
- BldPtLoadMesh(iWT)%RemapFlag = .FALSE.
-
- ! Temp mesh for load transfer
- CALL MeshCopy( SrcMesh = BldPtLoadMesh(iWT) ,&
- DestMesh = BldPtLoadMesh_tmp(iWT),&
- CtrlCode = MESH_COUSIN ,&
- IOS = COMPONENT_OUTPUT ,&
- ErrStat = ErrStat2 ,&
- ErrMess = ErrMsg2 ,&
- Force = .TRUE. ,&
- Moment = .TRUE. )
- if(Failed()) return
- BldPtLoadMesh_tmp(iWT)%RemapFlag = .FALSE.
-
-
- ! For checking the mesh
- ! note: CU is is output unit (platform dependent).
- if (debugverbose >= 4) call MeshPrintInfo( CU, BldPtLoadMesh(iWT), MeshName='BldPtLoadMesh'//trim(Num2LStr(iWT)) )
-
-
-! !-------------------------------------------------------------
-! ! Load mesh for nacelle
-! CALL MeshCopy( SrcMesh = NacMotionMesh(iWT) ,&
-! DestMesh = NacLoadMesh(iWT) ,&
-! CtrlCode = MESH_SIBLING ,&
-! IOS = COMPONENT_OUTPUT ,&
-! ErrStat = ErrStat2 ,&
-! ErrMess = ErrMsg2 ,&
-! Force = .TRUE. ,&
-! Moment = .TRUE. )
-! if(Failed()) return
-! NacLoadMesh(iWT)%RemapFlag = .FALSE.
-!
-! ! For checking the mesh, uncomment this.
-! ! note: CU is is output unit (platform dependent).
-! if (debugverbose >= 4) call MeshPrintInfo( CU, NacLoadMesh(iWT), MeshName='NacLoadMesh'//trim(Num2LStr(iWT)) )
-
-
- !-------------------------------------------------------------
- ! Set the mapping meshes
- ! blades
- do i=1,Sim%WT(iWT)%NumBlades
- call MeshMapCreate( BldPtMotionMesh(iWT), ADI%u(1)%AD%rotors(iWT)%BladeMotion(i), Map_BldPtMotion_2_AD_Blade(i,iWT), ErrStat2, ErrMsg2 ); if(Failed()) return
- call MeshMapCreate( ADI%y%AD%rotors(iWT)%BladeLoad(i), BldPtLoadMesh(iWT), Map_AD_BldLoad_P_2_BldPtLoad(i,iWT), ErrStat2, ErrMsg2 ); if(Failed()) return
- enddo
- ! nacelle -- TODO: add this mesh for nacelle load transfers
-! if ( ADI%y%AD%rotors(iWT)%NacelleLoad%Committed ) then
-! call MeshMapCreate( NacMotionMesh(iWT), ADI%u(1)%AD%rotors(iWT)%NacelleMotion, Map_NacPtMotion_2_AD_Nac(iWT), ErrStat2, ErrMsg2 ); if(Failed()) return
-! call MeshMapCreate( ADI%y%AD%rotors(iWT)%NacelleLoad, NacLoadMesh(iWT), Map_AD_Nac_2_NacPtLoad(iWT), ErrStat2, ErrMsg2 ); if(Failed()) return
-! endif
-
+ ! Step through all blades on this rotor
+ do iBlade=1,Sim%WT(iWT)%NumBlades
+ !-------------------------------------------------------------
+ ! Load mesh for blades
+ CALL MeshCopy( SrcMesh = BldPtMotionMesh(iWT)%Mesh(iBlade) ,&
+ DestMesh = BldPtLoadMesh(iWT)%Mesh(iBlade) ,&
+ CtrlCode = MESH_SIBLING ,&
+ IOS = COMPONENT_OUTPUT ,&
+ ErrStat = ErrStat2 ,&
+ ErrMess = ErrMsg2 ,&
+ Force = .TRUE. ,&
+ Moment = .TRUE. )
+ if(Failed()) return
+ BldPtMotionMesh(iWT)%Mesh(iBlade)%RemapFlag = .FALSE.
+
+ ! Temp mesh for load transfer
+ CALL MeshCopy( SrcMesh = BldPtLoadMesh(iWT)%Mesh(iBlade) ,&
+ DestMesh = BldPtLoadMesh_tmp(iWT)%Mesh(iBlade) ,&
+ CtrlCode = MESH_COUSIN ,&
+ IOS = COMPONENT_OUTPUT ,&
+ ErrStat = ErrStat2 ,&
+ ErrMess = ErrMsg2 ,&
+ Force = .TRUE. ,&
+ Moment = .TRUE. )
+ if(Failed()) return
+ BldPtLoadMesh_tmp(iWT)%Mesh(iBlade)%RemapFlag = .FALSE.
+
+ ! For checking the mesh
+ ! Note: CU is is output unit (platform dependent).
+ if (debugverbose >= 4) call MeshPrintInfo( CU, BldPtLoadMesh(iWT)%Mesh(iBlade), MeshName='BldPtLoadMesh'//trim(Num2LStr(iWT))//'_'//trim(Num2LStr(iBlade)) )
+
+ !-------------------------------------------------------------
+ ! Set the mapping meshes
+ ! blades
+ call MeshMapCreate( BldPtMotionMesh(iWT)%Mesh(iBlade), ADI%u(1)%AD%rotors(iWT)%BladeMotion(iBlade), Map_BldPtMotion_2_AD_Blade(iBlade, iWT), ErrStat2, ErrMsg2 ); if(Failed()) return
+ call MeshMapCreate( ADI%y%AD%rotors(iWT)%BladeLoad(iBlade), BldPtLoadMesh(iWT)%Mesh(iBlade), Map_AD_BldLoad_P_2_BldPtLoad(iBlade, iWT), ErrStat2, ErrMsg2 ); if(Failed()) return
+ enddo ! iBlade
enddo ! iWT
end subroutine SetupMotionLoadsInterfaceMeshes
@@ -1291,7 +1295,7 @@ subroutine ADI_C_SetupRotor(iWT_c, TurbineIsHAWT_c, TurbOrigin_C, &
NacPos_C, NacOri_C, &
NumBlades_C, BldRootPos_C, BldRootOri_C, &
NumMeshPts_C, InitMeshPos_C, InitMeshOri_C, &
- ErrStat_C, ErrMsg_C) BIND (C, NAME='ADI_C_SetupRotor')
+ MeshPtToBladeNum_C, ErrStat_C, ErrMsg_C) BIND (C, NAME='ADI_C_SetupRotor')
implicit none
#ifndef IMPLICIT_DLLEXPORT
!DEC$ ATTRIBUTES DLLEXPORT :: ADI_C_SetupRotor
@@ -1308,21 +1312,23 @@ subroutine ADI_C_SetupRotor(iWT_c, TurbineIsHAWT_c, TurbOrigin_C, &
integer(c_int), intent(in ) :: NumBlades_C !< Number of blades
real(c_float), intent(in ) :: BldRootPos_C( 3*NumBlades_C ) !< Blade root positions
real(c_double), intent(in ) :: BldRootOri_C( 9*NumBlades_C ) !< Blade root orientations
- ! Initial nodes
- integer(c_int), intent(in ) :: NumMeshPts_C !< Number of mesh points we are transfering motions to and output loads to
+ ! Initial nodes
+ integer(c_int), intent(in ) :: NumMeshPts_C !< Number of mesh points we are transferring motions and outputting loads to
real(c_float), intent(in ) :: InitMeshPos_C( 3*NumMeshPts_C ) !< A 3xNumMeshPts_C array [x,y,z]
real(c_double), intent(in ) :: InitMeshOri_C( 9*NumMeshPts_C ) !< A 9xNumMeshPts_C array [r11,r12,r13,r21,r22,r23,r31,r32,r33]
+ integer(c_int), intent(in ) :: MeshPtToBladeNum_C( NumMeshPts_C ) !< A NumMeshPts_C array of blade numbers associated with each mesh point
integer(c_int), intent( out) :: ErrStat_C !< Error status
character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) !< Error message (C_NULL_CHAR terminated)
- ! local vars
- integer(IntKi) :: iWT !< current turbine
+ ! Local variables
+ integer(IntKi) :: iWT !< current turbine
+ integer(IntKi) :: iBlade !< current blade
logical :: TurbineIsHAWT !< true for HAWT, false for VAWT
integer(IntKi) :: ErrStat !< aggregated error messagNumBlades_ee
character(ErrMsgLen) :: ErrMsg !< aggregated error message
integer(IntKi) :: ErrStat2 !< temporary error status from a call
character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call
- integer(IntKi) :: i,j,k !< generic counters
+ integer(IntKi) :: i,j,k !< generic index variables
character(*), parameter :: RoutineName = 'ADI_C_Init' !< for error handling
! Initialize error handling
@@ -1375,7 +1381,7 @@ subroutine ADI_C_SetupRotor(iWT_c, TurbineIsHAWT_c, TurbOrigin_C, &
enddo
endif
- ! Remap the orientation DCM just in case there is some issue with passed
+ ! Remap the orientation DCM just in case there is some issue with passed
call OrientRemap(InitInp%AD%rotors(iWT)%HubOrientation)
call OrientRemap(InitInp%AD%rotors(iWT)%NacelleOrientation)
do i=1,Sim%WT(iWT)%NumBlades
@@ -1472,56 +1478,112 @@ end subroutine ShowPassedData
subroutine SetupMotionMesh()
real(ReKi) :: InitPos(3)
real(R8Ki) :: Orient(3,3)
- integer(IntKi) :: iNode
- ! this may not be super efficient since we are allocating for every turbine, but this makes things a bit simpler
- if ( allocated(tmpBldPtMeshPos) ) deallocate(tmpBldPtMeshPos)
- if ( allocated(tmpBldPtMeshOri) ) deallocate(tmpBldPtMeshOri)
- call AllocAry( tmpBldPtMeshPos, 3, NumMeshPts(iWT), "tmpBldPtMeshPos", ErrStat2, ErrMsg2 ); if (Failed()) return
- call AllocAry( tmpBldPtMeshOri, 3, 3, NumMeshPts(iWT), "tmpBldPtMeshOri", ErrStat2, ErrMsg2 ); if (Failed()) return
+ integer(IntKi) :: count
- ! Reshape mesh position, orientation, velocity, acceleration
- tmpBldPtMeshPos(1:3,1:NumMeshPts(iWT)) = reshape( real(InitMeshPos_C(1:3*NumMeshPts(iWT)),ReKi), (/3, NumMeshPts(iWT)/) )
- tmpBldPtMeshOri(1:3,1:3,1:NumMeshPts(iWT)) = reshape( real(InitMeshOri_C(1:9*NumMeshPts(iWT)),R8Ki), (/3,3,NumMeshPts(iWT)/) )
+ !-------------------------------------------------------------
+ ! Allocate and define the components of StrucPts_2_Bld_Map
+ !-------------------------------------------------------------
+ StrucPts_2_Bld_Map(iWT)%NumBlades = Sim%WT(iWT)%NumBlades
+
+ call AllocAry(StrucPts_2_Bld_Map(iWT)%NumMeshPtsPerBlade, Sim%WT(iWT)%NumBlades, "NumMeshPtsPerBlade", ErrStat2, ErrMsg2 ); if (Failed()) return
+ call AllocAry( StrucPts_2_Bld_Map(iWT)%MeshPt_2_BladeNum, NumMeshPts(iWT), "MeshPt_2_BladeNum", ErrStat2, ErrMsg2 ); if (Failed()) return
+
+ allocate(StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt( Sim%WT(iWT)%NumBlades ), STAT=ErrStat2); if (Failed0('StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt' )) return
+
+ ! Calculate the number of mesh points per blade
+ do i=1,Sim%WT(iWT)%NumBlades
+ count = 0
+ do j=1,NumMeshPts(iWT)
+ if (MeshPtToBladeNum_C(j) == i) then
+ count = count + 1
+ endif
+ enddo
+ StrucPts_2_Bld_Map(iWT)%NumMeshPtsPerBlade(i) = count
+ enddo
+
+ StrucPts_2_Bld_Map(iWT)%MeshPt_2_BladeNum(1:NumMeshPts(iWT)) = MeshPtToBladeNum_C(1:NumMeshPts(iWT))
+
+ ! Allocate remaining components of StrucPts_2_Bld_Map based on the number of mesh points per blade
+ do i=1,Sim%WT(iWT)%NumBlades
+ call AllocAry(StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint, StrucPts_2_Bld_Map(iWT)%NumMeshPtsPerBlade(i), "BladeNodeToMeshPoint", ErrStat2, ErrMsg2); if (Failed()) return
+ enddo
+
+ do i=1,Sim%WT(iWT)%NumBlades
+ count = 0
+ do j=1,NumMeshPts(iWT)
+ if (MeshPtToBladeNum_C(j) == i) then
+ count = count + 1
+ StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint(count) = j
+ endif
+ enddo
+ enddo
+
+ ! Allocate and define the components of BladePtMeshCoords
+ allocate(StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(Sim%WT(iWT)%NumBlades), STAT=ErrStat2); if (Failed0('StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords')) return
+ do i=1,Sim%WT(iWT)%NumBlades
+ call AllocAry(StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(i)%Position, 3, StrucPts_2_Bld_Map(iWT)%NumMeshPtsPerBlade(i), "BladePtMeshCoords(i)%Position", ErrStat2, ErrMsg2 ); if (Failed()) return
+ call AllocAry(StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(i)%Orient, 3, 3, StrucPts_2_Bld_Map(iWT)%NumMeshPtsPerBlade(i), "BladePtMeshCoords(i)%Orient", ErrStat2, ErrMsg2 ); if (Failed()) return
+ call AllocAry(StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(i)%Velocity, 6, StrucPts_2_Bld_Map(iWT)%NumMeshPtsPerBlade(i), "BladePtMeshCoords(i)%Velocity", ErrStat2, ErrMsg2 ); if (Failed()) return
+ call AllocAry(StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(i)%Accln, 6, StrucPts_2_Bld_Map(iWT)%NumMeshPtsPerBlade(i), "BladePtMeshCoords(i)%Accln", ErrStat2, ErrMsg2 ); if (Failed()) return
+ call AllocAry(StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(i)%Force, 6, StrucPts_2_Bld_Map(iWT)%NumMeshPtsPerBlade(i), "BladePtMeshCoords(i)%Force", ErrStat2, ErrMsg2 ); if (Failed()) return
+ enddo
+
+ do i=1,Sim%WT(iWT)%NumBlades
+ do j=1,StrucPts_2_Bld_Map(iWT)%NumMeshPtsPerBlade(i)
+ StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(i)%Position(1:3,j) = reshape( real(InitMeshPos_C(3 * StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint(j) - 2 : 3 * StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint(j)),ReKi), (/3/) )
+ StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(i)%Orient(1:3,1:3,j) = reshape( real(InitMeshOri_C(9 * StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint(j) - 8 : 9 * StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint(j)),R8Ki), (/3,3/) )
+ enddo
+ enddo
+
+ ! Allocate the meshes
+ allocate(BldPtMotionMesh(iWT)%Mesh( Sim%WT(iWT)%NumBlades ), STAT=ErrStat2); if (Failed0('BldPtMotionMesh( iWT )%Mesh' )) return
+ allocate(BldPtLoadMesh(iWT)%Mesh( Sim%WT(iWT)%NumBlades ), STAT=ErrStat2); if (Failed0('BldPtLoadMesh( iWT )%Mesh' )) return
+ allocate(BldPtLoadMesh_tmp(iWT)%Mesh( Sim%WT(iWT)%NumBlades ), STAT=ErrStat2); if (Failed0('BldPtLoadMesh_tmp( iWT )%Mesh' )) return
!-------------------------------------------------------------
! Set the interface meshes for motion inputs and loads output
!-------------------------------------------------------------
! Motion mesh for blades
- call MeshCreate( BldPtMotionMesh(iWT) , &
- IOS = COMPONENT_INPUT , &
- Nnodes = NumMeshPts(iWT) , &
- ErrStat = ErrStat2 , &
- ErrMess = ErrMsg2 , &
- TranslationDisp = .TRUE., Orientation = .TRUE., &
- TranslationVel = .TRUE., RotationVel = .TRUE., &
- TranslationAcc = .TRUE., RotationAcc = .FALSE. )
- if(Failed()) return
-
- do iNode=1,NumMeshPts(iWT)
- ! initial position and orientation of node
- InitPos = tmpBldPtMeshPos(1:3,iNode) + Sim%WT(iWT)%OriginInit(1:3)
- if (TransposeDCM) then
- Orient = transpose(tmpBldPtMeshOri(1:3,1:3,iNode))
- else
- Orient = tmpBldPtMeshOri(1:3,1:3,iNode)
- endif
- call OrientRemap(Orient)
- call MeshPositionNode( BldPtMotionMesh(iWT) , &
- iNode , &
- InitPos , & ! position
- ErrStat2, ErrMsg2 , &
- Orient ) ! orientation
+ do iBlade=1,Sim%WT(iWT)%NumBlades
+ call MeshCreate( BldPtMotionMesh(iWT)%Mesh(iBlade) , &
+ IOS = COMPONENT_INPUT , &
+ Nnodes = StrucPts_2_Bld_Map(iWT)%NumMeshPtsPerBlade(iBlade) , &
+ ErrStat = ErrStat2 , &
+ ErrMess = ErrMsg2 , &
+ TranslationDisp = .TRUE., Orientation = .TRUE. , &
+ TranslationVel = .TRUE., RotationVel = .TRUE. , &
+ TranslationAcc = .TRUE., RotationAcc = .FALSE. )
if(Failed()) return
- call MeshConstructElement ( BldPtMotionMesh(iWT), ELEMENT_POINT, ErrStat2, ErrMsg2, iNode ); if(Failed()) return
enddo
- call MeshCommit ( BldPtMotionMesh(iWT), ErrStat2, ErrMsg2 ); if(Failed()) return
- BldPtMotionMesh(iWT)%RemapFlag = .FALSE.
+ do iBlade=1,Sim%WT(iWT)%NumBlades
+ do j=1,StrucPts_2_Bld_Map(iWT)%NumMeshPtsPerBlade(iBlade)
+ ! Initial position and orientation of node
+ InitPos = StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(iBlade)%Position(1:3,j) + Sim%WT(iWT)%OriginInit(1:3)
+ if (TransposeDCM) then
+ Orient = transpose(StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(iBlade)%Orient(1:3,1:3,j))
+ else
+ Orient = StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(iBlade)%Orient(1:3,1:3,j)
+ endif
+ call OrientRemap(Orient)
+ call MeshPositionNode( BldPtMotionMesh(iWT)%Mesh(iBlade) , &
+ j , &
+ InitPos , & ! position
+ ErrStat2, ErrMsg2 , &
+ Orient ) ! orientation
+ if(Failed()) return
+ call MeshConstructElement ( BldPtMotionMesh(iWT)%Mesh(iBlade), ELEMENT_POINT, ErrStat2, ErrMsg2, j ); if(Failed()) return
+ enddo
+ enddo
- ! For checking the mesh, uncomment this.
- ! note: CU is is output unit (platform dependent).
- if (debugverbose >= 4) call MeshPrintInfo( CU, BldPtMotionMesh(iWT), MeshName='BldPtMotionMesh'//trim(Num2LStr(iWT)) )
+ do iBlade=1,Sim%WT(iWT)%NumBlades
+ call MeshCommit ( BldPtMotionMesh(iWT)%Mesh(iBlade), ErrStat2, ErrMsg2 ); if(Failed()) return
+ BldPtMotionMesh(iWT)%Mesh(iBlade)%RemapFlag = .FALSE.
+ ! For checking the mesh
+ ! Note: CU is is output unit (platform dependent)
+ if (debugverbose >= 4) call MeshPrintInfo( CU, BldPtMotionMesh(iWT)%Mesh(iBlade), MeshName='BldPtMotionMesh'//trim(Num2LStr(iWT))//'_'//trim(Num2LStr(iBlade)) )
+ enddo
! !-------------------------------------------------------------
! ! Motion mesh for nacelle -- TODO: add this mesh for nacelle load transfers
@@ -1554,11 +1616,8 @@ subroutine SetupMotionMesh()
! ! note: CU is is output unit (platform dependent).
! if (debugverbose >= 4) call MeshPrintInfo( CU, NacMotionMesh(iWT), MeshName='NacMotionMesh'//trim(Num2LStr(iWT)) )
- ! clear the tmp stuff (will need to resize this later)
- if ( allocated(tmpBldPtMeshPos) ) deallocate(tmpBldPtMeshPos)
- if ( allocated(tmpBldPtMeshOri) ) deallocate(tmpBldPtMeshOri)
end subroutine SetupMotionMesh
-end subroutine ADI_C_SetupRotor
+end subroutine ADI_C_SetupRotor
!===============================================================================================================
!--------------------------------------------- AeroDyn SetRotorMotion ------------------------------------------
@@ -1576,7 +1635,7 @@ subroutine ADI_C_SetRotorMotion( iWT_c, &
!DEC$ ATTRIBUTES DLLEXPORT :: ADI_C_SetRotorMotion
!GCC$ ATTRIBUTES DLLEXPORT :: ADI_C_SetRotorMotion
#endif
- integer(c_int), intent(in ) :: iWT_c !< Wind turbine / rotor number
+ integer(c_int), intent(in ) :: iWT_c !< Wind turbine / rotor number
real(c_float), intent(in ) :: HubPos_C( 3 ) !< Hub position
real(c_double), intent(in ) :: HubOri_C( 9 ) !< Hub orientation
real(c_float), intent(in ) :: HubVel_C( 6 ) !< Hub velocity
@@ -1600,8 +1659,8 @@ subroutine ADI_C_SetRotorMotion( iWT_c, &
! Local variables
real(DbKi) :: Time
- integer(IntKi) :: iNode
integer(IntKi) :: iWT !< current wind turbine / rotor
+ integer(IntKi) :: i,j !< generic index variables
integer(IntKi) :: ErrStat !< aggregated error status
character(ErrMsgLen) :: ErrMsg !< aggregated error message
integer(IntKi) :: ErrStat2 !< temporary error status from a call
@@ -1628,11 +1687,14 @@ subroutine ADI_C_SetRotorMotion( iWT_c, &
endif
! Reshape mesh position, orientation, velocity, acceleration
- tmpBldPtMeshPos(1:3,1:NumMeshPts(iWT)) = reshape( real(MeshPos_C(1:3*NumMeshPts(iWT)),ReKi), (/3, NumMeshPts(iWT)/) )
- tmpBldPtMeshOri(1:3,1:3,1:NumMeshPts(iWT)) = reshape( real(MeshOri_C(1:9*NumMeshPts(iWT)),R8Ki), (/3,3,NumMeshPts(iWT)/) )
- tmpBldPtMeshVel(1:6,1:NumMeshPts(iWT)) = reshape( real(MeshVel_C(1:6*NumMeshPts(iWT)),ReKi), (/6, NumMeshPts(iWT)/) )
- tmpBldPtMeshAcc(1:6,1:NumMeshPts(iWT)) = reshape( real(MeshAcc_C(1:6*NumMeshPts(iWT)),ReKi), (/6, NumMeshPts(iWT)/) )
-
+ do i=1,Sim%WT(iWT)%NumBlades
+ do j=1,StrucPts_2_Bld_Map(iWT)%NumMeshPtsPerBlade(i)
+ StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(i)%Position( 1:3,j) = reshape( real(MeshPos_C(3 * StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint(j) - 2 : 3 * StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint(j)),ReKi), (/3/) )
+ StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(i)%Orient(1:3,1:3,j) = reshape( real(MeshOri_C(9 * StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint(j) - 8 : 9 * StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint(j)),R8Ki), (/3,3/) )
+ StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(i)%Velocity( 1:6,j) = reshape( real(MeshVel_C(6 * StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint(j) - 5 : 6 * StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint(j)),ReKi), (/6/) )
+ StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(i)%Accln( 1:6,j) = reshape( real(MeshAcc_C(6 * StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint(j) - 5 : 6 * StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint(j)),ReKi), (/6/) )
+ enddo
+ enddo
! Transfer motions to input meshes
do iWT=1,Sim%NumTurbines
@@ -1738,7 +1800,7 @@ end subroutine ADI_C_SetRotorMotion
!--------------------------------------------- AeroDyn GetRotorLoads -------------------------------------------
!===============================================================================================================
!> Get the loads from a single rotor. This must be called after ADI_C_CalcOutput
-subroutine ADI_C_GetRotorLoads(iWT_C, &
+subroutine ADI_C_GetRotorLoads(iWT_C, &
NumMeshPts_C, MeshFrc_C, &
ErrStat_C, ErrMsg_C) BIND (C, NAME='ADI_C_GetRotorLoads')
implicit none
@@ -1746,7 +1808,7 @@ subroutine ADI_C_GetRotorLoads(iWT_C, &
!DEC$ ATTRIBUTES DLLEXPORT :: ADI_C_GetRotorLoads
!GCC$ ATTRIBUTES DLLEXPORT :: ADI_C_GetRotorLoads
#endif
- integer(c_int), intent(in ) :: iWT_C !< Wind turbine / rotor number
+ integer(c_int), intent(in ) :: iWT_C !< Wind turbine / rotor number
integer(c_int), intent(in ) :: NumMeshPts_C !< Number of mesh points we are transfering motions to and output loads to
real(c_float), intent( out) :: MeshFrc_C( 6*NumMeshPts_C ) !< A 6xNumMeshPts_C array [Fx,Fy,Fz,Mx,My,Mz] -- forces and moments (global)
integer(c_int), intent( out) :: ErrStat_C
@@ -1759,6 +1821,7 @@ subroutine ADI_C_GetRotorLoads(iWT_C, &
integer(IntKi) :: ErrStat2 !< temporary error status from a call
character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call
character(*), parameter :: RoutineName = 'ADI_C_SetRotorMotion' !< for error handling
+ integer(IntKi) :: i,j !< generic index variables
! Initialize error handling
ErrStat = ErrID_None
@@ -1785,7 +1848,11 @@ subroutine ADI_C_GetRotorLoads(iWT_C, &
! Set output force/moment array
call Set_OutputLoadArray(iWT)
- MeshFrc_C(1:6*NumMeshPts(iWT)) = reshape( real(tmpBldPtMeshFrc(1:6,1:NumMeshPts(iWT)), c_float), (/6*NumMeshPts(iWT)/) )
+ do i=1,Sim%WT(iWT)%NumBlades
+ do j=1,StrucPts_2_Bld_Map(iWT)%NumMeshPtsPerBlade(i)
+ MeshFrc_C(6 * StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint(j) - 5 : 6 * StrucPts_2_Bld_Map(iWT)%BladeNode_2_MeshPt(i)%BladeNodeToMeshPoint(j)) = real(StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(i)%Force(1:6,j), c_float)
+ enddo
+ enddo
! Set error status
call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C)
@@ -1809,8 +1876,7 @@ subroutine ShowPassedData()
call WrScr(" NumMeshPts_C "//trim(Num2LStr( NumMeshPts_C )) )
call WrScr("-----------------------------------------------------------")
end subroutine ShowPassedData
-end subroutine ADI_C_GetRotorLoads
-
+end subroutine ADI_C_GetRotorLoads
!===================================================================================================================================
@@ -1820,24 +1886,27 @@ end subroutine ADI_C_GetRotorLoads
!> This routine is operating on module level data. Error handling here in case checks added
!! NOTE: the OriginInit is not included in the data passed in and must be added to the the position info here
subroutine Set_MotionMesh(iWT, ErrStat3, ErrMsg3)
- integer(IntKi), intent(in ) :: iWT !< current rotor/turbine
+ integer(IntKi), intent(in ) :: iWT !< current rotor/turbine
integer(IntKi), intent( out) :: ErrStat3
character(ErrMsgLen), intent( out) :: ErrMsg3
- integer(IntKi) :: iNode
+ integer(IntKi) :: iBlade !< current blade
+ integer(IntKi) :: j !< generic index variables
+
ErrStat3 = 0_IntKi
ErrMsg3 = ''
! Set mesh corresponding to input motions
- do iNode=1,NumMeshPts(iWT)
- BldPtMotionMesh(iWT)%TranslationDisp(1:3,iNode) = tmpBldPtMeshPos(1:3,iNode) + Sim%WT(iWT)%OriginInit(1:3) - real(BldPtMotionMesh(iWT)%Position(1:3,iNode), R8Ki)
- BldPtMotionMesh(iWT)%Orientation(1:3,1:3,iNode) = tmpBldPtMeshOri(1:3,1:3,iNode)
- BldPtMotionMesh(iWT)%TranslationVel( 1:3,iNode) = tmpBldPtMeshVel(1:3,iNode)
- BldPtMotionMesh(iWT)%RotationVel( 1:3,iNode) = tmpBldPtMeshVel(4:6,iNode)
- BldPtMotionMesh(iWT)%TranslationAcc( 1:3,iNode) = tmpBldPtMeshAcc(1:3,iNode)
- !BldPtMotionMesh(iWT)%RotationAcc( 1:3,iNode) = tmpBldPtMeshAcc(4:6,iNode) ! Rotational acc not included
- call OrientRemap(BldPtMotionMesh(iWT)%Orientation(1:3,1:3,iNode))
- if (TransposeDCM) then
- BldPtMotionMesh(iWT)%Orientation(1:3,1:3,iNode) = transpose(BldPtMotionMesh(iWT)%Orientation(1:3,1:3,iNode))
- endif
+ do iBlade=1,Sim%WT(iWT)%NumBlades
+ do j=1,StrucPts_2_Bld_Map(iWT)%NumMeshPtsPerBlade(iBlade)
+ BldPtMotionMesh(iWT)%Mesh(iBlade)%TranslationDisp(1:3,j) = StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(iBlade)%Position(1:3,j) + Sim%WT(iWT)%OriginInit(1:3) - real(BldPtMotionMesh(iWT)%Mesh(iBlade)%Position(1:3,j), R8Ki)
+ BldPtMotionMesh(iWT)%Mesh(iBlade)%Orientation(1:3,1:3,j) = StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(iBlade)%Orient(1:3,1:3,j)
+ BldPtMotionMesh(iWT)%Mesh(iBlade)%TranslationVel( 1:3,j) = StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(iBlade)%Velocity(1:3,j)
+ BldPtMotionMesh(iWT)%Mesh(iBlade)%RotationVel( 1:3,j) = StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(iBlade)%Velocity(4:6,j)
+ BldPtMotionMesh(iWT)%Mesh(iBlade)%TranslationAcc( 1:3,j) = StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(iBlade)%Accln(1:3,j)
+ call OrientRemap(BldPtMotionMesh(iWT)%Mesh(iBlade)%Orientation(1:3,1:3,j))
+ if (TransposeDCM) then
+ BldPtMotionMesh(iWT)%Mesh(iBlade)%Orientation(1:3,1:3,j) = transpose(BldPtMotionMesh(iWT)%Mesh(iBlade)%Orientation(1:3,1:3,j))
+ endif
+ enddo
enddo
end subroutine Set_MotionMesh
@@ -1849,7 +1918,7 @@ subroutine AD_SetInputMotion( iWT, u_local, &
NacPos_C, NacOri_C, NacVel_C, NacAcc_C, &
BldRootPos_C, BldRootOri_C, BldRootVel_C, BldRootAcc_C, &
ErrStat, ErrMsg )
- integer(IntKi), intent(in ) :: iWT !< this turbine
+ integer(IntKi), intent(in ) :: iWT !< this turbine
type(ADI_InputType), intent(inout) :: u_local ! Only one input (probably at T)
real(c_float), intent(in ) :: HubPos_C( 3 ) !< Hub position
real(c_double), intent(in ) :: HubOri_C( 9 ) !< Hub orientation
@@ -1865,9 +1934,12 @@ subroutine AD_SetInputMotion( iWT, u_local, &
real(c_float), intent(in ) :: BldRootAcc_C( 6*Sim%WT(iWT)%NumBlades ) !< Blade root accelerations
integer(IntKi), intent( out) :: ErrStat
character(ErrMsgLen), intent( out) :: ErrMsg
- integer(IntKi) :: i
+ integer(IntKi) :: iBlade !< current blade
+ integer(IntKi) :: i,j !< generic index variables
+ integer(IntKi) :: n_elems !< number of elements in the mesh
ErrStat = 0_IntKi
ErrMsg = ''
+
! Hub -- NOTE: RotationalAcc not present in the mesh
if ( u_local%AD%rotors(iWT)%HubMotion%Committed ) then
u_local%AD%rotors(iWT)%HubMotion%TranslationDisp(1:3,1) = real(HubPos_C(1:3),R8Ki) + Sim%WT(iWT)%OriginInit(1:3) - real(u_local%AD%rotors(iWT)%HubMotion%Position(1:3,1), R8Ki)
@@ -1880,6 +1952,7 @@ subroutine AD_SetInputMotion( iWT, u_local, &
u_local%AD%rotors(iWT)%HubMotion%Orientation(1:3,1:3,1) = transpose(u_local%AD%rotors(iWT)%HubMotion%Orientation(1:3,1:3,1))
endif
endif
+
! Nacelle -- NOTE: RotationalVel and RotationalAcc not present in the mesh
if ( u_local%AD%rotors(iWT)%NacelleMotion%Committed ) then
u_local%AD%rotors(iWT)%NacelleMotion%TranslationDisp(1:3,1) = real(NacPos_C(1:3),R8Ki) + Sim%WT(iWT)%OriginInit(1:3) - real(u_local%AD%rotors(iWT)%NacelleMotion%Position(1:3,1), R8Ki)
@@ -1891,6 +1964,7 @@ subroutine AD_SetInputMotion( iWT, u_local, &
u_local%AD%rotors(iWT)%NacelleMotion%Orientation(1:3,1:3,1) = transpose(u_local%AD%rotors(iWT)%NacelleMotion%Orientation(1:3,1:3,1))
endif
endif
+
! Blade root
do i=0,Sim%WT(iWT)%numBlades-1
if ( u_local%AD%rotors(iWT)%BladeRootMotion(i+1)%Committed ) then
@@ -1908,9 +1982,10 @@ subroutine AD_SetInputMotion( iWT, u_local, &
enddo
! Blade mesh
- do i=1,Sim%WT(iWT)%numBlades
- if ( u_local%AD%rotors(iWT)%BladeMotion(i)%Committed ) then
- call Transfer_Point_to_Line2( BldPtMotionMesh(iWT), u_local%AD%rotors(iWT)%BladeMotion(i), Map_BldPtMotion_2_AD_Blade(i,iWT), ErrStat, ErrMsg )
+ do iBlade=1,Sim%WT(iWT)%numBlades
+ n_elems = size(BldPtMotionMesh(iWT)%Mesh(iBlade)%Position, 2)
+ if (( u_local%AD%rotors(iWT)%BladeMotion(iBlade)%Committed ) .and. (n_elems > 0)) then
+ call Transfer_Point_to_Line2( BldPtMotionMesh(iWT)%Mesh(iBlade), u_local%AD%rotors(iWT)%BladeMotion(iBlade), Map_BldPtMotion_2_AD_Blade(i,iWT), ErrStat, ErrMsg )
if (ErrStat >= AbortErrLev) return
endif
enddo
@@ -1919,36 +1994,51 @@ end subroutine AD_SetInputMotion
!> Map the loads of the output mesh to the intermediate output mesh.
!! This routine is operating on module level data, hence few inputs
subroutine AD_TransferLoads( iWT, u_local, y_local, ErrStat3, ErrMsg3 )
- integer(IntKi), intent(in ) :: iWT !< this turbine
- type(ADI_InputType), intent(in ) :: u_local ! Only one input (probably at T)
- type(ADI_OutputType), intent(in ) :: y_local ! Only one input (probably at T)
+ integer(IntKi), intent(in ) :: iWT !< Current rotor/turbine
+ type(ADI_InputType), intent(in ) :: u_local !< Only one input (probably at T)
+ type(ADI_OutputType), intent(in ) :: y_local !< Only one input (probably at T)
integer(IntKi), intent( out) :: ErrStat3
character(ErrMsgLen), intent( out) :: ErrMsg3
- integer(IntKi) :: i
- BldPtLoadMesh(iWT)%Force = 0.0_ReKi
- BldPtLoadMesh(iWT)%Moment = 0.0_ReKi
- do i=1,Sim%WT(iWT)%NumBlades
- if ( y_local%AD%rotors(iWT)%BladeLoad(i)%Committed ) then
- if (debugverbose > 4) call MeshPrintInfo( CU, y_local%AD%rotors(iWT)%BladeLoad(i), MeshName='AD%rotors('//trim(Num2LStr(1))//')%BladeLoad('//trim(Num2LStr(i))//')' )
- call Transfer_Line2_to_Point( ADI%y%AD%rotors(iWT)%BladeLoad(i), BldPtLoadMesh_tmp(iWT), Map_AD_BldLoad_P_2_BldPtLoad(i,iWT), &
- ErrStat3, ErrMsg3, u_local%AD%rotors(iWT)%BladeMotion(i), BldPtMotionMesh(iWT) )
- if (ErrStat3 >= AbortErrLev) return
- BldPtLoadMesh(iWT)%Force = BldPtLoadMesh(iWT)%Force + BldPtLoadMesh_tmp(iWT)%Force
- BldPtLoadMesh(iWT)%Moment = BldPtLoadMesh(iWT)%Moment + BldPtLoadMesh_tmp(iWT)%Moment
+ integer(IntKi) :: iBlade !< current blade
+ integer(IntKi) :: n_elems !< number of elements in the mesh
+
+
+ do iBlade=1,Sim%WT(iWT)%NumBlades
+ n_elems = size(BldPtMotionMesh(iWT)%Mesh(iBlade)%Position, 2)
+ if (n_elems > 0) then
+ BldPtLoadMesh(iWT)%Mesh(iBlade)%Force = 0.0_ReKi
+ BldPtLoadMesh(iWT)%Mesh(iBlade)%Moment = 0.0_ReKi
endif
enddo
- if (debugverbose > 4) call MeshPrintInfo( CU, BldPtLoadMesh(iWT), MeshName='BldPtLoadMesh'//trim(Num2LStr(iWT)) )
+
+ do iBlade=1,Sim%WT(iWT)%NumBlades
+ if ( y_local%AD%rotors(iWT)%BladeLoad(iBlade)%Committed ) then
+ if (debugverbose > 4) call MeshPrintInfo( CU, y_local%AD%rotors(iWT)%BladeLoad(iBlade), MeshName='AD%rotors('//trim(Num2LStr(iWT))//')%BladeLoad('//trim(Num2LStr(iBlade))//')' )
+ n_elems = size(BldPtMotionMesh(iWT)%Mesh(iBlade)%Position, 2)
+ if (n_elems > 0) then
+ call Transfer_Line2_to_Point( ADI%y%AD%rotors(iWT)%BladeLoad(iBlade), BldPtLoadMesh_tmp(iWT)%Mesh(iBlade), Map_AD_BldLoad_P_2_BldPtLoad(iBlade,iWT), &
+ ErrStat3, ErrMsg3, u_local%AD%rotors(iWT)%BladeMotion(iBlade), BldPtMotionMesh(iWT)%Mesh(iBlade) )
+ if (ErrStat3 >= AbortErrLev) return
+ BldPtLoadMesh(iWT)%Mesh(iBlade)%Force = BldPtLoadMesh(iWT)%Mesh(iBlade)%Force + BldPtLoadMesh_tmp(iWT)%Mesh(iBlade)%Force
+ BldPtLoadMesh(iWT)%Mesh(iBlade)%Moment = BldPtLoadMesh(iWT)%Mesh(iBlade)%Moment + BldPtLoadMesh_tmp(iWT)%Mesh(iBlade)%Moment
+ endif
+ endif
+ if (debugverbose > 4) call MeshPrintInfo( CU, BldPtLoadMesh(iWT)%Mesh(iBlade), MeshName='BldPtLoadMesh'//trim(Num2LStr(iWT))//'_'//trim(Num2LStr(iBlade)) )
+ enddo
end subroutine AD_TransferLoads
!> Transfer the loads from the load mesh to the temporary array for output
!! This routine is operating on module level data, hence few inputs
subroutine Set_OutputLoadArray(iWT)
integer(IntKi), intent(in ) :: iWT !< current rotor/turbine
- integer(IntKi) :: iNode
+ integer(IntKi) :: iBlade !< current blade
+ integer(IntKi) :: j !< generic index variables
! Set mesh corresponding to input motions
- do iNode=1,NumMeshPts(iWT)
- tmpBldPtMeshFrc(1:3,iNode) = BldPtLoadMesh(iWT)%Force (1:3,iNode)
- tmpBldPtMeshFrc(4:6,iNode) = BldPtLoadMesh(iWT)%Moment(1:3,iNode)
+ do iBlade=1,Sim%WT(iWT)%NumBlades
+ do j=1,StrucPts_2_Bld_Map(iWT)%NumMeshPtsPerBlade(iBlade)
+ StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(iBlade)%Force(1:3,j) = BldPtLoadMesh(iWT)%Mesh(iBlade)%Force( 1:3,j)
+ StrucPts_2_Bld_Map(iWT)%BladePtMeshCoords(iBlade)%Force(4:6,j) = BldPtLoadMesh(iWT)%Mesh(iBlade)%Moment(1:3,j)
+ enddo
enddo
end subroutine Set_OutputLoadArray
@@ -1976,7 +2066,7 @@ subroutine WrVTK_refMeshes(rot_u, RefPoint, ErrStat, ErrMsg)
integer(IntKi), intent( out) :: ErrStat !< error status
character(ErrMsgLen), intent( out) :: ErrMsg !< error message
integer(IntKi) :: nBlades
- integer(IntKi) :: iWT, nWT, k
+ integer(IntKi) :: iWT, nWT, iBlade
character(*), parameter :: RoutineName = 'WrVTK_refMeshes' !< for error handling
integer(IntKi) :: ErrStat2 !< temporary error status
character(ErrMsgLen) :: ErrMsg2 !< temporary error message
@@ -1993,7 +2083,7 @@ subroutine WrVTK_refMeshes(rot_u, RefPoint, ErrStat, ErrMsg)
else
sWT = '.T'//trim(num2lstr(iWT))
endif
-
+
select case (WrOutputsData%WrVTK_Type)
case (1) ! surfaces -- don't write any surface references
call WrVTK_PointsRef( ErrStat2,ErrMsg2); if (Failed()) return;
@@ -2005,6 +2095,7 @@ subroutine WrVTK_refMeshes(rot_u, RefPoint, ErrStat, ErrMsg)
call WrVTK_LinesRef( ErrStat2,ErrMsg2); if (Failed()) return;
end select
enddo
+
contains
logical function Failed()
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
@@ -2019,20 +2110,22 @@ subroutine WrVTK_PointsRef(ErrStat3,ErrMsg3)
ErrMsg3 = ''
! Blade point motion (structural mesh from driver)
- call MeshWrVTKreference(RefPoint, BldPtMotionMesh(iWT), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.BldPtMotionMesh', ErrStat3, ErrMsg3)
- if (ErrStat3 >= AbortErrLev) return
+ do iBlade=1,Sim%WT(iWT)%NumBlades
+ call MeshWrVTKreference(RefPoint, BldPtMotionMesh(iWT)%Mesh(iBlade), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.BldPtMotionMesh', ErrStat3, ErrMsg3)
+ if (ErrStat3 >= AbortErrLev) return
+ enddo
! Blade root motion (point only)
if (allocated(rot_u(iWT)%BladeRootMotion)) then
- do k=1,Sim%WT(iWT)%NumBlades
- if (rot_u(iWT)%BladeRootMotion(k)%Committed) then
- call MeshWrVTKreference(RefPoint, rot_u(iWT)%BladeRootMotion(k), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.BladeRootMotion'//trim(num2lstr(k)), ErrStat3, ErrMsg3 )
+ do iBlade=1,Sim%WT(iWT)%NumBlades
+ if (rot_u(iWT)%BladeRootMotion(iBlade)%Committed) then
+ call MeshWrVTKreference(RefPoint, rot_u(iWT)%BladeRootMotion(iBlade), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.BladeRootMotion'//trim(num2lstr(iBlade)), ErrStat3, ErrMsg3 )
if (ErrStat3 >= AbortErrLev) return
endif
enddo
endif
- ! Nacelle (structural point input
+ ! Nacelle (structural point input)
if ( rot_u(iWT)%NacelleMotion%Committed ) call MeshWrVTKreference(RefPoint, rot_u(iWT)%NacelleMotion, trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.NacelleMotion', ErrStat3, ErrMsg3)
if (ErrStat3 >= AbortErrLev) return
end subroutine WrVTK_PointsRef
@@ -2058,9 +2151,9 @@ subroutine WrVTK_LinesRef(ErrStat3,ErrMsg3)
! Blades
if (allocated(rot_u(iWT)%BladeMotion)) then
- do k=1,Sim%WT(iWT)%NumBlades
- if (rot_u(iWT)%BladeMotion(k)%Committed) then
- call MeshWrVTKreference(RefPoint, rot_u(iWT)%BladeMotion(k), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.Blade'//trim(num2lstr(k)), ErrStat3, ErrMsg3 )
+ do iBlade=1,Sim%WT(iWT)%NumBlades
+ if (rot_u(iWT)%BladeMotion(iBlade)%Committed) then
+ call MeshWrVTKreference(RefPoint, rot_u(iWT)%BladeMotion(iBlade), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.Blade'//trim(num2lstr(iBlade)), ErrStat3, ErrMsg3 )
if (ErrStat3 >= AbortErrLev) return
endif
enddo
@@ -2076,7 +2169,7 @@ subroutine WrVTK_Meshes(rot_u, RefPoint, ErrStat, ErrMsg)
integer(IntKi), intent( out) :: ErrStat !< error status
character(ErrMsgLen), intent( out) :: ErrMsg !< error message
integer(IntKi) :: nBlades
- integer(IntKi) :: iWT, nWT, k
+ integer(IntKi) :: iWT, nWT, iBlade
character(IntfStrLen) :: TmpFileName
character(*), parameter :: RoutineName = 'WrVTK_Meshes' !< for error handling
integer(IntKi) :: ErrStat2 !< temporary error status
@@ -2094,7 +2187,7 @@ subroutine WrVTK_Meshes(rot_u, RefPoint, ErrStat, ErrMsg)
else
sWT = '.T'//trim(num2lstr(iWT))
endif
-
+
select case (WrOutputsData%WrVTK_Type)
case (1) ! surfaces
call WrVTK_Points( ErrStat2,ErrMsg2); if (Failed()) return;
@@ -2124,14 +2217,16 @@ subroutine WrVTK_Points(ErrStat3,ErrMsg3)
ErrMsg3 = ''
! Blade point motion (structural mesh from driver)
- call MeshWrVTK(RefPoint, BldPtMotionMesh(iWT), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.BldPtMotionMesh', n_Global, .true., ErrStat3, ErrMsg3, WrOutputsData%VTK_tWidth)
- if (ErrStat3 >= AbortErrLev) return
+ do iBlade=1,Sim%WT(iWT)%NumBlades
+ call MeshWrVTK(RefPoint, BldPtMotionMesh(iWT)%Mesh(iBlade), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.BldPtMotionMesh'//trim(num2lstr(iBlade)), n_Global, .true., ErrStat3, ErrMsg3, WrOutputsData%VTK_tWidth)
+ if (ErrStat3 >= AbortErrLev) return
+ enddo
! Blade root motion (point only)
if (allocated(rot_u(iWT)%BladeRootMotion)) then
- do k=1,Sim%WT(iWT)%NumBlades
- if (rot_u(iWT)%BladeRootMotion(k)%Committed) then
- call MeshWrVTK(RefPoint, rot_u(iWT)%BladeRootMotion(k), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.BladeRootMotion'//trim(num2lstr(k)), n_Global, .true., ErrStat3, ErrMsg3, WrOutputsData%VTK_tWidth)
+ do iBlade=1,Sim%WT(iWT)%NumBlades
+ if (rot_u(iWT)%BladeRootMotion(iBlade)%Committed) then
+ call MeshWrVTK(RefPoint, rot_u(iWT)%BladeRootMotion(iBlade), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.BladeRootMotion'//trim(num2lstr(iBlade)), n_Global, .true., ErrStat3, ErrMsg3, WrOutputsData%VTK_tWidth)
if (ErrStat3 >= AbortErrLev) return
endif
enddo
@@ -2159,8 +2254,8 @@ subroutine WrVTK_Surfaces(ErrStat3,ErrMsg3)
ErrStat3 = 0_IntKi
ErrMsg3 = ''
-!TODO: use this routine when it is moved out of the driver and into ADI
-! call AD_WrVTK_Surfaces(ADI%u(1)%AD, ADI%y%AD, RefPoint, ADI%m%VTK_Surfaces, n_Global, WrOutputsData%Root, WrOutputsData%VTK_tWidth, 25, WrOutputsData%VTKHubRad)
+ ! TODO: use this routine when it is moved out of the driver and into ADI
+ ! call AD_WrVTK_Surfaces(ADI%u(1)%AD, ADI%y%AD, RefPoint, ADI%m%VTK_Surfaces, n_Global, WrOutputsData%Root, WrOutputsData%VTK_tWidth, 25, WrOutputsData%VTKHubRad)
! Nacelle
if ( rot_u(iWT)%NacelleMotion%Committed ) then
@@ -2186,11 +2281,11 @@ subroutine WrVTK_Surfaces(ErrStat3,ErrMsg3)
! Blades
if (allocated(rot_u(iWT)%BladeMotion)) then
- do k=1,Sim%WT(iWT)%NumBlades
- if (rot_u(iWT)%BladeMotion(k)%Committed) then
- call MeshWrVTK_Ln2Surface (RefPoint, rot_u(iWT)%BladeMotion(k), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.Blade'//trim(num2lstr(k))//'Surface', &
- n_Global, OutputFields, errStat3, errMsg3, WrOutputsData%VTK_tWidth , verts=ADI%m%VTK_Surfaces(iWT)%BladeShape(k)%AirfoilCoords, &
- Sib=ADI%y%AD%rotors(iWT)%BladeLoad(k) )
+ do iBlade=1,Sim%WT(iWT)%NumBlades
+ if (rot_u(iWT)%BladeMotion(iBlade)%Committed) then
+ call MeshWrVTK_Ln2Surface (RefPoint, rot_u(iWT)%BladeMotion(iBlade), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.Blade'//trim(num2lstr(iBlade))//'Surface', &
+ n_Global, OutputFields, errStat3, errMsg3, WrOutputsData%VTK_tWidth , verts=ADI%m%VTK_Surfaces(iWT)%BladeShape(iBlade)%AirfoilCoords, &
+ Sib=ADI%y%AD%rotors(iWT)%BladeLoad(iBlade) )
if (ErrStat3 >= AbortErrLev) return
endif
enddo
@@ -2218,9 +2313,9 @@ subroutine WrVTK_Lines(ErrStat3,ErrMsg3)
! Blades
if (allocated(rot_u(iWT)%BladeMotion)) then
- do k=1,Sim%WT(iWT)%NumBlades
- if (rot_u(iWT)%BladeMotion(k)%Committed) then
- call MeshWrVTK(RefPoint, rot_u(iWT)%BladeMotion(k), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.Blade'//trim(num2lstr(k)), n_Global, .true., ErrStat3, ErrMsg3, WrOutputsData%VTK_tWidth)
+ do iBlade=1,Sim%WT(iWT)%NumBlades
+ if (rot_u(iWT)%BladeMotion(iBlade)%Committed) then
+ call MeshWrVTK(RefPoint, rot_u(iWT)%BladeMotion(iBlade), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.Blade'//trim(num2lstr(iBlade)), n_Global, .true., ErrStat3, ErrMsg3, WrOutputsData%VTK_tWidth)
if (ErrStat3 >= AbortErrLev) return
endif
enddo
@@ -2250,7 +2345,7 @@ subroutine WrVTK_Ground (RefPoint, HalfLengths, FileRootName, errStat, errMsg)
errStat = ErrID_None
errMsg = ""
FileName = TRIM(FileRootName)//'.vtp'
- call WrVTK_header( FileName, NumberOfPoints, NumberOfLines, NumberOfPolys, Un, errStat2, errMsg2 )
+ call WrVTK_header( FileName, NumberOfPoints, NumberOfLines, NumberOfPolys, Un, errStat2, errMsg2 )
call SetErrStat(errStat2,errMsg2,errStat,errMsg,'WrVTK_Ground'); if (errStat >= AbortErrLev) return
WRITE(Un,'(A)') ' '
WRITE(Un,'(A)') ' '
@@ -2260,16 +2355,16 @@ subroutine WrVTK_Ground (RefPoint, HalfLengths, FileRootName, errStat, errMsg)
WRITE(Un,VTK_AryFmt) RefPoint(1) - HalfLengths(1) , RefPoint(2) + HalfLengths(2), RefPoint(3)
WRITE(Un,'(A)') ' '
WRITE(Un,'(A)') ' '
- WRITE(Un,'(A)') ' '
- WRITE(Un,'(A)') ' '
- WRITE(Un,'('//trim(num2lstr(NumberOfPoints))//'(i7))') (ix, ix=0,NumberOfPoints-1)
- WRITE(Un,'(A)') ' '
-
- WRITE(Un,'(A)') ' '
+ WRITE(Un,'(A)') ' '
+ WRITE(Un,'(A)') ' '
+ WRITE(Un,'('//trim(num2lstr(NumberOfPoints))//'(i7))') (ix, ix=0,NumberOfPoints-1)
+ WRITE(Un,'(A)') ' '
+
+ WRITE(Un,'(A)') ' '
WRITE(Un,'(i7)') NumberOfPoints
WRITE(Un,'(A)') ' '
- WRITE(Un,'(A)') ' '
- call WrVTK_footer( Un )
+ WRITE(Un,'(A)') ' '
+ call WrVTK_footer( Un )
end subroutine WrVTK_Ground
@@ -2278,7 +2373,6 @@ end subroutine WrVTK_Ground
subroutine SetTempStorage(ErrStat,ErrMsg)
INTEGER(IntKi), intent(out) :: errStat !< Indicates whether an error occurred (see NWTC_Library)
character(*), intent(out) :: errMsg !< Error message associated with the errStat
- integer(IntKi) :: maxMeshPts, iWT
INTEGER(IntKi) :: errStat2
CHARACTER(ErrMsgLen) :: errMsg2
character(*), parameter :: RoutineName = 'SetTempStorage' !< for error handling
@@ -2288,28 +2382,12 @@ subroutine SetTempStorage(ErrStat,ErrMsg)
ErrStat = ErrID_Fatal
ErrMSg = "Pre-Init has not been called yet"
return
- endif
+ endif
if (minval(NumMeshPts) < 0) then
ErrStat = ErrID_Fatal
ErrMSg = "ADI_C_SetupRotor haven't been called for all rotors"
return
- endif
-
- ! Allocate temporary arrays to simplify data conversions
- maxMeshPts=0_IntKi
- do iWT=1,Sim%NumTurbines
- maxMeshPts=max(maxMeshPts,NumMeshPts(iWT))
- enddo
- if ( allocated(tmpBldPtMeshPos) ) deallocate(tmpBldPtMeshPos)
- if ( allocated(tmpBldPtMeshOri) ) deallocate(tmpBldPtMeshOri)
- if ( allocated(tmpBldPtMeshVel) ) deallocate(tmpBldPtMeshVel)
- if ( allocated(tmpBldPtMeshAcc) ) deallocate(tmpBldPtMeshAcc)
- if ( allocated(tmpBldPtMeshFrc) ) deallocate(tmpBldPtMeshFrc)
- call AllocAry( tmpBldPtMeshPos, 3, maxMeshPts, "tmpBldPtMeshPos", ErrStat2, ErrMsg2 ); if (Failed()) return
- call AllocAry( tmpBldPtMeshOri, 3, 3, maxMeshPts, "tmpBldPtMeshOri", ErrStat2, ErrMsg2 ); if (Failed()) return
- call AllocAry( tmpBldPtMeshVel, 6, maxMeshPts, "tmpBldPtMeshVel", ErrStat2, ErrMsg2 ); if (Failed()) return
- call AllocAry( tmpBldPtMeshAcc, 6, maxMeshPts, "tmpBldPtMeshAcc", ErrStat2, ErrMsg2 ); if (Failed()) return
- call AllocAry( tmpBldPtMeshFrc, 6, maxMeshPts, "tmpBldPtMeshFrc", ErrStat2, ErrMsg2 ); if (Failed()) return
+ endif
contains
logical function Failed()
@@ -2318,26 +2396,21 @@ logical function Failed()
end function Failed
end subroutine SetTempStorage
-
+!--------------------------------------------------------------------
+!> Don't leave junk in memory. So destroy meshes and mappings.
subroutine ClearTmpStorage()
- INTEGER(IntKi) :: errStat2
+ INTEGER(IntKi) :: errStat2, iWT
CHARACTER(ErrMsgLen) :: errMsg2
- ! Arrays
- if (allocated(tmpBldPtMeshPos)) deallocate(tmpBldPtMeshPos)
- if (allocated(tmpBldPtMeshOri)) deallocate(tmpBldPtMeshOri)
- if (allocated(tmpBldPtMeshVel)) deallocate(tmpBldPtMeshVel)
- if (allocated(tmpBldPtMeshAcc)) deallocate(tmpBldPtMeshAcc)
- if (allocated(tmpBldPtMeshFrc)) deallocate(tmpBldPtMeshFrc)
! Meshes
- if (allocated(BldPtMotionMesh )) call ClearMeshArr1(BldPtMotionMesh )
- if (allocated(BldPtLoadMesh )) call ClearMeshArr1(BldPtLoadMesh )
- if (allocated(BldPtLoadMesh_tmp)) call ClearMeshArr1(BldPtLoadMesh_tmp)
- !if (allocated(NacMotionMesh )) call ClearMeshArr1(NacMotionMesh )
- !if (allocated(NacLoadMesh )) call ClearMeshArr1(NacLoadMesh )
+ do iWT=1,Sim%NumTurbines
+ if (allocated(BldPtMotionMesh(iWT)%Mesh)) call ClearMeshArr1(BldPtMotionMesh(iWT)%Mesh)
+ if (allocated(BldPtLoadMesh(iWT)%Mesh)) call ClearMeshArr1(BldPtLoadMesh(iWT)%Mesh)
+ if (allocated(BldPtLoadMesh_tmp(iWT)%Mesh)) call ClearMeshArr1(BldPtLoadMesh_tmp(iWT)%Mesh)
+ enddo
+ ! if (allocated(NacMotionMesh )) call ClearMeshArr1(NacMotionMesh )
+ ! if (allocated(NacLoadMesh )) call ClearMeshArr1(NacLoadMesh )
if (allocated(Map_BldPtMotion_2_AD_Blade )) call ClearMeshMapArr2(Map_BldPtMotion_2_AD_Blade )
- if (allocated(Map_AD_BldLoad_P_2_BldPtLoad)) call ClearMeshMapArr2(Map_AD_BldLoad_P_2_BldPtLoad)
contains
- !> Don't leave junk in memory. So destroy meshes and mappings.
subroutine ClearMeshArr1(MeshName)
type(MeshType), allocatable :: MeshName(:)
integer :: i
diff --git a/modules/hydrodyn/src/WAMIT.f90 b/modules/hydrodyn/src/WAMIT.f90
index e256247382..e5049a9a8c 100644
--- a/modules/hydrodyn/src/WAMIT.f90
+++ b/modules/hydrodyn/src/WAMIT.f90
@@ -189,7 +189,7 @@ SUBROUTINE WAMIT_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, ErrS
real(ReKi) :: WaveNmbr ! Frequency-dependent wave number
COMPLEX(SiKi) :: Fxy ! Phase correction term for Wave excitation forces
real(ReKi) :: tmpAngle ! Frequency and heading and platform offset dependent phase shift angle for Euler's Equation e^(-j*tmpAngle)
- COMPLEX(SiKi), ALLOCATABLE :: HdroExctn_Local (:,:,:) ! Temporary Frequency- and direction-dependent complex hydrodynamic wave excitation force per unit wave amplitude vector (kg/s^2, kg-m/s^2)
+ REAL(ReKi) :: RotateZdegOffset ! PtfmRefZtRot converted to degrees
! Error handling
CHARACTER(ErrMsgLen) :: ErrMsg2 ! Temporary error message for calls
INTEGER(IntKi) :: ErrStat2 ! Temporary error status for calls
@@ -954,7 +954,13 @@ SUBROUTINE WAMIT_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, ErrS
! NOTE: we may end up inadvertantly aborting if the wave direction crosses
! the -Pi / Pi boundary (-180/180 degrees).
- IF ( ( p%WaveField%WaveDirMin < HdroWvDir(1) ) .OR. ( p%WaveField%WaveDirMax > HdroWvDir(NInpWvDir) ) ) THEN
+ ! Need to account for PtfmRefZtRot if NBodyMod=2 (NBody=1 in the case)
+ IF ( p%NBodyMod == 2 ) THEN
+ RotateZdegOffset = InitInp%PtfmRefztRot(1)*R2D
+ ELSE
+ RotateZdegOffset = 0.0
+ END IF
+ IF ( ( (p%WaveField%WaveDirMin-RotateZdegOffset) HdroWvDir(NInpWvDir) ) ) THEN
ErrMsg2 = 'All Wave directions must be within the wave heading angle range available in "' &
//TRIM(InitInp%WAMITFile)//'.3" (inclusive).'
CALL SetErrStat( ErrID_Fatal, ErrMsg2, ErrStat, ErrMsg, RoutineName)
@@ -1003,46 +1009,35 @@ SUBROUTINE WAMIT_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, ErrS
if ( p%NBodyMod == 2 ) then
! Since NBodyMod = 2, then NBody = 1 for this WAMIT object (this requirement is encoded at the HydroDyn module level)
-
- allocate ( HdroExctn_Local(NInpFreq, NInpWvDir, 6), STAT=ErrStat2 )
- IF ( ErrStat2 /= 0 ) THEN
- CALL SetErrStat( ErrID_Fatal, 'Error allocating memory for the HdroExctn_Local array.', ErrStat, ErrMsg, RoutineName)
- CALL Cleanup()
- RETURN
- END IF
-
- do K = 1,6 ! Loop through all wave excitation forces and moments
- do J = 1, NInpWvDir
- TmpCoord(2) = HdroWvDir(J) - InitInp%PtfmRefztRot(1)*R2D ! apply locale Z rotation to heading angle (degrees)
- do I = 1, NInpFreq
- TmpCoord(1) = HdroFreq(I)
- ! Iterpolate to find new coef
- call WAMIT_Interp2D_Cplx( TmpCoord, HdroExctn(:,:,K), HdroFreq, HdroWvDir, LastInd2, HdroExctn_Local(I,J,K), ErrStat2, ErrMsg2 )
- end do
- end do
- end do
-
- ! Now apply rotation and phase shift
+ ! HdroWvDir from WAMIT is effectively in the body-local coordinate system. Need to offset HdroWvDir to be back in the global frame
+ HdroWvDir = HdroWvDir + InitInp%PtfmRefztRot(1)*R2D
+ ! Now apply rotation and phase shift
do J = 1, NInpWvDir
do I = 1, NInpFreq
- ! Fxy = exp(-j * k(w) * ( X*cos(Beta(w)) + Y*sin(Beta(w)) )
+ ! Fxy = exp(-j * k(w) * ( X*cos(Beta(w)) + Y*sin(Beta(w)) )
WaveNmbr = WaveNumber ( HdroFreq(I), InitInp%Gravity, p%WaveField%EffWtrDpth )
tmpAngle = WaveNmbr * ( InitInp%PtfmRefxt(1)*cos(HdroWvDir(J)*D2R) + InitInp%PtfmRefyt(1)*sin(HdroWvDir(J)*D2R) )
TmpRe = cos(tmpAngle)
TmpIm = -sin(tmpAngle)
Fxy = CMPLX( TmpRe, TmpIm )
- HdroExctn(I,J,1) = Fxy*( HdroExctn_Local(I,J,1)*cos(InitInp%PtfmRefztRot(1)) - HdroExctn_Local(I,J,2)*sin(InitInp%PtfmRefztRot(1)) )
- HdroExctn(I,J,2) = Fxy*( HdroExctn_Local(I,J,1)*sin(InitInp%PtfmRefztRot(1)) + HdroExctn_Local(I,J,2)*cos(InitInp%PtfmRefztRot(1)) )
- HdroExctn(I,J,3) = Fxy*( HdroExctn_Local(I,J,3) )
- HdroExctn(I,J,4) = Fxy*( HdroExctn_Local(I,J,4)*cos(InitInp%PtfmRefztRot(1)) - HdroExctn_Local(I,J,5)*sin(InitInp%PtfmRefztRot(1)) )
- HdroExctn(I,J,5) = Fxy*( HdroExctn_Local(I,J,4)*sin(InitInp%PtfmRefztRot(1)) + HdroExctn_Local(I,J,5)*cos(InitInp%PtfmRefztRot(1)) )
- HdroExctn(I,J,6) = Fxy*( HdroExctn_Local(I,J,6) )
+ Ctmp1 = Fxy*( HdroExctn(I,J,1)*cos(InitInp%PtfmRefztRot(1)) - HdroExctn(I,J,2)*sin(InitInp%PtfmRefztRot(1)) )
+ Ctmp2 = Fxy*( HdroExctn(I,J,1)*sin(InitInp%PtfmRefztRot(1)) + HdroExctn(I,J,2)*cos(InitInp%PtfmRefztRot(1)) )
+ Ctmp4 = Fxy*( HdroExctn(I,J,4)*cos(InitInp%PtfmRefztRot(1)) - HdroExctn(I,J,5)*sin(InitInp%PtfmRefztRot(1)) )
+ Ctmp5 = Fxy*( HdroExctn(I,J,4)*sin(InitInp%PtfmRefztRot(1)) + HdroExctn(I,J,5)*cos(InitInp%PtfmRefztRot(1)) )
+
+ HdroExctn(I,J,1) = Ctmp1
+ HdroExctn(I,J,2) = Ctmp2
+ HdroExctn(I,J,4) = Ctmp4
+ HdroExctn(I,J,5) = Ctmp5
+
+ HdroExctn(I,J,3) = Fxy*( HdroExctn(I,J,3) )
+ HdroExctn(I,J,6) = Fxy*( HdroExctn(I,J,6) )
end do
end do
- deallocate(HdroExctn_Local)
+
else
! Apply rotation only for NBodyMod = 1,3
@@ -1204,23 +1199,18 @@ SUBROUTINE WAMIT_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, ErrS
! Now apply the phase shift in the frequency space
- do J = 1, NInpWvDir
- do I = 0,p%WaveField%NStepWave2 ! Loop through the positive frequency components (including zero) of the discrete Fourier transform
-
- ! Compute the frequency of this component:
-
- Omega = I*p%WaveField%WaveDOmega
- ! Fxy = exp(-j * k(w) * ( X*cos(Beta(w)) + Y*sin(Beta(w)) )
- WaveNmbr = WaveNumber ( Omega, InitInp%Gravity, p%WaveField%EffWtrDpth )
- tmpAngle = WaveNmbr * ( InitInp%PtfmRefxt(1)*cos(HdroWvDir(J)*D2R) + InitInp%PtfmRefyt(1)*sin(HdroWvDir(J)*D2R) )
- TmpRe = cos(tmpAngle)
- TmpIm = -sin(tmpAngle)
- Fxy = CMPLX( TmpRe, TmpIm )
+ do I = 0,p%WaveField%NStepWave2 ! Loop through the positive frequency components (including zero) of the discrete Fourier transform
- tmpComplexArr(I) = Fxy*CMPLX(p%WaveField%WaveElevC0(1,I), p%WaveField%WaveElevC0(2,I))
-
+ ! Compute the phase shift of this component, Fxy = exp(-j * k(w) * ( X*cos(Beta(w)) + Y*sin(Beta(w)) ):
+ Omega = I*p%WaveField%WaveDOmega
+ WaveNmbr = WaveNumber ( Omega, InitInp%Gravity, p%WaveField%EffWtrDpth )
+ tmpAngle = WaveNmbr * ( InitInp%PtfmRefxt(1)*cos(p%WaveField%WaveDirArr(I)*D2R) + InitInp%PtfmRefyt(1)*sin(p%WaveField%WaveDirArr(I)*D2R) )
+ TmpRe = cos(tmpAngle)
+ TmpIm = -sin(tmpAngle)
+ Fxy = CMPLX( TmpRe, TmpIm )
- end do
+ ! Apply phase shift
+ tmpComplexArr(I) = Fxy*CMPLX(p%WaveField%WaveElevC0(1,I), p%WaveField%WaveElevC0(2,I))
end do
! Compute the inverse discrete Fourier transforms to find the time-domain
diff --git a/modules/hydrodyn/src/WAMIT2.f90 b/modules/hydrodyn/src/WAMIT2.f90
index 21c4f18ffc..8cbb1da71e 100644
--- a/modules/hydrodyn/src/WAMIT2.f90
+++ b/modules/hydrodyn/src/WAMIT2.f90
@@ -180,6 +180,10 @@ MODULE WAMIT2
TYPE(W2_InitData4D_Type) :: Data4D !< The 4D type from above
END TYPE W2_SumData_Type
+ INTERFACE CheckWamit2WvHdg
+ MODULE PROCEDURE CheckWAMIT2WvHdgDiffData
+ MODULE PROCEDURE CheckWAMIT2WvHdgSumData
+ END INTERFACE CheckWamit2WvHdg
CONTAINS
!----------------------------------------------------------------------------------------------------------------------------------
@@ -819,188 +823,10 @@ SUBROUTINE MnDrift_InitCalc( InitInp, p, MnDriftData, MnDriftForce, ErrMsg, ErrS
RETURN
ENDIF
-
-
- !> 2. Check the data to see if the wave frequencies are present in the QTF data. Since the mean drift term only uses
- !! frequencies where \f$ \omega_1=\omega_2 \f$, the data read in from the files must contain the full range of frequencies
- !! present in the waves.
-
- IF ( MnDriftData%DataIs3D ) THEN
-
- ! Check the low frequency cutoff
- IF ( MINVAL( MnDriftData%Data3D%WvFreq1 ) > InitInp%WaveField%WvLowCOffD ) THEN
- CALL SetErrStat( ErrID_Fatal,' The lowest frequency ( '//TRIM(Num2LStr(MINVAL(MnDriftData%Data3D%WvFreq1)))// &
- ' rad/s for first wave period) data in '//TRIM(MnDriftData%Filename)// &
- ' is above the low frequency cutoff set by WvLowCOffD.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ! Check the high frequency cutoff -- using the Difference high frequency cutoff. The first order high frequency
- ! cutoff is typically too high for this in most cases.
- IF ( (MAXVAL(MnDriftData%Data3D%WvFreq1 ) < InitInp%WaveField%WvHiCOffD) ) THEN
- CALL SetErrStat( ErrID_Fatal,' The highest frequency ( '//TRIM(Num2LStr(MAXVAL(MnDriftData%Data3D%WvFreq1)))// &
- ' rad/s for first wave period) data in '//TRIM(MnDriftData%Filename)// &
- ' is below the high frequency cutoff set by WvHiCOffD.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ELSE IF ( MnDriftData%DataIs4D ) THEN ! only check if not 3D data. If there is 3D data, we default to using it for calculations
-
- ! Check the low frequency cutoff
- IF ( MINVAL( MnDriftData%Data4D%WvFreq1 ) > InitInp%WaveField%WvLowCOffD ) THEN
- CALL SetErrStat( ErrID_Fatal,' The lowest frequency ( '//TRIM(Num2LStr(MINVAL(MnDriftData%Data4D%WvFreq1)))// &
- ' rad/s first wave period) data in '//TRIM(MnDriftData%Filename)// &
- ' is above the low frequency cutoff set by WvLowCOffD.',ErrStat,ErrMsg,RoutineName)
- ENDIF
- IF ( MINVAL( MnDriftData%Data4D%WvFreq2 ) > InitInp%WaveField%WvLowCOffD ) THEN
- CALL SetErrStat( ErrID_Fatal,' The lowest frequency ( '//TRIM(Num2LStr(MINVAL(MnDriftData%Data4D%WvFreq2)))// &
- ' rad/s for second wave period) data in '//TRIM(MnDriftData%Filename)// &
- ' is above the low frequency cutoff set by WvLowCOffD.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ! Check the high frequency cutoff -- using the Difference high frequency cutoff. The first order high frequency
- ! cutoff is typically too high for this in most cases.
- IF ( (MAXVAL(MnDriftData%Data4D%WvFreq1) < InitInp%WaveField%WvHiCOffD) ) THEN
- CALL SetErrStat( ErrID_Fatal,' The highest frequency ( '//TRIM(Num2LStr(MAXVAL(MnDriftData%Data4D%WvFreq1)))// &
- ' rad/s for first wave period) data in '//TRIM(MnDriftData%Filename)// &
- ' is below the high frequency cutoff set by WvHiCOffD.',ErrStat,ErrMsg,RoutineName)
- ENDIF
- IF ( (MAXVAL(MnDriftData%Data4D%WvFreq2) < InitInp%WaveField%WvHiCOffD) ) THEN
- CALL SetErrStat( ErrID_Fatal,' The highest frequency ( '//TRIM(Num2LStr(MAXVAL(MnDriftData%Data4D%WvFreq1)))// &
- ' rad/s second wave period) data in '//TRIM(MnDriftData%Filename)// &
- ' is below the high frequency cutoff set by WvHiCOffD.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ELSE
- ! This is a catastrophic issue. We should not have called this routine without data that is usable for the MnDrift calculation
- CALL SetErrStat( ErrID_Fatal, ' Mean drift calculation called without data.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
+ CALL CheckWAMIT2WvHdg(InitInp,MnDriftData,ErrStatTmp,ErrMsgTmp)
+ CALL SetErrStat(ErrStatTmp,ErrMsgTmp,ErrStat,ErrMsg,RoutineName)
IF ( ErrStat >= AbortErrLev ) RETURN
-
-
- !> 3. Check the data to see if the wave directions are present. May need to adjust for the boundary at +/- PI
- IF ( MnDriftData%DataIs3D ) THEN
-
- ! If we are using multidirectional waves, then we should have more than 1 wave direction in the WAMIT file.
- IF ( InitInp%WaveField%WaveMultiDir .AND. (MnDriftData%Data3D%NumWvDir1 == 1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' WAMIT output file '//TRIM(MnDriftData%Filename)//' only contains one wave '// &
- 'direction at '//TRIM(Num2LStr(MnDriftData%Data3D%WvDir1(1)))//' degrees (first wave direction). '// &
- 'It cannot be used with multidirectional waves. Set WaveDirMod to 0 to use this file.', &
- ErrStat,ErrMsg,RoutineName)
- ELSE IF ( InitInp%WaveField%WaveMultiDir .AND. (MnDriftData%Data3D%NumWvDir2 == 1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' WAMIT output file '//TRIM(MnDriftData%Filename)//' only contains one wave '// &
- 'direction at '//TRIM(Num2LStr(MnDriftData%Data3D%WvDir2(1)))//' degrees (second wave direction). '// &
- 'It cannot be used with multidirectional waves. Set WaveDirMod to 0 to use this file.', &
- ErrStat,ErrMsg,RoutineName)
- ELSE
-
- ! See Known Issues #1 at the top of this file. There may be problems if the data spans the +/- Pi boundary. For
- ! now (since time is limited) we will issue a warning if any of the wave directions for multidirectional waves
- ! or data from the WAMIT file for the wavedirections is close to the +/-pi boundary (>150 degrees, <-150 degrees),
- ! we will issue a warning.
- IF ( (InitInp%WaveField%WaveDirMin > 150.0_SiKi) .OR. (InitInp%WaveField%WaveDirMax < -150.0_SiKi) .OR. &
- (minval(MnDriftData%data3d%WvDir1) > 150.0_SiKi) .OR. (maxval(MnDriftData%data3d%WvDir1) < -150.0_SiKi) .OR. &
- (minval(MnDriftData%data3d%WvDir2) > 150.0_SiKi) .OR. (maxval(MnDriftData%data3d%WvDir2) < -150.0_SiKi) ) THEN
- CALL SetErrStat( ErrID_Warn,' There may be issues with how the wave direction data is handled when the wave '// &
- 'direction of interest is near the +/- 180 direction. This is a known issue with '// &
- 'the WAMIT2 module that has not yet been addressed.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ! Now check the limits for the first wave direction
- ! --> FIXME: sometime fix this to handle the above case. See Known Issue #1 at top of file.
- IF ( InitInp%WaveField%WaveDirMin < MINVAL(MnDriftData%Data3D%WvDir1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Minimum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMin))//' is not'//&
- 'found in the WAMIT data file '//TRIM(MnDriftData%Filename)//' for the first wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
- IF ( InitInp%WaveField%WaveDirMax > MAXVAL(MnDriftData%Data3D%WvDir1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Maximum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMax))//' is not'//&
- 'found in the WAMIT data file '//TRIM(MnDriftData%Filename)//' for the first wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
-
-
- ! Now check the limits for the second wave direction
- ! --> FIXME: sometime fix this to handle the above case. See Known Issue #1 at top of file.
- IF ( InitInp%WaveField%WaveDirMin < MINVAL(MnDriftData%Data3D%WvDir2) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Minimum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMin))//' is not'//&
- 'found in the WAMIT data file '//TRIM(MnDriftData%Filename)//' for the second wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
- IF ( InitInp%WaveField%WaveDirMax > MAXVAL(MnDriftData%Data3D%WvDir2) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Maximum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMax))//' is not'//&
- 'found in the WAMIT data file '//TRIM(MnDriftData%Filename)//' for the second wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
-
- ENDIF
-
- ELSEIF ( MnDriftData%DataIs4D ) THEN
-
- ! If we are using multidirectional waves, then we should have more than 1 wave direction in the WAMIT file.
- IF ( InitInp%WaveField%WaveMultiDir .AND. (MnDriftData%Data4D%NumWvDir1 == 1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' WAMIT output file '//TRIM(MnDriftData%Filename)//' only contains one wave '// &
- 'direction at '//TRIM(Num2LStr(MnDriftData%Data4D%WvDir1(1)))//' degrees (first wave direction). '// &
- 'It cannot be used with multidirectional waves. Set WaveDirMod to 0 to use this file.', &
- ErrStat,ErrMsg,RoutineName)
- ELSE IF ( InitInp%WaveField%WaveMultiDir .AND. (MnDriftData%Data4D%NumWvDir2 == 1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' WAMIT output file '//TRIM(MnDriftData%Filename)//' only contains one wave '// &
- 'direction at '//TRIM(Num2LStr(MnDriftData%Data4D%WvDir2(1)))//' degrees (second wave direction). '// &
- 'It cannot be used with multidirectional waves. Set WaveDirMod to 0 to use this file.', &
- ErrStat,ErrMsg,RoutineName)
- ELSE
-
- ! See Known Issues #1 at the top of this file. There may be problems if the data spans the +/- Pi boundary. For
- ! now (since time is limited) we will issue a warning if any of the wave directions for multidirectional waves
- ! or data from the WAMIT file for the wavedirections is close to the +/-pi boundary (>150 degrees, <-150 degrees),
- ! we will issue a warning.
- IF ( (InitInp%WaveField%WaveDirMin > 150.0_SiKi) .OR. (InitInp%WaveField%WaveDirMax < -150.0_SiKi) .OR. &
- (MINVAL(MnDriftData%Data4D%WvDir1) > 150.0_SiKi) .OR. (MAXVAL(MnDriftData%Data4D%WvDir1) < -150.0_SiKi) .OR. &
- (MINVAL(MnDriftData%Data4D%WvDir2) > 150.0_SiKi) .OR. (MAXVAL(MnDriftData%Data4D%WvDir2) < -150.0_SiKi) ) THEN
- CALL SetErrStat( ErrID_Warn,' There may be issues with how the wave direction data is handled when the wave '// &
- 'direction of interest is near the +/- 180 direction. This is a known issue with '// &
- 'the WAMIT2 module that has not yet been addressed.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ! Now check the limits for the first wave direction
- ! --> FIXME: sometime fix this to handle the above case. See Known Issue #1 at top of file.
- ! --> FIXME: modify to allow shifting values by TwoPi before comparing
- IF ( InitInp%WaveField%WaveDirMin < MINVAL(MnDriftData%Data4D%WvDir1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Minimum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMin))//' is not'//&
- 'found in the WAMIT data file '//TRIM(MnDriftData%Filename)//' for the first wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
- IF ( InitInp%WaveField%WaveDirMax > MAXVAL(MnDriftData%Data4D%WvDir1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Maximum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMax))//' is not'//&
- 'found in the WAMIT data file '//TRIM(MnDriftData%Filename)//' for the first wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
-
-
- ! Now check the limits for the second wave direction
- ! --> FIXME: sometime fix this to handle the above case. See Known Issue #1 at top of file.
- IF ( InitInp%WaveField%WaveDirMin < MINVAL(MnDriftData%Data4D%WvDir2) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Minimum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMin))//' is not'//&
- 'found in the WAMIT data file '//TRIM(MnDriftData%Filename)//' for the second wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
- IF ( InitInp%WaveField%WaveDirMax > MAXVAL(MnDriftData%Data4D%WvDir2) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Maximum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMax))//' is not'//&
- 'found in the WAMIT data file '//TRIM(MnDriftData%Filename)//' for the second wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
-
- ENDIF
-
- ELSE
- ! No data. This is a catastrophic issue. We should not have called this routine without data that is usable for the MnDrift calculation
- CALL SetErrStat( ErrID_Fatal, ' Mean drift calculation called without data.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- IF ( ErrStat >= AbortErrLev ) RETURN
-
-
-
!> 4. Check the data to see if we need to convert to 3D arrays before continuing (4D is sparse in any dimension we want and
!! frequency diagonal is complete). Only check if we don't have 3D data.
@@ -1105,8 +931,8 @@ SUBROUTINE MnDrift_InitCalc( InitInp, p, MnDriftData, MnDriftForce, ErrMsg, ErrS
endif
! NOTE: RotateZMatrixT is the rotation from local to global.
- RotateZMatrixT(:,1) = (/ cos(InitInp%PtfmRefztRot(IBody)), -sin(InitInp%PtfmRefztRot(IBody)) /)
- RotateZMatrixT(:,2) = (/ sin(InitInp%PtfmRefztRot(IBody)), cos(InitInp%PtfmRefztRot(IBody)) /)
+ RotateZMatrixT(1,:) = (/ cos(InitInp%PtfmRefztRot(IBody)), -sin(InitInp%PtfmRefztRot(IBody)) /)
+ RotateZMatrixT(2,:) = (/ sin(InitInp%PtfmRefztRot(IBody)), cos(InitInp%PtfmRefztRot(IBody)) /)
DO ThisDim=1,6
@@ -1220,8 +1046,9 @@ SUBROUTINE MnDrift_InitCalc( InitInp, p, MnDriftData, MnDriftForce, ErrMsg, ErrS
! Now rotate the force components with platform orientation
- MnDriftForce(1:2) = MATMUL( RotateZMatrixT, MnDriftForce(1:2) ) ! Fx and Fy, rotation about z
- MnDriftForce(4:5) = MATMUL( RotateZMatrixT, MnDriftForce(4:5) ) ! Mx and My, rotation about z
+ Idx = (IBody-1)*6
+ MnDriftForce( (Idx+1):(Idx+2) ) = MATMUL( RotateZMatrixT, MnDriftForce( (Idx+1):(Idx+2) ) ) ! Fx and Fy, rotation about z
+ MnDriftForce( (Idx+4):(Idx+5) ) = MATMUL( RotateZMatrixT, MnDriftForce( (Idx+4):(Idx+5) ) ) ! Mx and My, rotation about z
ENDDO ! IBody
@@ -1344,193 +1171,10 @@ SUBROUTINE NewmanApp_InitCalc( InitInp, p, NewmanAppData, NewmanAppForce, ErrMsg
ErrStat = ErrID_None
ErrStatTmp = ErrID_None
-
-
- !> 1. Check the data to see if the wave frequencies are present in the QTF data. Since Newman's approximation only uses
- !! frequencies where \f$ \omega_1=\omega_2 \f$, the data read in from the files must contain the full range of frequencies
- !! present in the waves.
- IF ( NewmanAppData%DataIs3D ) THEN
-
- ! Check the low frequency cutoff
- IF ( MINVAL( NewmanAppData%Data3D%WvFreq1 ) > InitInp%WaveField%WvLowCOff ) THEN
- CALL SetErrStat( ErrID_Fatal,' The lowest frequency ( '//TRIM(Num2LStr(MINVAL(NewmanAppData%Data3D%WvFreq1)))// &
- ' rad/s for first wave period) data in '//TRIM(NewmanAppData%Filename)// &
- ' is above the low frequency cutoff set by WvLowCOff.', &
- ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ! Check the high frequency cutoff -- using the Difference high frequency cutoff. The first order high frequency
- ! cutoff is typically too high for this in most cases.
- IF ( MAXVAL(NewmanAppData%Data3D%WvFreq1 ) < InitInp%WaveField%WvHiCOff ) THEN
- CALL SetErrStat( ErrID_Fatal,' The highest frequency ( '//TRIM(Num2LStr(MAXVAL(NewmanAppData%Data3D%WvFreq1)))// &
- ' rad/s for first wave period) data in '//TRIM(NewmanAppData%Filename)// &
- ' is below the high frequency cutoff set by WvHiCOff.', &
- ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ELSE IF ( NewmanAppData%DataIs4D ) THEN ! only check if not 3D data. If there is 3D data, we default to using it for calculations
-
- ! Check the low frequency cutoff
- IF ( MINVAL( NewmanAppData%Data4D%WvFreq1 ) > InitInp%WaveField%WvLowCOff ) THEN
- CALL SetErrStat( ErrID_Fatal,' The lowest frequency ( '//TRIM(Num2LStr(MINVAL(NewmanAppData%Data4D%WvFreq1)))// &
- ' rad/s first wave period) data in '//TRIM(NewmanAppData%Filename)// &
- ' is above the low frequency cutoff set by WvLowCOff.', &
- ErrStat,ErrMsg,RoutineName)
- ENDIF
- IF ( MINVAL( NewmanAppData%Data4D%WvFreq2 ) > InitInp%WaveField%WvLowCOff ) THEN
- CALL SetErrStat( ErrID_Fatal,' The lowest frequency ( '//TRIM(Num2LStr(MINVAL(NewmanAppData%Data4D%WvFreq2)))// &
- ' rad/s for second wave period) data in '//TRIM(NewmanAppData%Filename)// &
- ' is above the low frequency cutoff set by WvLowCOff.', &
- ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ! Check the high frequency cutoff -- using the Difference high frequency cutoff. The first order high frequency
- ! cutoff is typically too high for this in most cases.
- IF ( MAXVAL(NewmanAppData%Data4D%WvFreq1) < InitInp%WaveField%WvHiCOff ) THEN
- CALL SetErrStat( ErrID_Fatal,' The highest frequency ( '//TRIM(Num2LStr(MAXVAL(NewmanAppData%Data4D%WvFreq1)))// &
- ' rad/s for first wave period) data in '//TRIM(NewmanAppData%Filename)// &
- ' is below the high frequency cutoff set by WvHiCOff.', &
- ErrStat,ErrMsg,RoutineName)
- ENDIF
- IF ( MAXVAL(NewmanAppData%Data4D%WvFreq2) < InitInp%WaveField%WvHiCOff ) THEN
- CALL SetErrStat( ErrID_Fatal,' The highest frequency ( '//TRIM(Num2LStr(MAXVAL(NewmanAppData%Data4D%WvFreq1)))// &
- ' rad/s second wave period) data in '//TRIM(NewmanAppData%Filename)// &
- ' is below the high frequency cutoff set by WvHiCOff.', &
- ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ELSE
- ! This is a catastrophic issue. We should not have called this routine without data that is usable for the NewmanApp calculation
- CALL SetErrStat( ErrID_Fatal, ' Newman approximation calculation called without data.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
+ CALL CheckWAMIT2WvHdg(InitInp,NewmanAppData,ErrStatTmp,ErrMsgTmp)
+ CALL SetErrStat(ErrStatTmp,ErrMsgTmp,ErrStat,ErrMsg,RoutineName)
IF ( ErrStat >= AbortErrLev ) RETURN
-
-
- !> 2. Check the data to see if the wave directions are present. May need to adjust for the boundary at +/- PI
- IF ( NewmanAppData%DataIs3D ) THEN
-
- ! If we are using multidirectional waves, then we should have more than 1 wave direction in the WAMIT file.
- IF ( InitInp%WaveField%WaveMultiDir .AND. (NewmanAppData%Data3D%NumWvDir1 == 1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' WAMIT output file '//TRIM(NewmanAppData%Filename)//' only contains one wave '// &
- 'direction at '//TRIM(Num2LStr(NewmanAppData%Data3D%WvDir1(1)))//' degrees (first wave direction). '// &
- 'It cannot be used with multidirectional waves. Set WaveDirMod to 0 to use this file.', &
- ErrStat,ErrMsg,RoutineName)
- ELSE IF ( InitInp%WaveField%WaveMultiDir .AND. (NewmanAppData%Data3D%NumWvDir2 == 1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' WAMIT output file '//TRIM(NewmanAppData%Filename)//' only contains one wave '// &
- 'direction at '//TRIM(Num2LStr(NewmanAppData%Data3D%WvDir2(1)))//' degrees (second wave direction). '// &
- 'It cannot be used with multidirectional waves. Set WaveDirMod to 0 to use this file.', &
- ErrStat,ErrMsg,RoutineName)
- ELSE
-
- ! See Known Issues #1 at the top of this file. There may be problems if the data spans the +/- Pi boundary. For
- ! now (since time is limited) we will issue a warning if any of the wave directions for multidirectional waves
- ! or data from the WAMIT file for the wavedirections is close to the +/-pi boundary (>150 degrees, <-150 degrees),
- ! we will issue a warning.
- IF ( (InitInp%WaveField%WaveDirMin > 150.0_SiKi) .OR. (InitInp%WaveField%WaveDirMax < -150.0_SiKi) .OR. &
- (minval(NewmanAppData%data3d%WvDir1) > 150.0_SiKi) .OR. (maxval(NewmanAppData%data3d%WvDir1) < -150.0_SiKi) .OR. &
- (minval(NewmanAppData%data3d%WvDir2) > 150.0_SiKi) .OR. (maxval(NewmanAppData%data3d%WvDir2) < -150.0_SiKi) ) THEN
- CALL SetErrStat( ErrID_Warn,' There may be issues with how the wave direction data is handled when the wave '// &
- 'direction of interest is near the +/- 180 direction. This is a known issue with '// &
- 'the WAMIT2 module that has not yet been addressed.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ! Now check the limits for the first wave direction
- ! --> FIXME: sometime fix this to handle the above case. See Known Issue #1 at top of file.
- IF ( InitInp%WaveField%WaveDirMin < MINVAL(NewmanAppData%Data3D%WvDir1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Minimum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMin))//' is not'//&
- 'found in the WAMIT data file '//TRIM(NewmanAppData%Filename)//' for the first wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
- IF ( InitInp%WaveField%WaveDirMax > MAXVAL(NewmanAppData%Data3D%WvDir1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Maximum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMax))//' is not'//&
- 'found in the WAMIT data file '//TRIM(NewmanAppData%Filename)//' for the first wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
-
-
- ! Now check the limits for the second wave direction
- ! --> FIXME: sometime fix this to handle the above case. See Known Issue #1 at top of file.
- IF ( InitInp%WaveField%WaveDirMin < MINVAL(NewmanAppData%Data3D%WvDir2) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Minimum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMin))//' is not'//&
- 'found in the WAMIT data file '//TRIM(NewmanAppData%Filename)//' for the second wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
- IF ( InitInp%WaveField%WaveDirMax > MAXVAL(NewmanAppData%Data3D%WvDir2) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Maximum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMax))//' is not'//&
- 'found in the WAMIT data file '//TRIM(NewmanAppData%Filename)//' for the second wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
-
- ENDIF
-
- ELSEIF ( NewmanAppData%DataIs4D ) THEN
-
- ! If we are using multidirectional waves, then we should have more than 1 wave direction in the WAMIT file.
- IF ( InitInp%WaveField%WaveMultiDir .AND. (NewmanAppData%Data4D%NumWvDir1 == 1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' WAMIT output file '//TRIM(NewmanAppData%Filename)//' only contains one wave '// &
- 'direction at '//TRIM(Num2LStr(NewmanAppData%Data4D%WvDir1(1)))//' degrees (first wave direction). '// &
- 'It cannot be used with multidirectional waves. Set WaveDirMod to 0 to use this file.', &
- ErrStat,ErrMsg,RoutineName)
- ELSE IF ( InitInp%WaveField%WaveMultiDir .AND. (NewmanAppData%Data4D%NumWvDir2 == 1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' WAMIT output file '//TRIM(NewmanAppData%Filename)//' only contains one wave '// &
- 'direction at '//TRIM(Num2LStr(NewmanAppData%Data4D%WvDir2(1)))//' degrees (second wave direction). '// &
- 'It cannot be used with multidirectional waves. Set WaveDirMod to 0 to use this file.', &
- ErrStat,ErrMsg,RoutineName)
- ELSE
-
- ! See Known Issues #1 at the top of this file. There may be problems if the data spans the +/- Pi boundary. For
- ! now (since time is limited) we will issue a warning if any of the wave directions for multidirectional waves
- ! or data from the WAMIT file for the wavedirections is close to the +/-pi boundary (>150 degrees, <-150 degrees),
- ! we will issue a warning.
- IF ( (InitInp%WaveField%WaveDirMin > 150.0_SiKi) .OR. (InitInp%WaveField%WaveDirMax < -150.0_SiKi) .OR. &
- (MINVAL(NewmanAppData%Data4D%WvDir1) > 150.0_SiKi) .OR. (MAXVAL(NewmanAppData%Data4D%WvDir1) < -150.0_SiKi) .OR. &
- (MINVAL(NewmanAppData%Data4D%WvDir2) > 150.0_SiKi) .OR. (MAXVAL(NewmanAppData%Data4D%WvDir2) < -150.0_SiKi) ) THEN
- CALL SetErrStat( ErrID_Warn,' There may be issues with how the wave direction data is handled when the wave '// &
- 'direction of interest is near the +/- 180 direction. This is a known issue with '// &
- 'the WAMIT2 module that has not yet been addressed.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ! Now check the limits for the first wave direction
- ! --> FIXME: sometime fix this to handle the above case. See Known Issue #1 at top of file.
- ! --> FIXME: modify to allow shifting values by TwoPi before comparing
- IF ( InitInp%WaveField%WaveDirMin < MINVAL(NewmanAppData%Data4D%WvDir1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Minimum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMin))//' is not'//&
- 'found in the WAMIT data file '//TRIM(NewmanAppData%Filename)//' for the first wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
- IF ( InitInp%WaveField%WaveDirMax > MAXVAL(NewmanAppData%Data4D%WvDir1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Maximum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMax))//' is not'//&
- 'found in the WAMIT data file '//TRIM(NewmanAppData%Filename)//' for the first wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
-
-
- ! Now check the limits for the second wave direction
- ! --> FIXME: sometime fix this to handle the above case. See Known Issue #1 at top of file.
- IF ( InitInp%WaveField%WaveDirMin < MINVAL(NewmanAppData%Data4D%WvDir2) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Minimum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMin))//' is not'//&
- 'found in the WAMIT data file '//TRIM(NewmanAppData%Filename)//' for the second wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
- IF ( InitInp%WaveField%WaveDirMax > MAXVAL(NewmanAppData%Data4D%WvDir2) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Maximum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMax))//' is not'//&
- 'found in the WAMIT data file '//TRIM(NewmanAppData%Filename)//' for the second wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
-
- ENDIF
-
- ELSE
- ! No data. This is a catastrophic issue. We should not have called this routine without data that is usable for the NewmanApp calculation
- CALL SetErrStat( ErrID_Fatal, ' Newman approximation calculation called without data.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- IF ( ErrStat >= AbortErrLev ) RETURN
-
-
-
!> 3. Check the data to see if we need to convert to 3D arrays before continuing (4D is sparse in any dimension we want and
!! frequency diagonal is complete). Only check if we don't have 3D data.
@@ -1804,8 +1448,8 @@ SUBROUTINE NewmanApp_InitCalc( InitInp, p, NewmanAppData, NewmanAppForce, ErrMsg
! Set rotation
! NOTE: RotateZMatrixT is the rotation from local to global.
- RotateZMatrixT(:,1) = (/ cos(InitInp%PtfmRefztRot(IBody)), -sin(InitInp%PtfmRefztRot(IBody)) /)
- RotateZMatrixT(:,2) = (/ sin(InitInp%PtfmRefztRot(IBody)), cos(InitInp%PtfmRefztRot(IBody)) /)
+ RotateZMatrixT(1,:) = (/ cos(InitInp%PtfmRefztRot(IBody)), -sin(InitInp%PtfmRefztRot(IBody)) /)
+ RotateZMatrixT(2,:) = (/ sin(InitInp%PtfmRefztRot(IBody)), cos(InitInp%PtfmRefztRot(IBody)) /)
! Loop through all the frequencies
DO J=1,InitInp%WaveField%NStepWave2
@@ -2017,110 +1661,10 @@ SUBROUTINE DiffQTF_InitCalc( InitInp, p, DiffQTFData, DiffQTFForce, ErrMsg, ErrS
ErrStat = ErrID_None
ErrStatTmp = ErrID_None
-
- !> 1. Check the data to see if the wave frequencies are present in the QTF data.
-
- IF ( DiffQTFData%DataIs4D ) THEN ! We must have a 4D data set
-
- ! Check the low frequency cutoff
- IF ( MINVAL( DiffQTFData%Data4D%WvFreq1 ) > InitInp%WaveField%WvLowCOffD ) THEN
- CALL SetErrStat( ErrID_Fatal,' The lowest frequency ( '//TRIM(Num2LStr(MINVAL(DiffQTFData%Data4D%WvFreq1)))// &
- ' rad/s first wave period) data in '//TRIM(DiffQTFData%Filename)// &
- ' is above the low frequency cutoff set by WvLowCOffD.', &
- ErrStat,ErrMsg,RoutineName)
- ENDIF
- IF ( MINVAL( DiffQTFData%Data4D%WvFreq2 ) > InitInp%WaveField%WvLowCOffD ) THEN
- CALL SetErrStat( ErrID_Fatal,' The lowest frequency ( '//TRIM(Num2LStr(MINVAL(DiffQTFData%Data4D%WvFreq2)))// &
- ' rad/s for second wave period) data in '//TRIM(DiffQTFData%Filename)// &
- ' is above the low frequency cutoff set by WvLowCOffD.', &
- ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ! Check the high frequency cutoff -- using the Difference high frequency cutoff. The first order high frequency
- ! cutoff is typically too high for this in most cases.
- IF ( MAXVAL(DiffQTFData%Data4D%WvFreq1) < InitInp%WaveField%WvHiCOffD ) THEN
- CALL SetErrStat( ErrID_Fatal,' The highest frequency ( '//TRIM(Num2LStr(MAXVAL(DiffQTFData%Data4D%WvFreq1)))// &
- ' rad/s for first wave period) data in '//TRIM(DiffQTFData%Filename)// &
- ' is below the high frequency cutoff set by WvHiCOffD.', &
- ErrStat,ErrMsg,RoutineName)
- ENDIF
- IF ( MAXVAL(DiffQTFData%Data4D%WvFreq2) < InitInp%WaveField%WvHiCOffD ) THEN
- CALL SetErrStat( ErrID_Fatal,' The highest frequency ( '//TRIM(Num2LStr(MAXVAL(DiffQTFData%Data4D%WvFreq1)))// &
- ' rad/s second wave period) data in '//TRIM(DiffQTFData%Filename)// &
- ' is below the high frequency cutoff set by WvHiCOffD.', &
- ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ELSE
- ! This is a catastrophic issue. We should not have called this routine without data that is usable for the DiffQTF calculation
- CALL SetErrStat( ErrID_Fatal, ' The full Difference QTF method requires 4D data, and was not passed any.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- IF ( ErrStat >= AbortErrLev ) RETURN
-
-
-
-
- ! If we are using multidirectional waves, then we should have more than 1 wave direction in the WAMIT file.
- IF ( InitInp%WaveField%WaveMultiDir .AND. (DiffQTFData%Data4D%NumWvDir1 == 1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' WAMIT output file '//TRIM(DiffQTFData%Filename)//' only contains one wave '// &
- 'direction at '//TRIM(Num2LStr(DiffQTFData%Data4D%WvDir1(1)))//' degrees (first wave direction). '// &
- 'It cannot be used with multidirectional waves. Set WaveDirMod to 0 to use this file.', &
- ErrStat,ErrMsg,RoutineName)
- ELSE IF ( InitInp%WaveField%WaveMultiDir .AND. (DiffQTFData%Data4D%NumWvDir2 == 1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' WAMIT output file '//TRIM(DiffQTFData%Filename)//' only contains one wave '// &
- 'direction at '//TRIM(Num2LStr(DiffQTFData%Data4D%WvDir2(1)))//' degrees (second wave direction). '// &
- 'It cannot be used with multidirectional waves. Set WaveDirMod to 0 to use this file.', &
- ErrStat,ErrMsg,RoutineName)
- ELSE
-
- ! See Known Issues #1 at the top of this file. There may be problems if the data spans the +/- Pi boundary. For
- ! now (since time is limited) we will issue a warning if any of the wave directions for multidirectional waves
- ! or data from the WAMIT file for the wavedirections is close to the +/-pi boundary (>150 degrees, <-150 degrees),
- ! we will issue a warning.
- IF ( (InitInp%WaveField%WaveDirMin > 150.0_SiKi) .OR. (InitInp%WaveField%WaveDirMax < -150.0_SiKi) .OR. &
- (MINVAL(DiffQTFData%Data4D%WvDir1) > 150.0_SiKi) .OR. (MAXVAL(DiffQTFData%Data4D%WvDir1) < -150.0_SiKi) .OR. &
- (MINVAL(DiffQTFData%Data4D%WvDir2) > 150.0_SiKi) .OR. (MAXVAL(DiffQTFData%Data4D%WvDir2) < -150.0_SiKi) ) THEN
- CALL SetErrStat( ErrID_Warn,' There may be issues with how the wave direction data is handled when the wave '// &
- 'direction of interest is near the +/- 180 direction. This is a known issue with '// &
- 'the WAMIT2 module that has not yet been addressed.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ! Now check the limits for the first wave direction
- ! --> FIXME: sometime fix this to handle the above case. See Known Issue #1 at top of file.
- ! --> FIXME: modify to allow shifting values by TwoPi before comparing
- IF ( InitInp%WaveField%WaveDirMin < MINVAL(DiffQTFData%Data4D%WvDir1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Minimum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMin))//' is not'//&
- 'found in the WAMIT data file '//TRIM(DiffQTFData%Filename)//' for the first wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
- IF ( InitInp%WaveField%WaveDirMax > MAXVAL(DiffQTFData%Data4D%WvDir1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Maximum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMax))//' is not'//&
- 'found in the WAMIT data file '//TRIM(DiffQTFData%Filename)//' for the first wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
-
-
- ! Now check the limits for the second wave direction
- ! --> FIXME: sometime fix this to handle the above case. See Known Issue #1 at top of file.
- IF ( InitInp%WaveField%WaveDirMin < MINVAL(DiffQTFData%Data4D%WvDir2) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Minimum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMin))//' is not'//&
- 'found in the WAMIT data file '//TRIM(DiffQTFData%Filename)//' for the second wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
- IF ( InitInp%WaveField%WaveDirMax > MAXVAL(DiffQTFData%Data4D%WvDir2) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Maximum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMax))//' is not'//&
- 'found in the WAMIT data file '//TRIM(DiffQTFData%Filename)//' for the second wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
-
- ENDIF
-
+ CALL CheckWAMIT2WvHdg(InitInp,DiffQTFData,ErrStatTmp,ErrMsgTmp)
+ CALL SetErrStat(ErrStatTmp,ErrMsgTmp,ErrStat,ErrMsg,RoutineName)
IF ( ErrStat >= AbortErrLev ) RETURN
-
-
-
!> 4. Now check to make sure we have data that will work. For the 4D data, it must not be sparse.
!! To check this, we have to check the load components that we will use. So, we will loop through them
!! and set the TmpFlag to true if there is a sparse matrix for one of them.
@@ -2331,8 +1875,8 @@ SUBROUTINE DiffQTF_InitCalc( InitInp, p, DiffQTFData, DiffQTFForce, ErrMsg, ErrS
! Set rotation
! NOTE: RotateZMatrixT is the rotation from local to global.
- RotateZMatrixT(:,1) = (/ cos(InitInp%PtfmRefztRot(IBody)), -sin(InitInp%PtfmRefztRot(IBody)) /)
- RotateZMatrixT(:,2) = (/ sin(InitInp%PtfmRefztRot(IBody)), cos(InitInp%PtfmRefztRot(IBody)) /)
+ RotateZMatrixT(1,:) = (/ cos(InitInp%PtfmRefztRot(IBody)), -sin(InitInp%PtfmRefztRot(IBody)) /)
+ RotateZMatrixT(2,:) = (/ sin(InitInp%PtfmRefztRot(IBody)), cos(InitInp%PtfmRefztRot(IBody)) /)
! Loop through all the frequencies
DO J=1,InitInp%WaveField%NStepWave2
@@ -2523,110 +2067,10 @@ SUBROUTINE SumQTF_InitCalc( InitInp, p, SumQTFData, SumQTFForce, ErrMsg, ErrStat
ErrStat = ErrID_None
ErrStatTmp = ErrID_None
-
- !> 1. Check the data to see if the wave frequencies are present in the QTF data.
-
- IF ( SumQTFData%DataIs4D ) THEN ! We must have a 4D data set
-
- ! Check the low frequency cutoff
- IF ( MINVAL( SumQTFData%Data4D%WvFreq1 ) > InitInp%WaveField%WvLowCOffS ) THEN
- CALL SetErrStat( ErrID_Fatal,' The lowest frequency ( '//TRIM(Num2LStr(MINVAL(SumQTFData%Data4D%WvFreq1)))// &
- ' rad/s first wave period) data in '//TRIM(SumQTFData%Filename)// &
- ' is above the low frequency cutoff set by WvLowCOffS.', &
- ErrStat,ErrMsg,RoutineName)
- ENDIF
- IF ( MINVAL( SumQTFData%Data4D%WvFreq2 ) > InitInp%WaveField%WvLowCOffS ) THEN
- CALL SetErrStat( ErrID_Fatal,' The lowest frequency ( '//TRIM(Num2LStr(MINVAL(SumQTFData%Data4D%WvFreq2)))// &
- ' rad/s for second wave period) data in '//TRIM(SumQTFData%Filename)// &
- ' is above the low frequency cutoff set by WvLowCOffS.', &
- ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ! Check the high frequency cutoff -- using the Difference high frequency cutoff. The first order high frequency
- ! cutoff is typically too high for this in most cases.
- IF ( MAXVAL(SumQTFData%Data4D%WvFreq1) < InitInp%WaveField%WvHiCOffS ) THEN
- CALL SetErrStat( ErrID_Fatal,' The highest frequency ( '//TRIM(Num2LStr(MAXVAL(SumQTFData%Data4D%WvFreq1)))// &
- ' rad/s for first wave period) data in '//TRIM(SumQTFData%Filename)// &
- ' is below the high frequency cutoff set by WvHiCOffS.', &
- ErrStat,ErrMsg,RoutineName)
- ENDIF
- IF ( MAXVAL(SumQTFData%Data4D%WvFreq2) < InitInp%WaveField%WvHiCOffS ) THEN
- CALL SetErrStat( ErrID_Fatal,' The highest frequency ( '//TRIM(Num2LStr(MAXVAL(SumQTFData%Data4D%WvFreq1)))// &
- ' rad/s second wave period) data in '//TRIM(SumQTFData%Filename)// &
- ' is below the high frequency cutoff set by WvHiCOffS.', &
- ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ELSE
- ! This is a catastrophic issue. We should not have called this routine without data that is usable for the SumQTF calculation
- CALL SetErrStat( ErrID_Fatal, ' The full Sum QTF method requires 4D data, and was not passed any.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- IF ( ErrStat >= AbortErrLev ) RETURN
-
-
-
-
- ! If we are using multidirectional waves, then we should have more than 1 wave direction in the WAMIT file.
- IF ( InitInp%WaveField%WaveMultiDir .AND. (SumQTFData%Data4D%NumWvDir1 == 1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' WAMIT output file '//TRIM(SumQTFData%Filename)//' only contains one wave '// &
- 'direction at '//TRIM(Num2LStr(SumQTFData%Data4D%WvDir1(1)))//' degrees (first wave direction). '// &
- 'It cannot be used with multidirectional waves. Set WaveDirMod to 0 to use this file.', &
- ErrStat,ErrMsg,RoutineName)
- ELSE IF ( InitInp%WaveField%WaveMultiDir .AND. (SumQTFData%Data4D%NumWvDir2 == 1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' WAMIT output file '//TRIM(SumQTFData%Filename)//' only contains one wave '// &
- 'direction at '//TRIM(Num2LStr(SumQTFData%Data4D%WvDir2(1)))//' degrees (second wave direction). '// &
- 'It cannot be used with multidirectional waves. Set WaveDirMod to 0 to use this file.', &
- ErrStat,ErrMsg,RoutineName)
- ELSE
-
- ! See Known Issues #1 at the top of this file. There may be problems if the data spans the +/- Pi boundary. For
- ! now (since time is limited) we will issue a warning if any of the wave directions for multidirectional waves
- ! or data from the WAMIT file for the wavedirections is close to the +/-pi boundary (>150 degrees, <-150 degrees),
- ! we will issue a warning.
- IF ( (InitInp%WaveField%WaveDirMin > 150.0_SiKi) .OR. (InitInp%WaveField%WaveDirMax < -150.0_SiKi) .OR. &
- (MINVAL(SumQTFData%Data4D%WvDir1) > 150.0_SiKi) .OR. (MAXVAL(SumQTFData%Data4D%WvDir1) < -150.0_SiKi) .OR. &
- (MINVAL(SumQTFData%Data4D%WvDir2) > 150.0_SiKi) .OR. (MAXVAL(SumQTFData%Data4D%WvDir2) < -150.0_SiKi) ) THEN
- CALL SetErrStat( ErrID_Warn,' There may be issues with how the wave direction data is handled when the wave '// &
- 'direction of interest is near the +/- 180 direction. This is a known issue with '// &
- 'the WAMIT2 module that has not yet been addressed.',ErrStat,ErrMsg,RoutineName)
- ENDIF
-
- ! Now check the limits for the first wave direction
- ! --> FIXME: sometime fix this to handle the above case. See Known Issue #1 at top of file.
- ! --> FIXME: modify to allow shifting values by TwoPi before comparing
- IF ( InitInp%WaveField%WaveDirMin < MINVAL(SumQTFData%Data4D%WvDir1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Minimum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMin))//' is not'//&
- 'found in the WAMIT data file '//TRIM(SumQTFData%Filename)//' for the first wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
- IF ( InitInp%WaveField%WaveDirMax > MAXVAL(SumQTFData%Data4D%WvDir1) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Maximum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMax))//' is not'//&
- 'found in the WAMIT data file '//TRIM(SumQTFData%Filename)//' for the first wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
-
-
- ! Now check the limits for the second wave direction
- ! --> FIXME: sometime fix this to handle the above case. See Known Issue #1 at top of file.
- IF ( InitInp%WaveField%WaveDirMin < MINVAL(SumQTFData%Data4D%WvDir2) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Minimum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMin))//' is not'//&
- 'found in the WAMIT data file '//TRIM(SumQTFData%Filename)//' for the second wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
- IF ( InitInp%WaveField%WaveDirMax > MAXVAL(SumQTFData%Data4D%WvDir2) ) THEN
- CALL SetErrStat( ErrID_Fatal,' Maximum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMax))//' is not'//&
- 'found in the WAMIT data file '//TRIM(SumQTFData%Filename)//' for the second wave direction.', &
- ErrStat, ErrMsg, RoutineName)
- ENDIF
-
- ENDIF
-
+ CALL CheckWAMIT2WvHdg(InitInp,SumQTFData,ErrStatTmp,ErrMsgTmp)
+ CALL SetErrStat(ErrStatTmp,ErrMsgTmp,ErrStat,ErrMsg,RoutineName)
IF ( ErrStat >= AbortErrLev ) RETURN
-
-
-
!> 4. Now check to make sure we have data that will work. For the 4D data, it must not be sparse.
!! To check this, we have to check the load components that we will use. So, we will loop through them
!! and set the TmpFlag to true if there is a sparse matrix for one of them.
@@ -2931,8 +2375,8 @@ SUBROUTINE SumQTF_InitCalc( InitInp, p, SumQTFData, SumQTFForce, ErrMsg, ErrStat
! Set rotation
! NOTE: RotateZMatrixT is the rotation from local to global.
- RotateZMatrixT(:,1) = (/ cos(InitInp%PtfmRefztRot(IBody)), -sin(InitInp%PtfmRefztRot(IBody)) /)
- RotateZMatrixT(:,2) = (/ sin(InitInp%PtfmRefztRot(IBody)), cos(InitInp%PtfmRefztRot(IBody)) /)
+ RotateZMatrixT(1,:) = (/ cos(InitInp%PtfmRefztRot(IBody)), -sin(InitInp%PtfmRefztRot(IBody)) /)
+ RotateZMatrixT(2,:) = (/ sin(InitInp%PtfmRefztRot(IBody)), cos(InitInp%PtfmRefztRot(IBody)) /)
! Loop through all the frequencies
DO J=1,InitInp%WaveField%NStepWave2
@@ -5587,6 +5031,117 @@ SUBROUTINE Destroy_InitData4D(Data4D)
END SUBROUTINE Destroy_InitData4D
+SUBROUTINE CheckWAMIT2WvHdgDiffData(InitInp,W2Data,ErrStat,ErrMsg)
+ TYPE(WAMIT2_InitInputType), INTENT(IN ) :: InitInp !< Input data for initialization routine
+ TYPE(W2_DiffData_Type), INTENT(IN ) :: W2Data
+ INTEGER(IntKi), INTENT( OUT) :: ErrStat
+ CHARACTER(*), INTENT( OUT) :: ErrMsg
+
+ CHARACTER(*), PARAMETER :: RoutineName = 'CheckWAMIT2WvHdgDiffData'
+ INTEGER(IntKi) :: ErrStatTmp
+ CHARACTER(ErrMsgLen) :: ErrMsgTmp
+
+ ErrStat = ErrID_None
+ ErrMsg = ""
+
+ IF ( W2Data%DataIs3D ) THEN
+ CALL CheckWvHdg(InitInp,W2Data%Data3D%NumWvDir1,W2Data%data3d%WvDir1,'first',ErrStatTmp,ErrMsgTmp)
+ CALL SetErrStat(ErrStatTmp,' WAMIT output file '//TRIM(W2Data%Filename)//ErrMsgTmp,ErrStat,ErrMsg,RoutineName)
+ CALL CheckWvHdg(InitInp,W2Data%Data3D%NumWvDir2,W2Data%data3d%WvDir2,'second',ErrStatTmp,ErrMsgTmp)
+ CALL SetErrStat(ErrStatTmp,' WAMIT output file '//TRIM(W2Data%Filename)//ErrMsgTmp,ErrStat,ErrMsg,RoutineName)
+ ELSEIF ( W2Data%DataIs4D ) THEN
+ CALL CheckWvHdg(InitInp,W2Data%Data4D%NumWvDir1,W2Data%data4D%WvDir1,'first',ErrStatTmp,ErrMsgTmp)
+ CALL SetErrStat(ErrStatTmp,' WAMIT output file '//TRIM(W2Data%Filename)//ErrMsgTmp,ErrStat,ErrMsg,RoutineName)
+ CALL CheckWvHdg(InitInp,W2Data%Data4D%NumWvDir2,W2Data%data4D%WvDir2,'second',ErrStatTmp,ErrMsgTmp)
+ CALL SetErrStat(ErrStatTmp,' WAMIT output file '//TRIM(W2Data%Filename)//ErrMsgTmp,ErrStat,ErrMsg,RoutineName)
+ ELSE
+ ! No data. This is a catastrophic issue. We should not have called this routine without data that is usable for the MnDrift calculation
+ CALL SetErrStat( ErrID_Fatal, ' Second-order wave-load calculation called without data.',ErrStat,ErrMsg,RoutineName)
+ ENDIF
+
+ RETURN
+
+END SUBROUTINE CheckWAMIT2WvHdgDiffData
+
+SUBROUTINE CheckWAMIT2WvHdgSumData(InitInp,W2Data,ErrStat,ErrMsg)
+ TYPE(WAMIT2_InitInputType), INTENT(IN ) :: InitInp !< Input data for initialization routine
+ TYPE(W2_SumData_Type), INTENT(IN ) :: W2Data
+ INTEGER(IntKi), INTENT( OUT) :: ErrStat
+ CHARACTER(*), INTENT( OUT) :: ErrMsg
+
+ CHARACTER(*), PARAMETER :: RoutineName = 'CheckWAMIT2WvHdgSumData'
+ INTEGER(IntKi) :: ErrStatTmp
+ CHARACTER(ErrMsgLen) :: ErrMsgTmp
+
+ ErrStat = ErrID_None
+ ErrMsg = ""
+
+ IF ( W2Data%DataIs4D ) THEN
+ CALL CheckWvHdg(InitInp,W2Data%Data4D%NumWvDir1,W2Data%data4D%WvDir1,'first',ErrStatTmp,ErrMsgTmp)
+ CALL SetErrStat(ErrStatTmp,' WAMIT output file '//TRIM(W2Data%Filename)//ErrMsgTmp,ErrStat,ErrMsg,RoutineName)
+ CALL CheckWvHdg(InitInp,W2Data%Data4D%NumWvDir2,W2Data%data4D%WvDir2,'second',ErrStatTmp,ErrMsgTmp)
+ CALL SetErrStat(ErrStatTmp,' WAMIT output file '//TRIM(W2Data%Filename)//ErrMsgTmp,ErrStat,ErrMsg,RoutineName)
+ ELSE
+ ! No data. This is a catastrophic issue. We should not have called this routine without data that is usable for the MnDrift calculation
+ CALL SetErrStat( ErrID_Fatal, ' Second-order wave-load calculation called without data.',ErrStat,ErrMsg,RoutineName)
+ ENDIF
+
+ RETURN
+
+END SUBROUTINE CheckWAMIT2WvHdgSumData
+
+SUBROUTINE CheckWvHdg(InitInp,NumWAMITWvDir,WAMITWvDir,WvDirName,ErrStat,ErrMsg)
+ TYPE(WAMIT2_InitInputType), INTENT(IN ) :: InitInp !< Input data for initialization routine
+ INTEGER(IntKi), INTENT(IN ) :: NumWAMITWvDir
+ REAL(SiKi), INTENT(IN ) :: WAMITWvDir(:)
+ CHARACTER(*), INTENT(IN ) :: WvDirName
+ INTEGER(IntKi), INTENT( OUT) :: ErrStat
+ CHARACTER(*), INTENT( OUT) :: ErrMsg
+
+ REAL(ReKi), PARAMETER :: WvDirTol = 0.001 ! deg
+ REAL(ReKi) :: RotateZdegOffset
+
+ ErrStat = ErrID_None
+ ErrMsg = ""
+
+ ! If we are using multidirectional waves, then we should have more than 1 wave direction in the WAMIT file.
+ IF ( InitInp%WaveField%WaveMultiDir .AND. (NumWAMITWvDir == 1) ) THEN
+ CALL SetErrStat( ErrID_Fatal,' only contains one '//WvDirName//' wave direction at '//TRIM(Num2LStr(WAMITWvDir(1)))//' degrees'// &
+ 'It cannot be used with multidirectional waves. Set WaveDirMod to 0 to use this file.', &
+ ErrStat,ErrMsg,'')
+ ELSE
+ ! See Known Issues #1 at the top of this file. There may be problems if the data spans the +/- Pi boundary. For
+ ! now (since time is limited) we will issue a warning if any of the wave directions for multidirectional waves
+ ! or data from the WAMIT file for the wavedirections is close to the +/-pi boundary (>150 degrees, <-150 degrees),
+ ! we will issue a warning.
+ IF ( (InitInp%WaveField%WaveDirMin > 150.0_SiKi) .OR. (InitInp%WaveField%WaveDirMax < -150.0_SiKi) .OR. &
+ (minval(WAMITWvDir) > 150.0_SiKi) .OR. (maxval(WAMITWvDir) < -150.0_SiKi)) THEN
+ CALL SetErrStat( ErrID_Warn,' There may be issues with how the wave direction data is handled when the wave '// &
+ 'direction of interest is near the +/- 180 direction. This is a known issue with '// &
+ 'the WAMIT2 module that has not yet been addressed.',ErrStat,ErrMsg,'')
+ ENDIF
+
+ if (InitInp%NBodyMod==2) then ! Need to account for PtfmRefztRot (the current WAMIT2 module can only contain one body in this case)
+ RotateZdegOffset = InitInp%PtfmRefztRot(1)*R2D
+ else
+ RotateZdegOffset = 0.0
+ end if
+
+ ! Now check the limits for the wave direction
+ ! --> FIXME: sometime fix this to handle the above case. See Known Issue #1 at top of file.
+ IF ( (InitInp%WaveField%WaveDirMin-RotateZdegOffset) < (MINVAL(WAMITWvDir)-WvDirTol) ) THEN
+ CALL SetErrStat( ErrID_Fatal,' does not contain the minimum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMin))//' for the '//WvDirName//' wave direction.', &
+ ErrStat, ErrMsg, '')
+ ENDIF
+ IF ( (InitInp%WaveField%WaveDirMax-RotateZdegOffset) > (MAXVAL(WAMITWvDir)+WvDirTol) ) THEN
+ CALL SetErrStat( ErrID_Fatal,' does not contain the maximum wave direction required of '//TRIM(Num2LStr(InitInp%WaveField%WaveDirMax))//' for the '//WvDirName//' wave direction.', &
+ ErrStat, ErrMsg, '')
+ ENDIF
+
+ ENDIF
+
+END SUBROUTINE CheckWvHdg
+
!----------------------------------------------------------------------------------------------------------------------------------
END MODULE WAMIT2
diff --git a/modules/seastate/src/SeaState_Output.f90 b/modules/seastate/src/SeaState_Output.f90
index 55f529f490..7e0a6d0ebd 100644
--- a/modules/seastate/src/SeaState_Output.f90
+++ b/modules/seastate/src/SeaState_Output.f90
@@ -428,9 +428,9 @@ subroutine SeaStOut_WriteWaveElev0( Rootname, NStepWave, NGrid, WaveElev1, WaveE
CHARACTER(*), INTENT( IN ) :: Rootname ! filename including full path, minus any file extension.
INTEGER, INTENT( IN ) :: NStepWave ! Number of time steps for the wave kinematics arrays
INTEGER, INTENT( IN ) :: NGrid(3) ! Number of grid points for the wave kinematics arrays
- REAL(SiKi), pointer, INTENT( IN ) :: WaveElev1 (:,:,: ) ! Instantaneous wave elevations at requested locations - 1st order
- REAL(SiKi), pointer, INTENT( IN ) :: WaveElev2 (:,:,: ) ! Instantaneous wave elevations at requested locations - 2nd order
- REAL(SiKi), pointer, INTENT( IN ) :: WaveTime (: ) ! The time values for the wave kinematics (time)
+ REAL(SiKi), allocatable, INTENT( IN ) :: WaveElev1 (:,:,: ) ! Instantaneous wave elevations at requested locations - 1st order
+ REAL(SiKi), allocatable, INTENT( IN ) :: WaveElev2 (:,:,: ) ! Instantaneous wave elevations at requested locations - 2nd order
+ REAL(SiKi), allocatable, INTENT( IN ) :: WaveTime (: ) ! The time values for the wave kinematics (time)
INTEGER, INTENT( OUT ) :: ErrStat ! returns a non-zero value when an error occurs
CHARACTER(*), INTENT( OUT ) :: ErrMsg ! Error message if ErrStat /= ErrID_None
@@ -459,12 +459,11 @@ subroutine SeaStOut_WriteWaveElev0( Rootname, NStepWave, NGrid, WaveElev1, WaveE
! Write the increments from [0, NStepWave] even though for OpenFAST data, NStepWave = 0, but for arbitrary user data this may not be true.
! As a result for WaveMod=5,6 we shouldn't assume periodic waves over the period WaveTMax
DO m= 0,NStepWave
-
- if ( associated(WaveElev2) ) then
- WRITE(UnWv,Frmt) WaveTime(m), WaveElev1(m,i,j) + WaveElev2(m,i,j)
- else
- WRITE(UnWv,Frmt) WaveTime(m), WaveElev1(m,i,j)
- end if
+ if ( allocated(WaveElev2) ) then
+ WRITE(UnWv,Frmt) WaveTime(m), WaveElev1(m,i,j) + WaveElev2(m,i,j)
+ else
+ WRITE(UnWv,Frmt) WaveTime(m), WaveElev1(m,i,j)
+ end if
END DO
CLOSE( UnWv, IOSTAT=ErrStat )
diff --git a/reg_tests/r-test b/reg_tests/r-test
index f3601cfabb..f4cd893e88 160000
--- a/reg_tests/r-test
+++ b/reg_tests/r-test
@@ -1 +1 @@
-Subproject commit f3601cfabbe8bc5692bd31d60367ebecc42d064e
+Subproject commit f4cd893e886cb026435a5cf0e529d715ec62a3a4