diff --git a/modules/hydrodyn/CMakeLists.txt b/modules/hydrodyn/CMakeLists.txt index 93e8e7f06b..feee1f8f7c 100644 --- a/modules/hydrodyn/CMakeLists.txt +++ b/modules/hydrodyn/CMakeLists.txt @@ -62,7 +62,7 @@ target_link_libraries(hydrodyn_driver hydrodyn_driver_subs versioninfolib) add_library(hydrodyn_c_binding SHARED src/HydroDyn_C_Binding.f90 ) -target_link_libraries(hydrodyn_c_binding hydrodynlib versioninfolib) +target_link_libraries(hydrodyn_c_binding hydrodynlib seastlib versioninfolib) if(APPLE OR UNIX) target_compile_definitions(hydrodyn_c_binding PRIVATE IMPLICIT_DLLEXPORT) endif() diff --git a/modules/hydrodyn/python-lib/hydrodyn_library.py b/modules/hydrodyn/python-lib/hydrodyn_library.py index f9abfff7b2..3cf6847bfa 100644 --- a/modules/hydrodyn/python-lib/hydrodyn_library.py +++ b/modules/hydrodyn/python-lib/hydrodyn_library.py @@ -162,8 +162,10 @@ def __init__(self, library_path): def _initialize_routines(self): self.HydroDyn_C_Init.argtypes = [ POINTER(c_char), # OutRootName - POINTER(c_char_p), # input file string - POINTER(c_int), # input file string length + POINTER(c_char_p), # SeaState input file string + POINTER(c_int), # SeaState input file string length + POINTER(c_char_p), # HydroDyn input file string + POINTER(c_int), # HydroDyn input file string length POINTER(c_float), # gravity POINTER(c_float), # defWtrDens POINTER(c_float), # defWtrDpth @@ -216,14 +218,17 @@ def _initialize_routines(self): self.HydroDyn_C_End.restype = c_int # hydrodyn_init ------------------------------------------------------------------------------------------------------------ - def hydrodyn_init(self, input_string_array): + def hydrodyn_init(self, seast_input_string_array, hd_input_string_array): # nodePositions -- N x 6 array -- position info as [x1,y1,z1,Rx1,Ry1,Rz1] # Primary input file will be passed as a single string joined by # C_NULL_CHAR. - input_string = '\x00'.join(input_string_array) - input_string = input_string.encode('utf-8') - input_string_length = len(input_string) + seast_input_string = '\x00'.join(seast_input_string_array) + seast_input_string = seast_input_string.encode('utf-8') + seast_input_string_length = len(seast_input_string) + hd_input_string = '\x00'.join(hd_input_string_array) + hd_input_string = hd_input_string.encode('utf-8') + hd_input_string_length = len(hd_input_string) self._numChannels_c = c_int(0) @@ -264,8 +269,10 @@ def hydrodyn_init(self, input_string_array): # call HydroDyn_C_Init self.HydroDyn_C_Init( _outRootName_c, # IN: rootname for HD file writing - c_char_p(input_string), # IN: input file string - byref(c_int(input_string_length)), # IN: input file string length + c_char_p(seast_input_string), # IN: SeaState input file string + byref(c_int(seast_input_string_length)),# IN: SeaState input file string length + c_char_p(hd_input_string), # IN: HydroDyn input file string + byref(c_int(hd_input_string_length)), # IN: HydroDyn input file string length byref(c_float(self.gravity)), # IN: gravity byref(c_float(self.defWtrDens)), # IN: default water density byref(c_float(self.defWtrDpth)), # IN: default water depth diff --git a/modules/hydrodyn/src/HydroDyn.txt b/modules/hydrodyn/src/HydroDyn.txt index ab893ad5e7..0f4137f018 100644 --- a/modules/hydrodyn/src/HydroDyn.txt +++ b/modules/hydrodyn/src/HydroDyn.txt @@ -86,8 +86,8 @@ typedef ^ ^ DbKi typedef ^ ^ ReKi PtfmLocationX - - - "Supplied by Driver: X coordinate of platform location in the wave field" "m" typedef ^ ^ ReKi PtfmLocationY - - - "Supplied by Driver: Y coordinate of platform location in the wave field" "m" # -typedef ^ ^ INTEGER NStepWave - - - "Total number of frequency components = total number of time steps in the incident wave" - -typedef ^ ^ INTEGER NStepWave2 - - - "NStepWave / 2" - +typedef ^ ^ INTEGER NStepWave - 0 - "Total number of frequency components = total number of time steps in the incident wave" - +typedef ^ ^ INTEGER NStepWave2 - 0 - "NStepWave / 2" - typedef ^ ^ SiKi RhoXg - - - "= WtrDens*Gravity" - typedef ^ ^ INTEGER WaveMod - - - "Incident wave kinematics model {0: none=still water, 1: plane progressive (regular), 2: JONSWAP/Pierson-Moskowitz spectrum (irregular), 3: white-noise spectrum, 4: user-defind spectrum from routine UserWaveSpctrm (irregular), 5: GH BLADED }" - typedef ^ ^ INTEGER WaveStMod - - - "Model for stretching incident wave kinematics to instantaneous free surface {0: none=no stretching, 1: vertical stretching, 2: extrapolation stretching, 3: Wheeler stretching}" - diff --git a/modules/hydrodyn/src/HydroDyn_C_Binding.f90 b/modules/hydrodyn/src/HydroDyn_C_Binding.f90 index 602962c362..37278bac6a 100644 --- a/modules/hydrodyn/src/HydroDyn_C_Binding.f90 +++ b/modules/hydrodyn/src/HydroDyn_C_Binding.f90 @@ -20,6 +20,8 @@ MODULE HydroDyn_C_BINDING USE ISO_C_BINDING + USE SeaState + USE SeaState_Types USE HydroDyn USE HydroDyn_Types USE NWTC_Library @@ -42,7 +44,7 @@ MODULE HydroDyn_C_BINDING !------------------------------------------------------------------------------------ ! Version info for display - type(ProgDesc), parameter :: version = ProgDesc( 'HydroDyn library', '', '' ) + type(ProgDesc), parameter :: version = ProgDesc( 'HydroDyn+SeaState library', '', '' ) !------------------------------------------------------------------------------------ ! Potential issues @@ -70,16 +72,40 @@ MODULE HydroDyn_C_BINDING integer(IntKi) :: InterpOrder !------------------------------ ! Primary HD derived data types - type(HydroDyn_InputType), allocatable :: u(:) !< Inputs at T, T-dt, T-2*dt (history kept for updating states) - type(HydroDyn_InitInputType) :: InitInp !< Initialization data - type(HydroDyn_InitOutputType) :: InitOutData !< Initial output data -- Names, units, and version info. - type(HydroDyn_ParameterType) :: p !< Parameters - type(HydroDyn_ContinuousStateType) :: x(0:2) !< continuous states at Time t and t+dt (predicted) - type(HydroDyn_DiscreteStateType) :: xd(0:2) !< discrete states at Time t and t+dt (predicted) - type(HydroDyn_ConstraintStateType) :: z(0:2) !< Constraint states at Time t and t+dt (predicted) - type(HydroDyn_OtherStateType) :: OtherStates(0:2) !< Initial other/optimization states - type(HydroDyn_OutputType) :: y !< Initial output (outputs are not calculated; only the output mesh is initialized) - type(HydroDyn_MiscVarType) :: m !< Misc variables for optimization (not copied in glue code) + type :: HD_data + type(HydroDyn_InputType), allocatable :: u(:) !< Inputs at T, T-dt, T-2*dt (history kept for updating states) + type(HydroDyn_InitInputType) :: InitInp !< Initialization data + type(HydroDyn_InitOutputType) :: InitOutData !< Initial output data -- Names, units, and version info. + type(HydroDyn_ParameterType) :: p !< Parameters + type(HydroDyn_ContinuousStateType) :: x(0:2) !< continuous states at Time t and t+dt (predicted) + type(HydroDyn_DiscreteStateType) :: xd(0:2) !< discrete states at Time t and t+dt (predicted) + type(HydroDyn_ConstraintStateType) :: z(0:2) !< Constraint states at Time t and t+dt (predicted) + type(HydroDyn_OtherStateType) :: OtherStates(0:2) !< Initial other/optimization states + type(HydroDyn_OutputType) :: y !< Initial output (outputs are not calculated; only the output mesh is initialized) + type(HydroDyn_MiscVarType) :: m !< Misc variables for optimization (not copied in glue code) + logical :: Initialized = .FALSE. + end type HD_data + + type(HD_data) :: HD + + ! Primary SeaState derived data types + ! NOTE: SeaSt does not contain states, so only using single instance of states. + type :: SeaSt_data + type(SeaSt_InputType) :: u !< Inputs at T -- since no states, only need single + type(SeaSt_InitInputType) :: InitInp !< Initialization data + type(SeaSt_InitOutputType) :: InitOutData !< Initial output data -- Names, units, and version info. + type(SeaSt_ParameterType) :: p !< Parameters + type(SeaSt_ContinuousStateType) :: x !< continuous states -- contains no data + type(SeaSt_DiscreteStateType) :: xd !< discrete states -- contains no data + type(SeaSt_ConstraintStateType) :: z !< Constraint states -- contains no data + type(SeaSt_OtherStateType) :: OtherStates !< Initial other/optimization states + type(SeaSt_OutputType) :: y !< Initial output (outputs are not calculated; only the output mesh is initialized) + type(SeaSt_MiscVarType) :: m !< Misc variables for optimization (not copied in glue code) + logical :: Initialized = .FALSE. + end type SeaSt_data + + type(SeaSt_data) :: SeaSt + !------------------------------ ! Time tracking ! When we are performing a correction step, time information of previous @@ -174,7 +200,9 @@ end subroutine SetErr !=============================================================================================================== !--------------------------------------------- HydroDyn Init---------------------------------------------------- !=============================================================================================================== -SUBROUTINE HydroDyn_C_Init( OutRootName_C, InputFileString_C, InputFileStringLength_C, & +SUBROUTINE HydroDyn_C_Init( OutRootName_C, & + SeaSt_InputFileString_C, SeaSt_InputFileStringLength_C, & + HD_InputFileString_C, HD_InputFileStringLength_C, & Gravity_C, defWtrDens_C, defWtrDpth_C, defMSL2SWL_C, & PtfmRefPtPositionX_C, PtfmRefPtPositionY_C, & NumNodePts_C, InitNodePositions_C, & @@ -189,8 +217,10 @@ SUBROUTINE HydroDyn_C_Init( OutRootName_C, InputFileString_C, InputFileStringLen #endif character(kind=c_char), intent(in ) :: OutRootName_C(IntfStrLen) !< Root name to use for echo files and other - type(c_ptr), intent(in ) :: InputFileString_C !< Input file as a single string with lines deliniated by C_NULL_CHAR - integer(c_int), intent(in ) :: InputFileStringLength_C !< lenght of the input file string + type(c_ptr), intent(in ) :: SeaSt_InputFileString_C !< SeaSt input file as a single string with lines deliniated by C_NULL_CHAR + integer(c_int), intent(in ) :: SeaSt_InputFileStringLength_C !< SeaSt length of the input file string + type(c_ptr), intent(in ) :: HD_InputFileString_C !< HD input file as a single string with lines deliniated by C_NULL_CHAR + integer(c_int), intent(in ) :: HD_InputFileStringLength_C !< HD length of the input file string real(c_float), intent(in ) :: Gravity_C !< Gravitational constant (set by calling code) real(c_float), intent(in ) :: defWtrDens_C !< Default value for water density (may be overridden by input file) real(c_float), intent(in ) :: defWtrDpth_C !< Default value for water density (may be overridden by input file) @@ -213,16 +243,17 @@ SUBROUTINE HydroDyn_C_Init( OutRootName_C, InputFileString_C, InputFileStringLen character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) !< Error message (C_NULL_CHAR terminated) ! Local Variables - character(IntfStrLen) :: OutRootName !< Root name to use for echo files and other - character(kind=C_char, len=InputFileStringLength_C), pointer :: InputFileString !< Input file as a single string with NULL chracter separating lines - - real(DbKi) :: TimeInterval !< timestep for HD - integer(IntKi) :: ErrStat !< aggregated error message - 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 - character(*), parameter :: RoutineName = 'HydroDyn_C_Init' !< for error handling + character(IntfStrLen) :: OutRootName !< Root name to use for echo files and other + character(kind=C_char, len=SeaSt_InputFileStringLength_C), pointer :: SeaSt_InputFileString !< Input file as a single string with NULL chracter separating lines + character(kind=C_char, len=HD_InputFileStringLength_C), pointer :: HD_InputFileString !< Input file as a single string with NULL chracter separating lines + + real(DbKi) :: TimeInterval !< timestep for HD + integer(IntKi) :: ErrStat !< aggregated error message + 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 + character(*), parameter :: RoutineName = 'HydroDyn_C_Init' !< for error handling ! Initialize error handling ErrStat = ErrID_None @@ -240,43 +271,15 @@ SUBROUTINE HydroDyn_C_Init( OutRootName_C, InputFileString_C, InputFileStringLen if (Failed()) return endif - ! Get fortran pointer to C_NULL_CHAR deliniated input file as a string - call C_F_pointer(InputFileString_C, InputFileString) - - ! Get the data to pass to HD_Init - call InitFileInfo(InputFileString, InitInp%PassedFileData, ErrStat2, ErrMsg2); if (Failed()) return - - ! For diagnostic purposes, the following can be used to display the contents - ! of the InFileInfo data structure. - ! CU is the screen -- system dependent. - !call Print_FileInfo_Struct( CU, InitInp%PassedFileData ) - ! Set other inputs for calling HydroDyn_Init - InitInp%InputFile = "passed_hd_file" ! dummy - InitInp%UseInputFile = .FALSE. ! this probably should be passed in -! InitInp%HasIce = .FALSE. ! Always keep at false unless interfacing to ice modules - ! Linearization - ! for now, set linearization to false. Pass this in later when interface supports it - ! Note: we may want to linearize at T=0 for added mass effects, but that might be - ! special case - InitInp%Linearize = .FALSE. - - ! RootName -- for output of echo or other files - OutRootName = TRANSFER( OutRootName_C, OutRootName ) - i = INDEX(OutRootName,C_NULL_CHAR) - 1 ! if this has a c null character at the end... - if ( i > 0 ) OutRootName = OutRootName(1:I) ! remove it - InitInp%OutRootName = trim(OutRootName) - - ! Values passed in - InitInp%Gravity = REAL(Gravity_C, ReKi) - InitInp%WtrDens = REAL(defWtrDens_C, ReKi) ! use values from SeaState - InitInp%WtrDpth = REAL(defWtrDpth_C, ReKi) ! use values from SeaState - InitInp% MSL2SWL = REAL(defMSL2SWL_C, ReKi) ! use values from SeaState + !-------------------------------------------------------------------------------------------------------------------------------- + ! Initialize wrapper variables + !-------------------------------------------------------------------------------------------------------------------------------- + ! Simulation time TimeInterval = REAL(DT_C, DbKi) dT_Global = TimeInterval ! Assume this DT is constant for all simulation N_Global = 0_IntKi ! Assume we are on timestep 0 at start t_initial = REAL(T_Initial_C, DbKi) - InitInp%TMax = REAL(TMax_C, DbKi) ! Number of bodies and initial positions ! - NumNodePts is the number of interface Mesh points we are expecting on the python @@ -294,31 +297,15 @@ SUBROUTINE HydroDyn_C_Init( OutRootName_C, InputFileString_C, InputFileStringLen call AllocAry( tmpNodeFrc, 6, NumNodePts, "tmpNodeFrc", ErrStat2, ErrMsg2 ); if (Failed()) return tmpNodePos(1:6,1:NumNodePts) = reshape( real(InitNodePositions_C(1:6*NumNodePts),ReKi), (/6,NumNodePts/) ) - ! Platform reference position - ! The HD model uses this for building the moddel. This is only specified as an (X,Y) - ! position (no Z). - InitInp%PtfmLocationX = REAL(PtfmRefPtPositionX_C, ReKi) - InitInp%PtfmLocationY = REAL(PtfmRefPtPositionY_C, ReKi) - - - ! Wave eleveation output - ! Wave elevations can be exported for a set of points (grid or any other layout). - ! This feature is used only in the driver codes for exporting for visualization - ! and could be added to this inteface. - ! Skipping this for now. Maybe add later. - !InitInp%WaveElevXY - - !---------------------------------------------------- - ! Allocate input array u and corresponding InputTimes - !---------------------------------------------------- - ! These inputs are used in the time stepping algorithm within HD_UpdateStates + ! Allocate input array u and corresponding InputTimes for SeaState and HD + ! These inputs are used in the time stepping algorithm within HD%UpdateStates ! For quadratic interpolation (InterpOrder==2), 3 timesteps are used. For ! linear (InterOrder==1), 2 timesteps (the HD code can handle either). ! u(1) inputs at t ! u(2) inputs at t - dt ! u(3) inputs at t - 2*dt ! quadratic only - allocate(u(InterpOrder+1), STAT=ErrStat2) + allocate(HD%u(InterpOrder+1), STAT=ErrStat2) if (ErrStat2 /= 0) then ErrStat2 = ErrID_Fatal ErrMsg2 = "Could not allocate inuput" @@ -327,28 +314,163 @@ SUBROUTINE HydroDyn_C_Init( OutRootName_C, InputFileString_C, InputFileStringLen call AllocAry( InputTimes, InterpOrder+1, "InputTimes", ErrStat2, ErrMsg2 ); if (Failed()) return + + !-------------------------------------------------------------------------------------------------------------------------------- + ! SeaState initialize + !-------------------------------------------------------------------------------------------------------------------------------- + + ! Get fortran pointer to C_NULL_CHAR deliniated input file as a string + call C_F_pointer(SeaSt_InputFileString_C, SeaSt_InputFileString) + + ! Get the data to pass to SeaSt%Init + call InitFileInfo(SeaSt_InputFileString, SeaSt%InitInp%PassedFileData, ErrStat2, ErrMsg2); if (Failed()) return + + ! For diagnostic purposes, the following can be used to display the contents + ! of the InFileInfo data structure. + ! CU is the screen -- system dependent. + !call Print_FileInfo_Struct( CU, SeaSt%InitInp%PassedFileData ) + + ! Set other inputs for calling SeaState_Init + SeaSt%InitInp%hasIce = .FALSE. ! Always keep at false unless interfacing to ice modules + SeaSt%InitInp%InputFile = "passed_SeaSt_file" ! dummy + SeaSt%InitInp%UseInputFile = .FALSE. ! this probably should be passed in + ! Linearization + ! for now, set linearization to false. Pass this in later when interface supports it + ! Note: we may want to linearize at T=0 for added mass effects, but that might be + ! special case + HD%InitInp%Linearize = .FALSE. + + ! RootName -- for output of echo or other files + OutRootName = TRANSFER( OutRootName_C, OutRootName ) + i = INDEX(OutRootName,C_NULL_CHAR) - 1 ! if this has a c null character at the end... + if ( i > 0 ) OutRootName = OutRootName(1:I) ! remove it + SeaSt%InitInp%OutRootName = trim(OutRootName)//".SEA" + + ! Values passed in + SeaSt%InitInp%Gravity = REAL(Gravity_C, ReKi) + SeaSt%InitInp%defWtrDens = REAL(defWtrDens_C, ReKi) ! use values from SeaState + SeaSt%InitInp%defWtrDpth = REAL(defWtrDpth_C, ReKi) ! use values from SeaState + SeaSt%InitInp%defMSL2SWL = REAL(defMSL2SWL_C, ReKi) ! use values from SeaState + SeaSt%InitInp%TMax = REAL(TMax_C, DbKi) + + ! Wave elevation output + ! Wave elevations can be exported for a set of points (grid or any other layout). + ! This feature is used only in the driver codes for exporting for visualization + ! and could be added to this inteface. + ! Skipping this for now. Maybe add later. + !SeaSt%InitInp%WaveElevXY + + call SeaSt_Init( SeaSt%InitInp, SeaSt%u, SeaSt%p, SeaSt%x, SeaSt%xd, SeaSt%z, SeaSt%OtherStates, SeaSt%y, SeaSt%m, TimeInterval, SeaSt%InitOutData, ErrStat, ErrMsg ) + if (Failed()) return + SeaSt%Initialized = .true. + + + + !-------------------------------------------------------------------------------------------------------------------------------- + ! HydroDyn initialize + !-------------------------------------------------------------------------------------------------------------------------------- + + ! Get fortran pointer to C_NULL_CHAR deliniated input file as a string + call C_F_pointer(HD_InputFileString_C, HD_InputFileString) + + ! Get the data to pass to HD%Init + call InitFileInfo(HD_InputFileString, HD%InitInp%PassedFileData, ErrStat2, ErrMsg2); if (Failed()) return + + ! For diagnostic purposes, the following can be used to display the contents + ! of the InFileInfo data structure. + ! CU is the screen -- system dependent. + !call Print_FileInfo_Struct( CU, HD%InitInp%PassedFileData ) + + ! Set other inputs for calling HydroDyn_Init + HD%InitInp%InputFile = "passed_hd_file" ! dummy + HD%InitInp%UseInputFile = .FALSE. ! this probably should be passed in + ! Linearization + ! for now, set linearization to false. Pass this in later when interface supports it + ! Note: we may want to linearize at T=0 for added mass effects, but that might be + ! special case + HD%InitInp%Linearize = .FALSE. + + ! RootName -- for output of echo or other files + OutRootName = TRANSFER( OutRootName_C, OutRootName ) + i = INDEX(OutRootName,C_NULL_CHAR) - 1 ! if this has a c null character at the end... + if ( i > 0 ) OutRootName = OutRootName(1:I) ! remove it + HD%InitInp%OutRootName = trim(OutRootName)//".HD" + + ! Values passed in + HD%InitInp%Gravity = REAL(Gravity_C, ReKi) + HD%InitInp%WtrDens = REAL(defWtrDens_C, ReKi) ! use values from SeaState + HD%InitInp%WtrDpth = REAL(defWtrDpth_C, ReKi) ! use values from SeaState + HD%InitInp% MSL2SWL = REAL(defMSL2SWL_C, ReKi) ! use values from SeaState + HD%InitInp%TMax = REAL(TMax_C, DbKi) + + ! Transfer data from SeaState + ! Need to set up other module's InitInput data here because we will also need to clean up SeaState data and would rather not defer that cleanup + HD%InitInp%NStepWave = SeaSt%InitOutData%NStepWave + HD%InitInp%NStepWave2 = SeaSt%InitOutData%NStepWave2 + HD%InitInp%RhoXg = SeaSt%InitOutData%RhoXg + HD%InitInp%WaveMod = SeaSt%InitOutData%WaveMod + HD%InitInp%WaveStMod = SeaSt%InitOutData%WaveStMod + HD%InitInp%WaveDirMod = SeaSt%InitOutData%WaveDirMod + HD%InitInp%WvLowCOff = SeaSt%InitOutData%WvLowCOff + HD%InitInp%WvHiCOff = SeaSt%InitOutData%WvHiCOff + HD%InitInp%WvLowCOffD = SeaSt%InitOutData%WvLowCOffD + HD%InitInp%WvHiCOffD = SeaSt%InitOutData%WvHiCOffD + HD%InitInp%WvLowCOffS = SeaSt%InitOutData%WvLowCOffS + HD%InitInp%WvHiCOffS = SeaSt%InitOutData%WvHiCOffS + HD%InitInp%InvalidWithSSExctn = SeaSt%InitOutData%InvalidWithSSExctn + + HD%InitInp%WaveDirMin = SeaSt%InitOutData%WaveDirMin + HD%InitInp%WaveDirMax = SeaSt%InitOutData%WaveDirMax + HD%InitInp%WaveDir = SeaSt%InitOutData%WaveDir + HD%InitInp%WaveMultiDir = SeaSt%InitOutData%WaveMultiDir + HD%InitInp%WaveDOmega = SeaSt%InitOutData%WaveDOmega + HD%InitInp%MCFD = SeaSt%InitOutData%MCFD + + CALL MOVE_ALLOC( SeaSt%InitOutData%WaveElev0, HD%InitInp%WaveElev0 ) + if(associated(SeaSt%InitOutData%WaveTime )) HD%InitInp%WaveTime => SeaSt%InitOutData%WaveTime + if(associated(SeaSt%InitOutData%WaveDynP )) HD%InitInp%WaveDynP => SeaSt%InitOutData%WaveDynP + if(associated(SeaSt%InitOutData%WaveAcc )) HD%InitInp%WaveAcc => SeaSt%InitOutData%WaveAcc + if(associated(SeaSt%InitOutData%WaveVel )) HD%InitInp%WaveVel => SeaSt%InitOutData%WaveVel + if(associated(SeaSt%InitOutData%PWaveDynP0)) HD%InitInp%PWaveDynP0 => SeaSt%InitOutData%PWaveDynP0 + if(associated(SeaSt%InitOutData%PWaveAcc0 )) HD%InitInp%PWaveAcc0 => SeaSt%InitOutData%PWaveAcc0 + if(associated(SeaSt%InitOutData%PWaveVel0 )) HD%InitInp%PWaveVel0 => SeaSt%InitOutData%PWaveVel0 + if(associated(SeaSt%InitOutData%WaveElevC0)) HD%InitInp%WaveElevC0 => SeaSt%InitOutData%WaveElevC0 + CALL MOVE_ALLOC( SeaSt%InitOutData%WaveElevC, HD%InitInp%WaveElevC ) + if(associated(SeaSt%InitOutData%WaveDirArr)) HD%InitInp%WaveDirArr => SeaSt%InitOutData%WaveDirArr + if(associated(SeaSt%InitOutData%WaveElev1 )) HD%InitInp%WaveElev1 => SeaSt%InitOutData%WaveElev1 + if(associated(SeaSt%InitOutData%WaveElev2 )) HD%InitInp%WaveElev2 => SeaSt%InitOutData%WaveElev2 + + HD%InitInp%WaveAccMCF => SeaSt%InitOutData%WaveAccMCF + HD%InitInp%PWaveAccMCF0 => SeaSt%InitOutData%PWaveAccMCF0 + + call SeaSt_Interp_CopyParam(SeaSt%InitOutData%SeaSt_Interp_p, HD%InitInp%SeaSt_Interp_p, MESH_NEWCOPY, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) + + + ! Platform reference position + ! The HD model uses this for building the moddel. This is only specified as an (X,Y) + ! position (no Z). + HD%InitInp%PtfmLocationX = REAL(PtfmRefPtPositionX_C, ReKi) + HD%InitInp%PtfmLocationY = REAL(PtfmRefPtPositionY_C, ReKi) + + + !------------------------------------------------------------- ! Call the main subroutine HydroDyn_Init - ! TimeInterval and InitInp are passed into HD_Init, all the rest are set by HD_Init + ! TimeInterval and HD%InitInp are passed into HD%Init, all the rest are set by HD%Init ! ! NOTE: Pass u(1) only (this is empty and will be set inside Init). We will copy ! this to u(2) and u(3) afterwards - call HydroDyn_Init( InitInp, u(1), p, x(STATE_CURR), xd(STATE_CURR), z(STATE_CURR), OtherStates(STATE_CURR), y, m, TimeInterval, InitOutData, ErrStat2, ErrMsg2 ) + call HydroDyn_Init( HD%InitInp, HD%u(1), HD%p, HD%x(STATE_CURR), HD%xd(STATE_CURR), HD%z(STATE_CURR), HD%OtherStates(STATE_CURR), HD%y, HD%m, TimeInterval, HD%InitOutData, ErrStat2, ErrMsg2 ) if (Failed()) return + HD%Initialized = .true. !------------------------------------------------------------- ! Sanity checks - !------------------------------------------------------------- call CheckDepth(ErrStat2,ErrMsg2); if (Failed()) return call CheckNodes(ErrStat2,ErrMsg2); if (Failed()) return - !------------------------------------------------------------- - ! Set the interface meshes for motion inputs and loads output - !------------------------------------------------------------- - call SetMotionLoadsInterfaceMeshes(ErrStat2,ErrMsg2); if (Failed()) return - - !------------------------------------------------------------- ! Setup other prior timesteps ! We fill InputTimes with negative times, but the Input values are identical for each of those times; this allows @@ -357,7 +479,7 @@ SUBROUTINE HydroDyn_C_Init( OutRootName_C, InputFileString_C, InputFileStringLen ! order = SIZE(Input) !------------------------------------------------------------- do i=2,InterpOrder+1 - call HydroDyn_CopyInput (u(1), u(i), MESH_NEWCOPY, Errstat2, ErrMsg2) + call HydroDyn_CopyInput (HD%u(1), HD%u(i), MESH_NEWCOPY, Errstat2, ErrMsg2) if (Failed()) return enddo do i = 1, InterpOrder + 1 @@ -368,40 +490,46 @@ SUBROUTINE HydroDyn_C_Init( OutRootName_C, InputFileString_C, InputFileStringLen !------------------------------------------------------------- ! Initial setup of other pieces of x,xd,z,OtherStates - !------------------------------------------------------------- - CALL HydroDyn_CopyContState ( x( STATE_CURR), x( STATE_PRED), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyDiscState ( xd( STATE_CURR), xd( STATE_PRED), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyConstrState( z( STATE_CURR), z( STATE_PRED), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyOtherState ( OtherStates(STATE_CURR), OtherStates(STATE_PRED), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyContState ( HD%x( STATE_CURR), HD%x( STATE_PRED), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyDiscState ( HD%xd( STATE_CURR), HD%xd( STATE_PRED), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyConstrState( HD%z( STATE_CURR), HD%z( STATE_PRED), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyOtherState ( HD%OtherStates(STATE_CURR), HD%OtherStates(STATE_PRED), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return !------------------------------------------------------------- ! Setup the previous timestep copies of states - !------------------------------------------------------------- - CALL HydroDyn_CopyContState ( x( STATE_CURR), x( STATE_LAST), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyDiscState ( xd( STATE_CURR), xd( STATE_LAST), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyConstrState( z( STATE_CURR), z( STATE_LAST), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyOtherState ( OtherStates(STATE_CURR), OtherStates(STATE_LAST), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyContState ( HD%x( STATE_CURR), HD%x( STATE_LAST), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyDiscState ( HD%xd( STATE_CURR), HD%xd( STATE_LAST), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyConstrState( HD%z( STATE_CURR), HD%z( STATE_LAST), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyOtherState ( HD%OtherStates(STATE_CURR), HD%OtherStates(STATE_LAST), MESH_NEWCOPY, Errstat2, ErrMsg2); if (Failed()) return -!TODO -! Is there any other InitOutData should be returned -! Any additional warnings or error handling necessary - !------------------------------------------------- - ! Set output channel information for driver code - !------------------------------------------------- + !-------------------------------------------------------------------------------------------------------------------------------- + ! Set the interface meshes and outputs + !-------------------------------------------------------------------------------------------------------------------------------- + call SetMotionLoadsInterfaceMeshes(ErrStat2,ErrMsg2); if (Failed()) return + + !------------------------------------------------------------- + ! Set output channel information for driver code ! Number of channels - NumChannels_C = size(InitOutData%WriteOutputHdr) + NumChannels_C = size(SeaSt%InitOutData%WriteOutputHdr) + size(HD%InitOutData%WriteOutputHdr) ! transfer the output channel names and units to c_char arrays for returning ! Upgrade idea: use C_NULL_CHAR as delimiters. Requires rework of Python ! side of code. k=1 - do i=1,NumChannels_C + do i=1,size(SeaSt%InitOutData%WriteOutputHdr) do j=1,ChanLen ! max length of channel name. Same for units - OutputChannelNames_C(k)=InitOutData%WriteOutputHdr(i)(j:j) - OutputChannelUnits_C(k)=InitOutData%WriteOutputUnt(i)(j:j) + OutputChannelNames_C(k)=SeaSt%InitOutData%WriteOutputHdr(i)(j:j) + OutputChannelUnits_C(k)=SeaSt%InitOutData%WriteOutputUnt(i)(j:j) + k=k+1 + enddo + enddo + do i=1,size(HD%InitOutData%WriteOutputHdr) + do j=1,ChanLen ! max length of channel name. Same for units + OutputChannelNames_C(k)=HD%InitOutData%WriteOutputHdr(i)(j:j) + OutputChannelUnits_C(k)=HD%InitOutData%WriteOutputUnt(i)(j:j) k=k+1 enddo enddo @@ -526,24 +654,24 @@ subroutine SetMotionLoadsInterfaceMeshes(ErrStat3,ErrMsg3) !------------------------------------------------------------- ! Set the mapping meshes ! PRP - principle reference point - call MeshMapCreate( HD_MotionMesh, u(1)%PRPMesh, Map_Motion_2_HD_PRP_P, ErrStat3, ErrMsg3 ) + call MeshMapCreate( HD_MotionMesh, HD%u(1)%PRPMesh, Map_Motion_2_HD_PRP_P, ErrStat3, ErrMsg3 ) if (ErrStat3 >= AbortErrLev) return ! WAMIT - floating bodies using potential flow - if ( u(1)%WAMITMesh%Committed ) then ! input motions - call MeshMapCreate( HD_MotionMesh, u(1)%WAMITMesh, Map_Motion_2_HD_WB_P, ErrStat3, ErrMsg3 ) + if ( HD%u(1)%WAMITMesh%Committed ) then ! input motions + call MeshMapCreate( HD_MotionMesh, HD%u(1)%WAMITMesh, Map_Motion_2_HD_WB_P, ErrStat3, ErrMsg3 ) if (ErrStat3 >= AbortErrLev) return endif - if ( y%WAMITMesh%Committed ) then ! output loads - call MeshMapCreate( y%WAMITMesh, HD_LoadMesh, Map_HD_WB_P_2_Load, ErrStat3, ErrMsg3 ) + if ( HD%y%WAMITMesh%Committed ) then ! output loads + call MeshMapCreate( HD%y%WAMITMesh, HD_LoadMesh, Map_HD_WB_P_2_Load, ErrStat3, ErrMsg3 ) if (ErrStat3 >= AbortErrLev) return endif ! Morison - nodes for strip theory - if ( u(1)%Morison%Mesh%Committed ) then ! input motions - call MeshMapCreate( HD_MotionMesh, u(1)%Morison%Mesh, Map_Motion_2_HD_Mo_P, ErrStat3, ErrMsg3 ) + if ( HD%u(1)%Morison%Mesh%Committed ) then ! input motions + call MeshMapCreate( HD_MotionMesh, HD%u(1)%Morison%Mesh, Map_Motion_2_HD_Mo_P, ErrStat3, ErrMsg3 ) if (ErrStat3 >= AbortErrLev) return endif - if ( y%Morison%Mesh%Committed ) then ! output loads - call MeshMapCreate( y%Morison%Mesh, HD_LoadMesh, Map_HD_Mo_P_2_Load, ErrStat3, ErrMsg3 ) + if ( HD%y%Morison%Mesh%Committed ) then ! output loads + call MeshMapCreate( HD%y%Morison%Mesh, HD_LoadMesh, Map_HD_Mo_P_2_Load, ErrStat3, ErrMsg3 ) if (ErrStat3 >= AbortErrLev) return endif @@ -560,18 +688,18 @@ subroutine CheckNodes(ErrStat3,ErrMsg3) ErrStat3 = ErrID_None ErrMsg3 = "" if ( NumNodePts > 1 ) then - if ( u(1)%Morison%Mesh%Committed .and. u(1)%WAMITMesh%Committed ) then - if ( (u(1)%Morison%Mesh%Nnodes + u(1)%WAMITMesh%Nnodes) < NumNodePts ) then + if ( HD%u(1)%Morison%Mesh%Committed .and. HD%u(1)%WAMITMesh%Committed ) then + if ( (HD%u(1)%Morison%Mesh%Nnodes + HD%u(1)%WAMITMesh%Nnodes) < NumNodePts ) then ErrStat3 = ErrID_Fatal ErrMsg3 = "More nodes passed into library than exist in HydroDyn model" endif - elseif ( u(1)%Morison%Mesh%Committed ) then ! No WAMIT - if ( u(1)%Morison%Mesh%Nnodes < NumNodePts ) then + elseif ( HD%u(1)%Morison%Mesh%Committed ) then ! No WAMIT + if ( HD%u(1)%Morison%Mesh%Nnodes < NumNodePts ) then ErrStat3 = ErrID_Fatal ErrMsg3 = "More nodes passed into library than exist in HydroDyn model Morison mesh" endif - elseif ( u(1)%WAMITMesh%Committed ) then ! No Morison - if ( u(1)%WAMITMesh%Nnodes < NumNodePts ) then + elseif ( HD%u(1)%WAMITMesh%Committed ) then ! No Morison + if ( HD%u(1)%WAMITMesh%Nnodes < NumNodePts ) then ErrStat3 = ErrID_Fatal ErrMsg3 = "More nodes passed into library than exist in HydroDyn model WAMIT mesh" endif @@ -591,15 +719,15 @@ subroutine CheckDepth(ErrStat3,ErrMsg3) real(ReKi) :: tmpZpos !< temporary z-position ErrStat3 = ErrID_None ErrMsg3 = "" - tmpZpos=-0.001_ReKi*abs(p%WtrDpth) ! Initial comparison value close to surface - if ( NumNodePts == 1 .and. u(1)%Morison%Mesh%Committed ) then - do i=1,u(1)%Morison%Mesh%Nnodes + tmpZpos=-0.001_ReKi*abs(HD%p%WtrDpth) ! Initial comparison value close to surface + if ( NumNodePts == 1 .and. HD%u(1)%Morison%Mesh%Committed ) then + do i=1,HD%u(1)%Morison%Mesh%Nnodes ! Find lowest Morison node - if (u(1)%Morison%Mesh%Position(3,i) < tmpZpos) then - tmpZpos = u(1)%Morison%Mesh%Position(3,i) + if (HD%u(1)%Morison%Mesh%Position(3,i) < tmpZpos) then + tmpZpos = HD%u(1)%Morison%Mesh%Position(3,i) endif enddo - if (tmpZpos < -abs(p%WtrDpth)*0.9_ReKi) then ! within 10% of the seafloor + if (tmpZpos < -abs(HD%p%WtrDpth)*0.9_ReKi) then ! within 10% of the seafloor ErrStat3 = ErrID_Severe ErrMsg3 = "Inconsistent model"//NewLine//" -- Single library input node for simulating rigid floating structure."// & NewLine//" -- Lowest Morison node is is in lowest 10% of water depth indicating fixed bottom structure from HydroDyn."// & @@ -628,13 +756,13 @@ SUBROUTINE HydroDyn_C_CalcOutput(Time_C, NumNodePts_C, NodePos_C, NodeVel_C, Nod real(c_float), intent(in ) :: NodeVel_C( 6*NumNodePts_C ) !< A 6xNumNodePts_C array [Vx,Vy,Vz,RVx,RVy,RVz] -- velocities (global) real(c_float), intent(in ) :: NodeAcc_C( 6*NumNodePts_C ) !< A 6xNumNodePts_C array [Ax,Ay,Az,RAx,RAy,RAz] -- accelerations (global) real(c_float), intent( out) :: NodeFrc_C( 6*NumNodePts_C ) !< A 6xNumNodePts_C array [Fx,Fy,Fz,Mx,My,Mz] -- forces and moments (global) - real(c_float), intent( out) :: OutputChannelValues_C(p%NumOuts) + real(c_float), intent( out) :: OutputChannelValues_C(SeaSt%p%NumOuts+HD%p%NumOuts) integer(c_int), intent( out) :: ErrStat_C character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) ! Local variables real(DbKi) :: Time - integer(IntKi) :: iNode + integer(IntKi) :: iNode,i,k integer(IntKi) :: ErrStat !< aggregated error status character(ErrMsgLen) :: ErrMsg !< aggregated error message integer(IntKi) :: ErrStat2 !< temporary error status from a call @@ -664,17 +792,17 @@ SUBROUTINE HydroDyn_C_CalcOutput(Time_C, NumNodePts_C, NodePos_C, NodeVel_C, Nod ! Transfer motions to input meshes call Set_MotionMesh( ErrStat2, ErrMsg2 ) ! update motion mesh with input motion arrays if (Failed()) return - call HD_SetInputMotion( u(1), ErrStat2, ErrMsg2 ) ! transfer input motion mesh to u(1) meshes + call HD_SetInputMotion( HD%u(1), ErrStat2, ErrMsg2 ) ! transfer input motion mesh to u(1) meshes if (Failed()) return ! Call the main subroutine HydroDyn_CalcOutput to get the resulting forces and moments at time T - CALL HydroDyn_CalcOutput( Time, u(1), p, x(STATE_CURR), xd(STATE_CURR), z(STATE_CURR), OtherStates(STATE_CURR), y, m, ErrStat2, ErrMsg2 ) + CALL HydroDyn_CalcOutput( Time, HD%u(1), HD%p, HD%x(STATE_CURR), HD%xd(STATE_CURR), HD%z(STATE_CURR), HD%OtherStates(STATE_CURR), HD%y, HD%m, ErrStat2, ErrMsg2 ) if (Failed()) return ! Transfer resulting load meshes to intermediate mesh - call HD_TransferLoads( u(1), y, ErrStat2, ErrMsg2 ) + call HD_TransferLoads( HD%u(1), HD%y, ErrStat2, ErrMsg2 ) if (Failed()) return @@ -683,8 +811,20 @@ SUBROUTINE HydroDyn_C_CalcOutput(Time_C, NumNodePts_C, NodePos_C, NodeVel_C, Nod ! Reshape for return NodeFrc_C(1:6*NumNodePts) = reshape( real(tmpNodeFrc(1:6,1:NumNodePts), c_float), (/6*NumNodePts/) ) + + ! call SeaState to get outputs of WaveElev, etc + call SeaSt_CalcOutput( Time, SeaSt%u, SeaSt%p, SeaSt%x, SeaSt%xd, SeaSt%z, SeaSt%OtherStates, SeaSt%y, SeaSt%m, ErrStat2, ErrMsg2 ) + ! Get the output channel info out of y - OutputChannelValues_C = REAL(y%WriteOutput, C_FLOAT) + k=1 + do i=1,size(SeaSt%y%WriteOutput) + OutputChannelValues_C(k) = REAL(SeaSt%y%WriteOutput(i), C_FLOAT) + k=k+1 + enddo + do i=1,size(HD%y%WriteOutput) + OutputChannelValues_C(k) = REAL(HD%y%WriteOutput(i), C_FLOAT) + k=k+1 + enddo ! Set error status call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) @@ -746,7 +886,7 @@ SUBROUTINE HydroDyn_C_UpdateStates( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, !------------------------------------------------------- ! Check the time for current timestep and next timestep !------------------------------------------------------- - ! These inputs are used in the time stepping algorithm within HD_UpdateStates + ! These inputs are used in the time stepping algorithm within HD%UpdateStates ! For quadratic interpolation (InterpOrder==2), 3 timesteps are used. For ! linear (InterOrder==1), 2 timesteps (the HD code can handle either). ! u(1) inputs at t + dt ! Next timestep @@ -781,17 +921,17 @@ SUBROUTINE HydroDyn_C_UpdateStates( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, ! Step back to previous state because we are doing a correction step ! -- repeating the T -> T+dt update with new inputs at T+dt ! -- the STATE_CURR contains states at T+dt from the previous call, so revert those - CALL HydroDyn_CopyContState (x( STATE_LAST), x( STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyDiscState (xd( STATE_LAST), xd( STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyConstrState (z( STATE_LAST), z( STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyOtherState (OtherStates(STATE_LAST), OtherStates(STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyContState (HD%x( STATE_LAST), HD%x( STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyDiscState (HD%xd( STATE_LAST), HD%xd( STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyConstrState (HD%z( STATE_LAST), HD%z( STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyOtherState (HD%OtherStates(STATE_LAST), HD%OtherStates(STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return else ! Cycle inputs back one timestep since we are moving forward in time. if (InterpOrder>1) then ! quadratic, so keep the old time - call HydroDyn_CopyInput( u(INPUT_CURR), u(INPUT_LAST), MESH_UPDATECOPY, ErrStat2, ErrMsg2); if (Failed()) return + call HydroDyn_CopyInput( HD%u(INPUT_CURR), HD%u(INPUT_LAST), MESH_UPDATECOPY, ErrStat2, ErrMsg2); if (Failed()) return endif ! Move inputs from previous t+dt (now t) to t - call HydroDyn_CopyInput( u(INPUT_PRED), u(INPUT_CURR), MESH_UPDATECOPY, ErrStat2, ErrMsg2); if (Failed()) return + call HydroDyn_CopyInput( HD%u(INPUT_PRED), HD%u(INPUT_CURR), MESH_UPDATECOPY, ErrStat2, ErrMsg2); if (Failed()) return endif !------------------------------------------------------- @@ -805,21 +945,21 @@ SUBROUTINE HydroDyn_C_UpdateStates( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, ! Transfer motions to input meshes call Set_MotionMesh( ErrStat2, ErrMsg2 ) ! update motion mesh with input motion arrays if (Failed()) return - call HD_SetInputMotion( u(INPUT_PRED), ErrStat2, ErrMsg2 ) ! transfer input motion mesh to u(1) meshes + call HD_SetInputMotion( HD%u(INPUT_PRED), ErrStat2, ErrMsg2 ) ! transfer input motion mesh to u(1) meshes if (Failed()) return ! Set copy the current state over to the predicted state for sending to UpdateStates ! -- The STATE_PREDicted will get updated in the call. ! -- The UpdateStates routine expects this to contain states at T at the start of the call (history not passed in) - CALL HydroDyn_CopyContState (x( STATE_CURR), x( STATE_PRED), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyDiscState (xd( STATE_CURR), xd( STATE_PRED), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyConstrState (z( STATE_CURR), z( STATE_PRED), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyOtherState (OtherStates(STATE_CURR), OtherStates(STATE_PRED), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyContState (HD%x( STATE_CURR), HD%x( STATE_PRED), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyDiscState (HD%xd( STATE_CURR), HD%xd( STATE_PRED), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyConstrState (HD%z( STATE_CURR), HD%z( STATE_PRED), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyOtherState (HD%OtherStates(STATE_CURR), HD%OtherStates(STATE_PRED), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return ! Call the main subroutine HydroDyn_UpdateStates to get the velocities - CALL HydroDyn_UpdateStates( InputTimes(INPUT_CURR), N_Global, u, InputTimes, p, x(STATE_PRED), xd(STATE_PRED), z(STATE_PRED), OtherStates(STATE_PRED), m, ErrStat2, ErrMsg2 ) + CALL HydroDyn_UpdateStates( InputTimes(INPUT_CURR), N_Global, HD%u, InputTimes, HD%p, HD%x(STATE_PRED), HD%xd(STATE_PRED), HD%z(STATE_PRED), HD%OtherStates(STATE_PRED), HD%m, ErrStat2, ErrMsg2 ) if (Failed()) return @@ -829,16 +969,16 @@ SUBROUTINE HydroDyn_C_UpdateStates( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, ! move current state at T to previous state at T-dt ! -- STATE_LAST now contains info at time T ! -- this allows repeating the T --> T+dt update - CALL HydroDyn_CopyContState (x( STATE_CURR), x( STATE_LAST), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyDiscState (xd( STATE_CURR), xd( STATE_LAST), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyConstrState (z( STATE_CURR), z( STATE_LAST), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyOtherState (OtherStates(STATE_CURR), OtherStates(STATE_LAST), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyContState (HD%x( STATE_CURR), HD%x( STATE_LAST), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyDiscState (HD%xd( STATE_CURR), HD%xd( STATE_LAST), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyConstrState (HD%z( STATE_CURR), HD%z( STATE_LAST), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyOtherState (HD%OtherStates(STATE_CURR), HD%OtherStates(STATE_LAST), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return ! Update the predicted state as the new current state ! -- we have now advanced from T to T+dt. This allows calling with CalcOuput to get the outputs at T+dt - CALL HydroDyn_CopyContState (x( STATE_PRED), x( STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyDiscState (xd( STATE_PRED), xd( STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyConstrState (z( STATE_PRED), z( STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return - CALL HydroDyn_CopyOtherState (OtherStates(STATE_PRED), OtherStates(STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyContState (HD%x( STATE_PRED), HD%x( STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyDiscState (HD%xd( STATE_PRED), HD%xd( STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyConstrState (HD%z( STATE_PRED), HD%z( STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + CALL HydroDyn_CopyOtherState (HD%OtherStates(STATE_PRED), HD%OtherStates(STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return @@ -887,10 +1027,10 @@ SUBROUTINE HydroDyn_C_End(ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_C_End') ! Call the main subroutine HydroDyn_End ! If u is not allocated, then we didn't get far at all in initialization, - ! or HD_C_End got called before Init. We don't want a segfault, so check + ! or HD%C_End got called before Init. We don't want a segfault, so check ! for allocation. - if (allocated(u)) then - call HydroDyn_End( u(1), p, x(STATE_CURR), xd(STATE_CURR), z(STATE_CURR), OtherStates(STATE_CURR), y, m, ErrStat2, ErrMsg2 ) + if (allocated(HD%u)) then + call HydroDyn_End( HD%u(1), HD%p, HD%x(STATE_CURR), HD%xd(STATE_CURR), HD%z(STATE_CURR), HD%OtherStates(STATE_CURR), HD%y, HD%m, ErrStat2, ErrMsg2 ) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) endif @@ -898,27 +1038,38 @@ SUBROUTINE HydroDyn_C_End(ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_C_End') ! logic is required here (this isn't necessary in the fortran driver ! or in openfast, but may be when this code is called from C, Python, ! or some other code using the c-bindings. - if (allocated(u)) then - do i=2,size(u) - call HydroDyn_DestroyInput( u(i), ErrStat2, ErrMsg2 ) + if (allocated(HD%u)) then + do i=2,size(HD%u) + call HydroDyn_DestroyInput( HD%u(i), ErrStat2, ErrMsg2 ) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) enddo - if (allocated(u)) deallocate(u) + if (allocated(HD%u)) deallocate(HD%u) endif ! Destroy any other copies of states (rerun on (STATE_CURR) is ok) - call HydroDyn_DestroyContState( x( STATE_LAST), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - call HydroDyn_DestroyDiscState( xd( STATE_LAST), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - call HydroDyn_DestroyConstrState( z( STATE_LAST), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - call HydroDyn_DestroyOtherState( OtherStates(STATE_LAST), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - call HydroDyn_DestroyContState( x( STATE_CURR), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - call HydroDyn_DestroyDiscState( xd( STATE_CURR), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - call HydroDyn_DestroyConstrState( z( STATE_CURR), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - call HydroDyn_DestroyOtherState( OtherStates(STATE_CURR), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - call HydroDyn_DestroyContState( x( STATE_PRED), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - call HydroDyn_DestroyDiscState( xd( STATE_PRED), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - call HydroDyn_DestroyConstrState( z( STATE_PRED), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - call HydroDyn_DestroyOtherState( OtherStates(STATE_PRED), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyContState( HD%x( STATE_LAST), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyDiscState( HD%xd( STATE_LAST), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyConstrState( HD%z( STATE_LAST), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyOtherState( HD%OtherStates(STATE_LAST), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyContState( HD%x( STATE_CURR), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyDiscState( HD%xd( STATE_CURR), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyConstrState( HD%z( STATE_CURR), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyOtherState( HD%OtherStates(STATE_CURR), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyContState( HD%x( STATE_PRED), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyDiscState( HD%xd( STATE_PRED), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyConstrState( HD%z( STATE_PRED), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyOtherState( HD%OtherStates(STATE_PRED), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + + + ! Call the main subroutine SeaSt_End + call SeaSt_End( SeaSt%u, SeaSt%p, SeaSt%x, SeaSt%xd, SeaSt%z, SeaSt%OtherStates, SeaSt%y, SeaSt%m, ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + + ! Destroy any other copies of states (rerun on (STATE_CURR) is ok) + call SeaSt_DestroyContState( SeaSt%x , ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call SeaSt_DestroyDiscState( SeaSt%xd , ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call SeaSt_DestroyConstrState( SeaSt%z , ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call SeaSt_DestroyOtherState( SeaSt%OtherStates, ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ! if deallocate other items now @@ -973,21 +1124,21 @@ end subroutine Set_MotionMesh !> Map the motion of the intermediate input mesh over to the input meshes !! This routine is operating on module level data, hence few inputs -subroutine HD_SetInputMotion( u_local, ErrStat3, ErrMsg3 ) - type(HydroDyn_InputType), intent(inout) :: u_local ! Only one input (probably at T) +subroutine HD_SetInputMotion( HD_u_local, ErrStat3, ErrMsg3 ) + type(HydroDyn_InputType), intent(inout) :: HD_u_local ! Only one input (probably at T) integer(IntKi), intent( out) :: ErrStat3 character(ErrMsgLen), intent( out) :: ErrMsg3 ! Principle reference point - CALL Transfer_Point_to_Point( HD_MotionMesh, u_local%PRPMesh, Map_Motion_2_HD_PRP_P, ErrStat3, ErrMsg3 ) + CALL Transfer_Point_to_Point( HD_MotionMesh, HD_u_local%PRPMesh, Map_Motion_2_HD_PRP_P, ErrStat3, ErrMsg3 ) if (ErrStat3 >= AbortErrLev) return ! WAMIT mesh - if ( u_local%WAMITMesh%Committed ) then - call Transfer_Point_to_Point( HD_MotionMesh, u_local%WAMITMesh, Map_Motion_2_HD_WB_P, ErrStat3, ErrMsg3 ) + if ( HD_u_local%WAMITMesh%Committed ) then + call Transfer_Point_to_Point( HD_MotionMesh, HD_u_local%WAMITMesh, Map_Motion_2_HD_WB_P, ErrStat3, ErrMsg3 ) if (ErrStat3 >= AbortErrLev) return endif ! Morison mesh - if ( u_local%Morison%Mesh%Committed ) then - call Transfer_Point_to_Point( HD_MotionMesh, u_local%Morison%Mesh, Map_Motion_2_HD_Mo_P, ErrStat3, ErrMsg3 ) + if ( HD_u_local%Morison%Mesh%Committed ) then + call Transfer_Point_to_Point( HD_MotionMesh, HD_u_local%Morison%Mesh, Map_Motion_2_HD_Mo_P, ErrStat3, ErrMsg3 ) if (ErrStat3 >= AbortErrLev) return endif end subroutine HD_SetInputMotion @@ -998,9 +1149,9 @@ end subroutine HD_SetInputMotion !! temporary mesh -- prevents accidental overwrite of WAMIT loads on HD_LoadMesh !! with the mapping of the Morison loads. !! This routine is operating on module level data, hence few inputs -subroutine HD_TransferLoads( u_local, y_local, ErrStat3, ErrMsg3 ) - type(HydroDyn_InputType), intent(in ) :: u_local ! Only one input (probably at T) - type(HydroDyn_OutputType), intent(in ) :: y_local ! Only one input (probably at T) +subroutine HD_TransferLoads( HD_u_local, HD_y_local, ErrStat3, ErrMsg3 ) + type(HydroDyn_InputType), intent(in ) :: HD_u_local ! Only one input (probably at T) + type(HydroDyn_OutputType), intent(in ) :: HD_y_local ! Only one input (probably at T) integer(IntKi), intent( out) :: ErrStat3 character(ErrMsgLen), intent( out) :: ErrMsg3 @@ -1008,19 +1159,19 @@ subroutine HD_TransferLoads( u_local, y_local, ErrStat3, ErrMsg3 ) HD_LoadMesh%Moment = 0.0_ReKi ! WAMIT mesh - if ( y_local%WAMITMesh%Committed ) then + if ( HD_y_local%WAMITMesh%Committed ) then HD_LoadMesh_tmp%Force = 0.0_ReKi HD_LoadMesh_tmp%Moment = 0.0_ReKi - call Transfer_Point_to_Point( y_local%WAMITMesh, HD_LoadMesh_tmp, Map_HD_WB_P_2_Load, ErrStat3, ErrMsg3, u_local%WAMITMesh, HD_MotionMesh ) + call Transfer_Point_to_Point( HD_y_local%WAMITMesh, HD_LoadMesh_tmp, Map_HD_WB_P_2_Load, ErrStat3, ErrMsg3, HD_u_local%WAMITMesh, HD_MotionMesh ) if (ErrStat3 >= AbortErrLev) return HD_LoadMesh%Force = HD_LoadMesh%Force + HD_LoadMesh_tmp%Force HD_LoadMesh%Moment = HD_LoadMesh%Moment + HD_LoadMesh_tmp%Moment endif ! Morison mesh - if ( y_local%Morison%Mesh%Committed ) then + if ( HD_y_local%Morison%Mesh%Committed ) then HD_LoadMesh_tmp%Force = 0.0_ReKi HD_LoadMesh_tmp%Moment = 0.0_ReKi - call Transfer_Point_to_Point( y_local%Morison%Mesh, HD_LoadMesh_tmp, Map_HD_Mo_P_2_Load, ErrStat3, ErrMsg3, u_local%Morison%Mesh, HD_MotionMesh ) + call Transfer_Point_to_Point( HD_y_local%Morison%Mesh, HD_LoadMesh_tmp, Map_HD_Mo_P_2_Load, ErrStat3, ErrMsg3, HD_u_local%Morison%Mesh, HD_MotionMesh ) if (ErrStat3 >= AbortErrLev) return HD_LoadMesh%Force = HD_LoadMesh%Force + HD_LoadMesh_tmp%Force HD_LoadMesh%Moment = HD_LoadMesh%Moment + HD_LoadMesh_tmp%Moment diff --git a/modules/hydrodyn/src/HydroDyn_Input.f90 b/modules/hydrodyn/src/HydroDyn_Input.f90 index c059a92cd3..1f39d2c7fe 100644 --- a/modules/hydrodyn/src/HydroDyn_Input.f90 +++ b/modules/hydrodyn/src/HydroDyn_Input.f90 @@ -1152,7 +1152,12 @@ SUBROUTINE HydroDynInput_ProcessInitData( InitInp, Interval, InputFileData, ErrS END IF - ! WaveMod - Wave kinematics model switch. + + ! WaveMod - Wave kinematics model switch. -- Check that actual data was passed in from SeaState. If none exists, then set WaveMod=0 and warn + if (.not. associated(InitInp%WaveTime) .or. InitInp%NStepWave == 0) then + call SetErrStat( ErrID_Fatal,' No SeaState wave information available. Setting WaveMod=0.',ErrStat,ErrMsg,RoutineName) + return + endif IF ( InputFileData%PotMod > 0 .and. InitInp%WaveMod == 6 ) THEN CALL SetErrStat( ErrID_Fatal,'WaveMod must be 0, 1, 1P#, 2, 3, 4, or 5 when PotMod is not 0',ErrStat,ErrMsg,RoutineName) diff --git a/modules/hydrodyn/src/HydroDyn_Types.f90 b/modules/hydrodyn/src/HydroDyn_Types.f90 index 41b2a62fbe..af79c5895d 100644 --- a/modules/hydrodyn/src/HydroDyn_Types.f90 +++ b/modules/hydrodyn/src/HydroDyn_Types.f90 @@ -92,8 +92,8 @@ MODULE HydroDyn_Types REAL(DbKi) :: TMax !< Supplied by Driver: The total simulation time [(sec)] REAL(ReKi) :: PtfmLocationX !< Supplied by Driver: X coordinate of platform location in the wave field [m] REAL(ReKi) :: PtfmLocationY !< Supplied by Driver: Y coordinate of platform location in the wave field [m] - INTEGER(IntKi) :: NStepWave !< Total number of frequency components = total number of time steps in the incident wave [-] - INTEGER(IntKi) :: NStepWave2 !< NStepWave / 2 [-] + INTEGER(IntKi) :: NStepWave = 0 !< Total number of frequency components = total number of time steps in the incident wave [-] + INTEGER(IntKi) :: NStepWave2 = 0 !< NStepWave / 2 [-] REAL(SiKi) :: RhoXg !< = WtrDens*Gravity [-] INTEGER(IntKi) :: WaveMod !< Incident wave kinematics model {0: none=still water, 1: plane progressive (regular), 2: JONSWAP/Pierson-Moskowitz spectrum (irregular), 3: white-noise spectrum, 4: user-defind spectrum from routine UserWaveSpctrm (irregular), 5: GH BLADED } [-] INTEGER(IntKi) :: WaveStMod !< Model for stretching incident wave kinematics to instantaneous free surface {0: none=no stretching, 1: vertical stretching, 2: extrapolation stretching, 3: Wheeler stretching} [-] diff --git a/reg_tests/executeHydrodynPyRegressionCase.py b/reg_tests/executeHydrodynPyRegressionCase.py index 77044846a4..241094a1ff 100644 --- a/reg_tests/executeHydrodynPyRegressionCase.py +++ b/reg_tests/executeHydrodynPyRegressionCase.py @@ -88,24 +88,15 @@ if not os.path.isdir(inputsDirectory): rtl.exitWithError("The test data inputs directory, {}, does not exist. Verify your local repository is up to date.".format(inputsDirectory)) -# create the local output directory if it does not already exist -# and initialize it with input files for all test cases -if not os.path.isdir(testBuildDirectory): - os.makedirs(testBuildDirectory) - for file in glob.glob(os.path.join(inputsDirectory,"*py")): - filename = file.split(os.path.sep)[-1] - shutil.copy(os.path.join(inputsDirectory,filename), os.path.join(testBuildDirectory,filename)) - for file in glob.glob(os.path.join(inputsDirectory,"hd_*inp")): - filename = file.split(os.path.sep)[-1] - shutil.copy(os.path.join(inputsDirectory,filename), os.path.join(testBuildDirectory,filename)) - for file in glob.glob(os.path.join(inputsDirectory,"*dat")): - filename = file.split(os.path.sep)[-1] - shutil.copy(os.path.join(inputsDirectory,filename), os.path.join(testBuildDirectory,filename)) +# create the local output directory and initialize it with input files +rtl.copyTree(inputsDirectory, testBuildDirectory) + # , excludeExt=['.out','.outb']) + ### Run HydroDyn on the test case if not noExec: caseInputFile = os.path.join(testBuildDirectory, "hydrodyn_driver.py") - returnCode = openfastDrivers.runHydrodynDriverCase(caseInputFile, executable) + returnCode = openfastDrivers.runHydrodynDriverCase(caseInputFile, executable, verbose=verbose) if returnCode != 0: sys.exit(returnCode*10) diff --git a/reg_tests/r-test b/reg_tests/r-test index 1ab766b7b1..b5bd1b547c 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 1ab766b7b1ce35217ae6092ae524753648b5790b +Subproject commit b5bd1b547c8057776796ba55a71ea917484a4b25