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 diff --git a/modules/hydrodyn/CMakeLists.txt b/modules/hydrodyn/CMakeLists.txt index 3c10213253..ccf0a80cad 100644 --- a/modules/hydrodyn/CMakeLists.txt +++ b/modules/hydrodyn/CMakeLists.txt @@ -59,6 +59,13 @@ set(HYDRODYN_SOURCES add_library(hydrodynlib ${HYDRODYN_SOURCES}) target_link_libraries(hydrodynlib nwtclibs) +# c-bindings interface library +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_binding PUBLIC -DIMPLICIT_DLLEXPORT) +endif() + add_executable(hydrodyn_driver src/HydroDyn_DriverCode.f90) target_link_libraries(hydrodyn_driver hydrodynlib nwtclibs versioninfolib) @@ -66,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 +install(TARGETS hydrodynlib hydrodyn_driver hydrodyn_c_binding 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..b52e2eb74f --- /dev/null +++ b/modules/hydrodyn/python-lib/hydrodyn_library.py @@ -0,0 +1,621 @@ +#********************************************************************************************************************************** +# LICENSING +# Copyright (C) 2021 National Renewable Energy Laboratory +# +# 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. +# +#********************************************************************************************************************************** +# +# 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. +# +# 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, + create_string_buffer, + byref, + c_int, + c_double, + c_float, + c_char, + c_char_p, +) +import numpy as np +import datetime + +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 + + # 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 + + self._initialize_routines() + self.ended = False # For error handling at end + + + # Create buffers for class data + 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 + # 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) + + # Initial environmental conditions + self.gravity = 9.80665 # Gravity (m/s^2) + self.defWtrDens = 1025.0 # Water density (kg/m^3) + 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) + self.InterpOrder = 1 # default of linear interpolation + + # 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 + self.numTimeSteps = 0 + + 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 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.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 + # root of the file name. + self.outRootName = "Output_HDlib_default" + + # _initialize_routines() ------------------------------------------------------------------------------------------------------------ + def _initialize_routines(self): + self.HydroDyn_C_Init.argtypes = [ + POINTER(c_char), # OutRootName + POINTER(c_char_p), # input file string + POINTER(c_int), # input file string length + POINTER(c_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), # 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 + POINTER(c_char), # output channel names + POINTER(c_char), # output channel units + POINTER(c_int), # ErrStat_C + POINTER(c_char) # ErrMsg_C + ] + self.HydroDyn_C_Init.restype = c_int + + 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 + 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 + ] + self.HydroDyn_C_CalcOutput.restype = c_int + + 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 + 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_int), # ErrStat_C + POINTER(c_char) # ErrMsg_C + ] + self.HydroDyn_C_UpdateStates.restype = c_int + + self.HydroDyn_C_End.argtypes = [ + POINTER(c_int), # ErrStat_C + POINTER(c_char) # ErrMsg_C + ] + self.HydroDyn_C_End.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. + input_string = '\x00'.join(input_string_array) + input_string = input_string.encode('utf-8') + input_string_length = len(input_string) + + self._numChannels_c = c_int(0) + + # 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 + + # 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.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 + # 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.") + + # 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.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) + + + # 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 + 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_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: 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 + 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 + ) + + self.check_error() + + # Initialize output channels + self.numChannels = self._numChannels_c.value + + + # hydrodyn_calcOutput ------------------------------------------------------------------------------------------------------------ + 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) + + # 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_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 + 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]) + + # 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_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) + nodePos_flat_c, # IN: positions - specified by user + nodeVel_flat_c, # IN: velocities at desired positions + nodeAcc_flat_c, # IN: accelerations at desired positions + 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: + self.ended = True + # Run HydroDyn_C_End + self.HydroDyn_C_End( + byref(self.error_status_c), + self.error_message_c + ) + + self.check_error() + + # 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"HydroDyn error status: {self.error_levels[self.error_status_c.value]}: {self.error_message_c.value.decode('ascii')}") + else: + 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.") + + + 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: + 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 + + @property + 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 + + +#=============================================================================== +# 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,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 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,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: + 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.f90 b/modules/hydrodyn/src/HydroDyn.f90 index 7ea184437a..1f8cedf0ee 100644 --- a/modules/hydrodyn/src/HydroDyn.f90 +++ b/modules/hydrodyn/src/HydroDyn.f90 @@ -1793,52 +1793,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_C_Binding.f90 b/modules/hydrodyn/src/HydroDyn_C_Binding.f90 new file mode 100644 index 0000000000..d4cda71634 --- /dev/null +++ b/modules/hydrodyn/src/HydroDyn_C_Binding.f90 @@ -0,0 +1,1032 @@ +!********************************************************************************************************************************** +! LICENSING +! Copyright (C) 2021 National Renewable Energy Lab +! +! 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 HydroDyn_C_BINDING + + USE ISO_C_BINDING + USE HydroDyn + USE HydroDyn_Types + USE NWTC_Library + + IMPLICIT NONE + + PUBLIC :: HydroDyn_C_Init + PUBLIC :: HydroDyn_C_CalcOutput + PUBLIC :: HydroDyn_C_UpdateStates + PUBLIC :: HydroDyn_C_End + + !------------------------------------------------------------------------------------ + ! 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 + 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 + ! 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. + type(HydroDyn_ParameterType) :: p !< Parameters + type(HydroDyn_ContinuousStateType) :: x(0:2) !< continuous states at Time t and t+dt (predicted) + type(HydroDyn_DiscreteStateType) :: xd(0:2) !< discrete states at Time t and t+dt (predicted) + type(HydroDyn_ConstraintStateType) :: z(0:2) !< Constraint states at Time t and t+dt (predicted) + type(HydroDyn_OtherStateType) :: OtherStates(0:2) !< Initial other/optimization states + type(HydroDyn_OutputType) :: y !< Initial output (outputs are not calculated; only the output mesh is initialized) + type(HydroDyn_MiscVarType) :: m !< Misc variables for optimization (not copied in glue code) + !------------------------------ + ! 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. + ! 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. + 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 + !------------------------------------------------------------------------------------ + + + + !------------------------------------------------------------------------------------ + ! 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 + 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 + 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 + ! 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. + !------------------------------------------------------------------------------------ + + +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) + 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 ) + else + ErrMsg_C = TRANSFER( trim(ErrMsg)//C_NULL_CHAR, ErrMsg_C ) + endif +end subroutine SetErr + + +!=============================================================================================================== +!--------------------------------------------- HydroDyn Init---------------------------------------------------- +!=============================================================================================================== +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_C_Init') + implicit none +#ifndef IMPLICIT_DLLEXPORT +!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 + 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] + 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) + 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) !< 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) + + ! Local Variables + character(IntfStrLen) :: OutRootName !< Root name to use for echo files and other + character(kind=C_char, len=InputFileStringLength_C), pointer :: InputFileString !< Input file as a single string with NULL chracter separating lines + + real(DbKi) :: TimeInterval !< timestep for HD + integer(IntKi) :: ErrStat !< aggregated error message + character(ErrMsgLen) :: ErrMsg !< aggregated error message + integer(IntKi) :: ErrStat2 !< temporary error status from a call + character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call + integer(IntKi) :: i,j,k !< generic counters + character(*), parameter :: RoutineName = 'HydroDyn_C_Init' !< for error handling + + ! 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) + + ! Get the data to pass to HD_Init + call InitFileInfo(InputFileString, InitInp%PassedFileData, ErrStat2, ErrMsg2); if (Failed()) return + + ! For diagnostic purposes, the following can be used to display the contents + ! of the InFileInfo data structure. + ! CU is the screen -- system dependent. + !call Print_FileInfo_Struct( CU, InitInp%PassedFileData ) + + ! Set other inputs for calling HydroDyn_Init + InitInp%InputFile = "passed_hd_file" ! dummy + InitInp%UseInputFile = .FALSE. ! this probably should be passed in + InitInp%HasIce = .FALSE. ! Always keep at false unless interfacing to ice modules + ! Linearization + ! for now, set linearization to false. Pass this in later when interface supports it + ! Note: we may want to linearize at T=0 for added mass effects, but that might be + ! special case + InitInp%Linearize = .FALSE. + + ! RootName -- for output of echo or other files + OutRootName = TRANSFER( OutRootName_C, OutRootName ) + i = INDEX(OutRootName,C_NULL_CHAR) - 1 ! if this has a c null character at the end... + if ( i > 0 ) OutRootName = OutRootName(1:I) ! remove it + InitInp%OutRootName = trim(OutRootName) + + ! Values passed in + InitInp%Gravity = REAL(Gravity_C, ReKi) + InitInp%defWtrDens = REAL(defWtrDens_C, ReKi) + 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 + ! - 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 + ! 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) + ! position (no Z). + InitInp%PtfmLocationX = REAL(PtfmRefPtPositionX_C, ReKi) + InitInp%PtfmLocationY = REAL(PtfmRefPtPositionY_C, ReKi) + + + ! Wave eleveation output + ! Wave elevations can be exported for a set of points (grid or any other layout). + ! This feature is used only in the driver codes for exporting for visualization + ! and could be added to this inteface. + ! Skipping this for now. Maybe add later. + !InitInp%WaveElevXY + + + !---------------------------------------------------- + ! Allocate input array u and corresponding InputTimes + !---------------------------------------------------- + ! These inputs are used in the time stepping algorithm within HD_UpdateStates + ! For quadratic interpolation (InterpOrder==2), 3 timesteps are used. For + ! linear (InterOrder==1), 2 timesteps (the HD code can handle either). + ! u(1) inputs at t + ! u(2) inputs at t - dt + ! u(3) inputs at t - 2*dt ! quadratic only + allocate(u(InterpOrder+1), STAT=ErrStat2) + 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 + ! 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(STATE_CURR), xd(STATE_CURR), z(STATE_CURR), OtherStates(STATE_CURR), y, m, TimeInterval, InitOutData, ErrStat2, ErrMsg2 ) + if (Failed()) return + + + !------------------------------------------------------------- + ! Sanity checks + !------------------------------------------------------------- + call CheckDepth(ErrStat2,ErrMsg2); if (Failed()) return + call CheckNodes(ErrStat2,ErrMsg2); if (Failed()) return + + + !------------------------------------------------------------- + ! Set the interface meshes for motion inputs and loads output + !------------------------------------------------------------- + call SetMotionLoadsInterfaceMeshes(ErrStat2,ErrMsg2); if (Failed()) return + + + !------------------------------------------------------------- + ! Setup other prior timesteps + ! We fill InputTimes with negative times, but the Input values are identical for each of those times; this allows + ! 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) - dT_Global ! Initialize 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 +! Any additional warnings or error handling necessary + + + !------------------------------------------------- + ! Set output channel information for driver code + !------------------------------------------------- + + ! Number of channels + 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 + 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 + + + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + +CONTAINS + logical function Failed() + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + if (Failed) then + call FailCleanup() + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + endif + end function Failed + + 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 + !! 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 = tmpNodePos(1:3,iNode) + theta = real(tmpNodePos(4:6,iNode),DbKi) ! convert ReKi to DbKi to avoid roundoff + 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 + + 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 ) + + !------------------------------------------------------------- + ! 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 + 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 = "" + 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. 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 = "" + 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 + +END SUBROUTINE HydroDyn_C_Init + + +!=============================================================================================================== +!--------------------------------------------- HydroDyn CalcOutput --------------------------------------------- +!=============================================================================================================== + +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_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 + 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) :: 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_C_CalcOutput' !< 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) + + ! 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( 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 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 + + + ! 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 + logical function Failed() + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + if (Failed) call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + end function Failed +END SUBROUTINE HydroDyn_C_CalcOutput + +!=============================================================================================================== +!--------------------------------------------- HydroDyn UpdateStates ------------------------------------------- +!=============================================================================================================== +!> 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_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_C_UpdateStates +!GCC$ ATTRIBUTES DLLEXPORT :: HydroDyn_C_UpdateStates +#endif + real(c_double), intent(in ) :: Time_C + real(c_double), intent(in ) :: TimeNext_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) + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_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 + 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_C_UpdateStates' !< for error handling + + ! Initialize error handling + ErrStat = ErrID_None + ErrMsg = "" + CorrectionStep = .false. + + ! 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 + + + !------------------------------------------------------- + ! 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) + ! + ! 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. + 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) = ( N_Global - 1 ) * dT_Global ! u(3) at T-dT + endif + 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 + + + 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( 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. + ! -- 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_CURR), 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 + + + + 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(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + end function Failed +END SUBROUTINE HydroDyn_C_UpdateStates + +!=============================================================================================================== +!--------------------------------------------------- HydroDyn End----------------------------------------------- +!=============================================================================================================== +! NOTE: the error handling in this routine is slightly different than the other routines + +SUBROUTINE HydroDyn_C_End(ErrStat_C,ErrMsg_C) BIND (C, NAME='HydroDyn_C_End') + implicit none +#ifndef IMPLICIT_DLLEXPORT +!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) + + ! Local variables + 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 = "" + + ! 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) + + + ! Call the main subroutine HydroDyn_End + ! If u is not allocated, then we didn't get far at all in initialization, + ! or HD_C_End got called before Init. We don't want a segfault, so check + ! 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) + call HydroDyn_DestroyInput( u(i), ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + enddo + if (allocated(u)) deallocate(u) + endif + + ! Destroy any other copies of states (rerun on (STATE_CURR) is ok) + call HydroDyn_DestroyContState( x( STATE_LAST), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyDiscState( xd( STATE_LAST), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyConstrState( z( STATE_LAST), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyOtherState( OtherStates(STATE_LAST), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyContState( x( STATE_CURR), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyDiscState( xd( STATE_CURR), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyConstrState( z( STATE_CURR), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyOtherState( OtherStates(STATE_CURR), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyContState( x( STATE_PRED), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyDiscState( xd( STATE_PRED), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyConstrState( z( STATE_PRED), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call HydroDyn_DestroyOtherState( OtherStates(STATE_PRED), ErrStat2, ErrMsg2 ); call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + + + ! if deallocate other items now + if (allocated(InputTimes)) deallocate(InputTimes) + + ! 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_C_End + + +!> This routine is operating on module level data, hence few inputs +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 + 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) + 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_BINDING diff --git a/modules/hydrodyn/src/HydroDyn_Output.f90 b/modules/hydrodyn/src/HydroDyn_Output.f90 index 465797a7ce..fc1c73e6bc 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 diff --git a/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index 32dd86b78b..7542433367 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -159,6 +159,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") @@ -301,6 +310,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;python") + # 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 11ec72ba9d..fb50bbfc34 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 11ec72ba9d17b3b7d66b181a5a23b5ad1917fca1 +Subproject commit fb50bbfc347c2942e5b00e91bd3c56234dab4416 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 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +