From e5762db97c544dda354c063c6860d14ce0fa99ba Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Mon, 17 May 2021 11:29:23 -0600 Subject: [PATCH 01/22] hd_c_lib: placeholder interface and rough-in successfully can call HD_Init and HD_End. Issues: - error handling is lacking - inflexible file passing (set number of lines and line length -- this must be generalized) - setup reading of input files to pass - setup driver code to use the existing driver input files - setup test cases that mirror existing driver test cases - add the calcoutput and updatestates for single point rigid -- extend to flexible later. --- modules/hydrodyn/CMakeLists.txt | 6 +- .../hydrodyn/python-lib/hydrodyn_library.py | 182 +++++++++++ modules/hydrodyn/src/HydroDyn_C.f90 | 293 ++++++++++++++++++ reg_tests/r-test | 2 +- 4 files changed, 481 insertions(+), 2 deletions(-) create mode 100644 modules/hydrodyn/python-lib/hydrodyn_library.py create mode 100644 modules/hydrodyn/src/HydroDyn_C.f90 diff --git a/modules/hydrodyn/CMakeLists.txt b/modules/hydrodyn/CMakeLists.txt index 3c10213253..36ba07451f 100644 --- a/modules/hydrodyn/CMakeLists.txt +++ b/modules/hydrodyn/CMakeLists.txt @@ -59,6 +59,10 @@ set(HYDRODYN_SOURCES add_library(hydrodynlib ${HYDRODYN_SOURCES}) target_link_libraries(hydrodynlib nwtclibs) +# c-bindings interface library +add_library(hydrodyn_c_lib SHARED src/HydroDyn_C.f90) +target_link_libraries(hydrodyn_c_lib hydrodynlib) + add_executable(hydrodyn_driver src/HydroDyn_DriverCode.f90) target_link_libraries(hydrodyn_driver hydrodynlib nwtclibs versioninfolib) @@ -66,7 +70,7 @@ target_link_libraries(hydrodyn_driver hydrodynlib nwtclibs versioninfolib) # src/SS_Radiation_DriverCode.f90) #target_link_libraries(ss_radiation hydrodynlib nwtclibs) -install(TARGETS hydrodynlib hydrodyn_driver +install(TARGETS hydrodynlib hydrodyn_driver hydrodyn_c_lib EXPORT "${CMAKE_PROJECT_NAME}Libraries" RUNTIME DESTINATION bin LIBRARY DESTINATION lib diff --git a/modules/hydrodyn/python-lib/hydrodyn_library.py b/modules/hydrodyn/python-lib/hydrodyn_library.py new file mode 100644 index 0000000000..61f52d4eb5 --- /dev/null +++ b/modules/hydrodyn/python-lib/hydrodyn_library.py @@ -0,0 +1,182 @@ +#********************************************************************************************************************************** +# LICENSING +# Copyright (C) 2021 NREL +# +# This file is part of InflowWind. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +#********************************************************************************************************************************** +# +# This is the Python-C interface library for HydroDyn +# Usage: THIS LIBRARY IS NOT TO BE CHANGED OR EDITED BY THE USER +from ctypes import ( + CDLL, + POINTER, + create_string_buffer, + byref, + c_byte, + c_int, + c_double, + c_float, + c_char, + c_char_p, + c_wchar, + c_wchar_p, + c_bool +) +import numpy as np + +class HydroDynLibAPI(CDLL): + def __init__(self, library_path): + super().__init__(library_path) + self.library_path = library_path + + self._initialize_routines() + + # Create buffers for class data + self.abort_error_level = c_int(4) + self.error_status = c_int(0) + self.error_message = create_string_buffer(1025) + +# This is not sufficient for HD + self._channel_names = create_string_buffer(20 * 4000) + self._channel_units = create_string_buffer(20 * 4000) + + self.dt = c_double(0) + self.total_time = c_double(0) + self.numTimeSteps = c_int(0) + + # _initialize_routines() ------------------------------------------------------------------------------------------------------------ + def _initialize_routines(self): + self.HydroDyn_Init_C.argtypes = [ + POINTER(c_char_p), # input file string + POINTER(c_int), # input file string length + POINTER(c_double), # dt + POINTER(c_int), # number of channels + POINTER(c_char), # output channel names + POINTER(c_char), # output channel units + POINTER(c_int), # ErrStat_C + POINTER(c_char) # ErrMsg_C + ] + self.HydroDyn_Init_C.restype = c_int + + self.HydroDyn_CalcOutput_C.argtypes = [ + POINTER(c_double), # Time_C + POINTER(c_float), # Output Channel Values + POINTER(c_int), # ErrStat_C + POINTER(c_char) # ErrMsg_C + ] + self.HydroDyn_CalcOutput_C.restype = c_int + + self.HydroDyn_End_C.argtypes = [ + POINTER(c_int), # ErrStat_C + POINTER(c_char) # ErrMsg_C + ] + self.HydroDyn_End_C.restype = c_int + + # hydrodyn_init ------------------------------------------------------------------------------------------------------------ + def hydrodyn_init(self, input_strings, input_string_length): + + #print('hydrodyn_library.py: Running HydroDyn_Init_C .....') + + # Set up inputs + input_string_array = (c_char_p * len(input_strings))() + for i, param in enumerate(input_strings): + input_string_array[i] = param.encode('utf-8') + + self._numChannels = c_int(0) + + # Run HydroDyn_Init_C + self.HydroDyn_Init_C( + input_string_array, # IN: input file string + byref(c_int(input_string_length)), # IN: input file string length + byref(c_double(self.dt)), # IN: time step (dt) + byref(self._numChannels), # OUT: number of channels + self._channel_names, # OUT: output channel names + self._channel_units, # OUT: output channel units + byref(self.error_status), # OUT: ErrStat_C + self.error_message # OUT: ErrMsg_C + ) + if self.fatal_error: + print(f"Error {self.error_status.value}: {self.error_message.value}") + return + + # Initialize output channels + self._channel_output_array = (c_double * self._numChannels.value)(0.0, ) + self._channel_output_values = np.empty( (self.numTimeSteps, self._numChannels.value) ) + + #print('hydrodyn_library.py: Completed HydroDyn_Init_C') + + # hydrodyn_calcOutput ------------------------------------------------------------------------------------------------------------ + def hydrodyn_calcOutput(self, time, positions, velocities, outputChannelValues): + + #print('hydrodyn_library.py: Running HydroDyn_CalcOutput_C .....') + + # Set up inputs +# positions_flat = [pp for p in positions for pp in p] # need to flatten to pass through to Fortran (to reshape) +# positions_flat_c = (c_float * (3 * self.numWindPts))(0.0, ) +# for i, p in enumerate(positions_flat): +# positions_flat_c[i] = c_float(p) +# +# velocities_flat_c = (c_float * (3 * self.numWindPts))(0.0, ) + + outputChannelValues_c = (c_float * self._numChannels.value)(0.0, ) + + # Run HydroDyn_CalcOutput_C + self.HydroDyn_CalcOutput_C( + byref(c_double(time)), # IN: time at which to calculate velocities +# positions_flat_c, # IN: positions - specified by user +# velocities_flat_c, # OUT: velocities at desired positions + outputChannelValues_c, # OUT: output channel values as described in input file + byref(self.error_status), # OUT: ErrStat_C + self.error_message # OUT: ErrMsg_C + ) + if self.fatal_error: + print(f"Error {self.error_status.value}: {self.error_message.value}") + return + + # Convert output channel values back into python + for k in range(0,self._numChannels.value): + outputChannelValues[k] = float(outputChannelValues_c[k]) + + # Reshape velocities into [N,3] + count = 0 + for j in range(0,self.numWindPts): + velocities[j,0] = velocities_flat_c[count] + velocities[j,1] = velocities_flat_c[count+1] + velocities[j,2] = velocities_flat_c[count+2] + count = count + 3 + + #print('hydrodyn_library.py: Completed HydroDyn_CalcOutput_C') + + # hydrodyn_end ------------------------------------------------------------------------------------------------------------ + def hydrodyn_end(self): + + #print('hydrodyn_library.py: Running HydroDyn_End_C .....') + + # Run HydroDyn_End_C + self.HydroDyn_End_C( + byref(self.error_status), + self.error_message + ) + if self.fatal_error: + print(f"Error {self.error_status.value}: {self.error_message.value}") + return + + #print('hydrodyn_library.py: Completed HydroDyn_End_C') + + # other functions ---------------------------------------------------------------------------------------------------------- + @property + def fatal_error(self): + return self.error_status.value >= self.abort_error_level.value diff --git a/modules/hydrodyn/src/HydroDyn_C.f90 b/modules/hydrodyn/src/HydroDyn_C.f90 new file mode 100644 index 0000000000..6c1b6d8024 --- /dev/null +++ b/modules/hydrodyn/src/HydroDyn_C.f90 @@ -0,0 +1,293 @@ +!********************************************************************************************************************************** +! LICENSING +! Copyright (C) 2021 Nicole Mendoza +! +! This file is part of HydroDyn. +! +! Licensed under the Apache License, Version 2.0 (the "License"); +! you may not use this file except in compliance with the License. +! You may obtain a copy of the License at +! +! http://www.apache.org/licenses/LICENSE-2.0 +! +! Unless required by applicable law or agreed to in writing, software +! distributed under the License is distributed on an "AS IS" BASIS, +! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +! See the License for the specific language governing permissions and +! limitations under the License. +! +!********************************************************************************************************************************** +MODULE HydroDynAPI + + USE ISO_C_BINDING + USE HydroDyn + USE HydroDyn_Types + USE NWTC_Library + +IMPLICIT NONE + +PUBLIC :: HydroDyn_Init_C +PUBLIC :: HydroDyn_CalcOutput_C +PUBLIC :: HydroDyn_UpdateStates_C +PUBLIC :: HydroDyn_End_C + +! Accessible to all routines inside module +TYPE(HydroDyn_InputType) :: InputGuess !< An initial guess for the input; the input mesh must be defined, returned by Init +TYPE(HydroDyn_InputType) :: InputData !< Created by IFW_CALCOUTPUT_C and used by IFW_END_C +TYPE(HydroDyn_InitInputType) :: InitInp +TYPE(HydroDyn_InitOutputType) :: InitOutData !< Initial output data -- Names, units, and version info. +TYPE(HydroDyn_ParameterType) :: p !< Parameters +TYPE(HydroDyn_ContinuousStateType) :: ContStates !< Initial continuous states +TYPE(HydroDyn_DiscreteStateType) :: DiscStates !< Initial discrete states +TYPE(HydroDyn_ConstraintStateType) :: ConstrStateGuess !< Initial guess of the constraint states +TYPE(HydroDyn_ConstraintStateType) :: ConstrStates !< Constraint states at Time +TYPE(HydroDyn_OtherStateType) :: OtherStates !< 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) + +!FIXME: we cannot have this fixed!!!!! +INTEGER, PARAMETER :: InputStringLength = 647 !< Fixed length for all lines of the string-based input file + +CONTAINS + +!=============================================================================================================== +!--------------------------------------------- HydroDyn Init---------------------------------------------------- +!=============================================================================================================== +SUBROUTINE HydroDyn_Init_C(InputFileStrings_C, InputFileStringLength_C, DT_C, NumChannels_C, OutputChannelNames_C, OutputChannelUnits_C, ErrStat_C, ErrMsg_C) BIND (C, NAME='HydroDyn_Init_C') + + TYPE(C_PTR) , INTENT(IN ) :: InputFileStrings_C + INTEGER(C_INT) , INTENT(IN ) :: InputFileStringLength_C + REAL(C_DOUBLE) , INTENT(IN ) :: DT_C + INTEGER(C_INT) , INTENT( OUT) :: NumChannels_C + TYPE(C_PTR) , INTENT( OUT) :: OutputChannelNames_C + TYPE(C_PTR) , INTENT( OUT) :: OutputChannelUnits_C + INTEGER(C_INT) , INTENT( OUT) :: ErrStat_C + CHARACTER(KIND=C_CHAR) , INTENT( OUT) :: ErrMsg_C(1025) + + ! Local Variables + CHARACTER(InputStringLength), DIMENSION(InputFileStringLength_C) :: InputFileStrings + CHARACTER(kind=C_char, len=1), DIMENSION(:), POINTER :: character_pointer + CHARACTER(kind=C_char, len=1), DIMENSION(:), POINTER :: character_pointer2 + CHARACTER, DIMENSION(InputStringLength) :: single_line_character_array + CHARACTER, DIMENSION(InputStringLength) :: single_line_character_array2 + CHARACTER(InputStringLength) :: single_line_chars + CHARACTER(InputStringLength) :: single_line_chars2 + + CHARACTER(CHANLEN+1), ALLOCATABLE, TARGET :: tmp_OutputChannelNames_C(:) + CHARACTER(CHANLEN+1), ALLOCATABLE, TARGET :: tmp_OutputChannelUnits_C(:) + REAL(DbKi) :: TimeInterval + INTEGER :: ErrStat !< aggregated error message + CHARACTER(ErrMsgLen) :: ErrMsg !< aggregated error message + INTEGER :: ErrStat2 !< temporary error status from a call + CHARACTER(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call + INTEGER :: I, J, K + character(*), parameter :: RoutineName = 'HydroDyn_Init_C' !< for error handling + + ! Initialize error handling + ErrStat = ErrID_None + ErrMsg = "" + + ! Convert the string-input from C-style character arrays (char**) to a Fortran-style array of characters. + ! TODO: Add error checking + CALL C_F_pointer(InputFileStrings_C, character_pointer, [InputFileStringLength_C * InputStringLength]) + DO i = 0, InputFileStringLength_C - 1 + single_line_character_array = character_pointer(i * InputStringLength + 1 : i * InputStringLength + InputStringLength) + DO j = 1, InputStringLength + single_line_chars(j:j) = single_line_character_array(j) + END DO + InputFileStrings(i + 1) = single_line_chars + END DO + + ! Store string-inputs as type FileInfoType within HydroDyn_InitInputType + CALL InitFileInfo(InputFileStrings, InitInp%PassedFileData, ErrStat2, ErrMsg2) + if (Failed()) return + + ! Set other inputs for calling HydroDyn_Init + InitInp%InputFile = "passed_hd_file" ! dummy + InitInp%OutRootName = "HDpassed" ! used for making echo files + InitInp%UseInputFile = .FALSE. + TimeInterval = REAL(DT_C, DbKi) + InitInp%Linearize = .FALSE. + InitInp%Gravity = 9.80665 ! Gravity (m/s^2) + InitInp%defWtrDens = 1025 ! Water density (kg/m^3) + InitInp%defWtrDpth = 200 ! Water depth (m) + InitInp%defMSL2SWL = 0 ! Offset between still-water level and mean sea level (m) [positive upward] + InitInp%TMax = 600 + InitInp%HasIce = .FALSE. +! InitInp%WaveElevXY ! Don't allocate this + InitInp%PtfmLocationX = 0.0_ReKi + InitInp%PtfmLocationY = 0.0_ReKi + + + ! Call the main subroutine HydroDyn_Init - only need InitInp and TimeInterval as inputs, the rest are set by HydroDyn_Init + CALL HydroDyn_Init( InitInp, InputGuess, p, ContStates, DiscStates, ConstrStateGuess, OtherStates, y, m, TimeInterval, InitOutData, ErrStat2, ErrMsg2 ) + if (Failed()) return + + ! Convert the outputs of HydroDyn_Init from Fortran to C + ALLOCATE(tmp_OutputChannelNames_C(size(InitOutData%WriteOutputHdr)),STAT=ErrStat2) + if (ErrStat2 /= 0) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "Could not allocate WriteOutputHdr array" + endif + if (Failed()) return + ALLOCATE(tmp_OutputChannelUnits_C(size(InitOutData%WriteOutputUnt)),STAT=ErrStat2) + if (ErrStat2 /= 0) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "Could not allocate WriteOutputUnt array" + endif + if (Failed()) return + NumChannels_C = size(InitOutData%WriteOutputHdr) + + DO I = 1,NumChannels_C + tmp_OutputChannelNames_C(I) = TRANSFER(InitOutData%WriteOutputHdr(I)//C_NULL_CHAR, tmp_OutputChannelNames_C(I)) + tmp_OutputChannelUnits_C(I) = TRANSFER(InitOutData%WriteOutputUnt(I)//C_NULL_CHAR, tmp_OutputChannelUnits_C(I)) + END DO + OutputChannelNames_C = C_LOC(tmp_OutputChannelNames_C) + OutputChannelUnits_C = C_LOC(tmp_OutputChannelUnits_C) + + ! Clean up variables and set up for HydroDyn_CalcOutput_C + CALL HydroDyn_CopyInput(InputGuess, InputData, MESH_NEWCOPY, ErrStat2, ErrMsg2 ) + if (Failed()) return + CALL HydroDyn_CopyConstrState(ConstrStateGuess, ConstrStates, MESH_NEWCOPY, ErrStat2, ErrMsg2 ) + if (Failed()) return + + call Cleanup() + call SetErr() + +CONTAINS + logical function Failed() + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + if (Failed) then + call Cleanup() + call SetErr() + endif + end function Failed + subroutine Cleanup() ! NOTE: we are ignoring any error reporting from here + CALL HydroDyn_DestroyInput(InputGuess, ErrStat2, ErrMsg2 ) + CALL HydroDyn_DestroyConstrState(ConstrStateGuess, ErrStat2, ErrMsg2 ) + end subroutine Cleanup + subroutine SetErr() + ErrStat_C = ErrStat ! We will send back the same error status that is used in OpenFAST + ErrMsg_C = TRANSFER( trim(ErrMsg)//C_NULL_CHAR, ErrMsg_C ) + end subroutine SetErr + +END SUBROUTINE HydroDyn_Init_C + +!=============================================================================================================== +!--------------------------------------------- HydroDyn CalcOutput --------------------------------------------- +!=============================================================================================================== + +SUBROUTINE HydroDyn_CalcOutput_C(Time_C,OutputChannelValues_C,ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_CalcOutput_C') + REAL(C_DOUBLE) , INTENT(IN ) :: Time_C + REAL(C_FLOAT) , INTENT( OUT) :: OutputChannelValues_C(p%NumOuts) + INTEGER(C_INT) , INTENT( OUT) :: ErrStat_C + CHARACTER(KIND=C_CHAR) , INTENT( OUT) :: ErrMsg_C + + ! Local variables + REAL(DbKi) :: Time + INTEGER :: ErrStat !< aggregated error message + 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 = 'HydroDyn_CalcOutput_C' !< for error handling + + ! Initialize error handling + ErrStat = ErrID_None + ErrMsg = "" + + ! Convert the inputs from C to Fortrn + Time = REAL(Time_C,DbKi) +! InputData%PositionXYZ = reshape( real(Positions_C,ReKi), (/3, InitInp%NumWindPoints/) ) + +! ! Call the main subroutine HydroDyn_CalcOutput to get the velocities +! CALL HydroDyn_CalcOutput( Time, InputData, p, ContStates, DiscStates, ConstrStates, OtherStates, y, m, ErrStat2, ErrMsg2 ) +! if (Failed()) return + +! ! Get velocities out of y and flatten them (still in same spot in memory) +! Velocities_C = reshape( REAL(y%VelocityUVW, C_FLOAT), (/3*InitInp%NumWindPoints/) ) ! VelocityUVW is 2D array of ReKi (might need reshape or make into pointer); size [3,N] + +! ! Get the output channel info out of y +! OutputChannelValues_C = REAL(y%WriteOutput, C_FLOAT) + + call SetErr() + +CONTAINS + logical function Failed() + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + if (Failed) call SetErr() + end function Failed + subroutine SetErr() + ErrStat_C = ErrStat ! We will send back the same error status that is used in OpenFAST + ErrMsg_C = TRANSFER( trim(ErrMsg)//C_NULL_CHAR, ErrMsg_C ) + end subroutine SetErr +END SUBROUTINE HydroDyn_CalcOutput_C + +!=============================================================================================================== +!--------------------------------------------- HydroDyn UpdateStates --------------------------------------------- +!=============================================================================================================== + +SUBROUTINE HydroDyn_UpdateStates_C(Time_C,ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_UpdateStates_C') + REAL(C_DOUBLE) , INTENT(IN ) :: Time_C + INTEGER(C_INT) , INTENT( OUT) :: ErrStat_C + CHARACTER(KIND=C_CHAR) , INTENT( OUT) :: ErrMsg_C + ! Local variables + REAL(DbKi) :: Time + INTEGER :: ErrStat !< aggregated error message + 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 = 'HydroDyn_UpdateStates_C' !< for error handling + + ! Initialize error handling + ErrStat = ErrID_None + ErrMsg = "" + + ! Convert the inputs from C to Fortran + Time = REAL(Time_C,DbKi) + +! ! Call the main subroutine HydroDyn_UpdateStates to get the velocities +! CALL HydroDyn_UpdateStates( Time, InputData, p, ContStates, DiscStates, ConstrStates, OtherStates, y, m, ErrStat2, ErrMsg2 ) +! if (Failed()) return + + call SetErr() + +CONTAINS + logical function Failed() + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + if (Failed) call SetErr() + end function Failed + subroutine SetErr() + ErrStat_C = ErrStat ! We will send back the same error status that is used in OpenFAST + ErrMsg_C = TRANSFER( trim(ErrMsg)//C_NULL_CHAR, ErrMsg_C ) + end subroutine SetErr +END SUBROUTINE HydroDyn_UpdateStates_C + +!=============================================================================================================== +!--------------------------------------------------- HydroDyn End----------------------------------------------- +!=============================================================================================================== + +SUBROUTINE HydroDyn_End_C(ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_End_C') + + INTEGER(C_INT) , INTENT( OUT) :: ErrStat_C + CHARACTER(KIND=C_CHAR) , INTENT( OUT) :: ErrMsg_C + + ! Local variables + INTEGER :: ErrStat + CHARACTER(ErrMsgLen) :: ErrMsg + + ! Call the main subroutine HydroDyn_End + CALL HydroDyn_End( InputData, p, ContStates, DiscStates, ConstrStates, OtherStates, y, m, ErrStat, ErrMsg ) + + call SetErr() + +CONTAINS + subroutine SetErr() + ErrStat_C = ErrStat ! We will send back the same error status that is used in OpenFAST + ErrMsg_C = TRANSFER( trim(ErrMsg)//C_NULL_CHAR, ErrMsg_C ) + end subroutine SetErr +END SUBROUTINE HydroDyn_End_C + +END MODULE HydroDynAPI diff --git a/reg_tests/r-test b/reg_tests/r-test index 72264a44ee..7224fbaafd 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 72264a44ee1da3afe55a98b7cc23e9bf1208cdfd +Subproject commit 7224fbaafd060df59ba3cc37e3599b35573e77b9 From 005f9095bd3963a4f8bb768681ef5549a6a42d6c Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Tue, 18 May 2021 14:12:15 -0600 Subject: [PATCH 02/22] NWTC lib: add InitFileInfo_FromNullCString subroutine This routine allows for passing an input file as one long string with C_NULL_CHAR deliations between lines of the input file Set interface for InitFileInfo. --- modules/nwtc-library/src/NWTC_Base.f90 | 1 + modules/nwtc-library/src/NWTC_IO.f90 | 96 ++++++++++++++++++- .../nwtc-library/src/NWTC_Library_Types.f90 | 4 +- .../src/Registry_NWTC_Library.txt | 4 +- 4 files changed, 96 insertions(+), 9 deletions(-) diff --git a/modules/nwtc-library/src/NWTC_Base.f90 b/modules/nwtc-library/src/NWTC_Base.f90 index dddf9b63fe..6690f0842c 100644 --- a/modules/nwtc-library/src/NWTC_Base.f90 +++ b/modules/nwtc-library/src/NWTC_Base.f90 @@ -38,6 +38,7 @@ MODULE NWTC_Base INTEGER(IntKi), PARAMETER :: ChanLen = 20 !< The maximum allowable length of channel names (i.e., width of output columns) in the FAST framework INTEGER(IntKi), PARAMETER :: MinChanLen = 10 !< The min allowable length of channel names (i.e., width of output columns), used because some modules (like Bladed DLL outputs) have excessively long names INTEGER(IntKi), PARAMETER :: LinChanLen = 200 !< The allowable length of row/column names in linearization files + INTEGER(IntKi), PARAMETER :: MaxFileInfoLineLen = 1024 !< The allowable length of an input line stored in FileInfoType%Lines INTEGER(IntKi), PARAMETER :: NWTC_Verbose = 10 !< The maximum level of verbosity INTEGER(IntKi), PARAMETER :: NWTC_VerboseLevel = 5 !< a number in [0, NWTC_Verbose]: 0 = no output; NWTC_Verbose=verbose; diff --git a/modules/nwtc-library/src/NWTC_IO.f90 b/modules/nwtc-library/src/NWTC_IO.f90 index 6a9e77019b..af6d3ccce3 100644 --- a/modules/nwtc-library/src/NWTC_IO.f90 +++ b/modules/nwtc-library/src/NWTC_IO.f90 @@ -83,6 +83,11 @@ MODULE NWTC_IO !======================================================================= + INTERFACE InitFileInfo + MODULE PROCEDURE InitFileInfo_FromNullCString + MODULE PROCEDURE InitFileInfo_FromStringArray + END INTERFACE + !> \copydoc nwtc_io::allcary1 INTERFACE AllocAry MODULE PROCEDURE AllCAry1 @@ -3356,6 +3361,10 @@ subroutine Print_FileInfo_Struct( U, FileInfo ) character(45) :: TmpStr45 write(U,*) '-------------- Print_FileInfo_Struct ------------' write(U,*) ' info stored in the FileInfo data type' + if (.not. allocated(FileInfo%FileLine) .and. .not. allocated(FileInfo%FileIndx) .and. .not. allocated(FileInfo%Lines)) then + write(U,*) ' Data not allocated' + return + endif write(U,*) ' %NumLines (integer): ',FileInfo%NumLines write(U,*) ' %NumFiles (integer): ',FileInfo%NumFiles write(U,*) ' %FileList (array of strings): ',size(FileInfo%FileList) @@ -4683,7 +4692,84 @@ SUBROUTINE PremEOF ( Fil , Variable, TrapErrors, ErrMsg ) RETURN END SUBROUTINE PremEOF !======================================================================= - SUBROUTINE InitFileInfo( StringArray, FileInfo, ErrStat, ErrMsg ) +!> The following takes an input file as a C_Char string with C_NULL_CHAR deliniating line endings + subroutine InitFileInfo_FromNullCString(FileString, FileInfo, ErrStat, ErrMsg) + CHARACTER(kind=C_char,len=*), intent(in ) :: FileString !< input file as single C string with C_NULL_CHAR separated lines + TYPE(FileInfoType), intent( out) :: FileInfo + INTEGER(IntKi), intent( out) :: ErrStat + CHARACTER(*), intent( out) :: ErrMsg + + integer :: ErrStat2 !< temporary error status from a call + character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call + character(MaxFileInfoLineLen), allocatable :: FileStringArray(:) + character(*), parameter :: RoutineName = 'InitFileInfo_FromNullCString' + integer :: idx, NumLines, MaxLineLen, NullLoc, Line + + ErrStat = ErrID_None + ErrMsg = "" + NumLines = 0 ! Initialize counter for lines + NullLoc = 0 + MaxLineLen = 0 + + ! Find how many non-comment lines we have + do idx=1,len(FileString) + if(FileString(idx:idx) == C_NULL_CHAR) then + MaxLineLen = max(MaxLineLen,idx-NullLoc) + NumLines = NumLines + 1 ! Increment line number + NullLoc = idx + endif + enddo + + if (NumLines == 0) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "Input string contains no C_NULL_CHAR characters. "// & + " Cannot separete passed file info into string array." + if (Failed()) return; + endif + if (MaxLineLen > MaxFileInfoLineLen) then + ErrStat2 = ErrID_Warn + ErrMsg2 = "Input string contains lines longer than "//trim(Num2LStr(MaxFileInfoLineLen))// & + " characters. Check that the flat input file string is properly C_NULL_CHAR delineated" + if (Failed()) return; + endif + + ! Now allocate a string array + call AllocAry( FileStringArray, NumLines, "FileStringArray", ErrStat2, ErrMsg2 ) + if (Failed()) return; + FileStringArray = "" + + ! Now step through the FileString and parse it into the array + idx = 1 ! Index of start + do Line=1,NumLines + ! Index into the next segment + NullLoc = index(FileString(idx:len(FileString)),C_NULL_CHAR) + ! started indexing at idx, so add that back in for location in FileString + NullLoc = NullLoc + idx - 1 + if (NullLoc > 0) then + FileStringArray(Line) = trim(FileString(idx:NullLoc-1)) + else + exit ! exit loop as we didn't find any more + endif + idx = min(NullLoc + 1,len(FileString)) ! Start next segment of file, but overstep end + enddo + + ! Pass through to the FileInfo initialize routine + call InitFileInfo(FileStringArray, FileInfo, ErrStat2, ErrMsg2) + if (Failed()) return; + contains + logical function Failed() + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + if (Failed) then + call Cleanup() + endif + end function Failed + subroutine Cleanup() + if (allocated(FileStringArray)) deallocate(FileStringArray) + end subroutine Cleanup + end subroutine InitFileInfo_FromNullCString +!======================================================================= + SUBROUTINE InitFileInfo_FromStringArray( StringArray, FileInfo, ErrStat, ErrMsg ) CHARACTER(*), DIMENSION(:), INTENT(IN ) :: StringArray TYPE(FileInfoType), INTENT( OUT) :: FileInfo @@ -4694,7 +4780,7 @@ SUBROUTINE InitFileInfo( StringArray, FileInfo, ErrStat, ErrMsg ) character(len=len(StringArray)) :: Line integer :: TmpFileLine(size(StringArray)) - CHARACTER(*), PARAMETER :: RoutineName = 'InitFileInfo' + CHARACTER(*), PARAMETER :: RoutineName = 'InitFileInfo_FromStringArray' INTEGER :: i, NumLines, IC, NumCommChars, LineLen, FirstComm, CommLoc ErrStat = ErrID_None @@ -4735,8 +4821,8 @@ SUBROUTINE InitFileInfo( StringArray, FileInfo, ErrStat, ErrMsg ) ALLOCATE( FileInfo%FileList(FileInfo%NumFiles) ) DO i = 1, FileInfo%NumLines - IF ( LEN(TmpStringArray(i)) > LEN(FileInfo%Lines(i)) ) THEN - CALL SetErrStat( ErrID_Fatal, 'Input string exceeds the bounds of FileInfoType.' , ErrStat, ErrMsg, RoutineName ) + IF ( LEN_TRIM(TmpStringArray(i)) > MaxFileInfoLineLen ) THEN + CALL SetErrStat( ErrID_Fatal, 'Input string '//trim(Num2LStr(i))//' exceeds the bounds of FileInfoType.' , ErrStat, ErrMsg, RoutineName ) RETURN END IF FileInfo%Lines(i) = TmpStringArray(i) @@ -4745,7 +4831,7 @@ SUBROUTINE InitFileInfo( StringArray, FileInfo, ErrStat, ErrMsg ) FileInfo%FileIndx = FileInfo%NumFiles FileInfo%FileList = (/ "passed file info" /) - END SUBROUTINE + END SUBROUTINE InitFileInfo_FromStringArray !======================================================================= !> This routine calls ScanComFile (nwtc_io::scancomfile) and ReadComFile (nwtc_io::readcomfile) !! to move non-comments in a set of nested files starting with TopFile into the FileInfo (nwtc_io::fileinfo) structure. diff --git a/modules/nwtc-library/src/NWTC_Library_Types.f90 b/modules/nwtc-library/src/NWTC_Library_Types.f90 index 61640b5329..47550fd321 100644 --- a/modules/nwtc-library/src/NWTC_Library_Types.f90 +++ b/modules/nwtc-library/src/NWTC_Library_Types.f90 @@ -65,8 +65,8 @@ MODULE NWTC_Library_Types INTEGER(IntKi) :: NumFiles INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: FileLine INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: FileIndx - CHARACTER(1024) , DIMENSION(:), ALLOCATABLE :: FileList - CHARACTER(1024) , DIMENSION(:), ALLOCATABLE :: Lines + CHARACTER(MaxFileInfoLineLen) , DIMENSION(:), ALLOCATABLE :: FileList + CHARACTER(MaxFileInfoLineLen) , DIMENSION(:), ALLOCATABLE :: Lines END TYPE FileInfoType ! ======================= ! ========= Quaternion ======= diff --git a/modules/nwtc-library/src/Registry_NWTC_Library.txt b/modules/nwtc-library/src/Registry_NWTC_Library.txt index a237d6de6a..a77c60c073 100644 --- a/modules/nwtc-library/src/Registry_NWTC_Library.txt +++ b/modules/nwtc-library/src/Registry_NWTC_Library.txt @@ -28,8 +28,8 @@ usefrom NWTC_Library FileInfoType IntKi NumLines usefrom ^ ^ IntKi NumFiles usefrom ^ ^ IntKi FileLine {:} usefrom ^ ^ IntKi FileIndx {:} -usefrom ^ ^ CHARACTER(1024) FileList {:} -usefrom ^ ^ CHARACTER(1024) Lines {:} +usefrom ^ ^ CHARACTER(MaxFileInfoLineLen) FileList {:} +usefrom ^ ^ CHARACTER(MaxFileInfoLineLen) Lines {:} usefrom NWTC_Library Quaternion ReKi q0 usefrom ^ ^ ReKi v {3} From e60f275ffcdf7bae5953bb2ba8741615d61d7860 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Tue, 18 May 2021 16:27:02 -0600 Subject: [PATCH 03/22] NWTC lib: add unit test for InitFileInfo_FromNullCString subroutine and fix broken test Also reorder part of the InitFileInfo_FromStringArray routine --- modules/nwtc-library/src/NWTC_IO.f90 | 10 +++-- .../tests/test_NWTC_IO_FileInfo.F90 | 43 +++++++++++++++++-- 2 files changed, 45 insertions(+), 8 deletions(-) diff --git a/modules/nwtc-library/src/NWTC_IO.f90 b/modules/nwtc-library/src/NWTC_IO.f90 index af6d3ccce3..3ba5602a30 100644 --- a/modules/nwtc-library/src/NWTC_IO.f90 +++ b/modules/nwtc-library/src/NWTC_IO.f90 @@ -4788,6 +4788,7 @@ SUBROUTINE InitFileInfo_FromStringArray( StringArray, FileInfo, ErrStat, ErrMsg NumLines = 0 ! Initialize counter for non-comment populated lines TmpFileLine = 0 ! Line number that was passed in NumCommChars = LEN_TRIM( CommChars ) ! Number of globally specified CommChars + TmpStringArray = "" ! Find how many non-comment lines we have do i=1,size(StringArray) @@ -4820,16 +4821,17 @@ SUBROUTINE InitFileInfo_FromStringArray( StringArray, FileInfo, ErrStat, ErrMsg ALLOCATE( FileInfo%FileIndx(FileInfo%NumLines) ) ALLOCATE( FileInfo%FileList(FileInfo%NumFiles) ) + FileInfo%FileIndx = FileInfo%NumFiles + FileInfo%FileList = (/ "passed file info" /) + FileInfo%Lines = "" ! initialize empty in case of error + FileInfo%FileLine = 0 ! initialize empyt in case of later error DO i = 1, FileInfo%NumLines IF ( LEN_TRIM(TmpStringArray(i)) > MaxFileInfoLineLen ) THEN CALL SetErrStat( ErrID_Fatal, 'Input string '//trim(Num2LStr(i))//' exceeds the bounds of FileInfoType.' , ErrStat, ErrMsg, RoutineName ) - RETURN END IF - FileInfo%Lines(i) = TmpStringArray(i) + FileInfo%Lines(i) = trim(TmpStringArray(i)) FileInfo%FileLine(i) = TmpFileLine(i) END DO - FileInfo%FileIndx = FileInfo%NumFiles - FileInfo%FileList = (/ "passed file info" /) END SUBROUTINE InitFileInfo_FromStringArray !======================================================================= diff --git a/modules/nwtc-library/tests/test_NWTC_IO_FileInfo.F90 b/modules/nwtc-library/tests/test_NWTC_IO_FileInfo.F90 index 3dd9a6ea37..ca205b2909 100644 --- a/modules/nwtc-library/tests/test_NWTC_IO_FileInfo.F90 +++ b/modules/nwtc-library/tests/test_NWTC_IO_FileInfo.F90 @@ -18,7 +18,7 @@ subroutine test_initfileinfo() integer(IntKi) :: error_status character(ErrMsgLen) :: error_message - character(1024) :: input_strings(num_lines) + character(MaxFileInfoLineLen) :: input_strings(num_lines) type(FileInfoType) :: file_info_type integer :: i @@ -62,7 +62,8 @@ subroutine test_initfileinfo2() integer(IntKi) :: error_status character(ErrMsgLen) :: error_message - character(1025) :: input_strings(num_lines) + character(MaxFileInfoLineLen*2) :: input_strings(num_lines) + character(MaxFileInfoLineLen*2) :: tmpstring type(FileInfoType) :: file_info_type integer :: i @@ -73,12 +74,46 @@ subroutine test_initfileinfo2() "line 3", & "line 4" & /) + ! make the last character a + so a trim does not reduce this string length + tmpstring=input_strings(5) + tmpstring(MaxFileInfoLineLen+1:MaxFileInfoLineLen+1) = 'a' + input_strings(5)=tmpstring call InitFileInfo( input_strings, file_info_type, error_status, error_message ) - @assertEqual(num_lines, file_info_type%NumLines) @assertEqual(num_files, file_info_type%NumFiles) @assertEqual( 4, error_status ) end subroutine -end module \ No newline at end of file +@test +subroutine test_initfileinfo3() + USE ISO_C_BINDING, only: C_CHAR, C_NULL_CHAR + + ! This case should result in zero error status. + ! It attempts to initialize FileInfoType with a C_NULL_CHAR delimited string + + integer, parameter :: num_lines = 7 + integer, parameter :: num_files = 1 + + integer(IntKi) :: error_status + character(ErrMsgLen) :: error_message + character(kind=C_CHAR,len=MaxFileInfoLineLen*2) :: input_string + type(FileInfoType) :: file_info_type + integer :: i + + ! Note: the rest of the input string is empty. That should pass ok + input_string = "line 0"//C_NULL_CHAR// & + "line 1"//C_NULL_CHAR// & + "line 2"//C_NULL_CHAR// & + "line 3"//C_NULL_CHAR// & + "line 4"//C_NULL_CHAR// & + "line 5"//C_NULL_CHAR// & + "line 6"//C_NULL_CHAR + call InitFileInfo( input_string, file_info_type, error_status, error_message ) + @assertEqual(num_lines, file_info_type%NumLines) + @assertEqual(num_files, file_info_type%NumFiles) + @assertEqual( 0, error_status ) + +end subroutine + +end module From 0feeb35dbc75dde0a7c1fa6531315e381e2dd898 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Fri, 28 May 2021 14:22:13 -0600 Subject: [PATCH 04/22] hd: update hd_end error handling a bit --- modules/hydrodyn/src/HydroDyn.f90 | 45 +++++++++++------------- modules/hydrodyn/src/HydroDyn_Output.f90 | 4 ++- 2 files changed, 23 insertions(+), 26 deletions(-) diff --git a/modules/hydrodyn/src/HydroDyn.f90 b/modules/hydrodyn/src/HydroDyn.f90 index 3bddcb07db..25bc084f98 100644 --- a/modules/hydrodyn/src/HydroDyn.f90 +++ b/modules/hydrodyn/src/HydroDyn.f90 @@ -1794,52 +1794,47 @@ SUBROUTINE HydroDyn_End( u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg ) CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None + INTEGER(IntKi) :: ErrStat2 ! local error status + CHARACTER(ErrMsgLen) :: ErrMsg2 ! local error message + CHARACTER(*), PARAMETER :: RoutineName = 'HydroDyn_End' ! Initialize ErrStat ErrStat = ErrID_None ErrMsg = "" - - ! Place any last minute operations or calculations here: - - ! Write the HydroDyn-level output file data if the user requested module-level output ! and the current time has advanced since the last stored time step. - IF ( p%OutSwtch == 1 .OR. p%OutSwtch == 3) THEN - CALL HDOut_WriteOutputs( m%LastOutTime, y, p, m%Decimate, ErrStat, ErrMsg ) + CALL HDOut_WriteOutputs( m%LastOutTime, y, p, m%Decimate, ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END IF ! Close files here: - CALL HDOut_CloseOutput( p, ErrStat, ErrMsg ) + CALL HDOut_CloseOutput( p, ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - ! Destroy the input data: - - CALL HydroDyn_DestroyInput( u, ErrStat, ErrMsg ) + ! Destroy the input data: (ignore errors) + CALL HydroDyn_DestroyInput( u, ErrStat2, ErrMsg2 ) - ! Destroy the parameter data: - - CALL HydroDyn_DestroyParam( p, ErrStat, ErrMsg ) + ! Destroy the parameter data: (ignore errors) + CALL HydroDyn_DestroyParam( p, ErrStat2, ErrMsg2 ) - ! Destroy the state data: - - CALL HydroDyn_DestroyContState( x, ErrStat, ErrMsg ) - CALL HydroDyn_DestroyDiscState( xd, ErrStat, ErrMsg ) - CALL HydroDyn_DestroyConstrState( z, ErrStat, ErrMsg ) - CALL HydroDyn_DestroyOtherState( OtherState, ErrStat, ErrMsg ) + ! Destroy the state data: (ignore errors) + CALL HydroDyn_DestroyContState( x, ErrStat2, ErrMsg2 ) + CALL HydroDyn_DestroyDiscState( xd, ErrStat2, ErrMsg2 ) + CALL HydroDyn_DestroyConstrState( z, ErrStat2, ErrMsg2 ) + CALL HydroDyn_DestroyOtherState( OtherState, ErrStat2, ErrMsg2 ) - ! Destroy misc variables: - - CALL HydroDyn_DestroyMisc( m, ErrStat, ErrMsg ) + ! Destroy misc variables: (ignore errors) + CALL HydroDyn_DestroyMisc( m, ErrStat2, ErrMsg2 ) - ! Destroy the output data: - - CALL HydroDyn_DestroyOutput( y, ErrStat, ErrMsg ) + ! Destroy the output data: (ignore errors) + CALL HydroDyn_DestroyOutput( y, ErrStat2, ErrMsg2 ) END SUBROUTINE HydroDyn_End diff --git a/modules/hydrodyn/src/HydroDyn_Output.f90 b/modules/hydrodyn/src/HydroDyn_Output.f90 index 590eeb1290..1951fad322 100644 --- a/modules/hydrodyn/src/HydroDyn_Output.f90 +++ b/modules/hydrodyn/src/HydroDyn_Output.f90 @@ -1231,6 +1231,8 @@ SUBROUTINE HDOut_WriteOutputs( Time, y, p, Decimate, ErrStat, ErrMsg ) ! Local variables INTEGER :: I ! Generic loop counter CHARACTER(200) :: Frmt ! a string to hold a format statement + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 IF (p%UnOutFile < 0 ) RETURN @@ -1268,7 +1270,7 @@ SUBROUTINE HDOut_WriteOutputs( Time, y, p, Decimate, ErrStat, ErrMsg ) WRITE(p%UnOutFile,Frmt,ADVANCE='no') ( p%Delim, y%WriteOutput(I) , I=1,p%NumTotalOuts ) END IF - WRITE (p%UnOutFile,'()', IOSTAT=ErrStat) ! write the line return + WRITE (p%UnOutFile,'()', IOSTAT=ErrStat2) ! write the line return -- ignore error status ELSE Decimate = Decimate + 1 From 45c967db2d912574cc4502a4634cc08f94a08a16 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Fri, 28 May 2021 17:29:38 -0600 Subject: [PATCH 05/22] hd_c_lib: some updates for data passing to init - revisions to hydrodyn_library.py - error handling - initialization vars - new routines for dealing with channel names and units - revisions to HydroDyn_C.f90 - notes on stuff to do - revised passing of channel units and names (old method didn't actually work) - initialization vars passing --- modules/hydrodyn/CMakeLists.txt | 3 + .../hydrodyn/python-lib/hydrodyn_library.py | 235 +++++---- modules/hydrodyn/src/HydroDyn_C.f90 | 445 +++++++++++------- reg_tests/r-test | 2 +- 4 files changed, 423 insertions(+), 262 deletions(-) diff --git a/modules/hydrodyn/CMakeLists.txt b/modules/hydrodyn/CMakeLists.txt index 36ba07451f..a368b16404 100644 --- a/modules/hydrodyn/CMakeLists.txt +++ b/modules/hydrodyn/CMakeLists.txt @@ -62,6 +62,9 @@ target_link_libraries(hydrodynlib nwtclibs) # c-bindings interface library add_library(hydrodyn_c_lib SHARED src/HydroDyn_C.f90) target_link_libraries(hydrodyn_c_lib hydrodynlib) +if(APPLE OR UNIX) + target_compile_definitions(hydrodyn_c_lib PUBLIC -DIMPLICIT_DLLEXPORT) +endif() add_executable(hydrodyn_driver src/HydroDyn_DriverCode.f90) target_link_libraries(hydrodyn_driver hydrodynlib nwtclibs versioninfolib) diff --git a/modules/hydrodyn/python-lib/hydrodyn_library.py b/modules/hydrodyn/python-lib/hydrodyn_library.py index 61f52d4eb5..7448aa6404 100644 --- a/modules/hydrodyn/python-lib/hydrodyn_library.py +++ b/modules/hydrodyn/python-lib/hydrodyn_library.py @@ -1,6 +1,6 @@ #********************************************************************************************************************************** # LICENSING -# Copyright (C) 2021 NREL +# Copyright (C) 2021 National Renewable Energy Laboratory # # This file is part of InflowWind. # @@ -18,10 +18,13 @@ # #********************************************************************************************************************************** # -# This is the Python-C interface library for HydroDyn -# Usage: THIS LIBRARY IS NOT TO BE CHANGED OR EDITED BY THE USER +# This is the Python-C interface library for HydroDyn. This may be used +# directly with Python based codes to call and run HydroDyn. An example of +# using this library from Python is given in the accompanying Python driver +# program. Additional notes and information on the interfacing is included +# there. from ctypes import ( - CDLL, + CDLL, POINTER, create_string_buffer, byref, @@ -30,153 +33,199 @@ c_double, c_float, c_char, - c_char_p, - c_wchar, + c_char_p, + c_wchar, c_wchar_p, c_bool ) import numpy as np -class HydroDynLibAPI(CDLL): +class HydroDynLib(CDLL): + # Human readable error levels from IfW. + error_levels = { + 0: "None", + 1: "Info", + 2: "Warning", + 3: "Severe Error", + 4: "Fatal Error" + } + + # NOTE: the error message length in Fortran is controlled by the + # ErrMsgLen variable in the NWTC_Base.f90 file. If that ever + # changes, it may be necessary to update the corresponding size + # here. + error_msg_c_len = 1025 + def __init__(self, library_path): super().__init__(library_path) self.library_path = library_path self._initialize_routines() + self.ended = False # For error handling at end + # Create buffers for class data - self.abort_error_level = c_int(4) - self.error_status = c_int(0) - self.error_message = create_string_buffer(1025) + self.abort_error_level = 4 + self.error_status_c = c_int(0) + self.error_message_c = create_string_buffer(self.error_msg_c_len) + + # This is not sufficient for HD + #FIXME: ChanLen may not always be 20 -- could be as much as 256 + # Number of channel names may exceeed 5000 + self._channel_names_c = create_string_buffer(20 * 4000) + self._channel_units_c = create_string_buffer(20 * 4000) -# This is not sufficient for HD - self._channel_names = create_string_buffer(20 * 4000) - self._channel_units = create_string_buffer(20 * 4000) + # Initial environmental conditions + self.gravity = 9.80665 # Gravity (m/s^2) + self.defWtrDens = 1025.0 # Water density (kg/m^3) + self.defWtrDpth = 120.0 # Water depth (m) + self.defMSL2SWL = 0.0 # Offset between still-water level and mean sea level (m) [positive upward] - self.dt = c_double(0) - self.total_time = c_double(0) - self.numTimeSteps = c_int(0) + + # Initial time related variables + self.dt = 0.1 # typical default for HD + self.tmax = 600.0 # typical default for HD waves FFT + #FIXME: check tmax/total_time and note exactly what is different between them. + self.total_time = 0.0 # may be longer than tmax + self.numTimeSteps = 0 + + self.numChannels = 0 # Number of channels returned # _initialize_routines() ------------------------------------------------------------------------------------------------------------ def _initialize_routines(self): - self.HydroDyn_Init_C.argtypes = [ + self.HydroDyn_Init_c.argtypes = [ POINTER(c_char_p), # input file string POINTER(c_int), # input file string length + POINTER(c_float), # gravity + POINTER(c_float), # defWtrDens + POINTER(c_float), # defWtrDpth + POINTER(c_float), # defMSL2SWL POINTER(c_double), # dt + POINTER(c_double), # tmax POINTER(c_int), # number of channels POINTER(c_char), # output channel names POINTER(c_char), # output channel units POINTER(c_int), # ErrStat_C POINTER(c_char) # ErrMsg_C ] - self.HydroDyn_Init_C.restype = c_int + self.HydroDyn_Init_c.restype = c_int - self.HydroDyn_CalcOutput_C.argtypes = [ + self.HydroDyn_CalcOutput_c.argtypes = [ POINTER(c_double), # Time_C - POINTER(c_float), # Output Channel Values + #POINTER(c_float), # Output Channel Values POINTER(c_int), # ErrStat_C POINTER(c_char) # ErrMsg_C ] - self.HydroDyn_CalcOutput_C.restype = c_int + self.HydroDyn_CalcOutput_c.restype = c_int - self.HydroDyn_End_C.argtypes = [ + self.HydroDyn_End_c.argtypes = [ POINTER(c_int), # ErrStat_C POINTER(c_char) # ErrMsg_C ] - self.HydroDyn_End_C.restype = c_int + self.HydroDyn_End_c.restype = c_int # hydrodyn_init ------------------------------------------------------------------------------------------------------------ - def hydrodyn_init(self, input_strings, input_string_length): - - #print('hydrodyn_library.py: Running HydroDyn_Init_C .....') + def hydrodyn_init(self, input_string_array): - # Set up inputs - input_string_array = (c_char_p * len(input_strings))() - for i, param in enumerate(input_strings): - input_string_array[i] = param.encode('utf-8') + # 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) - self._numChannels = c_int(0) - - # Run HydroDyn_Init_C - self.HydroDyn_Init_C( - input_string_array, # IN: input file string - byref(c_int(input_string_length)), # IN: input file string length - byref(c_double(self.dt)), # IN: time step (dt) - byref(self._numChannels), # OUT: number of channels - self._channel_names, # OUT: output channel names - self._channel_units, # OUT: output channel units - byref(self.error_status), # OUT: ErrStat_C - self.error_message # OUT: ErrMsg_C + self._numChannels_c = c_int(0) + + #FIXME: need to pass initial position and orientation info + # call HydroDyn_Init_c + self.HydroDyn_Init_c( + c_char_p(input_string), # IN: input file string + byref(c_int(input_string_length)), # IN: 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 + byref(c_float(self.defMSL2SWL)), # IN: default offset between still-water level and mean sea level + byref(c_double(self.dt)), # IN: time step (dt) + byref(c_double(self.tmax)), # IN: tmax + byref(self._numChannels_c), # OUT: number of channels + self._channel_names_c, # OUT: output channel names + self._channel_units_c, # OUT: output channel units + byref(self.error_status_c), # OUT: ErrStat_C + self.error_message_c # OUT: ErrMsg_C ) - if self.fatal_error: - print(f"Error {self.error_status.value}: {self.error_message.value}") - return + + self.check_error() + # Initialize output channels - self._channel_output_array = (c_double * self._numChannels.value)(0.0, ) - self._channel_output_values = np.empty( (self.numTimeSteps, self._numChannels.value) ) + self.numChannels = self._numChannels_c.value - #print('hydrodyn_library.py: Completed HydroDyn_Init_C') # hydrodyn_calcOutput ------------------------------------------------------------------------------------------------------------ - def hydrodyn_calcOutput(self, time, positions, velocities, outputChannelValues): + #def hydrodyn_calcOutput(self, time, positions, velocities, outputChannelValues): + def hydrodyn_calcOutput(self, time): - #print('hydrodyn_library.py: Running HydroDyn_CalcOutput_C .....') - - # Set up inputs -# positions_flat = [pp for p in positions for pp in p] # need to flatten to pass through to Fortran (to reshape) -# positions_flat_c = (c_float * (3 * self.numWindPts))(0.0, ) -# for i, p in enumerate(positions_flat): -# positions_flat_c[i] = c_float(p) -# -# velocities_flat_c = (c_float * (3 * self.numWindPts))(0.0, ) - outputChannelValues_c = (c_float * self._numChannels.value)(0.0, ) - - # Run HydroDyn_CalcOutput_C - self.HydroDyn_CalcOutput_C( + # Run HydroDyn_CalcOutput_c + self.HydroDyn_CalcOutput_c( byref(c_double(time)), # IN: time at which to calculate velocities -# positions_flat_c, # IN: positions - specified by user -# velocities_flat_c, # OUT: velocities at desired positions - outputChannelValues_c, # OUT: output channel values as described in input file - byref(self.error_status), # OUT: ErrStat_C - self.error_message # OUT: ErrMsg_C + #positions_flat_c, # IN: positions - specified by user + #velocities_flat_c, # IN: velocities at desired positions + #velocities_flat_c, # IN: accelerations at desired positions + #forces_flat_c, # OUT: resulting forces/moments array + #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 ) - if self.fatal_error: - print(f"Error {self.error_status.value}: {self.error_message.value}") - return + + self.check_error() # Convert output channel values back into python - for k in range(0,self._numChannels.value): + for k in range(0,self.numChannels): outputChannelValues[k] = float(outputChannelValues_c[k]) - - # Reshape velocities into [N,3] - count = 0 - for j in range(0,self.numWindPts): - velocities[j,0] = velocities_flat_c[count] - velocities[j,1] = velocities_flat_c[count+1] - velocities[j,2] = velocities_flat_c[count+2] - count = count + 3 - - #print('hydrodyn_library.py: Completed HydroDyn_CalcOutput_C') + + ## Reshape velocities into [N,3] + #count = 0 + #for j in range(0,self.numWindPts): + # velocities[j,0] = velocities_flat_c[count] + # velocities[j,1] = velocities_flat_c[count+1] + # velocities[j,2] = velocities_flat_c[count+2] + # count = count + 3 # hydrodyn_end ------------------------------------------------------------------------------------------------------------ def hydrodyn_end(self): + if not self.ended: + self.ended = True + # Run HydroDyn_End_c + self.HydroDyn_End_c( + byref(self.error_status_c), + self.error_message_c + ) - #print('hydrodyn_library.py: Running HydroDyn_End_C .....') + self.check_error() - # Run HydroDyn_End_C - self.HydroDyn_End_C( - byref(self.error_status), - self.error_message - ) - if self.fatal_error: - print(f"Error {self.error_status.value}: {self.error_message.value}") + # other functions ---------------------------------------------------------------------------------------------------------- + def check_error(self): + if self.error_status_c.value == 0: return + elif self.error_status_c.value < self.abort_error_level: + print(f"{self.error_levels[self.error_status_c.value]}: {self.error_message_c.value.decode('ascii')}") + else: + print(f"{self.error_levels[self.error_status_c.value]}: {self.error_message_c.value.decode('ascii')}") + self.hydrodyn_end() + raise Exception("\nHydroDyn terminated prematurely.") - #print('hydrodyn_library.py: Completed HydroDyn_End_C') + @property + def output_channel_names(self): + if len(self._channel_names_c.value.split()) == 0: + return [] + output_channel_names = self._channel_names_c.value.split() + output_channel_names = [n.decode('UTF-8') for n in output_channel_names] + return output_channel_names - # other functions ---------------------------------------------------------------------------------------------------------- @property - def fatal_error(self): - return self.error_status.value >= self.abort_error_level.value + def output_channel_units(self): + if len(self._channel_units_c.value.split()) == 0: + return [] + output_channel_units = self._channel_units_c.value.split() + output_channel_units = [n.decode('UTF-8') for n in output_channel_units] + return output_channel_units diff --git a/modules/hydrodyn/src/HydroDyn_C.f90 b/modules/hydrodyn/src/HydroDyn_C.f90 index 6c1b6d8024..a9171fbe72 100644 --- a/modules/hydrodyn/src/HydroDyn_C.f90 +++ b/modules/hydrodyn/src/HydroDyn_C.f90 @@ -1,6 +1,6 @@ !********************************************************************************************************************************** ! LICENSING -! Copyright (C) 2021 Nicole Mendoza +! Copyright (C) 2021 National Renewable Energy Lab ! ! This file is part of HydroDyn. ! @@ -17,142 +17,216 @@ ! limitations under the License. ! !********************************************************************************************************************************** -MODULE HydroDynAPI +MODULE HydroDyn_C USE ISO_C_BINDING USE HydroDyn USE HydroDyn_Types USE NWTC_Library -IMPLICIT NONE - -PUBLIC :: HydroDyn_Init_C -PUBLIC :: HydroDyn_CalcOutput_C -PUBLIC :: HydroDyn_UpdateStates_C -PUBLIC :: HydroDyn_End_C - -! Accessible to all routines inside module -TYPE(HydroDyn_InputType) :: InputGuess !< An initial guess for the input; the input mesh must be defined, returned by Init -TYPE(HydroDyn_InputType) :: InputData !< Created by IFW_CALCOUTPUT_C and used by IFW_END_C -TYPE(HydroDyn_InitInputType) :: InitInp -TYPE(HydroDyn_InitOutputType) :: InitOutData !< Initial output data -- Names, units, and version info. -TYPE(HydroDyn_ParameterType) :: p !< Parameters -TYPE(HydroDyn_ContinuousStateType) :: ContStates !< Initial continuous states -TYPE(HydroDyn_DiscreteStateType) :: DiscStates !< Initial discrete states -TYPE(HydroDyn_ConstraintStateType) :: ConstrStateGuess !< Initial guess of the constraint states -TYPE(HydroDyn_ConstraintStateType) :: ConstrStates !< Constraint states at Time -TYPE(HydroDyn_OtherStateType) :: OtherStates !< 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) - -!FIXME: we cannot have this fixed!!!!! -INTEGER, PARAMETER :: InputStringLength = 647 !< Fixed length for all lines of the string-based input file + IMPLICIT NONE + + PUBLIC :: HydroDyn_Init_C + PUBLIC :: HydroDyn_CalcOutput_C + PUBLIC :: HydroDyn_UpdateStates_C + PUBLIC :: HydroDyn_End_C + + integer(IntKi), parameter :: InterpOrder = 2 !< quadratic interpolation + + ! All HydroDyn data is stored within the following data structures inside this module. No data is + ! stored within HydroDyn itself, but is instead passed in from this module. This data is not available + ! to the calling code unless explicitly passed through the interface (derived types such as these are + ! non-trivial to pass through the c-bindings). + 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 !< continuous states at Time t + type(HydroDyn_DiscreteStateType) :: xd !< discrete states at Time t + type(HydroDyn_ConstraintStateType) :: z !< Constraint states at Time t + type(HydroDyn_OtherStateType) :: OtherStates !< 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) + + ! If we are doing a correction step, the previous state information will be needed so we can compute with the new inputs + type(HydroDyn_ContinuousStateType) :: x_prev !< continuous states at Time t of previous call + type(HydroDyn_DiscreteStateType) :: xd_prev !< discrete states at Time t of previous call + type(HydroDyn_ConstraintStateType) :: z_prev !< Constraint states at Time t of previous call + + ! Mesh mapping + type(MeshType) :: RefPtMesh ! 1-node Point mesh located at (0,0,0) in global system where all PRP-related driver inputs are set + type(MeshMapType) :: HD_Ref_2_WB_P ! Mesh mapping between Reference pt mesh and WAMIT body(ies) mesh + type(MeshMapType) :: HD_Ref_2_M_P ! Mesh mapping between Reference pt mesh and Morison mesh + real(R8Ki) :: theta(3) ! mesh creation helper data + + ! Time tracking for when we are repeating a timestep + integer(IntKi) :: T_Global ! Time of this call + integer(IntKi) :: T_Global_prev ! time of the previous call + integer(IntKi) :: N_T_Global ! count of which timestep we are on -- calculated internally based on the time passed into UpdateStates + integer(IntKi) :: N_T_Global_prev ! timestep of the previous call + logical :: CorrectionStep ! if we are repeating a timestep in UpdateStates, + + ! This must exactly match the value in the python-lib. If ErrMsgLen changes + ! at some point in the nwtc-library, this should be updated, but the logic + ! exists to correctly handle different lengths of the strings + integer(IntKi), parameter :: ErrMsgLen_C=1025 CONTAINS +!> This routine sets the error status in C_CHAR for export to calling code. +!! Make absolutely certain that we do not overrun the end of ErrMsg_C. That is hard coded to 1025, +!! but ErrMsgLen is set in the nwtc_library, and could change without updates here. We don't want an +!! inadvertant buffer overrun -- that can lead to bad things. +subroutine SetErr(ErrStat, ErrMsg, ErrStat_C, ErrMsg_C) + integer, intent(in ) :: ErrStat !< aggregated error message (fortran type) + character(ErrMsgLen), intent(in ) :: ErrMsg !< aggregated error message (fortran type) + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + ErrStat_C = ErrStat ! We will send back the same error status that is used in OpenFAST + if (ErrMsgLen > ErrMsgLen_C-1) then ! If ErrMsgLen is > the space in ErrMsg_C, do not copy everything over + ErrMsg_C = TRANSFER( trim(ErrMsg(1:ErrMsgLen_C-1))//C_NULL_CHAR, ErrMsg_C ) + else + ErrMsg_C = TRANSFER( trim(ErrMsg)//C_NULL_CHAR, ErrMsg_C ) + endif +end subroutine SetErr + + !=============================================================================================================== !--------------------------------------------- HydroDyn Init---------------------------------------------------- !=============================================================================================================== -SUBROUTINE HydroDyn_Init_C(InputFileStrings_C, InputFileStringLength_C, DT_C, NumChannels_C, OutputChannelNames_C, OutputChannelUnits_C, ErrStat_C, ErrMsg_C) BIND (C, NAME='HydroDyn_Init_C') - - TYPE(C_PTR) , INTENT(IN ) :: InputFileStrings_C - INTEGER(C_INT) , INTENT(IN ) :: InputFileStringLength_C - REAL(C_DOUBLE) , INTENT(IN ) :: DT_C - INTEGER(C_INT) , INTENT( OUT) :: NumChannels_C - TYPE(C_PTR) , INTENT( OUT) :: OutputChannelNames_C - TYPE(C_PTR) , INTENT( OUT) :: OutputChannelUnits_C - INTEGER(C_INT) , INTENT( OUT) :: ErrStat_C - CHARACTER(KIND=C_CHAR) , INTENT( OUT) :: ErrMsg_C(1025) - - ! Local Variables - CHARACTER(InputStringLength), DIMENSION(InputFileStringLength_C) :: InputFileStrings - CHARACTER(kind=C_char, len=1), DIMENSION(:), POINTER :: character_pointer - CHARACTER(kind=C_char, len=1), DIMENSION(:), POINTER :: character_pointer2 - CHARACTER, DIMENSION(InputStringLength) :: single_line_character_array - CHARACTER, DIMENSION(InputStringLength) :: single_line_character_array2 - CHARACTER(InputStringLength) :: single_line_chars - CHARACTER(InputStringLength) :: single_line_chars2 - - CHARACTER(CHANLEN+1), ALLOCATABLE, TARGET :: tmp_OutputChannelNames_C(:) - CHARACTER(CHANLEN+1), ALLOCATABLE, TARGET :: tmp_OutputChannelUnits_C(:) - REAL(DbKi) :: TimeInterval - INTEGER :: ErrStat !< aggregated error message - CHARACTER(ErrMsgLen) :: ErrMsg !< aggregated error message - INTEGER :: ErrStat2 !< temporary error status from a call - CHARACTER(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call - INTEGER :: I, J, K - character(*), parameter :: RoutineName = 'HydroDyn_Init_C' !< for error handling +!FIXME: pass in interporder +SUBROUTINE HydroDyn_Init_c( InputFileString_C, InputFileStringLength_C, & + Gravity_C, defWtrDens_C, defWtrDpth_C, defMSL2SWL_C, & + DT_C, TMax_C, & + NumChannels_C, OutputChannelNames_C, OutputChannelUnits_C, & + ErrStat_C, ErrMsg_C) BIND (C, NAME='HydroDyn_Init_c') + implicit none +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: HydroDyn_Init_c +!GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_Init_c +#endif + + 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 + 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) + real(c_float), intent(in ) :: defMSL2SWL_C !< Default Offset between still-water level and mean sea level (m) [positive upward] (may be overridden by input file) + real(c_double), intent(in ) :: DT_C !< Timestep used with HD for stepping forward from t to t+dt. Must be constant. + real(c_double), intent(in ) :: TMax_C !< Maximum time for simulation (used to set arrays for wave kinematics) + 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*MaxHDOutputs+1) + character(kind=c_char), intent( out) :: OutputChannelUnits_C(ChanLen*MaxHDOutputs+1) + 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 Variables + 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 + character(*), parameter :: RoutineName = 'HydroDyn_Init_C' !< for error handling ! Initialize error handling ErrStat = ErrID_None ErrMsg = "" - ! Convert the string-input from C-style character arrays (char**) to a Fortran-style array of characters. - ! TODO: Add error checking - CALL C_F_pointer(InputFileStrings_C, character_pointer, [InputFileStringLength_C * InputStringLength]) - DO i = 0, InputFileStringLength_C - 1 - single_line_character_array = character_pointer(i * InputStringLength + 1 : i * InputStringLength + InputStringLength) - DO j = 1, InputStringLength - single_line_chars(j:j) = single_line_character_array(j) - END DO - InputFileStrings(i + 1) = single_line_chars - END DO - - ! Store string-inputs as type FileInfoType within HydroDyn_InitInputType - CALL InitFileInfo(InputFileStrings, InitInp%PassedFileData, ErrStat2, ErrMsg2) - if (Failed()) return + ! 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. +! call Print_FileInfo_Struct( CU, InitInp%PassedFileData ) ! CU is the screen -- different number on different systems. + ! Set other inputs for calling HydroDyn_Init InitInp%InputFile = "passed_hd_file" ! dummy InitInp%OutRootName = "HDpassed" ! used for making echo files InitInp%UseInputFile = .FALSE. - TimeInterval = REAL(DT_C, DbKi) InitInp%Linearize = .FALSE. - InitInp%Gravity = 9.80665 ! Gravity (m/s^2) - InitInp%defWtrDens = 1025 ! Water density (kg/m^3) - InitInp%defWtrDpth = 200 ! Water depth (m) - InitInp%defMSL2SWL = 0 ! Offset between still-water level and mean sea level (m) [positive upward] - InitInp%TMax = 600 InitInp%HasIce = .FALSE. -! InitInp%WaveElevXY ! Don't allocate this + InitInp%Linearize = .FALSE. ! we may want to pass this flag to linearize at T=0 for added mass effects + +!FIXME: allocate this based on driver input file +! InitInp%WaveElevXY InitInp%PtfmLocationX = 0.0_ReKi InitInp%PtfmLocationY = 0.0_ReKi - - ! Call the main subroutine HydroDyn_Init - only need InitInp and TimeInterval as inputs, the rest are set by HydroDyn_Init - CALL HydroDyn_Init( InitInp, InputGuess, p, ContStates, DiscStates, ConstrStateGuess, OtherStates, y, m, TimeInterval, InitOutData, ErrStat2, ErrMsg2 ) - if (Failed()) return - - ! Convert the outputs of HydroDyn_Init from Fortran to C - ALLOCATE(tmp_OutputChannelNames_C(size(InitOutData%WriteOutputHdr)),STAT=ErrStat2) + ! Values passed in + InitInp%Gravity = REAL(Gravity_C, ReKi) + InitInp%defWtrDens = REAL(defWtrDens_C, ReKi) + InitInp%defWtrDpth = REAL(defWtrDpth_C, ReKi) + InitInp%defMSL2SWL = REAL(defMSL2SWL_C, ReKi) + TimeInterval = REAL(DT_C, DbKi) + InitInp%TMax = REAL(TMax_C, DbKi) + + + !----------------------- + ! Allocate input array u + !----------------------- + ! These inputs are used in the time stepping algorithm within HD_UpdateStates + ! For quadratic interpolation, 3 timesteps are used. For linear, 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 + allocate(u(3), STAT=ErrStat2) if (ErrStat2 /= 0) then ErrStat2 = ErrID_Fatal ErrMsg2 = "Could not allocate WriteOutputHdr array" + if (Failed()) return endif + +!FIXME: add handling of input mesh info +print*,'FIXME: add input platform location information' + + ! Call the main subroutine HydroDyn_Init + ! TimeInterval and 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, xd, z, OtherStates, y, m, TimeInterval, InitOutData, ErrStat2, ErrMsg2 ) if (Failed()) return - ALLOCATE(tmp_OutputChannelUnits_C(size(InitOutData%WriteOutputUnt)),STAT=ErrStat2) - if (ErrStat2 /= 0) then - ErrStat2 = ErrID_Fatal - ErrMsg2 = "Could not allocate WriteOutputUnt array" - endif - if (Failed()) return + +print*,'FIXME: add mesh handling and mapping' + +!FIXME: copy u(1) to others +!FIXME: what from InitOutData should be returned +!FIXME: what error handling is necessary + + + + !-------------------------------------------------- + ! Set output channel informatoion for driver code. + !-------------------------------------------------- + + ! Number of channels NumChannels_C = size(InitOutData%WriteOutputHdr) - DO I = 1,NumChannels_C - tmp_OutputChannelNames_C(I) = TRANSFER(InitOutData%WriteOutputHdr(I)//C_NULL_CHAR, tmp_OutputChannelNames_C(I)) - tmp_OutputChannelUnits_C(I) = TRANSFER(InitOutData%WriteOutputUnt(I)//C_NULL_CHAR, tmp_OutputChannelUnits_C(I)) - END DO - OutputChannelNames_C = C_LOC(tmp_OutputChannelNames_C) - OutputChannelUnits_C = C_LOC(tmp_OutputChannelUnits_C) + ! transfer the output channel names and units to c_char arrays for returning + k=1 + do i=1,NumChannels_C + 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) + k=k+1 + enddo + enddo + + ! null terminate the string + OutputChannelNames_C(k) = C_NULL_CHAR + OutputChannelUnits_C(k) = C_NULL_CHAR - ! Clean up variables and set up for HydroDyn_CalcOutput_C - CALL HydroDyn_CopyInput(InputGuess, InputData, MESH_NEWCOPY, ErrStat2, ErrMsg2 ) - if (Failed()) return - CALL HydroDyn_CopyConstrState(ConstrStateGuess, ConstrStates, MESH_NEWCOPY, ErrStat2, ErrMsg2 ) - if (Failed()) return call Cleanup() - call SetErr() + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) CONTAINS logical function Failed() @@ -160,37 +234,36 @@ logical function Failed() Failed = ErrStat >= AbortErrLev if (Failed) then call Cleanup() - call SetErr() + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) endif end function Failed subroutine Cleanup() ! NOTE: we are ignoring any error reporting from here - CALL HydroDyn_DestroyInput(InputGuess, ErrStat2, ErrMsg2 ) - CALL HydroDyn_DestroyConstrState(ConstrStateGuess, ErrStat2, ErrMsg2 ) end subroutine Cleanup - subroutine SetErr() - ErrStat_C = ErrStat ! We will send back the same error status that is used in OpenFAST - ErrMsg_C = TRANSFER( trim(ErrMsg)//C_NULL_CHAR, ErrMsg_C ) - end subroutine SetErr +END SUBROUTINE HydroDyn_Init_c -END SUBROUTINE HydroDyn_Init_C !=============================================================================================================== !--------------------------------------------- HydroDyn CalcOutput --------------------------------------------- !=============================================================================================================== -SUBROUTINE HydroDyn_CalcOutput_C(Time_C,OutputChannelValues_C,ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_CalcOutput_C') - REAL(C_DOUBLE) , INTENT(IN ) :: Time_C - REAL(C_FLOAT) , INTENT( OUT) :: OutputChannelValues_C(p%NumOuts) - INTEGER(C_INT) , INTENT( OUT) :: ErrStat_C - CHARACTER(KIND=C_CHAR) , INTENT( OUT) :: ErrMsg_C +SUBROUTINE HydroDyn_CalcOutput_c(Time_C,OutputChannelValues_C,ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_CalcOutput_c') + implicit none +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: HydroDyn_CalcOutput_c +!GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_CalcOutput_c +#endif + real(c_double), intent(in ) :: Time_C + real(c_float), intent( out) :: OutputChannelValues_C(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 :: ErrStat !< aggregated error message - 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 = 'HydroDyn_CalcOutput_C' !< for error handling + real(DbKi) :: Time + integer(IntKi) :: ErrStat !< aggregated error status + character(ErrMsgLen) :: ErrMsg !< aggregated error message + integer(IntKi) :: ErrStat2 !< temporary error status from a call + character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call + character(*), parameter :: RoutineName = 'HydroDyn_CalcOutput_c' !< for error handling ! Initialize error handling ErrStat = ErrID_None @@ -198,96 +271,132 @@ SUBROUTINE HydroDyn_CalcOutput_C(Time_C,OutputChannelValues_C,ErrStat_C,ErrMsg_C ! Convert the inputs from C to Fortrn Time = REAL(Time_C,DbKi) -! InputData%PositionXYZ = reshape( real(Positions_C,ReKi), (/3, InitInp%NumWindPoints/) ) -! ! Call the main subroutine HydroDyn_CalcOutput to get the velocities -! CALL HydroDyn_CalcOutput( Time, InputData, p, ContStates, DiscStates, ConstrStates, OtherStates, y, m, ErrStat2, ErrMsg2 ) -! if (Failed()) return +!FIXME: Reshape position, velocity, acceleration and set mesh accordingly. -! ! Get velocities out of y and flatten them (still in same spot in memory) -! Velocities_C = reshape( REAL(y%VelocityUVW, C_FLOAT), (/3*InitInp%NumWindPoints/) ) ! VelocityUVW is 2D array of ReKi (might need reshape or make into pointer); size [3,N] +!FIXME:Set inputs array and extrap/interp necessary -! ! Get the output channel info out of y -! OutputChannelValues_C = REAL(y%WriteOutput, C_FLOAT) + ! Call the main subroutine HydroDyn_CalcOutput to get the resulting forces and moments at time T + CALL HydroDyn_CalcOutput( Time, u, p, x, xd, z, OtherStates, y, m, ErrStat2, ErrMsg2 ) + if (Failed()) return + + ! Get the output channel info out of y + OutputChannelValues_C = REAL(y%WriteOutput, C_FLOAT) - call SetErr() + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) CONTAINS logical function Failed() CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) Failed = ErrStat >= AbortErrLev - if (Failed) call SetErr() + if (Failed) call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) end function Failed - subroutine SetErr() - ErrStat_C = ErrStat ! We will send back the same error status that is used in OpenFAST - ErrMsg_C = TRANSFER( trim(ErrMsg)//C_NULL_CHAR, ErrMsg_C ) - end subroutine SetErr -END SUBROUTINE HydroDyn_CalcOutput_C +END SUBROUTINE HydroDyn_CalcOutput_c !=============================================================================================================== !--------------------------------------------- HydroDyn UpdateStates --------------------------------------------- !=============================================================================================================== +SUBROUTINE HydroDyn_UpdateStates_c(Time_C,ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_UpdateStates_c') + implicit none +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: HydroDyn_UpdateStates_c +!GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_UpdateStates_c +#endif + real(c_double), intent(in ) :: Time_C + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) -SUBROUTINE HydroDyn_UpdateStates_C(Time_C,ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_UpdateStates_C') - REAL(C_DOUBLE) , INTENT(IN ) :: Time_C - INTEGER(C_INT) , INTENT( OUT) :: ErrStat_C - CHARACTER(KIND=C_CHAR) , INTENT( OUT) :: ErrMsg_C ! Local variables - REAL(DbKi) :: Time - INTEGER :: ErrStat !< aggregated error message - 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 = 'HydroDyn_UpdateStates_C' !< for error handling + real(DbKi) :: Time + integer :: ErrStat !< aggregated error status + 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 = 'HydroDyn_UpdateStates_c' !< for error handling ! Initialize error handling ErrStat = ErrID_None ErrMsg = "" +!FIXME: set time as array ! Convert the inputs from C to Fortran Time = REAL(Time_C,DbKi) +!FIXME: do I need to have a flag for stepping forward in time? What if we are doing a correction step???? +!Rotate values for inputs (do extrap interp if we are not passed new inputs here) + +!FIXME: Reshape position, velocity, acceleration and set mesh inputs + +!FIXME:Set inputs array and extrap/interp necessary + ! ! Call the main subroutine HydroDyn_UpdateStates to get the velocities -! CALL HydroDyn_UpdateStates( Time, InputData, p, ContStates, DiscStates, ConstrStates, OtherStates, y, m, ErrStat2, ErrMsg2 ) +! CALL HydroDyn_UpdateStates( Time, u, p, x, xd, z, OtherStates, y, m, ErrStat2, ErrMsg2 ) ! if (Failed()) return - call SetErr() +!FIXME: what do we need to update after the call? +! states are handled internally. If we are doing correction steps, we may need the old states -CONTAINS + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + +contains logical function Failed() CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) Failed = ErrStat >= AbortErrLev - if (Failed) call SetErr() + if (Failed) call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) end function Failed - subroutine SetErr() - ErrStat_C = ErrStat ! We will send back the same error status that is used in OpenFAST - ErrMsg_C = TRANSFER( trim(ErrMsg)//C_NULL_CHAR, ErrMsg_C ) - end subroutine SetErr -END SUBROUTINE HydroDyn_UpdateStates_C +END SUBROUTINE HydroDyn_UpdateStates_c !=============================================================================================================== !--------------------------------------------------- HydroDyn End----------------------------------------------- !=============================================================================================================== +! NOTE: the error handling in this routine is slightly different than the other routines -SUBROUTINE HydroDyn_End_C(ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_End_C') - - INTEGER(C_INT) , INTENT( OUT) :: ErrStat_C - CHARACTER(KIND=C_CHAR) , INTENT( OUT) :: ErrMsg_C +SUBROUTINE HydroDyn_End_C(ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_End_c') + implicit none +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: HydroDyn_End_c +!GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_End_c +#endif + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) ! Local variables - INTEGER :: ErrStat - CHARACTER(ErrMsgLen) :: ErrMsg + integer(IntKi) :: i !< generic loop counter + integer :: ErrStat !< aggregated error status + 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 = 'HydroDyn_End_c' !< for error handling + + ! Initialize error handling + ErrStat = ErrID_None + ErrMsg = "" + + ! NOTE: HydroDyn_End only takes 1 instance of u, not the array. So extra + ! 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) ! leave first one for passing to HD_End for destruction + call HydroDyn_DestroyInput( u, ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + enddo + endif ! Call the main subroutine HydroDyn_End - CALL HydroDyn_End( InputData, p, ContStates, DiscStates, ConstrStates, OtherStates, y, m, ErrStat, ErrMsg ) + ! If u is not allocated, then we didn't get far at all in initialization, + ! or HD_End_C 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, xd, z, OtherStates, y, m, ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + endif - call SetErr() + ! if u is still allocated, deallocate it now + if (allocated(u)) deallocate(u) -CONTAINS - subroutine SetErr() - ErrStat_C = ErrStat ! We will send back the same error status that is used in OpenFAST - ErrMsg_C = TRANSFER( trim(ErrMsg)//C_NULL_CHAR, ErrMsg_C ) - end subroutine SetErr -END SUBROUTINE HydroDyn_End_C + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + +END SUBROUTINE HydroDyn_End_c -END MODULE HydroDynAPI +END MODULE HydroDyn_C diff --git a/reg_tests/r-test b/reg_tests/r-test index 7224fbaafd..3dcf5751d0 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 7224fbaafd060df59ba3cc37e3599b35573e77b9 +Subproject commit 3dcf5751d03c4a049a2d0f58063ff88ebe239a14 From 55f00209fa684eb9d44d932acf7bdf38a06b39cf Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Fri, 11 Jun 2021 18:22:36 -0600 Subject: [PATCH 06/22] ifw_c_lib: update init inputs, mesh+mapping, move stuff from driver to lib - passing initial node location and orientation info from python to fortran - setup intermediate meshes - setup mesh mapping --- .../hydrodyn/python-lib/hydrodyn_library.py | 225 +++++++++- modules/hydrodyn/src/HydroDyn_C.f90 | 414 +++++++++++++++--- reg_tests/r-test | 2 +- 3 files changed, 548 insertions(+), 93 deletions(-) diff --git a/modules/hydrodyn/python-lib/hydrodyn_library.py b/modules/hydrodyn/python-lib/hydrodyn_library.py index 7448aa6404..5e0fc436bf 100644 --- a/modules/hydrodyn/python-lib/hydrodyn_library.py +++ b/modules/hydrodyn/python-lib/hydrodyn_library.py @@ -23,6 +23,45 @@ # using this library from Python is given in the accompanying Python driver # program. Additional notes and information on the interfacing is included # there. +# +# Note on angles: +# All angles passed in are assumed to be given in radians as an Euler +# angle sequence R(z)*R(y)*R(x) (notice the order as this important when +# passing angles in). Written in matrix form as (doxygen formatted): +# +# \f{eqnarray*}{ +# M & = & R(\theta_z) R(\theta_y) R(\theta_x) \\ +# & = & \begin{bmatrix} \cos(\theta_z) & \sin(\theta_z) & 0 \\ +# -\sin(\theta_z) & \cos(\theta_z) & 0 \\ +# 0 & 0 & 1 \end{bmatrix} +# \begin{bmatrix} \cos(\theta_y) & 0 & -\sin(\theta_y) \\ +# 0 & 1 & 0 \\ +# \sin(\theta_y) & 0 & \cos(\theta_y) \end{bmatrix} +# \begin{bmatrix} 1 & 0 & 0 \\ +# 0 & \cos(\theta_x) & \sin(\theta_x) \\ +# 0 & -\sin(\theta_x) & \cos(\theta_x) \end{bmatrix} \\ +# & = & \begin{bmatrix} +# \cos(\theta_y)\cos(\theta_z) & \cos(\theta_x)\sin(\theta_z)+\sin(\theta_x)\sin(\theta_y)\cos(\theta_z) & +# \sin(\theta_x)\sin(\theta_z)-\cos(\theta_x)\sin(\theta_y)\cos(\theta_z) \\ +# -\cos(\theta_y)\sin(\theta_z) & \cos(\theta_x)\cos(\theta_z)-\sin(\theta_x)\sin(\theta_y)\sin(\theta_z) & +# \sin(\theta_x)\cos(\theta_z)+\cos(\theta_x)\sin(\theta_y)\sin(\theta_z) \\ +# \sin(\theta_y) & -\sin(\theta_x)\cos(\theta_y) & \cos(\theta_x)\cos(\theta_y) \\ +# \end{bmatrix} +# \f} +# +# When passed into the Fortran library, this Euler angle set is converted +# into a DCM (direction cosine matrix) and stored on the input mesh. All +# calculations internally in HD are then performed using the DCM form of +# the input angles. +# +# It should be noted that a small angle assumption when returning the +# outputs for the platform roll, pitch, and yaw. These angles are +# assumed to be small enough that treating them as independent angles +# does not introduce significant error in the output channels. This may +# yield a small discrepency between the values passed in and the returned +# output channel values. These output channels should not be directly +# used -- they are only for reporting purposes. +# from ctypes import ( CDLL, POINTER, @@ -39,6 +78,7 @@ c_bool ) import numpy as np +import datetime class HydroDynLib(CDLL): # Human readable error levels from IfW. @@ -81,6 +121,8 @@ def __init__(self, library_path): self.defWtrDpth = 120.0 # Water depth (m) self.defMSL2SWL = 0.0 # Offset between still-water level and mean sea level (m) [positive upward] + # Interpolation order (must be 1: linear, or 2: quadratic) + self.InterpOrder = 1 # default of linear interpolation # Initial time related variables self.dt = 0.1 # typical default for HD @@ -91,41 +133,63 @@ def __init__(self, library_path): self.numChannels = 0 # Number of channels returned + # Number of bodies and initial reference point + # The initial position is only set as (X,Y). The Z value and + # orientation is set by HD and will be returned along with the full + # set of numBodies where it is expecting loads inputs. + self.ptfmRefPt_x = 0.0 + self.ptfmRefPt_y = 0.0 + + # Nodes + # The number of nodes must be constant throughout simulation. The + # initial position is given in the initNodeLocations array (resize as + # needed, should be Nx6). + # Rotations are given in radians assuming small angles. See note at + # top of this file. + self.numNodePts = 1 # Single ptfm attachment point for floating rigid + self.initNodeLocations = np.zeros((self.numNodePts,6)) # N x 6 array [x,y,z,Rx,Ry,Rz] + # _initialize_routines() ------------------------------------------------------------------------------------------------------------ def _initialize_routines(self): self.HydroDyn_Init_c.argtypes = [ - POINTER(c_char_p), # input file string - POINTER(c_int), # input file string length - POINTER(c_float), # gravity - POINTER(c_float), # defWtrDens - POINTER(c_float), # defWtrDpth - POINTER(c_float), # defMSL2SWL - POINTER(c_double), # dt - POINTER(c_double), # tmax - POINTER(c_int), # number of channels - POINTER(c_char), # output channel names - POINTER(c_char), # output channel units - POINTER(c_int), # ErrStat_C - POINTER(c_char) # ErrMsg_C + POINTER(c_char_p), # input file string + POINTER(c_int), # input file string length + POINTER(c_float), # gravity + POINTER(c_float), # defWtrDens + POINTER(c_float), # defWtrDpth + POINTER(c_float), # defMSL2SWL + POINTER(c_float), # PtfmRefPt_X + POINTER(c_float), # PtfmRefPt_Y + POINTER(c_int), # numNodePts -- number of points expecting motions/loads + POINTER(c_float), # initNodeLocations -- initial node positions in flat array of 6*numNodePts + POINTER(c_int), # InterpOrder + POINTER(c_double), # dt + POINTER(c_double), # tmax + POINTER(c_int), # number of channels + POINTER(c_char), # output channel names + POINTER(c_char), # output channel units + POINTER(c_int), # ErrStat_C + POINTER(c_char) # ErrMsg_C ] self.HydroDyn_Init_c.restype = c_int self.HydroDyn_CalcOutput_c.argtypes = [ - POINTER(c_double), # Time_C - #POINTER(c_float), # Output Channel Values - POINTER(c_int), # ErrStat_C - POINTER(c_char) # ErrMsg_C + POINTER(c_double), # Time_C + POINTER(c_float), # Output Channel Values + POINTER(c_int), # ErrStat_C + POINTER(c_char) # ErrMsg_C ] self.HydroDyn_CalcOutput_c.restype = c_int self.HydroDyn_End_c.argtypes = [ - POINTER(c_int), # ErrStat_C - POINTER(c_char) # ErrMsg_C + POINTER(c_int), # ErrStat_C + POINTER(c_char) # ErrMsg_C ] self.HydroDyn_End_c.restype = c_int # hydrodyn_init ------------------------------------------------------------------------------------------------------------ def hydrodyn_init(self, 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. @@ -135,6 +199,31 @@ def hydrodyn_init(self, input_string_array): self._numChannels_c = c_int(0) + #FIXME: for now we only allow for a single rigid motion platform. + if self.initNodeLocations.shape[1] != 6: + print("Expecting a Nx6 array of initial node locations (initNodeLocations) with second index for [x,y,z,Rx,Ry,Rz]") + self.hydrodyn_end() + raise Exception("\nHydroDyn terminated prematurely.") + if self.initNodeLocations.shape[0] != self.numNodePts: + print("Expecting a Nx6 array of initial node locations (initNodeLocations) with first index for number of nodes.") + self.hydrodyn_end() + raise Exception("\nHydroDyn terminated prematurely.") + # Placeholders exist for multiple input points for flexible + # platforms, but are for now not enabled. + if self.numNodePts != 1: + print(f"HydroDyn C interface does not currently support flexible structures ({self.numNodePts} input motion points were requested, but only 1 is currently supported).") + self.hydrodyn_end() + raise Exception("\nHydroDyn terminated prematurely.") + + # initNodeLocations + # Verify that the shape of initNodeLocations is correct + # Make a flat 1D array of position info: + # [x2,y1,z1,Rx1,Ry1,Rz1, x2,y2,z2,Rx2,Ry2,Rz2 ...] + nodeInitLoc_flat = [pp for p in self.initNodeLocations for pp in p] + nodeInitLoc_flat_c = (c_float * (6 * self.numNodePts))(0.0,) + for i, p in enumerate(nodeInitLoc_flat): + nodeInitLoc_flat_c[i] = c_float(p) + #FIXME: need to pass initial position and orientation info # call HydroDyn_Init_c self.HydroDyn_Init_c( @@ -144,6 +233,11 @@ def hydrodyn_init(self, input_string_array): byref(c_float(self.defWtrDens)), # IN: default water density byref(c_float(self.defWtrDpth)), # IN: default water depth byref(c_float(self.defMSL2SWL)), # IN: default offset between still-water level and mean sea level + byref(c_float(self.ptfmRefPt_x)), # IN: Platform initial position (X) + byref(c_float(self.ptfmRefPt_y)), # IN: Platform initial position (Y) + byref(c_int(self.numNodePts)), # IN: number of attachment points expected (where motions are transferred into HD) + nodeInitLoc_flat_c, # IN: initNodeLocations -- initial node positions in flat array of 6*numNodePts + byref(c_int(self.InterpOrder)), # IN: InterpOrder (1: linear, 2: quadratic) byref(c_double(self.dt)), # IN: time step (dt) byref(c_double(self.tmax)), # IN: tmax byref(self._numChannels_c), # OUT: number of channels @@ -154,7 +248,6 @@ def hydrodyn_init(self, input_string_array): ) self.check_error() - # Initialize output channels self.numChannels = self._numChannels_c.value @@ -162,8 +255,12 @@ def hydrodyn_init(self, input_string_array): # hydrodyn_calcOutput ------------------------------------------------------------------------------------------------------------ #def hydrodyn_calcOutput(self, time, positions, velocities, outputChannelValues): - def hydrodyn_calcOutput(self, time): + def hydrodyn_calcOutput(self, time, outputChannelValues): + # Set up inputs + + # Set up output channels + outputChannelValues_c = (c_float * self.numChannels)(0.0,) # Run HydroDyn_CalcOutput_c self.HydroDyn_CalcOutput_c( @@ -172,7 +269,7 @@ def hydrodyn_calcOutput(self, time): #velocities_flat_c, # IN: velocities at desired positions #velocities_flat_c, # IN: accelerations at desired positions #forces_flat_c, # OUT: resulting forces/moments array - #outputChannelValues_c, # OUT: output channel values as described in input file + 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 ) @@ -229,3 +326,87 @@ def output_channel_units(self): output_channel_units = self._channel_units_c.value.split() output_channel_units = [n.decode('UTF-8') for n in output_channel_units] return output_channel_units + + +#=============================================================================== +# Helper class for debugging the interface. This will write out all the +# input position/orientation, velocities, accelerations, and the resulting +# forces and moments at each input node. If all is functioning correctly, +# this will be identical to the corresponding values in the HydroDyn output +# channels. + +class DriverDbg(): + """ + This is only for debugging purposes only. The input motions and resulting + forces can be written to file with this class to verify the data I/O to the + Fortran library. + When coupled to another code, the force/moment array would be passed back + to the calling code for use in the structural solver. + """ + def __init__(self,filename): + self.DbgFile=open(filename,'wt') # open output file and write header info + # write file header + t_string=datetime.datetime.now() + dt_string=datetime.date.today() +# self.DbgFile.write(f"## This file was generated by InflowWind_Driver on {dt_string.strftime('%b-%d-%Y')} at {t_string.strftime('%H:%M:%S')}\n") +# self.DbgFile.write(f"## This file contains the wind velocity at the {numWindPts} points specified in the file ") +# self.DbgFile.write(f"{filename}\n") +# self.DbgFile.write("#\n") +# self.DbgFile.write("# T X Y Z U V W\n") +# self.DbgFile.write("# (s) (m) (m) (m) (m/s) (m/s) (m/s)\n") + self.opened = True + +# def write(self,t,positions,velocities): +# for p, v in zip(positions,velocities): +# self.DbgFile.write(' %14.7f %14.7f %14.7f %14.7f %14.7f %14.7f %14.7f\n' % (t,p[0],p[1],p[2],v[0],v[1],v[2])) + + def end(self): + if self.opened: + self.DbgFile.close() + self.opened = False + + +#=============================================================================== +# Helper class for writing channels to file. +# for the regression testing to mirror the output from the InfowWind Fortran +# driver. This may also have value for debugging the interfacing to IfW. + +class WriteOutChans(): + """ + This is only for testing purposes. Since we are not returning the + output channels to anything, we will write them to file. When coupled to + another code, this data would be passed back for inclusion the any output + file there. + """ + def __init__(self,filename,chan_names,chan_units): + chan_names.insert(0,'Time') # add time index header + chan_units.insert(0,'(s)') # add time index unit + self.OutFile=open(filename,'wt') # open output file and write header info + # write file header + t_string=datetime.datetime.now() + dt_string=datetime.date.today() + self.OutFile.write(f"## This file was generated by InflowWind_Driver on {dt_string.strftime('%b-%d-%Y')} at {t_string.strftime('%H:%M:%S')}\n") + self.OutFile.write(f"## This file contains output channels requested from the OutList section of the input file") + self.OutFile.write(f"{filename}\n") + self.OutFile.write("#\n") + self.OutFile.write("#\n") + self.OutFile.write("#\n") + self.OutFile.write("#\n") + l = len(chan_names) + f_string = "{:^15s}"+" {:^20s} "*(l-1) + self.OutFile.write(f_string.format(*chan_names) + '\n') + self.OutFile.write(f_string.format(*chan_units) + '\n') + self.opened = True + + def write(self,chan_data): + l = chan_data.shape[1] + f_string = "{:10.4f}"+"{:25.7f}"*(l-1) + for i in range(0,chan_data.shape[0]): + self.OutFile.write(f_string.format(*chan_data[i,:]) + '\n') + if i==0: + print(f"{chan_data[i,:]}") + + def end(self): + if self.opened: + self.OutFile.close() + self.opened = False diff --git a/modules/hydrodyn/src/HydroDyn_C.f90 b/modules/hydrodyn/src/HydroDyn_C.f90 index a9171fbe72..be57af1f45 100644 --- a/modules/hydrodyn/src/HydroDyn_C.f90 +++ b/modules/hydrodyn/src/HydroDyn_C.f90 @@ -31,12 +31,32 @@ MODULE HydroDyn_C PUBLIC :: HydroDyn_UpdateStates_C PUBLIC :: HydroDyn_End_C - integer(IntKi), parameter :: InterpOrder = 2 !< quadratic interpolation + !------------------------------------------------------------------------------------ + ! Error handling + ! This must exactly match the value in the python-lib. If ErrMsgLen changes at + ! some point in the nwtc-library, this should be updated, but the logic exists + ! to correctly handle different lengths of the strings + integer(IntKi), parameter :: ErrMsgLen_C=1025 + - ! All HydroDyn data is stored within the following data structures inside this module. No data is - ! stored within HydroDyn itself, but is instead passed in from this module. This data is not available - ! to the calling code unless explicitly passed through the interface (derived types such as these are - ! non-trivial to pass through the c-bindings). + !------------------------------------------------------------------------------------ + ! Data storage + ! All HydroDyn data is stored within the following data structures inside this + ! module. No data is stored within HydroDyn itself, but is instead passed in + ! from this module. This data is not available to the calling code unless + ! explicitly passed through the interface (derived types such as these are + ! non-trivial to pass through the c-bindings). + !------------------------------ + ! Extrapolation and interpolation + ! For the solver in HD, previous timesteps input must be stored for extrapolation + ! to the t+dt timestep. This can be either linear (1) quadratic (2). The + ! InterpOrder variable tracks what this is and sets the size of the inputs `u` + ! passed into HD. Inputs `u` will be sized as follows: + ! linear interp u(2) with inputs at T,T-dt + ! quadratic interp u(3) with inputs at T,T-dt,T-2*dt + 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. @@ -47,29 +67,66 @@ MODULE HydroDyn_C type(HydroDyn_OtherStateType) :: OtherStates !< 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) - - ! If we are doing a correction step, the previous state information will be needed so we can compute with the new inputs + !------------------------------ + ! Correction steps + ! OpenFAST has the ability to perform correction steps. During a correction + ! step, new input values are passed in but the timestep remains the same. + ! When this occurs the new input data at time t is used with the state + ! information from the previous timestep (t) to calculate new state values + ! time t+dt in the UpdateStates routine. In OpenFAST this is all handled by + ! the glue code. However, here we do not pass state information through the + ! interface and therefore must store it here analogously to how it is handled + ! in the OpenFAST glue code. type(HydroDyn_ContinuousStateType) :: x_prev !< continuous states at Time t of previous call type(HydroDyn_DiscreteStateType) :: xd_prev !< discrete states at Time t of previous call type(HydroDyn_ConstraintStateType) :: z_prev !< Constraint states at Time t of previous call - - ! Mesh mapping - type(MeshType) :: RefPtMesh ! 1-node Point mesh located at (0,0,0) in global system where all PRP-related driver inputs are set - type(MeshMapType) :: HD_Ref_2_WB_P ! Mesh mapping between Reference pt mesh and WAMIT body(ies) mesh - type(MeshMapType) :: HD_Ref_2_M_P ! Mesh mapping between Reference pt mesh and Morison mesh - real(R8Ki) :: theta(3) ! mesh creation helper data - - ! Time tracking for when we are repeating a timestep + !------------------------------ + ! Time tracking + ! When we are performing a correction step, time information of previous + ! calls is needed to decide how to apply correction logic or cycle the inputs + ! and resave the previous timestep states. + integer(IntKi) :: dT_Global ! dT of the code calling this module integer(IntKi) :: T_Global ! Time of this call integer(IntKi) :: T_Global_prev ! time of the previous call integer(IntKi) :: N_T_Global ! count of which timestep we are on -- calculated internally based on the time passed into UpdateStates integer(IntKi) :: N_T_Global_prev ! timestep of the previous call logical :: CorrectionStep ! if we are repeating a timestep in UpdateStates, - - ! This must exactly match the value in the python-lib. If ErrMsgLen changes - ! at some point in the nwtc-library, this should be updated, but the logic - ! exists to correctly handle different lengths of the strings - integer(IntKi), parameter :: ErrMsgLen_C=1025 + !------------------------------------------------------------------------------------ + + + + !------------------------------------------------------------------------------------ + ! Meshes for motions and loads + ! Meshes are used within HD to handle all motions and loads. Rather than directly + ! map to those nodes, we will create a mapping to go between the array of node + ! positions passed into this module and what is used inside HD. This is done + ! through a pair of meshes for the motion and loads corresponding to the node + ! positions passed in. + !------------------------------ + ! Meshes for external nodes + ! These point meshes are merely used to simplify the mapping of motions/loads + ! to/from HD using the library mesh mapping routines. These meshes may contain + ! one or multiple points. + ! - 1 point -- rigid floating body assumption + ! - N points -- flexible structure (either floating or fixed bottom) + integer(IntKi) :: NumNodePts ! Number of mesh points we are interfacing motions/loads to/from HD + type(MeshType) :: HD_MotionMesh ! mesh for motions of external nodes + type(MeshType) :: HD_LoadMesh ! mesh for loads for external nodes + !------------------------------ + ! Mesh mapping: motions + ! The mapping of motions from the nodes passed in to the corresponding HD meshes + type(MeshMapType) :: Map_Motion_2_HD_PRP_P ! Mesh mapping between input motion mesh and PRP + type(MeshMapType) :: Map_Motion_2_HD_WB_P ! Mesh mapping between input motion mesh and WAMIT body(ies) mesh + type(MeshMapType) :: Map_Motion_2_HD_Mo_P ! Mesh mapping between input motion mesh and Morison mesh + !------------------------------ + ! Mesh mapping: loads + ! The mapping of loads from the HD meshes to the corresponding external nodes + type(MeshMapType) :: Map_HD_WB_P_2_Load ! Mesh mapping between HD output WAMIT body loads mesh and external nodes mesh + type(MeshMapType) :: Map_HD_Mo_P_2_Load ! Mesh mapping between HD output Morison loads mesh and external nodes mesh + ! Meshes -- helper stuff + real(R8Ki) :: theta(3) ! mesh creation helper data + !------------------------------------------------------------------------------------ + CONTAINS @@ -94,10 +151,12 @@ end subroutine SetErr !=============================================================================================================== !--------------------------------------------- HydroDyn Init---------------------------------------------------- !=============================================================================================================== -!FIXME: pass in interporder SUBROUTINE HydroDyn_Init_c( InputFileString_C, InputFileStringLength_C, & Gravity_C, defWtrDens_C, defWtrDpth_C, defMSL2SWL_C, & - DT_C, TMax_C, & + PtfmRefPtPositionX_C, PtfmRefPtPositionY_C, & + NumNodePts_C, InitNodePositions_C, & + !NumWaveElev_C, WaveElevXY_C & !Placeholder for later + InterpOrder_C, DT_C, TMax_C, & NumChannels_C, OutputChannelNames_C, OutputChannelUnits_C, & ErrStat_C, ErrMsg_C) BIND (C, NAME='HydroDyn_Init_c') implicit none @@ -106,35 +165,52 @@ SUBROUTINE HydroDyn_Init_c( InputFileString_C, InputFileStringLength_C, !GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_Init_c #endif - 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 - 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) - real(c_float), intent(in ) :: defMSL2SWL_C !< Default Offset between still-water level and mean sea level (m) [positive upward] (may be overridden by input file) - real(c_double), intent(in ) :: DT_C !< Timestep used with HD for stepping forward from t to t+dt. Must be constant. - real(c_double), intent(in ) :: TMax_C !< Maximum time for simulation (used to set arrays for wave kinematics) - integer(c_int), intent( out) :: NumChannels_C !< Number of output channels requested from the input file + 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 + 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) + real(c_float), intent(in ) :: defMSL2SWL_C !< Default Offset between still-water level and mean sea level (m) [positive upward] (may be overridden by input file) + real(c_float), intent(in ) :: PtfmRefPtPositionX_C !< Initial position in wave field + real(c_float), intent(in ) :: PtfmRefPtPositionY_C !< Initial position in wave field + integer(c_int), intent(in ) :: NumNodePts_C !< Number of mesh points we are transfering motions to and output loads to + real(c_float), intent(in ) :: InitNodePositions_C( 6, NumNodePts_C ) !< A 6xNumNodePts_C array [x,y,z,Rx,Ry,Rz] + !NOTE: not setting up the WaveElev at this point. Leaving placeholder for future + !integer(c_int), intent(in ) :: NumWaveElev_C !< Number of mesh points we are transfering motions to and output loads to + !real(c_float), intent(in ) :: WaveElevXY_C !< A 2xNumWaveElev_C array [x,y] + integer(c_int), intent(in ) :: InterpOrder_C !< Interpolation order to use (must be 1 or 2) + real(c_double), intent(in ) :: DT_C !< Timestep used with HD for stepping forward from t to t+dt. Must be constant. + real(c_double), intent(in ) :: TMax_C !< Maximum time for simulation (used to set arrays for wave kinematics) + 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*MaxHDOutputs+1) character(kind=c_char), intent( out) :: OutputChannelUnits_C(ChanLen*MaxHDOutputs+1) - 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) + 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 Variables - 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 + character(kind=C_char, len=InputFileStringLength_C), pointer :: InputFileString !< Input file as a single string with NULL chracter separating lines + real(ReKi), allocatable :: InitNodePositions(:,:) !< temp array. Probably don't need this, but makes conversion from C clearer. + + 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_Init_C' !< for error handling ! Initialize error handling ErrStat = ErrID_None ErrMsg = "" + ! Sanity checks on values passed + InterpOrder = int(InterpOrder_C, IntKi) + if ( InterpOrder < 1_IntKi .or. InterpOrder > 2_IntKi ) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "InterpOrder passed into HydroDyn_C must be 1 (linear) or 2 (quadratic)" + if (Failed()) return + endif + ! Get fortran pointer to C_NULL_CHAR deliniated input file as a string call C_F_pointer(InputFileString_C, InputFileString) @@ -143,49 +219,76 @@ SUBROUTINE HydroDyn_Init_c( InputFileString_C, InputFileStringLength_C, ! For diagnostic purposes, the following can be used to display the contents ! of the InFileInfo data structure. -! call Print_FileInfo_Struct( CU, InitInp%PassedFileData ) ! CU is the screen -- different number on different systems. + ! 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%OutRootName = "HDpassed" ! used for making echo files - InitInp%UseInputFile = .FALSE. + InitInp%OutRootName = "HDpassed" ! FIXME: pass in a name to use here (For summary files etc) + 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. - InitInp%HasIce = .FALSE. - InitInp%Linearize = .FALSE. ! we may want to pass this flag to linearize at T=0 for added mass effects - -!FIXME: allocate this based on driver input file -! InitInp%WaveElevXY - InitInp%PtfmLocationX = 0.0_ReKi - InitInp%PtfmLocationY = 0.0_ReKi ! Values passed in - InitInp%Gravity = REAL(Gravity_C, ReKi) + InitInp%Gravity = REAL(Gravity_C, ReKi) InitInp%defWtrDens = REAL(defWtrDens_C, ReKi) InitInp%defWtrDpth = REAL(defWtrDpth_C, ReKi) InitInp%defMSL2SWL = REAL(defMSL2SWL_C, ReKi) - TimeInterval = REAL(DT_C, DbKi) - InitInp%TMax = REAL(TMax_C, DbKi) + TimeInterval = REAL(DT_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 + ! side. Will validate this against what HD reads from the initialization info. + NumNodePts = int(NumNodePts_C, IntKi) + if (NumNodePts < 1) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "At least one node point must be specified" + if (Failed()) return + endif + call AllocAry( InitNodePositions, 6, NumNodePts, "InitNodePositions", ErrStat2, ErrMsg2 ) + if (Failed()) return + do i=1,NumNodePts + InitNodePositions(1:6,i) = real(InitNodePositions_C(1:6,i),ReKi) + enddo + + ! 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. For now, skipping this feature. + ! Skipping this for now. Maybe add later. + !InitInp%WaveElevXY + !----------------------- ! Allocate input array u !----------------------- ! These inputs are used in the time stepping algorithm within HD_UpdateStates - ! For quadratic interpolation, 3 timesteps are used. For linear, 2 timesteps - ! (the HD code can handle either). + ! 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 - allocate(u(3), STAT=ErrStat2) + allocate(u(InterpOrder+1), STAT=ErrStat2) if (ErrStat2 /= 0) then ErrStat2 = ErrID_Fatal - ErrMsg2 = "Could not allocate WriteOutputHdr array" + ErrMsg2 = "Could not allocate inuput" if (Failed()) return endif -!FIXME: add handling of input mesh info -print*,'FIXME: add input platform location information' ! Call the main subroutine HydroDyn_Init ! TimeInterval and InitInp are passed into HD_Init, all the rest are set by HD_Init @@ -195,13 +298,27 @@ SUBROUTINE HydroDyn_Init_c( InputFileString_C, InputFileStringLength_C, call HydroDyn_Init( InitInp, u(1), p, x, xd, z, OtherStates, y, m, TimeInterval, InitOutData, ErrStat2, ErrMsg2 ) if (Failed()) return -print*,'FIXME: add mesh handling and mapping' -!FIXME: copy u(1) to others -!FIXME: what from InitOutData should be returned -!FIXME: what error handling is necessary + !------------------------------------------------------------- + ! 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 pieces of u(:) + !------------------------------------------------------------- +! copy u(1) to others +! what from InitOutData should be returned +! what error handling is necessary + !-------------------------------------------------- ! Set output channel informatoion for driver code. @@ -237,8 +354,145 @@ logical function Failed() call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) endif end function Failed - subroutine Cleanup() ! NOTE: we are ignoring any error reporting from here + + subroutine Cleanup() + if (allocated(InitNodePositions)) deallocate(InitNodePositions) end subroutine Cleanup + + !> This subroutine sets the interface meshes to map to the input motions to the HD + !! meshes for WAMIT and Morison. This subroutine also sets the meshes for the + !! output loads. + subroutine SetMotionLoadsInterfaceMeshes(ErrStat3,ErrMsg3) + integer(IntKi), intent( out) :: ErrStat3 !< temporary error status + character(ErrMsgLen), intent( out) :: ErrMsg3 !< temporary error message + integer(IntKi) :: iNode + real(ReKi) :: InitPos(3) + real(R8Ki) :: theta(3) + real(R8Ki) :: Orient(3,3) + !------------------------------------------------------------- + ! Set the interface meshes for motion inputs and loads output + !------------------------------------------------------------- + ! Motion mesh + ! This point mesh may contain more than one point. Mapping will be used to map + ! this to the input meshes for WAMIT and Morison. + call MeshCreate( HD_MotionMesh , & + IOS = COMPONENT_INPUT , & + Nnodes = NumNodePts , & + ErrStat = ErrStat3 , & + ErrMess = ErrMsg3 , & + TranslationDisp = .TRUE., Orientation = .TRUE., & + TranslationVel = .TRUE., RotationVel = .TRUE., & + TranslationAcc = .TRUE., RotationAcc = .TRUE. ) + if (ErrStat3 >= AbortErrLev) return + + do iNode=1,NumNodePts + ! initial position and orientation of node + InitPos = InitNodePositions(1:3,iNode) + theta = real(InitNodePositions(4:6,iNode),DbKi) ! convert ReKi to DbKi to avoid roundoff + Orient = EulerConstruct( theta ) + call MeshPositionNode( HD_MotionMesh , & + iNode , & + InitPos , & ! position + ErrStat3, ErrMsg3 , & + Orient ) ! orientation + if (ErrStat3 >= AbortErrLev) return + + call MeshConstructElement ( HD_MotionMesh, ELEMENT_POINT, ErrStat3, ErrMsg3, iNode ) + if (ErrStat3 >= AbortErrLev) return + enddo + + call MeshCommit ( HD_MotionMesh, ErrStat3, ErrMsg3 ) + if (ErrStat3 >= AbortErrLev) return + + HD_MotionMesh%RemapFlag = .TRUE. + + ! For checking the mesh, uncomment this. + ! note: CU is is output unit (platform dependent). + !call MeshPrintInfo( CU, HD_MotionMesh ) + + + !------------------------------------------------------------- + ! Loads mesh + ! This point mesh may contain more than one point. Mapping will be used to map + ! the loads from output meshes for WAMIT and Morison. + ! Output mesh for loads at each WAMIT body + CALL MeshCopy( SrcMesh = HD_MotionMesh ,& + DestMesh = HD_LoadMesh ,& + CtrlCode = MESH_SIBLING ,& + IOS = COMPONENT_OUTPUT ,& + ErrStat = ErrStat3 ,& + ErrMess = ErrMsg3 ,& + Force = .TRUE. ,& + Moment = .TRUE. ) + if (ErrStat3 >= AbortErrLev) return + + HD_LoadMesh%RemapFlag = .TRUE. + + ! For checking the mesh, uncomment this. + ! note: CU is is output unit (platform dependent). + !call MeshPrintInfo( CU, HD_LoadMesh ) + + !------------------------------------------------------------- + ! Set the mapping meshes + ! PRP - principle reference point + call MeshMapCreate( HD_MotionMesh, 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 (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 (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 (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 (ErrStat3 >= AbortErrLev) return + endif + + end subroutine SetMotionLoadsInterfaceMeshes + + !------------------------------------------------------------- + !> Sanity check the nodes + !! If more than one input node was passed in, but only a single HD node + !! exits (single Morison or single WAMIT), then give error that too many + !! nodes passed. + subroutine CheckNodes(ErrStat3,ErrMsg3) + integer(IntKi), intent( out) :: ErrStat3 !< temporary error status + character(ErrMsgLen), intent( out) :: ErrMsg3 !< temporary error message + ErrStat3 = ErrID_None + ErrMsg3 = "" +!TODO +! NumNodePts > y%WAMIT nodes +! or if no WAMIT, NumNodePts > y%Morison nodes +!I'm not exactly sure on the logic yet. + end subroutine CheckNodes + + !------------------------------------------------------------- + !> Sanity check the Morison mesh + !! If the Morison mesh has points near the bottom of the waterdepth, + !! then we could have some very strange results with only a single mesh + !! point for the HD_MotionMesh. + subroutine CheckDepth(ErrStat3,ErrMsg3) + integer(IntKi), intent( out) :: ErrStat3 !< temporary error status + character(ErrMsgLen), intent( out) :: ErrMsg3 !< temporary error message + ErrStat3 = ErrID_None + ErrMsg3 = "" + if ( NumNodePts == 1 ) then +!TODO: +!check if committed +! loop over nodes and find lowest +! check if lowest is really deep (close to wtrdepth) +! give severe warning if is case + endif + end subroutine CheckDepth + END SUBROUTINE HydroDyn_Init_c @@ -274,10 +528,8 @@ SUBROUTINE HydroDyn_CalcOutput_c(Time_C,OutputChannelValues_C,ErrStat_C,ErrMsg_C !FIXME: Reshape position, velocity, acceleration and set mesh accordingly. -!FIXME:Set inputs array and extrap/interp necessary - ! Call the main subroutine HydroDyn_CalcOutput to get the resulting forces and moments at time T - CALL HydroDyn_CalcOutput( Time, u, p, x, xd, z, OtherStates, y, m, ErrStat2, ErrMsg2 ) + CALL HydroDyn_CalcOutput( Time, u(1), p, x, xd, z, OtherStates, y, m, ErrStat2, ErrMsg2 ) if (Failed()) return ! Get the output channel info out of y @@ -378,7 +630,7 @@ SUBROUTINE HydroDyn_End_C(ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_End_c') ! or some other code using the c-bindings. if (allocated(u)) then do i=2,size(u) ! leave first one for passing to HD_End for destruction - call HydroDyn_DestroyInput( u, ErrStat2, ErrMsg2 ) + call HydroDyn_DestroyInput( u(i), ErrStat2, ErrMsg2 ) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) enddo endif @@ -395,8 +647,30 @@ SUBROUTINE HydroDyn_End_C(ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_End_c') ! if u is still allocated, deallocate it now if (allocated(u)) deallocate(u) - call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + ! Clear out mesh related data storage + call ClearMesh() + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) +CONTAINS + !> Don't leave junk in memory. So destroy meshes and mappings. + subroutine ClearMesh() + ! Destroy connection meshes + call MeshDestroy( HD_MotionMesh, ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call MeshDestroy( HD_LoadMesh, ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + ! Destroy mesh mappings + call NWTC_Library_Destroymeshmaptype( Map_Motion_2_HD_PRP_P, ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call NWTC_Library_Destroymeshmaptype( Map_Motion_2_HD_WB_P , ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call NWTC_Library_Destroymeshmaptype( Map_Motion_2_HD_Mo_P , ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call NWTC_Library_Destroymeshmaptype( Map_HD_WB_P_2_Load , ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call NWTC_Library_Destroymeshmaptype( Map_HD_Mo_P_2_Load , ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + end subroutine ClearMesh END SUBROUTINE HydroDyn_End_c END MODULE HydroDyn_C diff --git a/reg_tests/r-test b/reg_tests/r-test index 3dcf5751d0..8ea661f075 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 3dcf5751d03c4a049a2d0f58063ff88ebe239a14 +Subproject commit 8ea661f07534ee9fa8a8b243b069af81aae8a243 From 6f38412dcc7f83fec9777855bf3877f196986258 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Mon, 14 Jun 2021 14:28:39 -0600 Subject: [PATCH 07/22] hd_c_lib: add error checks at init, pass rootname --- .../hydrodyn/python-lib/hydrodyn_library.py | 29 ++++-- modules/hydrodyn/src/HydroDyn_C.f90 | 96 +++++++++++++------ 2 files changed, 88 insertions(+), 37 deletions(-) diff --git a/modules/hydrodyn/python-lib/hydrodyn_library.py b/modules/hydrodyn/python-lib/hydrodyn_library.py index 5e0fc436bf..07b4650721 100644 --- a/modules/hydrodyn/python-lib/hydrodyn_library.py +++ b/modules/hydrodyn/python-lib/hydrodyn_library.py @@ -96,6 +96,10 @@ class HydroDynLib(CDLL): # here. error_msg_c_len = 1025 + # NOTE: the length of the name used for any output file written by the + # HD Fortran code is 1025. + default_str_c_len = 1025 + def __init__(self, library_path): super().__init__(library_path) self.library_path = library_path @@ -149,9 +153,15 @@ def __init__(self, library_path): self.numNodePts = 1 # Single ptfm attachment point for floating rigid self.initNodeLocations = np.zeros((self.numNodePts,6)) # N x 6 array [x,y,z,Rx,Ry,Rz] + # OutRootName + # If HD writes a file (echo, summary, or other), use this for the + # root of the file name. + self.outRootName = "Output_HDlib_default" + # _initialize_routines() ------------------------------------------------------------------------------------------------------------ def _initialize_routines(self): self.HydroDyn_Init_c.argtypes = [ + POINTER(c_char), # OutRootName POINTER(c_char_p), # input file string POINTER(c_int), # input file string length POINTER(c_float), # gravity @@ -199,7 +209,12 @@ def hydrodyn_init(self, input_string_array): self._numChannels_c = c_int(0) - #FIXME: for now we only allow for a single rigid motion platform. + # Rootname for HD output files (echo etc). + _outRootName_c = create_string_buffer((self.outRootName.ljust(self.default_str_c_len)).encode('utf-8')) + + + # initNodeLocations + # Verify that the shape of initNodeLocations is correct if self.initNodeLocations.shape[1] != 6: print("Expecting a Nx6 array of initial node locations (initNodeLocations) with second index for [x,y,z,Rx,Ry,Rz]") self.hydrodyn_end() @@ -208,15 +223,16 @@ def hydrodyn_init(self, input_string_array): print("Expecting a Nx6 array of initial node locations (initNodeLocations) with first index for number of nodes.") self.hydrodyn_end() raise Exception("\nHydroDyn terminated prematurely.") - # Placeholders exist for multiple input points for flexible - # platforms, but are for now not enabled. + #FIXME: for now we only allow for a single rigid motion platform. If + # multiple nodes are allowed for substructure coupling or + # flexible floating platforms, then remove this check and verify + # the fortran side works as expected (theoretically should, but + # it is untested). if self.numNodePts != 1: print(f"HydroDyn C interface does not currently support flexible structures ({self.numNodePts} input motion points were requested, but only 1 is currently supported).") self.hydrodyn_end() raise Exception("\nHydroDyn terminated prematurely.") - # initNodeLocations - # Verify that the shape of initNodeLocations is correct # Make a flat 1D array of position info: # [x2,y1,z1,Rx1,Ry1,Rz1, x2,y2,z2,Rx2,Ry2,Rz2 ...] nodeInitLoc_flat = [pp for p in self.initNodeLocations for pp in p] @@ -224,9 +240,10 @@ def hydrodyn_init(self, input_string_array): for i, p in enumerate(nodeInitLoc_flat): nodeInitLoc_flat_c[i] = c_float(p) - #FIXME: need to pass initial position and orientation info + # call HydroDyn_Init_c self.HydroDyn_Init_c( + _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 byref(c_float(self.gravity)), # IN: gravity diff --git a/modules/hydrodyn/src/HydroDyn_C.f90 b/modules/hydrodyn/src/HydroDyn_C.f90 index be57af1f45..6f2d88fde4 100644 --- a/modules/hydrodyn/src/HydroDyn_C.f90 +++ b/modules/hydrodyn/src/HydroDyn_C.f90 @@ -36,7 +36,8 @@ MODULE HydroDyn_C ! This must exactly match the value in the python-lib. If ErrMsgLen changes at ! some point in the nwtc-library, this should be updated, but the logic exists ! to correctly handle different lengths of the strings - integer(IntKi), parameter :: ErrMsgLen_C=1025 + integer(IntKi), parameter :: ErrMsgLen_C = 1025 + integer(IntKi), parameter :: IntfStrLen = 1025 ! length of other strings through the C interface !------------------------------------------------------------------------------------ @@ -151,13 +152,13 @@ end subroutine SetErr !=============================================================================================================== !--------------------------------------------- HydroDyn Init---------------------------------------------------- !=============================================================================================================== -SUBROUTINE HydroDyn_Init_c( InputFileString_C, InputFileStringLength_C, & - Gravity_C, defWtrDens_C, defWtrDpth_C, defMSL2SWL_C, & - PtfmRefPtPositionX_C, PtfmRefPtPositionY_C, & - NumNodePts_C, InitNodePositions_C, & - !NumWaveElev_C, WaveElevXY_C & !Placeholder for later - InterpOrder_C, DT_C, TMax_C, & - NumChannels_C, OutputChannelNames_C, OutputChannelUnits_C, & +SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLength_C, & + Gravity_C, defWtrDens_C, defWtrDpth_C, defMSL2SWL_C, & + PtfmRefPtPositionX_C, PtfmRefPtPositionY_C, & + NumNodePts_C, InitNodePositions_C, & + !NumWaveElev_C, WaveElevXY_C & !Placeholder for later + InterpOrder_C, DT_C, TMax_C, & + NumChannels_C, OutputChannelNames_C, OutputChannelUnits_C, & ErrStat_C, ErrMsg_C) BIND (C, NAME='HydroDyn_Init_c') implicit none #ifndef IMPLICIT_DLLEXPORT @@ -165,6 +166,7 @@ SUBROUTINE HydroDyn_Init_c( InputFileString_C, InputFileStringLength_C, !GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_Init_c #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 real(c_float), intent(in ) :: Gravity_C !< Gravitational constant (set by calling code) @@ -188,6 +190,7 @@ SUBROUTINE HydroDyn_Init_c( InputFileString_C, InputFileStringLength_C, 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(ReKi), allocatable :: InitNodePositions(:,:) !< temp array. Probably don't need this, but makes conversion from C clearer. @@ -222,10 +225,8 @@ SUBROUTINE HydroDyn_Init_c( InputFileString_C, InputFileStringLength_C, ! 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%OutRootName = "HDpassed" ! FIXME: pass in a name to use here (For summary files etc) InitInp%UseInputFile = .FALSE. ! this probably should be passed in InitInp%HasIce = .FALSE. ! Always keep at false unless interfacing to ice modules ! Linearization @@ -234,6 +235,12 @@ SUBROUTINE HydroDyn_Init_c( InputFileString_C, InputFileStringLength_C, ! 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%defWtrDens = REAL(defWtrDens_C, ReKi) @@ -267,12 +274,11 @@ SUBROUTINE HydroDyn_Init_c( InputFileString_C, InputFileStringLength_C, ! 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. For now, skipping this feature. + ! and could be added to this inteface. ! Skipping this for now. Maybe add later. !InitInp%WaveElevXY - !----------------------- ! Allocate input array u !----------------------- @@ -305,24 +311,29 @@ SUBROUTINE HydroDyn_Init_c( InputFileString_C, InputFileStringLength_C, 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 pieces of u(:) !------------------------------------------------------------- -! copy u(1) to others -! what from InitOutData should be returned -! what error handling is necessary + do i=2,InterpOrder+1 + call HydroDyn_CopyInput (u(1), u(i), MESH_NEWCOPY, Errstat2, ErrMsg2) + if (Failed()) return + enddo + +!TODO +! Is there any other InitOutData should be returned +! Any additional warnings or error handling necessary - !-------------------------------------------------- - ! Set output channel informatoion for driver code. - !-------------------------------------------------- + !------------------------------------------------- + ! Set output channel information for driver code + !------------------------------------------------- ! Number of channels NumChannels_C = size(InitOutData%WriteOutputHdr) @@ -410,7 +421,6 @@ subroutine SetMotionLoadsInterfaceMeshes(ErrStat3,ErrMsg3) ! note: CU is is output unit (platform dependent). !call MeshPrintInfo( CU, HD_MotionMesh ) - !------------------------------------------------------------- ! Loads mesh ! This point mesh may contain more than one point. Mapping will be used to map @@ -468,28 +478,52 @@ subroutine CheckNodes(ErrStat3,ErrMsg3) character(ErrMsgLen), intent( out) :: ErrMsg3 !< temporary error message ErrStat3 = ErrID_None ErrMsg3 = "" -!TODO -! NumNodePts > y%WAMIT nodes -! or if no WAMIT, NumNodePts > y%Morison nodes -!I'm not exactly sure on the logic yet. + 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 + 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 + 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 + ErrStat3 = ErrID_Fatal + ErrMsg3 = "More nodes passed into library than exist in HydroDyn model WAMIT mesh" + endif + endif + endif end subroutine CheckNodes !------------------------------------------------------------- !> Sanity check the Morison mesh !! If the Morison mesh has points near the bottom of the waterdepth, !! then we could have some very strange results with only a single mesh - !! point for the HD_MotionMesh. + !! point for the HD_MotionMesh. All WAMIT mesh points should be near + !! the surface, so no checking is necessary. subroutine CheckDepth(ErrStat3,ErrMsg3) integer(IntKi), intent( out) :: ErrStat3 !< temporary error status character(ErrMsgLen), intent( out) :: ErrMsg3 !< temporary error message + real(ReKi) :: tmpZpos !< temporary z-position ErrStat3 = ErrID_None ErrMsg3 = "" - if ( NumNodePts == 1 ) then -!TODO: -!check if committed -! loop over nodes and find lowest -! check if lowest is really deep (close to wtrdepth) -! give severe warning if is case + 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 + ! Find lowest Morison node + if (u(1)%Morison%Mesh%Position(3,i) < tmpZpos) then + tmpZpos = u(1)%Morison%Mesh%Position(3,i) + endif + enddo + if (tmpZpos < -abs(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."// & + NewLine//" ---- Results may not make physical sense ----" + endif endif end subroutine CheckDepth From 8735e14bd8bfe5f750e95e6bef53828ade75d2cc Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Tue, 15 Jun 2021 12:27:09 -0600 Subject: [PATCH 08/22] hd_c_lib: add mesh transfers, CalcOutput interface, DBGoutput writing --- .../hydrodyn/python-lib/hydrodyn_library.py | 216 ++++++++++++++---- modules/hydrodyn/src/HydroDyn_C.f90 | 210 ++++++++++++++--- reg_tests/r-test | 2 +- 3 files changed, 361 insertions(+), 67 deletions(-) diff --git a/modules/hydrodyn/python-lib/hydrodyn_library.py b/modules/hydrodyn/python-lib/hydrodyn_library.py index 07b4650721..34219936f4 100644 --- a/modules/hydrodyn/python-lib/hydrodyn_library.py +++ b/modules/hydrodyn/python-lib/hydrodyn_library.py @@ -146,12 +146,12 @@ def __init__(self, library_path): # Nodes # The number of nodes must be constant throughout simulation. The - # initial position is given in the initNodeLocations array (resize as + # initial position is given in the initNodePos array (resize as # needed, should be Nx6). # Rotations are given in radians assuming small angles. See note at # top of this file. self.numNodePts = 1 # Single ptfm attachment point for floating rigid - self.initNodeLocations = np.zeros((self.numNodePts,6)) # N x 6 array [x,y,z,Rx,Ry,Rz] + self.initNodePos = np.zeros((self.numNodePts,6)) # N x 6 array [x,y,z,Rx,Ry,Rz] # OutRootName # If HD writes a file (echo, summary, or other), use this for the @@ -171,7 +171,7 @@ def _initialize_routines(self): POINTER(c_float), # PtfmRefPt_X POINTER(c_float), # PtfmRefPt_Y POINTER(c_int), # numNodePts -- number of points expecting motions/loads - POINTER(c_float), # initNodeLocations -- initial node positions in flat array of 6*numNodePts + POINTER(c_float), # initNodePos -- initial node positions in flat array of 6*numNodePts POINTER(c_int), # InterpOrder POINTER(c_double), # dt POINTER(c_double), # tmax @@ -185,6 +185,11 @@ def _initialize_routines(self): self.HydroDyn_CalcOutput_c.argtypes = [ POINTER(c_double), # Time_C + POINTER(c_int), # numNodePts -- number of points expecting motions/loads + POINTER(c_float), # nodePos -- node positions in flat array of 6*numNodePts + POINTER(c_float), # nodeVel -- node velocities in flat array of 6*numNodePts + POINTER(c_float), # nodeAcc -- node accelerations in flat array of 6*numNodePts + POINTER(c_float), # nodeFrc -- node forces/moments in flat array of 6*numNodePts POINTER(c_float), # Output Channel Values POINTER(c_int), # ErrStat_C POINTER(c_char) # ErrMsg_C @@ -212,15 +217,17 @@ def hydrodyn_init(self, input_string_array): # Rootname for HD output files (echo etc). _outRootName_c = create_string_buffer((self.outRootName.ljust(self.default_str_c_len)).encode('utf-8')) + # store initial number of node points for error handling at calls + self._initNumNodePts = self.numNodePts - # initNodeLocations - # Verify that the shape of initNodeLocations is correct - if self.initNodeLocations.shape[1] != 6: - print("Expecting a Nx6 array of initial node locations (initNodeLocations) with second index for [x,y,z,Rx,Ry,Rz]") + # initNodePos + # Verify that the shape of initNodePos is correct + if self.initNodePos.shape[1] != 6: + print("Expecting a Nx6 array of initial node locations (initNodePos) with second index for [x,y,z,Rx,Ry,Rz]") self.hydrodyn_end() raise Exception("\nHydroDyn terminated prematurely.") - if self.initNodeLocations.shape[0] != self.numNodePts: - print("Expecting a Nx6 array of initial node locations (initNodeLocations) with first index for number of nodes.") + if self.initNodePos.shape[0] != self.numNodePts: + print("Expecting a Nx6 array of initial node locations (initNodePos) with first index for node number.") self.hydrodyn_end() raise Exception("\nHydroDyn terminated prematurely.") #FIXME: for now we only allow for a single rigid motion platform. If @@ -235,7 +242,7 @@ def hydrodyn_init(self, input_string_array): # Make a flat 1D array of position info: # [x2,y1,z1,Rx1,Ry1,Rz1, x2,y2,z2,Rx2,Ry2,Rz2 ...] - nodeInitLoc_flat = [pp for p in self.initNodeLocations for pp in p] + nodeInitLoc_flat = [pp for p in self.initNodePos for pp in p] nodeInitLoc_flat_c = (c_float * (6 * self.numNodePts))(0.0,) for i, p in enumerate(nodeInitLoc_flat): nodeInitLoc_flat_c[i] = c_float(p) @@ -253,7 +260,7 @@ def hydrodyn_init(self, input_string_array): byref(c_float(self.ptfmRefPt_x)), # IN: Platform initial position (X) byref(c_float(self.ptfmRefPt_y)), # IN: Platform initial position (Y) byref(c_int(self.numNodePts)), # IN: number of attachment points expected (where motions are transferred into HD) - nodeInitLoc_flat_c, # IN: initNodeLocations -- initial node positions in flat array of 6*numNodePts + nodeInitLoc_flat_c, # IN: initNodePos -- initial node positions in flat array of 6*numNodePts byref(c_int(self.InterpOrder)), # IN: InterpOrder (1: linear, 2: quadratic) byref(c_double(self.dt)), # IN: time step (dt) byref(c_double(self.tmax)), # IN: tmax @@ -271,39 +278,65 @@ def hydrodyn_init(self, input_string_array): # hydrodyn_calcOutput ------------------------------------------------------------------------------------------------------------ - #def hydrodyn_calcOutput(self, time, positions, velocities, outputChannelValues): - def hydrodyn_calcOutput(self, time, outputChannelValues): + def hydrodyn_calcOutput(self, time, nodePos, nodeVel, nodeAcc, nodeFrcMom, outputChannelValues): + + # Check input motion info + self.check_input_motions(nodePos,nodeVel,nodeAcc) + + # set flat arrays for inputs of motion + # Position -- [x2,y1,z1,Rx1,Ry1,Rz1, x2,y2,z2,Rx2,Ry2,Rz2 ...] + nodePos_flat = [pp for p in nodePos for pp in p] + nodePos_flat_c = (c_float * (6 * self.numNodePts))(0.0,) + for i, p in enumerate(nodePos_flat): + nodePos_flat_c[i] = c_float(p) + + # Velocity -- [Vx2,Vy1,Vz1,RVx1,RVy1,RVz1, Vx2,Vy2,Vz2,RVx2,RVy2,RVz2 ...] + nodeVel_flat = [pp for p in nodeVel for pp in p] + nodeVel_flat_c = (c_float * (6 * self.numNodePts))(0.0,) + for i, p in enumerate(nodeVel_flat): + nodeVel_flat_c[i] = c_float(p) - # Set up inputs + # Acceleration -- [Ax1,Ay1,Az1,RAx1,RAy1,RAz1, Ax2,Ay2,Az2,RAx2,RAy2,RAz2 ...] + nodeAcc_flat = [pp for p in nodeAcc for pp in p] + nodeAcc_flat_c = (c_float * (6 * self.numNodePts))(0.0,) + for i, p in enumerate(nodeAcc_flat): + nodeAcc_flat_c[i] = c_float(p) + + # Resulting Forces/moments -- [Fx1,Fy1,Fz1,Mx1,My1,Mz1, Fx2,Fy2,Fz2,Mx2,My2,Mz2 ...] + nodeFrc_flat_c = (c_float * (6 * self.numNodePts))(0.0,) # Set up output channels outputChannelValues_c = (c_float * self.numChannels)(0.0,) # Run HydroDyn_CalcOutput_c self.HydroDyn_CalcOutput_c( - byref(c_double(time)), # IN: time at which to calculate velocities - #positions_flat_c, # IN: positions - specified by user - #velocities_flat_c, # IN: velocities at desired positions - #velocities_flat_c, # IN: accelerations at desired positions - #forces_flat_c, # OUT: resulting forces/moments array - 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 + byref(c_double(time)), # IN: time at which to calculate velocities + byref(c_int(self.numNodePts)), # IN: number of attachment points expected (where motions are transferred into HD) + nodePos_flat_c, # IN: positions - specified by user + nodeVel_flat_c, # IN: velocities at desired positions + nodeAcc_flat_c, # IN: accelerations at desired positions + nodeFrc_flat_c, # OUT: resulting forces/moments array + 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 ) self.check_error() + ## Reshape Force/Moment into [N,6] + count = 0 + for j in range(0,self.numNodePts): + nodeFrcMom[j,0] = nodeFrc_flat_c[count] + nodeFrcMom[j,1] = nodeFrc_flat_c[count+1] + nodeFrcMom[j,2] = nodeFrc_flat_c[count+2] + nodeFrcMom[j,3] = nodeFrc_flat_c[count+3] + nodeFrcMom[j,4] = nodeFrc_flat_c[count+4] + nodeFrcMom[j,5] = nodeFrc_flat_c[count+5] + count = count + 6 + # Convert output channel values back into python for k in range(0,self.numChannels): outputChannelValues[k] = float(outputChannelValues_c[k]) - - ## Reshape velocities into [N,3] - #count = 0 - #for j in range(0,self.numWindPts): - # velocities[j,0] = velocities_flat_c[count] - # velocities[j,1] = velocities_flat_c[count+1] - # velocities[j,2] = velocities_flat_c[count+2] - # count = count + 3 # hydrodyn_end ------------------------------------------------------------------------------------------------------------ def hydrodyn_end(self): @@ -328,6 +361,48 @@ def check_error(self): self.hydrodyn_end() raise Exception("\nHydroDyn terminated prematurely.") + + def check_input_motions(self,nodePos,nodeVel,nodeAcc): + # make sure number of nodes didn't change for some reason + if self._initNumNodePts != self.numNodePts: + print(f"At time {time}, the number of node points changed from initial value of {self._initNumNodePts}. This is not permitted during the simulation.") + self.hydrodyn_end() + raise Exception("\nError in calling HydroDyn library.") + + # Verify that the shape of positions array is correct + if nodePos.shape[1] != 6: + print("Expecting a Nx6 array of node positions (nodePos) with second index for [x,y,z,Rx,Ry,Rz]") + self.hydrodyn_end() + raise Exception("\nHydroDyn terminated prematurely.") + if nodePos.shape[0] != self.numNodePts: + print("Expecting a Nx6 array of node positions (nodePos) with first index for node number.") + self.hydrodyn_end() + raise Exception("\nHydroDyn terminated prematurely.") + + + # Verify that the shape of velocities array is correct + if nodeVel.shape[1] != 6: + print("Expecting a Nx6 array of node velocities (nodeVel) with second index for [x,y,z,Rx,Ry,Rz]") + self.hydrodyn_end() + raise Exception("\nHydroDyn terminated prematurely.") + if nodeVel.shape[0] != self.numNodePts: + print("Expecting a Nx6 array of node velocities (nodeVel) with first index for node number.") + self.hydrodyn_end() + raise Exception("\nHydroDyn terminated prematurely.") + + + # Verify that the shape of accelerations array is correct + if nodeAcc.shape[1] != 6: + print("Expecting a Nx6 array of node accelerations (nodeAcc) with second index for [x,y,z,Rx,Ry,Rz]") + self.hydrodyn_end() + raise Exception("\nHydroDyn terminated prematurely.") + if nodeAcc.shape[0] != self.numNodePts: + print("Expecting a Nx6 array of node accelerations (nodeAcc) with first index for node number.") + self.hydrodyn_end() + raise Exception("\nHydroDyn terminated prematurely.") + + + @property def output_channel_names(self): if len(self._channel_names_c.value.split()) == 0: @@ -360,22 +435,83 @@ class DriverDbg(): When coupled to another code, the force/moment array would be passed back to the calling code for use in the structural solver. """ - def __init__(self,filename): + def __init__(self,filename,numNodePts): self.DbgFile=open(filename,'wt') # open output file and write header info + self.numNodePts=numNodePts # write file header t_string=datetime.datetime.now() dt_string=datetime.date.today() -# self.DbgFile.write(f"## This file was generated by InflowWind_Driver on {dt_string.strftime('%b-%d-%Y')} at {t_string.strftime('%H:%M:%S')}\n") -# self.DbgFile.write(f"## This file contains the wind velocity at the {numWindPts} points specified in the file ") -# self.DbgFile.write(f"{filename}\n") -# self.DbgFile.write("#\n") -# self.DbgFile.write("# T X Y Z U V W\n") -# self.DbgFile.write("# (s) (m) (m) (m) (m/s) (m/s) (m/s)\n") + self.DbgFile.write(f"## This file was generated by hydrodyn_c_lib on {dt_string.strftime('%b-%d-%Y')} at {t_string.strftime('%H:%M:%S')}\n") + self.DbgFile.write(f"## This file contains the resulting forces/moments at each of {self.numNodePts} node(s) passed into the hydrodyn_c_lib\n") + self.DbgFile.write("#\n") + f_string = "{:^25s}" + self.DbgFile.write("# T ") + for i in range(1,self.numNodePts+1): + f_num = "N{0:04d}_".format(i) + self.DbgFile.write(f_string.format(f_num+"x" )) + self.DbgFile.write(f_string.format(f_num+"y" )) + self.DbgFile.write(f_string.format(f_num+"z" )) + self.DbgFile.write(f_string.format(f_num+"Rx" )) + self.DbgFile.write(f_string.format(f_num+"Ry" )) + self.DbgFile.write(f_string.format(f_num+"Rz" )) + self.DbgFile.write(f_string.format(f_num+"Vx" )) + self.DbgFile.write(f_string.format(f_num+"Vy" )) + self.DbgFile.write(f_string.format(f_num+"Vz" )) + self.DbgFile.write(f_string.format(f_num+"RVx")) + self.DbgFile.write(f_string.format(f_num+"RVy")) + self.DbgFile.write(f_string.format(f_num+"RVz")) + self.DbgFile.write(f_string.format(f_num+"Ax" )) + self.DbgFile.write(f_string.format(f_num+"Ay" )) + self.DbgFile.write(f_string.format(f_num+"Az" )) + self.DbgFile.write(f_string.format(f_num+"RAx")) + self.DbgFile.write(f_string.format(f_num+"RAy")) + self.DbgFile.write(f_string.format(f_num+"RAz")) + self.DbgFile.write(f_string.format(f_num+"Fx" )) + self.DbgFile.write(f_string.format(f_num+"Fy" )) + self.DbgFile.write(f_string.format(f_num+"Fz" )) + self.DbgFile.write(f_string.format(f_num+"Mx" )) + self.DbgFile.write(f_string.format(f_num+"My" )) + self.DbgFile.write(f_string.format(f_num+"Mz" )) + self.DbgFile.write("\n") + self.DbgFile.write("# (s) ") + for i in range(1,self.numNodePts+1): + self.DbgFile.write(f_string.format("(m)" )) + self.DbgFile.write(f_string.format("(m)" )) + self.DbgFile.write(f_string.format("(m)" )) + self.DbgFile.write(f_string.format("(rad)" )) + self.DbgFile.write(f_string.format("(rad)" )) + self.DbgFile.write(f_string.format("(rad)" )) + self.DbgFile.write(f_string.format("(m/s)" )) + self.DbgFile.write(f_string.format("(m/s)" )) + self.DbgFile.write(f_string.format("(m/s)" )) + self.DbgFile.write(f_string.format("(rad/s)" )) + self.DbgFile.write(f_string.format("(rad/s)" )) + self.DbgFile.write(f_string.format("(rad/s)" )) + self.DbgFile.write(f_string.format("(m/s^2)" )) + self.DbgFile.write(f_string.format("(m/s^2)" )) + self.DbgFile.write(f_string.format("(m/s^2)" )) + self.DbgFile.write(f_string.format("(rad/s^2)")) + self.DbgFile.write(f_string.format("(rad/s^2)")) + self.DbgFile.write(f_string.format("(rad/s^2)")) + self.DbgFile.write(f_string.format("(N)" )) + self.DbgFile.write(f_string.format("(N)" )) + self.DbgFile.write(f_string.format("(N)" )) + self.DbgFile.write(f_string.format("(N-m)" )) + self.DbgFile.write(f_string.format("(N-m)" )) + self.DbgFile.write(f_string.format("(N-m)" )) + self.DbgFile.write("\n") self.opened = True -# def write(self,t,positions,velocities): -# for p, v in zip(positions,velocities): -# self.DbgFile.write(' %14.7f %14.7f %14.7f %14.7f %14.7f %14.7f %14.7f\n' % (t,p[0],p[1],p[2],v[0],v[1],v[2])) + def write(self,t,nodePos,nodeVel,nodeAcc,nodeFrc): + t_string = "{:10.4f}" + f_string = "{:25.7f}"*6 + self.DbgFile.write(t_string.format(t)) + for i in range(0,self.numNodePts): + self.DbgFile.write(f_string.format(*nodePos[i,:])) + self.DbgFile.write(f_string.format(*nodeVel[i,:])) + self.DbgFile.write(f_string.format(*nodeAcc[i,:])) + self.DbgFile.write(f_string.format(*nodeFrc[i,:])) + self.DbgFile.write("\n") def end(self): if self.opened: diff --git a/modules/hydrodyn/src/HydroDyn_C.f90 b/modules/hydrodyn/src/HydroDyn_C.f90 index 6f2d88fde4..68d81a2d49 100644 --- a/modules/hydrodyn/src/HydroDyn_C.f90 +++ b/modules/hydrodyn/src/HydroDyn_C.f90 @@ -113,6 +113,7 @@ MODULE HydroDyn_C integer(IntKi) :: NumNodePts ! Number of mesh points we are interfacing motions/loads to/from HD type(MeshType) :: HD_MotionMesh ! mesh for motions of external nodes type(MeshType) :: HD_LoadMesh ! mesh for loads for external nodes + type(MeshType) :: HD_LoadMesh_tmp ! mesh for loads for external nodes -- temporary !------------------------------ ! Mesh mapping: motions ! The mapping of motions from the nodes passed in to the corresponding HD meshes @@ -126,6 +127,11 @@ MODULE HydroDyn_C type(MeshMapType) :: Map_HD_Mo_P_2_Load ! Mesh mapping between HD output Morison loads mesh and external nodes mesh ! Meshes -- helper stuff real(R8Ki) :: theta(3) ! mesh creation helper data + ! Motions input (so we don't have to reallocate all the time + real(ReKi), allocatable :: tmpNodePos(:,:) ! temp array. Probably don't need this, but makes conversion from C clearer. + real(ReKi), allocatable :: tmpNodeVel(:,:) ! temp array. Probably don't need this, but makes conversion from C clearer. + real(ReKi), allocatable :: tmpNodeAcc(:,:) ! temp array. Probably don't need this, but makes conversion from C clearer. + real(ReKi), allocatable :: tmpNodeFrc(:,:) ! temp array. Probably don't need this, but makes conversion to C clearer. !------------------------------------------------------------------------------------ @@ -176,7 +182,7 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen real(c_float), intent(in ) :: PtfmRefPtPositionX_C !< Initial position in wave field real(c_float), intent(in ) :: PtfmRefPtPositionY_C !< Initial position in wave field integer(c_int), intent(in ) :: NumNodePts_C !< Number of mesh points we are transfering motions to and output loads to - real(c_float), intent(in ) :: InitNodePositions_C( 6, NumNodePts_C ) !< A 6xNumNodePts_C array [x,y,z,Rx,Ry,Rz] + real(c_float), intent(in ) :: InitNodePositions_C( 6*NumNodePts_C ) !< A 6xNumNodePts_C array [x,y,z,Rx,Ry,Rz] !NOTE: not setting up the WaveElev at this point. Leaving placeholder for future !integer(c_int), intent(in ) :: NumWaveElev_C !< Number of mesh points we are transfering motions to and output loads to !real(c_float), intent(in ) :: WaveElevXY_C !< A 2xNumWaveElev_C array [x,y] @@ -192,7 +198,6 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen ! 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(ReKi), allocatable :: InitNodePositions(:,:) !< temp array. Probably don't need this, but makes conversion from C clearer. real(DbKi) :: TimeInterval !< timestep for HD integer(IntKi) :: ErrStat !< aggregated error message @@ -258,11 +263,12 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen ErrMsg2 = "At least one node point must be specified" if (Failed()) return endif - call AllocAry( InitNodePositions, 6, NumNodePts, "InitNodePositions", ErrStat2, ErrMsg2 ) - if (Failed()) return - do i=1,NumNodePts - InitNodePositions(1:6,i) = real(InitNodePositions_C(1:6,i),ReKi) - enddo + ! Allocate temporary arrays to simplify data conversions + call AllocAry( tmpNodePos, 6, NumNodePts, "tmpNodePos", ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( tmpNodeVel, 6, NumNodePts, "tmpNodeVel", ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( tmpNodeAcc, 6, NumNodePts, "tmpNodeAcc", ErrStat2, ErrMsg2 ); if (Failed()) return + 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) @@ -353,7 +359,6 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen OutputChannelUnits_C(k) = C_NULL_CHAR - call Cleanup() call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) CONTAINS @@ -361,14 +366,17 @@ logical function Failed() CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) Failed = ErrStat >= AbortErrLev if (Failed) then - call Cleanup() + call FailCleanup() call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) endif end function Failed - subroutine Cleanup() - if (allocated(InitNodePositions)) deallocate(InitNodePositions) - end subroutine Cleanup + subroutine FailCleanup() + if (allocated(tmpNodePos)) deallocate(tmpNodePos) + if (allocated(tmpNodeVel)) deallocate(tmpNodeVel) + if (allocated(tmpNodeAcc)) deallocate(tmpNodeAcc) + if (allocated(tmpNodeFrc)) deallocate(tmpNodeFrc) + end subroutine FailCleanup !> This subroutine sets the interface meshes to map to the input motions to the HD !! meshes for WAMIT and Morison. This subroutine also sets the meshes for the @@ -398,8 +406,8 @@ subroutine SetMotionLoadsInterfaceMeshes(ErrStat3,ErrMsg3) do iNode=1,NumNodePts ! initial position and orientation of node - InitPos = InitNodePositions(1:3,iNode) - theta = real(InitNodePositions(4:6,iNode),DbKi) ! convert ReKi to DbKi to avoid roundoff + InitPos = tmpNodePos(1:3,iNode) + theta = real(tmpNodePos(4:6,iNode),DbKi) ! convert ReKi to DbKi to avoid roundoff Orient = EulerConstruct( theta ) call MeshPositionNode( HD_MotionMesh , & iNode , & @@ -442,6 +450,27 @@ subroutine SetMotionLoadsInterfaceMeshes(ErrStat3,ErrMsg3) ! note: CU is is output unit (platform dependent). !call MeshPrintInfo( CU, HD_LoadMesh ) + !------------------------------------------------------------- + ! Loads mesh + ! This point mesh may contain more than one point. Mapping will be used to map + ! the loads from output meshes for WAMIT and Morison. + ! Output mesh for loads at each WAMIT body + CALL MeshCopy( SrcMesh = HD_LoadMesh ,& + DestMesh = HD_LoadMesh_tmp ,& + CtrlCode = MESH_COUSIN ,& + IOS = COMPONENT_OUTPUT ,& + ErrStat = ErrStat3 ,& + ErrMess = ErrMsg3 ,& + Force = .TRUE. ,& + Moment = .TRUE. ) + if (ErrStat3 >= AbortErrLev) return + + HD_LoadMesh_tmp%RemapFlag = .TRUE. + + ! For checking the mesh, uncomment this. + ! note: CU is is output unit (platform dependent). + !call MeshPrintInfo( CU, HD_LoadMesh_tmp ) + !------------------------------------------------------------- ! Set the mapping meshes ! PRP - principle reference point @@ -534,41 +563,77 @@ END SUBROUTINE HydroDyn_Init_c !--------------------------------------------- HydroDyn CalcOutput --------------------------------------------- !=============================================================================================================== -SUBROUTINE HydroDyn_CalcOutput_c(Time_C,OutputChannelValues_C,ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_CalcOutput_c') +SUBROUTINE HydroDyn_CalcOutput_c(Time_C, NumNodePts_C, NodePos_C, NodeVel_C, NodeAcc_C, & + NodeFrc_C, OutputChannelValues_C, ErrStat_C, ErrMsg_C) BIND (C, NAME='HydroDyn_CalcOutput_c') implicit none #ifndef IMPLICIT_DLLEXPORT !DEC$ ATTRIBUTES DLLEXPORT :: HydroDyn_CalcOutput_c !GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_CalcOutput_c #endif - real(c_double), intent(in ) :: Time_C - real(c_float), intent( out) :: OutputChannelValues_C(p%NumOuts) - integer(c_int), intent( out) :: ErrStat_C - character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + real(c_double), intent(in ) :: Time_C + integer(c_int), intent(in ) :: NumNodePts_C !< Number of mesh points we are transfering motions to and output loads to + real(c_float), intent(in ) :: NodePos_C( 6*NumNodePts_C ) !< A 6xNumNodePts_C array [x,y,z,Rx,Ry,Rz] -- positions (global) + 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) + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) ! Local variables - real(DbKi) :: Time - integer(IntKi) :: ErrStat !< aggregated error status - character(ErrMsgLen) :: ErrMsg !< aggregated error message - integer(IntKi) :: ErrStat2 !< temporary error status from a call - character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call - character(*), parameter :: RoutineName = 'HydroDyn_CalcOutput_c' !< for error handling + real(DbKi) :: Time + integer(IntKi) :: iNode + integer(IntKi) :: ErrStat !< aggregated error status + character(ErrMsgLen) :: ErrMsg !< aggregated error message + integer(IntKi) :: ErrStat2 !< temporary error status from a call + character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call + character(*), parameter :: RoutineName = 'HydroDyn_CalcOutput_c' !< for error handling ! Initialize error handling ErrStat = ErrID_None ErrMsg = "" + ! Sanity check -- number of node points cannot change + if ( NumNodePts /= int(NumNodePts_C, IntKi) ) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "Number of node points passed in changed. This must be constant throughout simulation" + if (Failed()) return + endif + ! Convert the inputs from C to Fortrn Time = REAL(Time_C,DbKi) -!FIXME: Reshape position, velocity, acceleration and set mesh accordingly. + ! Reshape position, velocity, acceleration + tmpNodePos(1:6,1:NumNodePts) = reshape( real(NodePos_C(1:6*NumNodePts),ReKi), (/6,NumNodePts/) ) + tmpNodeVel(1:6,1:NumNodePts) = reshape( real(NodeVel_C(1:6*NumNodePts),ReKi), (/6,NumNodePts/) ) + tmpNodeAcc(1:6,1:NumNodePts) = reshape( real(NodeAcc_C(1:6*NumNodePts),ReKi), (/6,NumNodePts/) ) + + ! Transfer motions to input meshes + call Set_MotionMesh() ! update motion mesh with input motion arrays + call HD_SetInputMotion( 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, xd, z, OtherStates, y, m, ErrStat2, ErrMsg2 ) if (Failed()) return + + ! Transfer resulting load meshes to intermediate mesh + call HD_TransferLoads( u(1), y, ErrStat2, ErrMsg2 ) + if (Failed()) return + + + ! Set output force/moment array + call Set_OutputLoadArray( ) + ! Reshape for return + NodeFrc_C(1:6*NumNodePts) = reshape( real(tmpNodeFrc(1:6,1:NumNodePts), c_float), (/6*NumNodePts/) ) + ! Get the output channel info out of y OutputChannelValues_C = REAL(y%WriteOutput, C_FLOAT) + ! Set error status call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) CONTAINS @@ -658,6 +723,13 @@ SUBROUTINE HydroDyn_End_C(ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_End_c') ErrStat = ErrID_None ErrMsg = "" + ! clear out any globably allocated helper arrays + if (allocated(tmpNodePos)) deallocate(tmpNodePos) + if (allocated(tmpNodeVel)) deallocate(tmpNodeVel) + if (allocated(tmpNodeAcc)) deallocate(tmpNodeAcc) + if (allocated(tmpNodeFrc)) deallocate(tmpNodeFrc) + + ! NOTE: HydroDyn_End only takes 1 instance of u, not the array. So extra ! 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, @@ -707,4 +779,90 @@ subroutine ClearMesh() end subroutine ClearMesh END SUBROUTINE HydroDyn_End_c + +!> This routine is operating on module level data, hence few inputs +subroutine Set_MotionMesh() + integer(IntKi) :: iNode + real(R8Ki) :: theta(3) + real(R8Ki) :: Orient(3,3) + ! Set mesh corresponding to input motions + do iNode=1,NumNodePts + theta = real(tmpNodePos(4:6,iNode),DbKi) ! convert ReKi to DbKi to avoid roundoff + Orient = EulerConstruct( theta ) + HD_MotionMesh%TranslationDisp(1:3,iNode) = tmpNodePos(1:3,iNode) - HD_MotionMesh%Position(1:3,iNode) ! relative displacement only + HD_MotionMesh%Orientation(1:3,1:3,iNode) = Orient + HD_MotionMesh%TranslationVel( 1:3,iNode) = tmpNodeVel(1:3,iNode) + HD_MotionMesh%RotationVel( 1:3,iNode) = tmpNodeVel(4:6,iNode) + HD_MotionMesh%TranslationAcc( 1:3,iNode) = tmpNodeAcc(1:3,iNode) + HD_MotionMesh%RotationAcc( 1:3,iNode) = tmpNodeAcc(4:6,iNode) + enddo +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) + 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 ) + 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 (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 (ErrStat3 >= AbortErrLev) return + endif +end subroutine HD_SetInputMotion + + +!> Map the loads of the output meshes to the intermediate output mesh. Since +!! we are mapping two meshes over to a single one, we use an intermediate +!! 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) + integer(IntKi), intent( out) :: ErrStat3 + character(ErrMsgLen), intent( out) :: ErrMsg3 + + HD_LoadMesh%Force = 0.0_ReKi + HD_LoadMesh%Moment = 0.0_ReKi + + ! WAMIT mesh + if ( 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 ) + 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 + 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 ) + 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 +end subroutine HD_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() + integer(IntKi) :: iNode + ! Set mesh corresponding to input motions + do iNode=1,NumNodePts + tmpNodeFrc(1:3,iNode) = HD_LoadMesh%Force (1:3,iNode) + tmpNodeFrc(4:6,iNode) = HD_LoadMesh%Moment(1:3,iNode) + enddo +end subroutine Set_OutputLoadArray + END MODULE HydroDyn_C diff --git a/reg_tests/r-test b/reg_tests/r-test index 8ea661f075..61ff3ab369 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 8ea661f07534ee9fa8a8b243b069af81aae8a243 +Subproject commit 61ff3ab3697ff9be357708195d3230492deadeea From 00c90838a667c5d664ba81f36c24719775112607 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Fri, 25 Jun 2021 09:53:45 -0600 Subject: [PATCH 09/22] hd_c_lib: add updatestates interface, refine state time handling --- .../hydrodyn/python-lib/hydrodyn_library.py | 61 +++- modules/hydrodyn/src/HydroDyn_C.f90 | 264 ++++++++++++++---- reg_tests/r-test | 2 +- 3 files changed, 263 insertions(+), 64 deletions(-) diff --git a/modules/hydrodyn/python-lib/hydrodyn_library.py b/modules/hydrodyn/python-lib/hydrodyn_library.py index 34219936f4..5443187c0e 100644 --- a/modules/hydrodyn/python-lib/hydrodyn_library.py +++ b/modules/hydrodyn/python-lib/hydrodyn_library.py @@ -130,6 +130,7 @@ def __init__(self, library_path): # Initial time related variables self.dt = 0.1 # typical default for HD + self.t_start = 0.0 # initial time self.tmax = 600.0 # typical default for HD waves FFT #FIXME: check tmax/total_time and note exactly what is different between them. self.total_time = 0.0 # may be longer than tmax @@ -173,6 +174,7 @@ def _initialize_routines(self): POINTER(c_int), # numNodePts -- number of points expecting motions/loads POINTER(c_float), # initNodePos -- initial node positions in flat array of 6*numNodePts POINTER(c_int), # InterpOrder + POINTER(c_double), # t_initial POINTER(c_double), # dt POINTER(c_double), # tmax POINTER(c_int), # number of channels @@ -196,6 +198,19 @@ def _initialize_routines(self): ] self.HydroDyn_CalcOutput_c.restype = c_int + self.HydroDyn_UpdateStates_c.argtypes = [ + POINTER(c_double), # Time_C + POINTER(c_double), # TimeNext_C + POINTER(c_int), # numNodePts -- number of points expecting motions/loads + POINTER(c_float), # nodePos -- node positions in flat array of 6*numNodePts + POINTER(c_float), # nodeVel -- node velocities in flat array of 6*numNodePts + POINTER(c_float), # nodeAcc -- node accelerations in flat array of 6*numNodePts + POINTER(c_float), # nodeFrc -- node forces/moments in flat array of 6*numNodePts + POINTER(c_int), # ErrStat_C + POINTER(c_char) # ErrMsg_C + ] + self.HydroDyn_CalcOutput_c.restype = c_int + self.HydroDyn_End_c.argtypes = [ POINTER(c_int), # ErrStat_C POINTER(c_char) # ErrMsg_C @@ -262,6 +277,7 @@ def hydrodyn_init(self, input_string_array): byref(c_int(self.numNodePts)), # IN: number of attachment points expected (where motions are transferred into HD) nodeInitLoc_flat_c, # IN: initNodePos -- initial node positions in flat array of 6*numNodePts byref(c_int(self.InterpOrder)), # IN: InterpOrder (1: linear, 2: quadratic) + byref(c_double(self.t_start)), # IN: time initial byref(c_double(self.dt)), # IN: time step (dt) byref(c_double(self.tmax)), # IN: tmax byref(self._numChannels_c), # OUT: number of channels @@ -310,7 +326,7 @@ def hydrodyn_calcOutput(self, time, nodePos, nodeVel, nodeAcc, nodeFrcMom, outpu # Run HydroDyn_CalcOutput_c self.HydroDyn_CalcOutput_c( - byref(c_double(time)), # IN: time at which to calculate velocities + byref(c_double(time)), # IN: time at which to calculate output forces byref(c_int(self.numNodePts)), # IN: number of attachment points expected (where motions are transferred into HD) nodePos_flat_c, # IN: positions - specified by user nodeVel_flat_c, # IN: velocities at desired positions @@ -338,6 +354,49 @@ def hydrodyn_calcOutput(self, time, nodePos, nodeVel, nodeAcc, nodeFrcMom, outpu for k in range(0,self.numChannels): outputChannelValues[k] = float(outputChannelValues_c[k]) + # hydrodyn_updateStates ------------------------------------------------------------------------------------------------------------ + def hydrodyn_updateStates(self, time, timeNext, nodePos, nodeVel, nodeAcc, nodeFrcMom): + + # Check input motion info + self.check_input_motions(nodePos,nodeVel,nodeAcc) + + # set flat arrays for inputs of motion + # Position -- [x2,y1,z1,Rx1,Ry1,Rz1, x2,y2,z2,Rx2,Ry2,Rz2 ...] + nodePos_flat = [pp for p in nodePos for pp in p] + nodePos_flat_c = (c_float * (6 * self.numNodePts))(0.0,) + for i, p in enumerate(nodePos_flat): + nodePos_flat_c[i] = c_float(p) + + # Velocity -- [Vx2,Vy1,Vz1,RVx1,RVy1,RVz1, Vx2,Vy2,Vz2,RVx2,RVy2,RVz2 ...] + nodeVel_flat = [pp for p in nodeVel for pp in p] + nodeVel_flat_c = (c_float * (6 * self.numNodePts))(0.0,) + for i, p in enumerate(nodeVel_flat): + nodeVel_flat_c[i] = c_float(p) + + # Acceleration -- [Ax1,Ay1,Az1,RAx1,RAy1,RAz1, Ax2,Ay2,Az2,RAx2,RAy2,RAz2 ...] + nodeAcc_flat = [pp for p in nodeAcc for pp in p] + nodeAcc_flat_c = (c_float * (6 * self.numNodePts))(0.0,) + for i, p in enumerate(nodeAcc_flat): + nodeAcc_flat_c[i] = c_float(p) + + # Resulting Forces/moments -- [Fx1,Fy1,Fz1,Mx1,My1,Mz1, Fx2,Fy2,Fz2,Mx2,My2,Mz2 ...] + nodeFrc_flat_c = (c_float * (6 * self.numNodePts))(0.0,) + + # Run HydroDyn_UpdateStates_c + self.HydroDyn_UpdateStates_c( + 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(c_int(self.numNodePts)), # IN: number of attachment points expected (where motions are transferred into HD) + nodePos_flat_c, # IN: positions - specified by user + nodeVel_flat_c, # IN: velocities at desired positions + nodeAcc_flat_c, # IN: accelerations at desired positions + nodeFrc_flat_c, # OUT: resulting forces/moments array + byref(self.error_status_c), # OUT: ErrStat_C + self.error_message_c # OUT: ErrMsg_C + ) + + self.check_error() + # hydrodyn_end ------------------------------------------------------------------------------------------------------------ def hydrodyn_end(self): if not self.ended: diff --git a/modules/hydrodyn/src/HydroDyn_C.f90 b/modules/hydrodyn/src/HydroDyn_C.f90 index 68d81a2d49..ff017c17b5 100644 --- a/modules/hydrodyn/src/HydroDyn_C.f90 +++ b/modules/hydrodyn/src/HydroDyn_C.f90 @@ -58,17 +58,21 @@ MODULE HydroDyn_C 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_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 !< continuous states at Time t - type(HydroDyn_DiscreteStateType) :: xd !< discrete states at Time t - type(HydroDyn_ConstraintStateType) :: z !< Constraint states at Time t - type(HydroDyn_OtherStateType) :: OtherStates !< Initial other/optimization states + 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) !------------------------------ + ! Time tracking + ! When we are performing a correction step, time information of previous + ! calls is needed to decide how to apply correction logic or cycle the inputs + ! and resave the previous timestep states. ! Correction steps ! OpenFAST has the ability to perform correction steps. During a correction ! step, new input values are passed in but the timestep remains the same. @@ -78,20 +82,19 @@ MODULE HydroDyn_C ! the glue code. However, here we do not pass state information through the ! interface and therefore must store it here analogously to how it is handled ! in the OpenFAST glue code. - type(HydroDyn_ContinuousStateType) :: x_prev !< continuous states at Time t of previous call - type(HydroDyn_DiscreteStateType) :: xd_prev !< discrete states at Time t of previous call - type(HydroDyn_ConstraintStateType) :: z_prev !< Constraint states at Time t of previous call - !------------------------------ - ! Time tracking - ! When we are performing a correction step, time information of previous - ! calls is needed to decide how to apply correction logic or cycle the inputs - ! and resave the previous timestep states. - integer(IntKi) :: dT_Global ! dT of the code calling this module - integer(IntKi) :: T_Global ! Time of this call - integer(IntKi) :: T_Global_prev ! time of the previous call - integer(IntKi) :: N_T_Global ! count of which timestep we are on -- calculated internally based on the time passed into UpdateStates - integer(IntKi) :: N_T_Global_prev ! timestep of the previous call - logical :: CorrectionStep ! if we are repeating a timestep in UpdateStates, + real(DbKi) :: dT_Global ! dT of the code calling this module + integer(IntKi) :: N_Global ! global timestep + real(DbKi) :: T_Initial ! initial Time of simulation + real(DbKi), allocatable :: InputTimes(:) ! input times corresponding to u(:) array + 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) + integer(IntKi), parameter :: STATE_CURR = 1 ! Index for current state + integer(IntKi), parameter :: STATE_PRED = 2 ! Index for predicted state + ! Note the indexing is different on inputs (no clue why, but thats how OF handles it) + integer(IntKi), parameter :: INPUT_LAST = 3 ! Index for previous input at t-dt + integer(IntKi), parameter :: INPUT_CURR = 2 ! Index for current input at t + integer(IntKi), parameter :: INPUT_PRED = 1 ! Index for predicted input at t+dt !------------------------------------------------------------------------------------ @@ -163,7 +166,7 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen PtfmRefPtPositionX_C, PtfmRefPtPositionY_C, & NumNodePts_C, InitNodePositions_C, & !NumWaveElev_C, WaveElevXY_C & !Placeholder for later - InterpOrder_C, DT_C, TMax_C, & + InterpOrder_C, T_initial_C, DT_C, TMax_C, & NumChannels_C, OutputChannelNames_C, OutputChannelUnits_C, & ErrStat_C, ErrMsg_C) BIND (C, NAME='HydroDyn_Init_c') implicit none @@ -186,6 +189,7 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen !NOTE: not setting up the WaveElev at this point. Leaving placeholder for future !integer(c_int), intent(in ) :: NumWaveElev_C !< Number of mesh points we are transfering motions to and output loads to !real(c_float), intent(in ) :: WaveElevXY_C !< A 2xNumWaveElev_C array [x,y] + real(c_double), intent(in ) :: T_initial_C integer(c_int), intent(in ) :: InterpOrder_C !< Interpolation order to use (must be 1 or 2) real(c_double), intent(in ) :: DT_C !< Timestep used with HD for stepping forward from t to t+dt. Must be constant. real(c_double), intent(in ) :: TMax_C !< Maximum time for simulation (used to set arrays for wave kinematics) @@ -252,6 +256,9 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen InitInp%defWtrDpth = REAL(defWtrDpth_C, ReKi) InitInp%defMSL2SWL = REAL(defMSL2SWL_C, ReKi) 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 @@ -285,21 +292,22 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen !InitInp%WaveElevXY - !----------------------- - ! Allocate input array u - !----------------------- + !---------------------------------------------------- + ! Allocate input array u and corresponding InputTimes + !---------------------------------------------------- ! 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 + ! u(3) inputs at t - 2*dt ! quadratic only allocate(u(InterpOrder+1), STAT=ErrStat2) if (ErrStat2 /= 0) then ErrStat2 = ErrID_Fatal ErrMsg2 = "Could not allocate inuput" if (Failed()) return endif + call AllocAry( InputTimes, InterpOrder+1, "InputTimes", ErrStat2, ErrMsg2 ); if (Failed()) return ! Call the main subroutine HydroDyn_Init @@ -307,7 +315,7 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen ! ! 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, xd, z, OtherStates, y, m, TimeInterval, InitOutData, ErrStat2, ErrMsg2 ) + 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 ) if (Failed()) return @@ -325,12 +333,37 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen !------------------------------------------------------------- - ! Setup other pieces of u(:) + ! Setup other prior timesteps + ! We fill InputTimes with negative times, but the Input values are identical for each of those times; this allows + ! us to use, e.g., quadratic interpolation that effectively acts as a zeroth-order extrapolation and first-order extrapolation + ! for the first and second time steps. (The interpolation order in the ExtrapInput routines are determined as + ! order = SIZE(Input) !------------------------------------------------------------- do i=2,InterpOrder+1 call HydroDyn_CopyInput (u(1), u(i), MESH_NEWCOPY, Errstat2, ErrMsg2) if (Failed()) return enddo + do i = 1, InterpOrder + 1 + InputTimes(i) = t_initial - (i - 1) * dT_Global + enddo + InputTimePrev = InputTimes(1) ! Set this as the last call time for UpdateStates + + + !------------------------------------------------------------- + ! 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 + + !------------------------------------------------------------- + ! 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 !TODO ! Is there any other InitOutData should be returned @@ -616,7 +649,7 @@ SUBROUTINE HydroDyn_CalcOutput_c(Time_C, NumNodePts_C, NodePos_C, NodeVel_C, Nod ! Call the main subroutine HydroDyn_CalcOutput to get the resulting forces and moments at time T - CALL HydroDyn_CalcOutput( Time, u(1), p, x, xd, z, OtherStates, y, m, ErrStat2, ErrMsg2 ) + CALL HydroDyn_CalcOutput( Time, u(1), p, x(STATE_CURR), xd(STATE_CURR), z(STATE_CURR), OtherStates(STATE_CURR), y, m, ErrStat2, ErrMsg2 ) if (Failed()) return @@ -645,47 +678,138 @@ end function Failed END SUBROUTINE HydroDyn_CalcOutput_c !=============================================================================================================== -!--------------------------------------------- HydroDyn UpdateStates --------------------------------------------- +!--------------------------------------------- HydroDyn UpdateStates ------------------------------------------- !=============================================================================================================== -SUBROUTINE HydroDyn_UpdateStates_c(Time_C,ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_UpdateStates_c') +!> This routine updates the states from Time_C to TimeNext_C. It is assumed that the inputs are given for +!! TimeNext_C, but will be checked against the previous timestep values. +!! Since we don't really know if we are doing correction steps or not, we will track the previous state and +!! reset to those if we are repeating a timestep (normally this would be handled by the OF glue code, but since +!! the states are not passed across the interface, we must handle them here). +SUBROUTINE HydroDyn_UpdateStates_c( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, NodeVel_C, NodeAcc_C, & + ErrStat_C, ErrMsg_C) BIND (C, NAME='HydroDyn_UpdateStates_c') implicit none #ifndef IMPLICIT_DLLEXPORT !DEC$ ATTRIBUTES DLLEXPORT :: HydroDyn_UpdateStates_c !GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_UpdateStates_c #endif - real(c_double), intent(in ) :: Time_C - integer(c_int), intent( out) :: ErrStat_C - character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + real(c_double), intent(in ) :: Time_C + real(c_double), intent(in ) :: TimeNext_C +!FIXME: estimated inputs at t+dt, or inputs at t???? + integer(c_int), intent(in ) :: NumNodePts_C !< Number of mesh points we are transfering motions to and output loads to + real(c_float), intent(in ) :: NodePos_C( 6*NumNodePts_C ) !< A 6xNumNodePts_C array [x,y,z,Rx,Ry,Rz] -- positions (global) + 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) + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) ! Local variables - real(DbKi) :: Time - integer :: ErrStat !< aggregated error status - 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 = 'HydroDyn_UpdateStates_c' !< for error handling + logical :: CorrectionStep ! if we are repeating a timestep in UpdateStates, don't update the inputs array + integer(IntKi) :: iNode + integer(IntKi) :: ErrStat !< aggregated error status + character(ErrMsgLen) :: ErrMsg !< aggregated error message + integer(IntKi) :: ErrStat2 !< temporary error status from a call + character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call + character(*), parameter :: RoutineName = 'HydroDyn_UpdateStates_c' !< for error handling ! Initialize error handling ErrStat = ErrID_None ErrMsg = "" + CorrectionStep = .false. -!FIXME: set time as array - ! Convert the inputs from C to Fortran - Time = REAL(Time_C,DbKi) + ! Sanity check -- number of node points cannot change + if ( NumNodePts /= int(NumNodePts_C, IntKi) ) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "Number of node points passed in changed. This must be constant throughout simulation" + if (Failed()) return + endif -!FIXME: do I need to have a flag for stepping forward in time? What if we are doing a correction step???? -!Rotate values for inputs (do extrap interp if we are not passed new inputs here) -!FIXME: Reshape position, velocity, acceleration and set mesh inputs + !------------------------------------------------------- + ! Check the time for current timestep and next timestep + !------------------------------------------------------- + ! 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 + ! u(2) inputs at t ! This timestep + ! u(3) inputs at t - dt ! previous timestep (quadratic only) + + ! Check if we are repeating an UpdateStates call (for example in a predictor/corrector loop) + if ( EqualRealNos( real(Time_C,DbKi), InputTimePrev ) ) then + CorrectionStep = .true. + N_Global = N_Global + 1_IntKi + else ! Setup time input times array + if (InterpOrder>1) then ! quadratic, so keep the old time + InputTimes(INPUT_LAST) = InputTimes(INPUT_CURR) ! u(3) + endif + InputTimes(INPUT_CURR) = REAL(Time_C,DbKi) ! u(2) + InputTimes(INPUT_PRED) = REAL(TimeNext_C,DbKi) ! u(1) + endif -!FIXME:Set inputs array and extrap/interp necessary -! ! Call the main subroutine HydroDyn_UpdateStates to get the velocities -! CALL HydroDyn_UpdateStates( Time, u, p, x, xd, z, OtherStates, y, m, ErrStat2, ErrMsg2 ) -! if (Failed()) return + if (CorrectionStep) then + ! 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 + 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 + 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 + endif + + !------------------------------------------------------- + ! Set inputs for time T+dt -- u(1) + !------------------------------------------------------- + ! Reshape position, velocity, acceleration + tmpNodePos(1:6,1:NumNodePts) = reshape( real(NodePos_C(1:6*NumNodePts),ReKi), (/6,NumNodePts/) ) + tmpNodeVel(1:6,1:NumNodePts) = reshape( real(NodeVel_C(1:6*NumNodePts),ReKi), (/6,NumNodePts/) ) + tmpNodeAcc(1:6,1:NumNodePts) = reshape( real(NodeAcc_C(1:6*NumNodePts),ReKi), (/6,NumNodePts/) ) + + ! Transfer motions to input meshes + call Set_MotionMesh() ! update motion mesh with input motion arrays + call HD_SetInputMotion( 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 the main subroutine HydroDyn_UpdateStates to get the velocities + CALL HydroDyn_UpdateStates( InputTimes(INPUT_PRED), N_Global, u, InputTimes, p, x(STATE_PRED), xd(STATE_PRED), z(STATE_PRED), OtherStates(STATE_PRED), m, ErrStat2, ErrMsg2 ) + if (Failed()) return + + + !------------------------------------------------------- + ! cycle the states + !------------------------------------------------------- + ! 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 + ! 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 + -!FIXME: what do we need to update after the call? -! states are handled internally. If we are doing correction steps, we may need the old states call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) @@ -730,28 +854,44 @@ SUBROUTINE HydroDyn_End_C(ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_End_c') if (allocated(tmpNodeFrc)) deallocate(tmpNodeFrc) + ! Call the main subroutine HydroDyn_End + ! If u is not allocated, then we didn't get far at all in initialization, + ! or HD_End_C 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 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + endif + ! NOTE: HydroDyn_End only takes 1 instance of u, not the array. So extra ! 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) ! leave first one for passing to HD_End for destruction + do i=2,size(u) call HydroDyn_DestroyInput( u(i), ErrStat2, ErrMsg2 ) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) enddo + if (allocated(u)) deallocate(u) endif - ! Call the main subroutine HydroDyn_End - ! If u is not allocated, then we didn't get far at all in initialization, - ! or HD_End_C 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, xd, z, OtherStates, y, m, ErrStat2, ErrMsg2 ) - call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - endif - - ! if u is still allocated, deallocate it now - if (allocated(u)) deallocate(u) + ! 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 ) + + + ! if deallocate other items now + if (allocated(InputTimes)) deallocate(InputTimes) ! Clear out mesh related data storage call ClearMesh() diff --git a/reg_tests/r-test b/reg_tests/r-test index 3dcd5b07f5..1c89781e91 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 3dcd5b07f5a86ef4d2659dc57bafa3d58b44a96f +Subproject commit 1c89781e9147ea01416a813f84674c3df7c95ddb From ba5ec1b3f2908e2a5d8951f1709aae1e2df233df Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Tue, 27 Jul 2021 11:08:51 -0600 Subject: [PATCH 10/22] HD_c_lib: add VS solution --- .../HydroDyn_c_lib/HydroDynDriver_c_lib.sln | 61 ++++ .../HydroDynDriver_c_lib.vfproj | 319 ++++++++++++++++++ 2 files changed, 380 insertions(+) create mode 100644 vs-build/HydroDyn_c_lib/HydroDynDriver_c_lib.sln create mode 100644 vs-build/HydroDyn_c_lib/HydroDynDriver_c_lib.vfproj diff --git a/vs-build/HydroDyn_c_lib/HydroDynDriver_c_lib.sln b/vs-build/HydroDyn_c_lib/HydroDynDriver_c_lib.sln new file mode 100644 index 0000000000..2c9e9c94f2 --- /dev/null +++ b/vs-build/HydroDyn_c_lib/HydroDynDriver_c_lib.sln @@ -0,0 +1,61 @@ + +Microsoft Visual Studio Solution File, Format Version 12.00 +# Visual Studio Version 16 +VisualStudioVersion = 16.0.30503.244 +MinimumVisualStudioVersion = 10.0.40219.1 +Project("{A2215CCC-531E-454C-A1F0-E502FB697697}") = "HydroDynDriver_c_lib.dll", "HydroDynDriver_c_lib.vfproj", "{FDA4A02B-B3A7-4D06-847C-941BE44E76FB}" + ProjectSection(ProjectDependencies) = postProject + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16} = {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16} + EndProjectSection +EndProject +Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "FAST_Registry", "..\Registry\FAST_Registry.vcxproj", "{DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}" +EndProject +Global + GlobalSection(SolutionConfigurationPlatforms) = preSolution + Debug_Double|Win32 = Debug_Double|Win32 + Debug_Double|x64 = Debug_Double|x64 + Debug|Win32 = Debug|Win32 + Debug|x64 = Debug|x64 + Release_Double|Win32 = Release_Double|Win32 + Release_Double|x64 = Release_Double|x64 + Release|Win32 = Release|Win32 + Release|x64 = Release|x64 + EndGlobalSection + GlobalSection(ProjectConfigurationPlatforms) = postSolution + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Debug_Double|Win32.ActiveCfg = Debug_Double|Win32 + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Debug_Double|Win32.Build.0 = Debug_Double|Win32 + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Debug_Double|x64.ActiveCfg = Debug_Double|x64 + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Debug_Double|x64.Build.0 = Debug_Double|x64 + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Debug|Win32.ActiveCfg = Debug|Win32 + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Debug|Win32.Build.0 = Debug|Win32 + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Debug|x64.ActiveCfg = Debug|x64 + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Debug|x64.Build.0 = Debug|x64 + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Release_Double|Win32.ActiveCfg = Release_Double|Win32 + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Release_Double|Win32.Build.0 = Release_Double|Win32 + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Release_Double|x64.ActiveCfg = Release_Double|x64 + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Release_Double|x64.Build.0 = Release_Double|x64 + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Release|Win32.ActiveCfg = Release|Win32 + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Release|Win32.Build.0 = Release|Win32 + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Release|x64.ActiveCfg = Release|x64 + {FDA4A02B-B3A7-4D06-847C-941BE44E76FB}.Release|x64.Build.0 = Release|x64 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Debug_Double|Win32.ActiveCfg = Release|Win32 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Debug_Double|Win32.Build.0 = Release|Win32 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Debug_Double|x64.ActiveCfg = Release|Win32 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Debug_Double|x64.Build.0 = Release|Win32 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Debug|Win32.ActiveCfg = Release|Win32 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Debug|Win32.Build.0 = Release|Win32 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Debug|x64.ActiveCfg = Release|Win32 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Debug|x64.Build.0 = Release|Win32 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Release_Double|Win32.ActiveCfg = Release|Win32 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Release_Double|Win32.Build.0 = Release|Win32 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Release_Double|x64.ActiveCfg = Release|Win32 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Release_Double|x64.Build.0 = Release|Win32 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Release|Win32.ActiveCfg = Release|Win32 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Release|Win32.Build.0 = Release|Win32 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Release|x64.ActiveCfg = Release|Win32 + {DA16A3A6-3297-4628-9E46-C6FA0E3C4D16}.Release|x64.Build.0 = Release|Win32 + EndGlobalSection + GlobalSection(SolutionProperties) = preSolution + HideSolutionNode = FALSE + EndGlobalSection +EndGlobal diff --git a/vs-build/HydroDyn_c_lib/HydroDynDriver_c_lib.vfproj b/vs-build/HydroDyn_c_lib/HydroDynDriver_c_lib.vfproj new file mode 100644 index 0000000000..ae2211356e --- /dev/null +++ b/vs-build/HydroDyn_c_lib/HydroDynDriver_c_lib.vfproj @@ -0,0 +1,319 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 2c6ca6dab49260eba90ecd9b6b1692fcce7ad78a Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Thu, 29 Jul 2021 18:51:50 -0600 Subject: [PATCH 11/22] HD_c_lib: remove extra arguments from UpdateStates_C interface Extra arguments existed on the python side that did not match arguments on fortran side. Nasty segfault to debug --- modules/hydrodyn/python-lib/hydrodyn_library.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/modules/hydrodyn/python-lib/hydrodyn_library.py b/modules/hydrodyn/python-lib/hydrodyn_library.py index 5443187c0e..9d676262e8 100644 --- a/modules/hydrodyn/python-lib/hydrodyn_library.py +++ b/modules/hydrodyn/python-lib/hydrodyn_library.py @@ -122,7 +122,7 @@ def __init__(self, library_path): # Initial environmental conditions self.gravity = 9.80665 # Gravity (m/s^2) self.defWtrDens = 1025.0 # Water density (kg/m^3) - self.defWtrDpth = 120.0 # Water depth (m) + self.defWtrDpth = 200.0 # Water depth (m) self.defMSL2SWL = 0.0 # Offset between still-water level and mean sea level (m) [positive upward] # Interpolation order (must be 1: linear, or 2: quadratic) @@ -205,7 +205,6 @@ def _initialize_routines(self): POINTER(c_float), # nodePos -- node positions in flat array of 6*numNodePts POINTER(c_float), # nodeVel -- node velocities in flat array of 6*numNodePts POINTER(c_float), # nodeAcc -- node accelerations in flat array of 6*numNodePts - POINTER(c_float), # nodeFrc -- node forces/moments in flat array of 6*numNodePts POINTER(c_int), # ErrStat_C POINTER(c_char) # ErrMsg_C ] @@ -390,7 +389,6 @@ def hydrodyn_updateStates(self, time, timeNext, nodePos, nodeVel, nodeAcc, nodeF nodePos_flat_c, # IN: positions - specified by user nodeVel_flat_c, # IN: velocities at desired positions nodeAcc_flat_c, # IN: accelerations at desired positions - nodeFrc_flat_c, # OUT: resulting forces/moments array byref(self.error_status_c), # OUT: ErrStat_C self.error_message_c # OUT: ErrMsg_C ) @@ -414,9 +412,9 @@ def check_error(self): if self.error_status_c.value == 0: return elif self.error_status_c.value < self.abort_error_level: - print(f"{self.error_levels[self.error_status_c.value]}: {self.error_message_c.value.decode('ascii')}") + print(f"HydroDyn error status: {self.error_levels[self.error_status_c.value]}: {self.error_message_c.value.decode('ascii')}") else: - print(f"{self.error_levels[self.error_status_c.value]}: {self.error_message_c.value.decode('ascii')}") + print(f"HydroDyn error status: {self.error_levels[self.error_status_c.value]}: {self.error_message_c.value.decode('ascii')}") self.hydrodyn_end() raise Exception("\nHydroDyn terminated prematurely.") @@ -615,8 +613,8 @@ def write(self,chan_data): f_string = "{:10.4f}"+"{:25.7f}"*(l-1) for i in range(0,chan_data.shape[0]): self.OutFile.write(f_string.format(*chan_data[i,:]) + '\n') - if i==0: - print(f"{chan_data[i,:]}") + #if i==0: + # print(f"{chan_data[i,:]}") def end(self): if self.opened: From 57d7e873e6ebc28369b7cdc954ae3fdccd1faa4c Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Thu, 29 Jul 2021 18:53:15 -0600 Subject: [PATCH 12/22] HD_c_lib: update timestep handling UpdateStates_C and fix angle handling - Angles had been handled using Euler angle orders, but that did not match the way HD_driver handled them, or how the output angles are set. This could still be a precision issue for coupling to other codes later. - Double precision handling of time: - time is passed in as double precision. Double precision compiles promote time internally to quad precision resulting in a mismatch. The radiation module checks the dT between input times and will throw an error due to the precision mismatch. So converted to using the global timestep value in calculating InputTimes to pass. --- modules/hydrodyn/src/HydroDyn_C.f90 | 96 +++++++++++++++++------------ 1 file changed, 57 insertions(+), 39 deletions(-) diff --git a/modules/hydrodyn/src/HydroDyn_C.f90 b/modules/hydrodyn/src/HydroDyn_C.f90 index ff017c17b5..e7bbac27fe 100644 --- a/modules/hydrodyn/src/HydroDyn_C.f90 +++ b/modules/hydrodyn/src/HydroDyn_C.f90 @@ -1,6 +1,6 @@ !********************************************************************************************************************************** ! LICENSING -! Copyright (C) 2021 National Renewable Energy Lab +! Copyright (C) 2021 National Renewable Energy Lab ! ! This file is part of HydroDyn. ! @@ -82,7 +82,7 @@ MODULE HydroDyn_C ! the glue code. However, here we do not pass state information through the ! interface and therefore must store it here analogously to how it is handled ! in the OpenFAST glue code. - real(DbKi) :: dT_Global ! dT of the code calling this module + real(DbKi) :: dT_Global ! dT of the code calling this module integer(IntKi) :: N_Global ! global timestep real(DbKi) :: T_Initial ! initial Time of simulation real(DbKi), allocatable :: InputTimes(:) ! input times corresponding to u(:) array @@ -115,11 +115,11 @@ MODULE HydroDyn_C ! - N points -- flexible structure (either floating or fixed bottom) integer(IntKi) :: NumNodePts ! Number of mesh points we are interfacing motions/loads to/from HD type(MeshType) :: HD_MotionMesh ! mesh for motions of external nodes - type(MeshType) :: HD_LoadMesh ! mesh for loads for external nodes + type(MeshType) :: HD_LoadMesh ! mesh for loads for external nodes type(MeshType) :: HD_LoadMesh_tmp ! mesh for loads for external nodes -- temporary !------------------------------ ! Mesh mapping: motions - ! The mapping of motions from the nodes passed in to the corresponding HD meshes + ! The mapping of motions from the nodes passed in to the corresponding HD meshes type(MeshMapType) :: Map_Motion_2_HD_PRP_P ! Mesh mapping between input motion mesh and PRP type(MeshMapType) :: Map_Motion_2_HD_WB_P ! Mesh mapping between input motion mesh and WAMIT body(ies) mesh type(MeshMapType) :: Map_Motion_2_HD_Mo_P ! Mesh mapping between input motion mesh and Morison mesh @@ -136,7 +136,7 @@ MODULE HydroDyn_C real(ReKi), allocatable :: tmpNodeAcc(:,:) ! temp array. Probably don't need this, but makes conversion from C clearer. real(ReKi), allocatable :: tmpNodeFrc(:,:) ! temp array. Probably don't need this, but makes conversion to C clearer. !------------------------------------------------------------------------------------ - + CONTAINS @@ -149,6 +149,7 @@ subroutine SetErr(ErrStat, ErrMsg, ErrStat_C, ErrMsg_C) character(ErrMsgLen), intent(in ) :: ErrMsg !< aggregated error message (fortran type) integer(c_int), intent( out) :: ErrStat_C character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + integer :: i ErrStat_C = ErrStat ! We will send back the same error status that is used in OpenFAST if (ErrMsgLen > ErrMsgLen_C-1) then ! If ErrMsgLen is > the space in ErrMsg_C, do not copy everything over ErrMsg_C = TRANSFER( trim(ErrMsg(1:ErrMsgLen_C-1))//C_NULL_CHAR, ErrMsg_C ) @@ -182,8 +183,8 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen 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) real(c_float), intent(in ) :: defMSL2SWL_C !< Default Offset between still-water level and mean sea level (m) [positive upward] (may be overridden by input file) - real(c_float), intent(in ) :: PtfmRefPtPositionX_C !< Initial position in wave field - real(c_float), intent(in ) :: PtfmRefPtPositionY_C !< Initial position in wave field + real(c_float), intent(in ) :: PtfmRefPtPositionX_C !< Initial position in wave field + real(c_float), intent(in ) :: PtfmRefPtPositionY_C !< Initial position in wave field integer(c_int), intent(in ) :: NumNodePts_C !< Number of mesh points we are transfering motions to and output loads to real(c_float), intent(in ) :: InitNodePositions_C( 6*NumNodePts_C ) !< A 6xNumNodePts_C array [x,y,z,Rx,Ry,Rz] !NOTE: not setting up the WaveElev at this point. Leaving placeholder for future @@ -203,7 +204,7 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen 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 + 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 @@ -223,7 +224,7 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen if (Failed()) return endif - ! Get fortran pointer to C_NULL_CHAR deliniated input file as a string + ! 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 @@ -257,7 +258,7 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen InitInp%defMSL2SWL = REAL(defMSL2SWL_C, ReKi) 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 + 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) @@ -346,7 +347,7 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen do i = 1, InterpOrder + 1 InputTimes(i) = t_initial - (i - 1) * dT_Global enddo - InputTimePrev = InputTimes(1) ! Set this as the last call time for UpdateStates + InputTimePrev = InputTimes(1) - dT_Global ! Initialize for UpdateStates !------------------------------------------------------------- @@ -367,7 +368,7 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen !TODO ! Is there any other InitOutData should be returned -! Any additional warnings or error handling necessary +! Any additional warnings or error handling necessary !------------------------------------------------- @@ -441,14 +442,14 @@ subroutine SetMotionLoadsInterfaceMeshes(ErrStat3,ErrMsg3) ! initial position and orientation of node InitPos = tmpNodePos(1:3,iNode) theta = real(tmpNodePos(4:6,iNode),DbKi) ! convert ReKi to DbKi to avoid roundoff - Orient = EulerConstruct( theta ) + CALL SmllRotTrans( 'InputRotation', theta(1), theta(2), theta(3), Orient, 'Orient', ErrStat, ErrMsg ) call MeshPositionNode( HD_MotionMesh , & iNode , & InitPos , & ! position ErrStat3, ErrMsg3 , & Orient ) ! orientation if (ErrStat3 >= AbortErrLev) return - + call MeshConstructElement ( HD_MotionMesh, ELEMENT_POINT, ErrStat3, ErrMsg3, iNode ) if (ErrStat3 >= AbortErrLev) return enddo @@ -461,7 +462,7 @@ subroutine SetMotionLoadsInterfaceMeshes(ErrStat3,ErrMsg3) ! For checking the mesh, uncomment this. ! note: CU is is output unit (platform dependent). !call MeshPrintInfo( CU, HD_MotionMesh ) - + !------------------------------------------------------------- ! Loads mesh ! This point mesh may contain more than one point. Mapping will be used to map @@ -476,9 +477,9 @@ subroutine SetMotionLoadsInterfaceMeshes(ErrStat3,ErrMsg3) Force = .TRUE. ,& Moment = .TRUE. ) if (ErrStat3 >= AbortErrLev) return - + HD_LoadMesh%RemapFlag = .TRUE. - + ! For checking the mesh, uncomment this. ! note: CU is is output unit (platform dependent). !call MeshPrintInfo( CU, HD_LoadMesh ) @@ -497,9 +498,9 @@ subroutine SetMotionLoadsInterfaceMeshes(ErrStat3,ErrMsg3) Force = .TRUE. ,& Moment = .TRUE. ) if (ErrStat3 >= AbortErrLev) return - + HD_LoadMesh_tmp%RemapFlag = .TRUE. - + ! For checking the mesh, uncomment this. ! note: CU is is output unit (platform dependent). !call MeshPrintInfo( CU, HD_LoadMesh_tmp ) @@ -534,7 +535,7 @@ end subroutine SetMotionLoadsInterfaceMeshes !> Sanity check the nodes !! If more than one input node was passed in, but only a single HD node !! exits (single Morison or single WAMIT), then give error that too many - !! nodes passed. + !! nodes passed. subroutine CheckNodes(ErrStat3,ErrMsg3) integer(IntKi), intent( out) :: ErrStat3 !< temporary error status character(ErrMsgLen), intent( out) :: ErrMsg3 !< temporary error message @@ -551,7 +552,7 @@ subroutine CheckNodes(ErrStat3,ErrMsg3) 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 + elseif ( u(1)%WAMITMesh%Committed ) then ! No Morison if ( u(1)%WAMITMesh%Nnodes < NumNodePts ) then ErrStat3 = ErrID_Fatal ErrMsg3 = "More nodes passed into library than exist in HydroDyn model WAMIT mesh" @@ -616,7 +617,7 @@ SUBROUTINE HydroDyn_CalcOutput_c(Time_C, NumNodePts_C, NodePos_C, NodeVel_C, Nod ! Local variables real(DbKi) :: Time integer(IntKi) :: iNode - integer(IntKi) :: ErrStat !< aggregated error status + integer(IntKi) :: ErrStat !< aggregated error status character(ErrMsgLen) :: ErrMsg !< aggregated error message integer(IntKi) :: ErrStat2 !< temporary error status from a call character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call @@ -641,14 +642,17 @@ SUBROUTINE HydroDyn_CalcOutput_c(Time_C, NumNodePts_C, NodePos_C, NodeVel_C, Nod tmpNodeVel(1:6,1:NumNodePts) = reshape( real(NodeVel_C(1:6*NumNodePts),ReKi), (/6,NumNodePts/) ) tmpNodeAcc(1:6,1:NumNodePts) = reshape( real(NodeAcc_C(1:6*NumNodePts),ReKi), (/6,NumNodePts/) ) +write(20,*) Time,NodePos_C,NodeVel_C,NodeAcc_C +write(21,*) Time,tmpNodePos,tmpNodeVel,tmpNodeAcc ! Transfer motions to input meshes - call Set_MotionMesh() ! update motion mesh with input motion arrays + 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 if (Failed()) return - - ! Call the main subroutine HydroDyn_CalcOutput to get the resulting forces and moments at time T + + ! 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 ) if (Failed()) return @@ -705,7 +709,7 @@ SUBROUTINE HydroDyn_UpdateStates_c( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, ! Local variables logical :: CorrectionStep ! if we are repeating a timestep in UpdateStates, don't update the inputs array integer(IntKi) :: iNode - integer(IntKi) :: ErrStat !< aggregated error status + integer(IntKi) :: ErrStat !< aggregated error status character(ErrMsgLen) :: ErrMsg !< aggregated error message integer(IntKi) :: ErrStat2 !< temporary error status from a call character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call @@ -725,7 +729,7 @@ SUBROUTINE HydroDyn_UpdateStates_c( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, !------------------------------------------------------- - ! Check the time for current timestep and next timestep + ! Check the time for current timestep and next timestep !------------------------------------------------------- ! These inputs are used in the time stepping algorithm within HD_UpdateStates ! For quadratic interpolation (InterpOrder==2), 3 timesteps are used. For @@ -733,17 +737,28 @@ SUBROUTINE HydroDyn_UpdateStates_c( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, ! u(1) inputs at t + dt ! Next timestep ! u(2) inputs at t ! This timestep ! u(3) inputs at t - dt ! previous timestep (quadratic only) + ! + ! NOTE: Within HD, the Radiation calculations can be done at an integer multiple of the + ! timestep. This is checked at each UpdateStates call. However, if we compile + ! in double precision, the values of Time_C and TimeNext_C are in double precison, + ! but InputTimes is in DbKi (which is promoted quad precision when compiling in + ! double precision) and the check may fail. So we are going to set the times we + ! we pass over to UpdateStates using the global timestep and the stored DbKi value + ! for the timestep rather than the lower precision (when compiled double) time + ! values passed in. It is a bit of a clumsy workaround for this precision loss, + ! but should not affect any results. ! Check if we are repeating an UpdateStates call (for example in a predictor/corrector loop) if ( EqualRealNos( real(Time_C,DbKi), InputTimePrev ) ) then CorrectionStep = .true. - N_Global = N_Global + 1_IntKi else ! Setup time input times array + InputTimePrev = real(Time_C,DbKi) ! Store for check next time if (InterpOrder>1) then ! quadratic, so keep the old time - InputTimes(INPUT_LAST) = InputTimes(INPUT_CURR) ! u(3) + InputTimes(INPUT_LAST) = ( N_Global - 1 ) * dT_Global ! u(3) at T-dT endif - InputTimes(INPUT_CURR) = REAL(Time_C,DbKi) ! u(2) - InputTimes(INPUT_PRED) = REAL(TimeNext_C,DbKi) ! u(1) + InputTimes(INPUT_CURR) = N_Global * dT_Global ! u(2) at T + InputTimes(INPUT_PRED) = ( N_Global + 1 ) * dT_Global ! u(1) at T+dT + N_Global = N_Global + 1_IntKi ! increment counter to T+dT endif @@ -758,7 +773,7 @@ SUBROUTINE HydroDyn_UpdateStates_c( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, 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( u(INPUT_CURR), 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 @@ -773,10 +788,11 @@ SUBROUTINE HydroDyn_UpdateStates_c( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, tmpNodeAcc(1:6,1:NumNodePts) = reshape( real(NodeAcc_C(1:6*NumNodePts),ReKi), (/6,NumNodePts/) ) ! Transfer motions to input meshes - call Set_MotionMesh() ! update motion mesh with input motion arrays + 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 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. @@ -785,7 +801,7 @@ SUBROUTINE HydroDyn_UpdateStates_c( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, 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 the main subroutine HydroDyn_UpdateStates to get the velocities CALL HydroDyn_UpdateStates( InputTimes(INPUT_PRED), N_Global, u, InputTimes, p, x(STATE_PRED), xd(STATE_PRED), z(STATE_PRED), OtherStates(STATE_PRED), m, ErrStat2, ErrMsg2 ) @@ -808,7 +824,7 @@ SUBROUTINE HydroDyn_UpdateStates_c( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, 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 SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) @@ -837,7 +853,7 @@ SUBROUTINE HydroDyn_End_C(ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_End_c') ! Local variables integer(IntKi) :: i !< generic loop counter - integer :: ErrStat !< aggregated error status + integer :: ErrStat !< aggregated error status character(ErrMsgLen) :: ErrMsg !< aggregated error message integer :: ErrStat2 !< temporary error status from a call character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call @@ -921,14 +937,16 @@ END SUBROUTINE HydroDyn_End_c !> This routine is operating on module level data, hence few inputs -subroutine Set_MotionMesh() +subroutine Set_MotionMesh(ErrStat3, ErrMsg3) + integer(IntKi), intent( out) :: ErrStat3 + character(ErrMsgLen), intent( out) :: ErrMsg3 integer(IntKi) :: iNode real(R8Ki) :: theta(3) real(R8Ki) :: Orient(3,3) ! Set mesh corresponding to input motions do iNode=1,NumNodePts theta = real(tmpNodePos(4:6,iNode),DbKi) ! convert ReKi to DbKi to avoid roundoff - Orient = EulerConstruct( theta ) + CALL SmllRotTrans( 'InputRotation', theta(1), theta(2), theta(3), Orient, 'Orient', ErrStat3, ErrMsg3 ) HD_MotionMesh%TranslationDisp(1:3,iNode) = tmpNodePos(1:3,iNode) - HD_MotionMesh%Position(1:3,iNode) ! relative displacement only HD_MotionMesh%Orientation(1:3,1:3,iNode) = Orient HD_MotionMesh%TranslationVel( 1:3,iNode) = tmpNodeVel(1:3,iNode) From fad2f0e7fe6bc9d20664b865baf079d867230ec6 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Thu, 29 Jul 2021 18:57:27 -0600 Subject: [PATCH 13/22] HD_c_lib: updates in r-test example driver case --- modules/hydrodyn/src/HydroDyn_C.f90 | 2 -- reg_tests/r-test | 2 +- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/modules/hydrodyn/src/HydroDyn_C.f90 b/modules/hydrodyn/src/HydroDyn_C.f90 index e7bbac27fe..70b1bd42d6 100644 --- a/modules/hydrodyn/src/HydroDyn_C.f90 +++ b/modules/hydrodyn/src/HydroDyn_C.f90 @@ -642,8 +642,6 @@ SUBROUTINE HydroDyn_CalcOutput_c(Time_C, NumNodePts_C, NodePos_C, NodeVel_C, Nod tmpNodeVel(1:6,1:NumNodePts) = reshape( real(NodeVel_C(1:6*NumNodePts),ReKi), (/6,NumNodePts/) ) tmpNodeAcc(1:6,1:NumNodePts) = reshape( real(NodeAcc_C(1:6*NumNodePts),ReKi), (/6,NumNodePts/) ) -write(20,*) Time,NodePos_C,NodeVel_C,NodeAcc_C -write(21,*) Time,tmpNodePos,tmpNodeVel,tmpNodeAcc ! Transfer motions to input meshes call Set_MotionMesh( ErrStat2, ErrMsg2 ) ! update motion mesh with input motion arrays diff --git a/reg_tests/r-test b/reg_tests/r-test index 1c89781e91..9d9e8f539a 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 1c89781e9147ea01416a813f84674c3df7c95ddb +Subproject commit 9d9e8f539a55b6d7e649ae0c291f892543aac3f9 From 9590ff165c7f7e5000e10e4a495b7eff4c861da6 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Fri, 6 Aug 2021 14:47:50 -0600 Subject: [PATCH 14/22] HD_c_lib: add notes about potential channel name buffer overruns and possible fixes --- modules/hydrodyn/python-lib/hydrodyn_library.py | 3 +++ modules/hydrodyn/src/HydroDyn_C.f90 | 12 ++++++++++-- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/modules/hydrodyn/python-lib/hydrodyn_library.py b/modules/hydrodyn/python-lib/hydrodyn_library.py index 9d676262e8..b949bed9c6 100644 --- a/modules/hydrodyn/python-lib/hydrodyn_library.py +++ b/modules/hydrodyn/python-lib/hydrodyn_library.py @@ -115,6 +115,9 @@ def __init__(self, library_path): # This is not sufficient for HD #FIXME: ChanLen may not always be 20 -- could be as much as 256 + # Possible fix is to pass this length over to Fortran side. + # Also may want to convert this at some point to C_NULL_CHAR + # delimeter instead of fixed width. Future problem though. # Number of channel names may exceeed 5000 self._channel_names_c = create_string_buffer(20 * 4000) self._channel_units_c = create_string_buffer(20 * 4000) diff --git a/modules/hydrodyn/src/HydroDyn_C.f90 b/modules/hydrodyn/src/HydroDyn_C.f90 index 70b1bd42d6..cb84cfd651 100644 --- a/modules/hydrodyn/src/HydroDyn_C.f90 +++ b/modules/hydrodyn/src/HydroDyn_C.f90 @@ -39,6 +39,13 @@ MODULE HydroDyn_C integer(IntKi), parameter :: ErrMsgLen_C = 1025 integer(IntKi), parameter :: IntfStrLen = 1025 ! length of other strings through the C interface + !------------------------------------------------------------------------------------ + ! Potential issues + ! - if MaxHDOutputs is sufficiently large, we may overrun the buffer on the Python + ! side (OutputChannelNames_C,OutputChannelUnits_C). Don't have a good method to + ! check this in code yet. Might be best to pass the max length over to Init and + ! do some checks here. May also want to convert this to C_NULL_CHAR delimiter + ! instead of fixed width. !------------------------------------------------------------------------------------ ! Data storage @@ -195,7 +202,7 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen real(c_double), intent(in ) :: DT_C !< Timestep used with HD for stepping forward from t to t+dt. Must be constant. real(c_double), intent(in ) :: TMax_C !< Maximum time for simulation (used to set arrays for wave kinematics) 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*MaxHDOutputs+1) + character(kind=c_char), intent( out) :: OutputChannelNames_C(ChanLen*MaxHDOutputs+1) !< NOTE: if MaxHDOutputs is sufficiently large, we may overrun the buffer on the Python side. character(kind=c_char), intent( out) :: OutputChannelUnits_C(ChanLen*MaxHDOutputs+1) 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) @@ -379,6 +386,8 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen NumChannels_C = size(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 j=1,ChanLen ! max length of channel name. Same for units @@ -696,7 +705,6 @@ SUBROUTINE HydroDyn_UpdateStates_c( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, #endif real(c_double), intent(in ) :: Time_C real(c_double), intent(in ) :: TimeNext_C -!FIXME: estimated inputs at t+dt, or inputs at t???? integer(c_int), intent(in ) :: NumNodePts_C !< Number of mesh points we are transfering motions to and output loads to real(c_float), intent(in ) :: NodePos_C( 6*NumNodePts_C ) !< A 6xNumNodePts_C array [x,y,z,Rx,Ry,Rz] -- positions (global) real(c_float), intent(in ) :: NodeVel_C( 6*NumNodePts_C ) !< A 6xNumNodePts_C array [Vx,Vy,Vz,RVx,RVy,RVz] -- velocities (global) From 9f8a537434135996aef9420a20e052b4e8f9ad0f Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Fri, 6 Aug 2021 16:32:56 -0600 Subject: [PATCH 15/22] HD_c_lib: add regression test for HD with python --- reg_tests/CTestList.cmake | 12 ++ reg_tests/executeHydrodynPyRegressionCase.py | 144 +++++++++++++++++++ reg_tests/r-test | 2 +- 3 files changed, 157 insertions(+), 1 deletion(-) create mode 100644 reg_tests/executeHydrodynPyRegressionCase.py diff --git a/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index d617fdbf15..5e2966bb25 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -135,6 +135,15 @@ function(hd_regression TESTNAME LABEL) regression(${TEST_SCRIPT} ${HYDRODYN_EXECUTABLE} ${SOURCE_DIRECTORY} ${BUILD_DIRECTORY} ${TESTNAME} "${LABEL}") endfunction(hd_regression) +# hydrodyn-Py +function(hd_py_regression TESTNAME LABEL) + set(TEST_SCRIPT "${CMAKE_CURRENT_LIST_DIR}/executeHydrodynPyRegressionCase.py") + set(HYDRODYN_EXECUTABLE "${PYTHON_EXECUTABLE}") + set(SOURCE_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/..") + set(BUILD_DIRECTORY "${CTEST_BINARY_DIR}/modules/hydrodyn") + regression(${TEST_SCRIPT} ${HYDRODYN_EXECUTABLE} ${SOURCE_DIRECTORY} ${BUILD_DIRECTORY} ${TESTNAME} "${LABEL}") +endfunction(hd_py_regression) + # subdyn function(sd_regression TESTNAME LABEL) set(TEST_SCRIPT "${CMAKE_CURRENT_LIST_DIR}/executeSubdynRegressionCase.py") @@ -230,6 +239,9 @@ hd_regression("hd_5MW_OC4Semi_WSt_WavesWN" "hydrodyn;offshore") hd_regression("hd_5MW_TLP_DLL_WTurb_WavesIrr_WavesMulti" "hydrodyn;offshore") hd_regression("hd_TaperCylinderPitchMoment" "hydrodyn;offshore") +# HydroDyn-Py regression tests +hd_py_regression("hd_py_5MW_OC4Semi_WSt_WavesWN" "hydrodyn;offshore") + # SubDyn regression tests sd_regression("SD_Cable_5Joints" "subdyn;offshore") sd_regression("SD_PendulumDamp" "subdyn;offshore") diff --git a/reg_tests/executeHydrodynPyRegressionCase.py b/reg_tests/executeHydrodynPyRegressionCase.py new file mode 100644 index 0000000000..6f8107cd4a --- /dev/null +++ b/reg_tests/executeHydrodynPyRegressionCase.py @@ -0,0 +1,144 @@ +# +# Copyright 2017 National Renewable Energy Laboratory +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +""" + This program executes HydroDyn and a regression test for a single test case. + The test data is contained in a git submodule, r-test, which must be initialized + prior to running. See the r-test README or OpenFAST documentation for more info. + + Get usage with: `executeHydrodynRegressionCase.py -h` +""" + +import os +import sys +basepath = os.path.dirname(__file__) +sys.path.insert(0, os.path.sep.join([basepath, "lib"])) +import argparse +import shutil +import glob +import subprocess +import rtestlib as rtl +import openfastDrivers +import pass_fail +from errorPlotting import exportCaseSummary + +##### Main program + +### Store the python executable for future python calls +pythonCommand = sys.executable + +### Verify input arguments +parser = argparse.ArgumentParser(description="Executes HydroDyn and a regression test for a single test case.") +parser.add_argument("caseName", metavar="Case-Name", type=str, nargs=1, help="The name of the test case.") +parser.add_argument("executable", metavar="HydroDyn-Python", type=str, nargs=1, help="The path to the Python executable.") +parser.add_argument("sourceDirectory", metavar="path/to/openfast_repo", type=str, nargs=1, help="The path to the OpenFAST repository.") +parser.add_argument("buildDirectory", metavar="path/to/openfast_repo/build", type=str, nargs=1, help="The path to the OpenFAST repository build directory.") +parser.add_argument("tolerance", metavar="Test-Tolerance", type=float, nargs=1, help="Tolerance defining pass or failure in the regression test.") +parser.add_argument("systemName", metavar="System-Name", type=str, nargs=1, help="The current system\'s name: [Darwin,Linux,Windows]") +parser.add_argument("compilerId", metavar="Compiler-Id", type=str, nargs=1, help="The compiler\'s id: [Intel,GNU]") +parser.add_argument("-p", "-plot", dest="plot", default=False, metavar="Plotting-Flag", type=bool, nargs="?", help="bool to include matplotlib plots in failed cases") +parser.add_argument("-n", "-no-exec", dest="noExec", default=False, metavar="No-Execution", type=bool, nargs="?", help="bool to prevent execution of the test cases") +parser.add_argument("-v", "-verbose", dest="verbose", default=False, metavar="Verbose-Flag", type=bool, nargs="?", help="bool to include verbose system output") + +args = parser.parse_args() + +caseName = args.caseName[0] +executable = args.executable[0] +sourceDirectory = args.sourceDirectory[0] +buildDirectory = args.buildDirectory[0] +tolerance = args.tolerance[0] +plotError = args.plot if args.plot is False else True +noExec = args.noExec if args.noExec is False else True +verbose = args.verbose if args.verbose is False else True + +# validate inputs +rtl.validateExeOrExit(executable) +rtl.validateDirOrExit(sourceDirectory) +if not os.path.isdir(buildDirectory): + os.makedirs(buildDirectory) + +### Build the filesystem navigation variables for running the test case +regtests = os.path.join(sourceDirectory, "reg_tests") +lib = os.path.join(regtests, "lib") +rtest = os.path.join(regtests, "r-test") +moduleDirectory = os.path.join(rtest, "modules", "hydrodyn") +inputsDirectory = os.path.join(moduleDirectory, caseName) +targetOutputDirectory = os.path.join(inputsDirectory) +testBuildDirectory = os.path.join(buildDirectory, caseName) + +# verify all the required directories exist +if not os.path.isdir(rtest): + rtl.exitWithError("The test data directory, {}, does not exist. If you haven't already, run `git submodule update --init --recursive`".format(rtest)) +if not os.path.isdir(targetOutputDirectory): + rtl.exitWithError("The test data outputs directory, {}, does not exist. Try running `git submodule update`".format(targetOutputDirectory)) +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)) + +### Run HydroDyn on the test case +if not noExec: + caseInputFile = os.path.join(testBuildDirectory, "hydrodyn_driver.py") + returnCode = openfastDrivers.runHydrodynDriverCase(caseInputFile, executable) + if returnCode != 0: + rtl.exitWithError("") + +### Build the filesystem navigation variables for running the regression test +localOutFile = os.path.join(testBuildDirectory, "hd_py.out") +baselineOutFile = os.path.join(targetOutputDirectory, "hd_py.out") +rtl.validateFileOrExit(localOutFile) +rtl.validateFileOrExit(baselineOutFile) + +testData, testInfo, testPack = pass_fail.readFASTOut(localOutFile) +baselineData, baselineInfo, _ = pass_fail.readFASTOut(baselineOutFile) +performance = pass_fail.calculateNorms(testData, baselineData) +normalizedNorm = performance[:, 1] + +# export all case summaries +results = list(zip(testInfo["attribute_names"], [*performance])) +results_max = performance.max(axis=0) +exportCaseSummary(testBuildDirectory, caseName, results, results_max, tolerance) + +# failing case +if not pass_fail.passRegressionTest(normalizedNorm, tolerance): + if plotError: + from errorPlotting import finalizePlotDirectory, plotOpenfastError + ixFailChannels = [i for i in range(len(testInfo["attribute_names"])) if normalizedNorm[i] > tolerance] + failChannels = [channel for i, channel in enumerate(testInfo["attribute_names"]) if i in ixFailChannels] + failResults = [res for i, res in enumerate(results) if i in ixFailChannels] + for channel in failChannels: + try: + plotOpenfastError(localOutFile, baselineOutFile, channel) + except: + error = sys.exc_info()[1] + print("Error generating plots: {}".format(error.msg)) + finalizePlotDirectory(localOutFile, failChannels, caseName) + sys.exit(1) + +# passing case +sys.exit(0) diff --git a/reg_tests/r-test b/reg_tests/r-test index 9d9e8f539a..f415e17908 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 9d9e8f539a55b6d7e649ae0c291f892543aac3f9 +Subproject commit f415e17908a23fc381cb2e3d5f273878091387af From 08aba1a1043438d550fbeb2c8fb2dd963e68054c Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Tue, 31 Aug 2021 11:19:30 -0600 Subject: [PATCH 16/22] hd_c_lib: update interp order in test case to match corresponding OF case --- reg_tests/r-test | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reg_tests/r-test b/reg_tests/r-test index f415e17908..d9a2684569 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit f415e17908a23fc381cb2e3d5f273878091387af +Subproject commit d9a26845690370a06c335230645956ba30803c26 From 5bbc63e8309d899fe1f981373253c015a138d892 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Wed, 8 Dec 2021 08:54:31 -0700 Subject: [PATCH 17/22] HD_C_lib: renaming HydroDyn_C.f90 to HydroDyn_C_Binding.f90 --- modules/hydrodyn/CMakeLists.txt | 8 ++++---- .../src/{HydroDyn_C.f90 => HydroDyn_C_Binding.f90} | 0 2 files changed, 4 insertions(+), 4 deletions(-) rename modules/hydrodyn/src/{HydroDyn_C.f90 => HydroDyn_C_Binding.f90} (100%) diff --git a/modules/hydrodyn/CMakeLists.txt b/modules/hydrodyn/CMakeLists.txt index a368b16404..ccf0a80cad 100644 --- a/modules/hydrodyn/CMakeLists.txt +++ b/modules/hydrodyn/CMakeLists.txt @@ -60,10 +60,10 @@ add_library(hydrodynlib ${HYDRODYN_SOURCES}) target_link_libraries(hydrodynlib nwtclibs) # c-bindings interface library -add_library(hydrodyn_c_lib SHARED src/HydroDyn_C.f90) -target_link_libraries(hydrodyn_c_lib hydrodynlib) +add_library(hydrodyn_c_binding SHARED src/HydroDyn_C_Binding.f90) +target_link_libraries(hydrodyn_c_binding hydrodynlib) if(APPLE OR UNIX) - target_compile_definitions(hydrodyn_c_lib PUBLIC -DIMPLICIT_DLLEXPORT) + target_compile_definitions(hydrodyn_c_binding PUBLIC -DIMPLICIT_DLLEXPORT) endif() add_executable(hydrodyn_driver src/HydroDyn_DriverCode.f90) @@ -73,7 +73,7 @@ target_link_libraries(hydrodyn_driver hydrodynlib nwtclibs versioninfolib) # src/SS_Radiation_DriverCode.f90) #target_link_libraries(ss_radiation hydrodynlib nwtclibs) -install(TARGETS hydrodynlib hydrodyn_driver hydrodyn_c_lib +install(TARGETS hydrodynlib hydrodyn_driver hydrodyn_c_binding EXPORT "${CMAKE_PROJECT_NAME}Libraries" RUNTIME DESTINATION bin LIBRARY DESTINATION lib diff --git a/modules/hydrodyn/src/HydroDyn_C.f90 b/modules/hydrodyn/src/HydroDyn_C_Binding.f90 similarity index 100% rename from modules/hydrodyn/src/HydroDyn_C.f90 rename to modules/hydrodyn/src/HydroDyn_C_Binding.f90 From f1f4dfc42afc6229e889ffff4e5fcc6a6a4f1dd7 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Wed, 8 Dec 2021 09:25:20 -0700 Subject: [PATCH 18/22] hd_c_lib: Update routine names to match how was done with ifw_c_lib --- .../hydrodyn/python-lib/hydrodyn_library.py | 30 +++++----- modules/hydrodyn/src/HydroDyn_C_Binding.f90 | 58 +++++++++---------- reg_tests/r-test | 2 +- 3 files changed, 45 insertions(+), 45 deletions(-) diff --git a/modules/hydrodyn/python-lib/hydrodyn_library.py b/modules/hydrodyn/python-lib/hydrodyn_library.py index b949bed9c6..ec4926d36f 100644 --- a/modules/hydrodyn/python-lib/hydrodyn_library.py +++ b/modules/hydrodyn/python-lib/hydrodyn_library.py @@ -164,7 +164,7 @@ def __init__(self, library_path): # _initialize_routines() ------------------------------------------------------------------------------------------------------------ def _initialize_routines(self): - self.HydroDyn_Init_c.argtypes = [ + self.HydroDyn_C_Init.argtypes = [ POINTER(c_char), # OutRootName POINTER(c_char_p), # input file string POINTER(c_int), # input file string length @@ -186,9 +186,9 @@ def _initialize_routines(self): POINTER(c_int), # ErrStat_C POINTER(c_char) # ErrMsg_C ] - self.HydroDyn_Init_c.restype = c_int + self.HydroDyn_C_Init.restype = c_int - self.HydroDyn_CalcOutput_c.argtypes = [ + self.HydroDyn_C_CalcOutput.argtypes = [ POINTER(c_double), # Time_C POINTER(c_int), # numNodePts -- number of points expecting motions/loads POINTER(c_float), # nodePos -- node positions in flat array of 6*numNodePts @@ -199,9 +199,9 @@ def _initialize_routines(self): POINTER(c_int), # ErrStat_C POINTER(c_char) # ErrMsg_C ] - self.HydroDyn_CalcOutput_c.restype = c_int + self.HydroDyn_C_CalcOutput.restype = c_int - self.HydroDyn_UpdateStates_c.argtypes = [ + self.HydroDyn_C_UpdateStates.argtypes = [ POINTER(c_double), # Time_C POINTER(c_double), # TimeNext_C POINTER(c_int), # numNodePts -- number of points expecting motions/loads @@ -211,13 +211,13 @@ def _initialize_routines(self): POINTER(c_int), # ErrStat_C POINTER(c_char) # ErrMsg_C ] - self.HydroDyn_CalcOutput_c.restype = c_int + self.HydroDyn_C_UpdateStates.restype = c_int - self.HydroDyn_End_c.argtypes = [ + self.HydroDyn_C_End.argtypes = [ POINTER(c_int), # ErrStat_C POINTER(c_char) # ErrMsg_C ] - self.HydroDyn_End_c.restype = c_int + self.HydroDyn_C_End.restype = c_int # hydrodyn_init ------------------------------------------------------------------------------------------------------------ def hydrodyn_init(self, input_string_array): @@ -265,8 +265,8 @@ def hydrodyn_init(self, input_string_array): nodeInitLoc_flat_c[i] = c_float(p) - # call HydroDyn_Init_c - self.HydroDyn_Init_c( + # 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 @@ -326,8 +326,8 @@ def hydrodyn_calcOutput(self, time, nodePos, nodeVel, nodeAcc, nodeFrcMom, outpu # Set up output channels outputChannelValues_c = (c_float * self.numChannels)(0.0,) - # Run HydroDyn_CalcOutput_c - self.HydroDyn_CalcOutput_c( + # Run HydroDyn_C_CalcOutput + self.HydroDyn_C_CalcOutput( byref(c_double(time)), # IN: time at which to calculate output forces byref(c_int(self.numNodePts)), # IN: number of attachment points expected (where motions are transferred into HD) nodePos_flat_c, # IN: positions - specified by user @@ -385,7 +385,7 @@ def hydrodyn_updateStates(self, time, timeNext, nodePos, nodeVel, nodeAcc, nodeF nodeFrc_flat_c = (c_float * (6 * self.numNodePts))(0.0,) # Run HydroDyn_UpdateStates_c - self.HydroDyn_UpdateStates_c( + self.HydroDyn_C_UpdateStates( 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(c_int(self.numNodePts)), # IN: number of attachment points expected (where motions are transferred into HD) @@ -402,8 +402,8 @@ def hydrodyn_updateStates(self, time, timeNext, nodePos, nodeVel, nodeAcc, nodeF def hydrodyn_end(self): if not self.ended: self.ended = True - # Run HydroDyn_End_c - self.HydroDyn_End_c( + # Run HydroDyn_C_End + self.HydroDyn_C_End( byref(self.error_status_c), self.error_message_c ) diff --git a/modules/hydrodyn/src/HydroDyn_C_Binding.f90 b/modules/hydrodyn/src/HydroDyn_C_Binding.f90 index cb84cfd651..dede85911c 100644 --- a/modules/hydrodyn/src/HydroDyn_C_Binding.f90 +++ b/modules/hydrodyn/src/HydroDyn_C_Binding.f90 @@ -17,7 +17,7 @@ ! limitations under the License. ! !********************************************************************************************************************************** -MODULE HydroDyn_C +MODULE HydroDyn_C_BINDING USE ISO_C_BINDING USE HydroDyn @@ -26,10 +26,10 @@ MODULE HydroDyn_C IMPLICIT NONE - PUBLIC :: HydroDyn_Init_C - PUBLIC :: HydroDyn_CalcOutput_C - PUBLIC :: HydroDyn_UpdateStates_C - PUBLIC :: HydroDyn_End_C + PUBLIC :: HydroDyn_C_Init + PUBLIC :: HydroDyn_C_CalcOutput + PUBLIC :: HydroDyn_C_UpdateStates + PUBLIC :: HydroDyn_C_End !------------------------------------------------------------------------------------ ! Error handling @@ -169,18 +169,18 @@ end subroutine SetErr !=============================================================================================================== !--------------------------------------------- HydroDyn Init---------------------------------------------------- !=============================================================================================================== -SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLength_C, & +SUBROUTINE HydroDyn_C_Init( OutRootName_C, InputFileString_C, InputFileStringLength_C, & Gravity_C, defWtrDens_C, defWtrDpth_C, defMSL2SWL_C, & PtfmRefPtPositionX_C, PtfmRefPtPositionY_C, & NumNodePts_C, InitNodePositions_C, & !NumWaveElev_C, WaveElevXY_C & !Placeholder for later InterpOrder_C, T_initial_C, DT_C, TMax_C, & NumChannels_C, OutputChannelNames_C, OutputChannelUnits_C, & - ErrStat_C, ErrMsg_C) BIND (C, NAME='HydroDyn_Init_c') + ErrStat_C, ErrMsg_C) BIND (C, NAME='HydroDyn_C_Init') implicit none #ifndef IMPLICIT_DLLEXPORT -!DEC$ ATTRIBUTES DLLEXPORT :: HydroDyn_Init_c -!GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_Init_c +!DEC$ ATTRIBUTES DLLEXPORT :: HydroDyn_C_Init +!GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_C_Init #endif character(kind=c_char), intent(in ) :: OutRootName_C(IntfStrLen) !< Root name to use for echo files and other @@ -217,7 +217,7 @@ SUBROUTINE HydroDyn_Init_c( OutRootName_C, InputFileString_C, InputFileStringLen 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_Init_C' !< for error handling + character(*), parameter :: RoutineName = 'HydroDyn_C_Init' !< for error handling ! Initialize error handling ErrStat = ErrID_None @@ -599,19 +599,19 @@ subroutine CheckDepth(ErrStat3,ErrMsg3) endif end subroutine CheckDepth -END SUBROUTINE HydroDyn_Init_c +END SUBROUTINE HydroDyn_C_Init !=============================================================================================================== !--------------------------------------------- HydroDyn CalcOutput --------------------------------------------- !=============================================================================================================== -SUBROUTINE HydroDyn_CalcOutput_c(Time_C, NumNodePts_C, NodePos_C, NodeVel_C, NodeAcc_C, & - NodeFrc_C, OutputChannelValues_C, ErrStat_C, ErrMsg_C) BIND (C, NAME='HydroDyn_CalcOutput_c') +SUBROUTINE HydroDyn_C_CalcOutput(Time_C, NumNodePts_C, NodePos_C, NodeVel_C, NodeAcc_C, & + NodeFrc_C, OutputChannelValues_C, ErrStat_C, ErrMsg_C) BIND (C, NAME='HydroDyn_C_CalcOutput') implicit none #ifndef IMPLICIT_DLLEXPORT -!DEC$ ATTRIBUTES DLLEXPORT :: HydroDyn_CalcOutput_c -!GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_CalcOutput_c +!DEC$ ATTRIBUTES DLLEXPORT :: HydroDyn_C_CalcOutput +!GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_C_CalcOutput #endif real(c_double), intent(in ) :: Time_C integer(c_int), intent(in ) :: NumNodePts_C !< Number of mesh points we are transfering motions to and output loads to @@ -630,7 +630,7 @@ SUBROUTINE HydroDyn_CalcOutput_c(Time_C, NumNodePts_C, NodePos_C, NodeVel_C, Nod character(ErrMsgLen) :: ErrMsg !< aggregated error message integer(IntKi) :: ErrStat2 !< temporary error status from a call character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call - character(*), parameter :: RoutineName = 'HydroDyn_CalcOutput_c' !< for error handling + character(*), parameter :: RoutineName = 'HydroDyn_C_CalcOutput' !< for error handling ! Initialize error handling ErrStat = ErrID_None @@ -686,7 +686,7 @@ logical function Failed() Failed = ErrStat >= AbortErrLev if (Failed) call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) end function Failed -END SUBROUTINE HydroDyn_CalcOutput_c +END SUBROUTINE HydroDyn_C_CalcOutput !=============================================================================================================== !--------------------------------------------- HydroDyn UpdateStates ------------------------------------------- @@ -696,12 +696,12 @@ END SUBROUTINE HydroDyn_CalcOutput_c !! Since we don't really know if we are doing correction steps or not, we will track the previous state and !! reset to those if we are repeating a timestep (normally this would be handled by the OF glue code, but since !! the states are not passed across the interface, we must handle them here). -SUBROUTINE HydroDyn_UpdateStates_c( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, NodeVel_C, NodeAcc_C, & - ErrStat_C, ErrMsg_C) BIND (C, NAME='HydroDyn_UpdateStates_c') +SUBROUTINE HydroDyn_C_UpdateStates( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, NodeVel_C, NodeAcc_C, & + ErrStat_C, ErrMsg_C) BIND (C, NAME='HydroDyn_C_UpdateStates') implicit none #ifndef IMPLICIT_DLLEXPORT -!DEC$ ATTRIBUTES DLLEXPORT :: HydroDyn_UpdateStates_c -!GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_UpdateStates_c +!DEC$ ATTRIBUTES DLLEXPORT :: HydroDyn_C_UpdateStates +!GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_C_UpdateStates #endif real(c_double), intent(in ) :: Time_C real(c_double), intent(in ) :: TimeNext_C @@ -719,7 +719,7 @@ SUBROUTINE HydroDyn_UpdateStates_c( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, character(ErrMsgLen) :: ErrMsg !< aggregated error message integer(IntKi) :: ErrStat2 !< temporary error status from a call character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call - character(*), parameter :: RoutineName = 'HydroDyn_UpdateStates_c' !< for error handling + character(*), parameter :: RoutineName = 'HydroDyn_C_UpdateStates' !< for error handling ! Initialize error handling ErrStat = ErrID_None @@ -841,18 +841,18 @@ logical function Failed() Failed = ErrStat >= AbortErrLev if (Failed) call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) end function Failed -END SUBROUTINE HydroDyn_UpdateStates_c +END SUBROUTINE HydroDyn_C_UpdateStates !=============================================================================================================== !--------------------------------------------------- HydroDyn End----------------------------------------------- !=============================================================================================================== ! NOTE: the error handling in this routine is slightly different than the other routines -SUBROUTINE HydroDyn_End_C(ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_End_c') +SUBROUTINE HydroDyn_C_End(ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_C_End') implicit none #ifndef IMPLICIT_DLLEXPORT -!DEC$ ATTRIBUTES DLLEXPORT :: HydroDyn_End_c -!GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_End_c +!DEC$ ATTRIBUTES DLLEXPORT :: HydroDyn_C_End +!GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_C_End #endif integer(c_int), intent( out) :: ErrStat_C character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) @@ -878,7 +878,7 @@ SUBROUTINE HydroDyn_End_C(ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_End_c') ! Call the main subroutine HydroDyn_End ! If u is not allocated, then we didn't get far at all in initialization, - ! or HD_End_C 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 ) @@ -939,7 +939,7 @@ subroutine ClearMesh() call NWTC_Library_Destroymeshmaptype( Map_HD_Mo_P_2_Load , ErrStat2, ErrMsg2 ) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) end subroutine ClearMesh -END SUBROUTINE HydroDyn_End_c +END SUBROUTINE HydroDyn_C_End !> This routine is operating on module level data, hence few inputs @@ -1029,4 +1029,4 @@ subroutine Set_OutputLoadArray() enddo end subroutine Set_OutputLoadArray -END MODULE HydroDyn_C +END MODULE HydroDyn_C_BINDING diff --git a/reg_tests/r-test b/reg_tests/r-test index 59f7124681..9879f8931d 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 59f712468109e06a99b3b76b14de7b93131fb82e +Subproject commit 9879f8931d75e44ab6ad00be82bdbb15020191f2 From 43415b7d59a8d3b914fe1d361a05dbfa51d12106 Mon Sep 17 00:00:00 2001 From: Rafael M Mudafort Date: Tue, 3 May 2022 16:06:23 -0500 Subject: [PATCH 19/22] Fix format and typo in comment --- modules/hydrodyn/python-lib/hydrodyn_library.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/modules/hydrodyn/python-lib/hydrodyn_library.py b/modules/hydrodyn/python-lib/hydrodyn_library.py index ec4926d36f..b52e2eb74f 100644 --- a/modules/hydrodyn/python-lib/hydrodyn_library.py +++ b/modules/hydrodyn/python-lib/hydrodyn_library.py @@ -2,7 +2,7 @@ # LICENSING # Copyright (C) 2021 National Renewable Energy Laboratory # -# This file is part of InflowWind. +# This file is part of HydroDyn. # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. @@ -67,15 +67,11 @@ POINTER, create_string_buffer, byref, - c_byte, c_int, c_double, c_float, c_char, c_char_p, - c_wchar, - c_wchar_p, - c_bool ) import numpy as np import datetime @@ -466,7 +462,7 @@ def check_input_motions(self,nodePos,nodeVel,nodeAcc): @property def output_channel_names(self): if len(self._channel_names_c.value.split()) == 0: - return [] + return [] output_channel_names = self._channel_names_c.value.split() output_channel_names = [n.decode('UTF-8') for n in output_channel_names] return output_channel_names From de120164ec4725cf4288445703c475fc0843ef9b Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Fri, 6 May 2022 08:17:00 -0600 Subject: [PATCH 20/22] hd_c_lib: send current time to UpdateStates (was incorrectly sending next timestep time) --- modules/hydrodyn/src/HydroDyn_C_Binding.f90 | 2 +- reg_tests/r-test | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/modules/hydrodyn/src/HydroDyn_C_Binding.f90 b/modules/hydrodyn/src/HydroDyn_C_Binding.f90 index dede85911c..d4cda71634 100644 --- a/modules/hydrodyn/src/HydroDyn_C_Binding.f90 +++ b/modules/hydrodyn/src/HydroDyn_C_Binding.f90 @@ -810,7 +810,7 @@ SUBROUTINE HydroDyn_C_UpdateStates( Time_C, TimeNext_C, NumNodePts_C, NodePos_C, ! Call the main subroutine HydroDyn_UpdateStates to get the velocities - CALL HydroDyn_UpdateStates( InputTimes(INPUT_PRED), 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, u, InputTimes, p, x(STATE_PRED), xd(STATE_PRED), z(STATE_PRED), OtherStates(STATE_PRED), m, ErrStat2, ErrMsg2 ) if (Failed()) return diff --git a/reg_tests/r-test b/reg_tests/r-test index 03ac75f90a..37b2a6d79c 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 03ac75f90a43b3a15409a3620587b9a0ab646154 +Subproject commit 37b2a6d79cf797fe78f165e67b0a10b52069f860 From 7b6f870ac6ca2008c12c53ffb3de5defedbb4356 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Wed, 29 Jun 2022 15:48:12 -0600 Subject: [PATCH 21/22] HD_c_lib: update GH actions to run the python test as an interface test --- .github/actions/tests-module-hydrodyn/action.yml | 2 +- .github/workflows/automated-dev-tests.yml | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/actions/tests-module-hydrodyn/action.yml b/.github/actions/tests-module-hydrodyn/action.yml index 3e94eba00d..eb453ec203 100644 --- a/.github/actions/tests-module-hydrodyn/action.yml +++ b/.github/actions/tests-module-hydrodyn/action.yml @@ -4,6 +4,6 @@ author: 'Rafael Mudafort https://github.com/rafmudaf' runs: using: "composite" steps: - - run: ctest -VV -j7 -R hd_ + - run: ctest -VV -j7 -R hd_ -LE python working-directory: ${{runner.workspace}}/openfast/build shell: bash diff --git a/.github/workflows/automated-dev-tests.yml b/.github/workflows/automated-dev-tests.yml index 204b95fa45..007bc4e6c8 100644 --- a/.github/workflows/automated-dev-tests.yml +++ b/.github/workflows/automated-dev-tests.yml @@ -455,6 +455,7 @@ jobs: cmake --build . --target openfastlib -- -j ${{env.NUM_PROCS}} cmake --build . --target openfast_cpp -- -j ${{env.NUM_PROCS}} cmake --build . --target openfastcpp -- -j ${{env.NUM_PROCS}} + cmake --build . --target hydrodyn_c_binding -- -j ${{env.NUM_PROCS}} cmake --build . --target ifw_c_binding -- -j ${{env.NUM_PROCS}} cmake --build . --target regression_tests -- -j ${{env.NUM_PROCS}} @@ -493,5 +494,6 @@ jobs: name: c-interface-reg-tests path: | ${{runner.workspace}}/openfast/build/reg_tests/glue-codes/openfast-cpp + ${{runner.workspace}}/openfast/build/reg_tests/modules/hydrodyn ${{runner.workspace}}/openfast/build/reg_tests/modules/inflowwind !${{runner.workspace}}/openfast/build/reg_tests/glue-codes/openfast-cpp/5MW_Baseline From 6426609c5ac28053fe8b7b5d22f914c849248593 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Wed, 29 Jun 2022 16:06:43 -0600 Subject: [PATCH 22/22] HD_c_lib: hd_py_* regression tests missing `pthon` label --- reg_tests/CTestList.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index 915fd66ee4..7542433367 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -311,7 +311,7 @@ hd_regression("hd_5MW_TLP_DLL_WTurb_WavesIrr_WavesMulti" "hydrodyn;offshore") hd_regression("hd_TaperCylinderPitchMoment" "hydrodyn;offshore") # HydroDyn-Py regression tests -hd_py_regression("hd_py_5MW_OC4Semi_WSt_WavesWN" "hydrodyn;offshore") +hd_py_regression("hd_py_5MW_OC4Semi_WSt_WavesWN" "hydrodyn;offshore;python") # SubDyn regression tests sd_regression("SD_Cable_5Joints" "subdyn;offshore")