diff --git a/.github/workflows/automated-dev-tests.yml b/.github/workflows/automated-dev-tests.yml index 401d57bd23..52ed9f8cb2 100644 --- a/.github/workflows/automated-dev-tests.yml +++ b/.github/workflows/automated-dev-tests.yml @@ -372,6 +372,7 @@ jobs: ${{github.workspace}}/build/reg_tests/modules/moordyn ${{github.workspace}}/build/reg_tests/modules/inflowwind ${{github.workspace}}/build/reg_tests/modules/hydrodyn + ${{github.workspace}}/build/reg_tests/modules/seastate !${{github.workspace}}/build/reg_tests/glue-codes/openfast-cpp/5MW_Baseline diff --git a/glue-codes/python/examples/SeaState.py b/glue-codes/python/examples/SeaState.py deleted file mode 100644 index fff806b8ea..0000000000 --- a/glue-codes/python/examples/SeaState.py +++ /dev/null @@ -1,38 +0,0 @@ - - -from pyOpenFAST.seastate import SeaStateLib -import matplotlib.pyplot as plt -import numpy as np - -project_root = '/Users/rmudafor/Development/openfast' -library_path = project_root + '/build/modules/seastate/libseastate_c_binding.dylib' - -dt = 1.0 -time_steps = 10 - -seastatelib = SeaStateLib( - library_path, - "NRELOffshrBsline5MW_OC4DeepCwindSemi_SeaState_WaveMod5.dat" -) - -seastatelib.init( - time_interval=dt, - n_steps=time_steps, -) - -seastate_outputs = np.zeros((time_steps, seastatelib.num_outs)) -for i in range(time_steps): - seastatelib.calc_output(i) - seastate_outputs[i] = seastatelib.output_values - print(i, [f"{value:3.4f} - " for value in seastate_outputs[i]]) -seastatelib.end() - -# Plot the results -# plt.figure(figsize=(10, 6)) -# plt.plot(seastate_outputs[:, 0]) -# plt.plot(seastate_outputs[:, 1]) -# plt.xlabel('Time Step') -# plt.ylabel('Value') -# plt.title('Sea State Outputs') -# plt.legend() -# plt.show() diff --git a/glue-codes/python/pyOpenFAST/aerodyn_inflow.py b/glue-codes/python/pyOpenFAST/aerodyn_inflow.py index da0afd7aa7..5dd46cc794 100644 --- a/glue-codes/python/pyOpenFAST/aerodyn_inflow.py +++ b/glue-codes/python/pyOpenFAST/aerodyn_inflow.py @@ -346,7 +346,7 @@ def adi_preinit(self) -> None: byref(c_float(self.water_depth)), # IN -> water depth byref(c_float(self.mean_sea_level_offset)), # IN -> MSL to SWL offset byref(c_int(self.mhk)), # IN -> mhk flag (0=not MHK, 1=fixed bottom, 2=floating) - byref(c_int(self.debug_level)), # IN -> debug level (0=None to 4=Fatal) + byref(c_int(self.debug_level)), # IN -> debug level (0=None to 4=all meshes) byref(self.error_status_c), # OUT <- error status code self.error_message_c # OUT <- error message buffer ) diff --git a/glue-codes/python/pyOpenFAST/interface_abc.py b/glue-codes/python/pyOpenFAST/interface_abc.py index 3758deff92..d3bf9d737e 100644 --- a/glue-codes/python/pyOpenFAST/interface_abc.py +++ b/glue-codes/python/pyOpenFAST/interface_abc.py @@ -33,7 +33,7 @@ class OpenFASTInterfaceType(CDLL): ERROR_MSG_C_LEN = 8197 # NOTE: the length of the name used for any output file written by the - # HD Fortran code is 1025. + # Fortran code is 1025. default_str_c_len = 1025 abort_error_level = c_int(4) diff --git a/glue-codes/python/pyOpenFAST/seastate.py b/glue-codes/python/pyOpenFAST/seastate.py index d10cb0126f..4ae6ef0420 100644 --- a/glue-codes/python/pyOpenFAST/seastate.py +++ b/glue-codes/python/pyOpenFAST/seastate.py @@ -27,13 +27,31 @@ c_float, c_char, c_char_p, - c_bool + c_bool, + c_void_p ) import numpy as np +import numpy.typing as npt from pathlib import Path +import datetime +import os +from dataclasses import dataclass from .interface_abc import OpenFASTInterfaceType +@dataclass +class MotionData: + """ + POD-style container for motion-related data i.e. state of a node. Only + position information for SeaState + """ + position: npt.NDArray[np.float32] + #orientation: npt.NDArray[np.float64] + #velocity: npt.NDArray[np.float32] + #acceleration: npt.NDArray[np.float32] + + + class SeaStateLib(OpenFASTInterfaceType): """ This is the Python interface to the OpenFAST SeaState module. @@ -46,41 +64,68 @@ class SeaStateLib(OpenFASTInterfaceType): from the last call to calc_output. """ - def __init__(self, library_path: str, input_file_name: str): + def __init__(self, library_path): super().__init__(library_path) - - self.input_file_name = str( Path(input_file_name).absolute() ).encode('utf-8') + self.library_path = library_path self._initialize_routines() - - self.ended = False # For error handling at end + self.ended = False # For error handling at end # Create buffers for class data - # These will generally be overwritten by the Fortran code - self.num_outs_c = c_int(0) - self.output_channel_names = [] - self.output_channel_units = [] - self.output_values = None + #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 buffer for the channel names and units is set arbitrarily large + # to start. Channel name and unit lengths are currently hard + # coded to 20 (this must match ChanLen in NWTC_Base.f90). + self._channel_names_c = create_string_buffer(20 * 4000 + 1) + self._channel_units_c = create_string_buffer(20 * 4000 + 1) + + self.numResPts = 0 # Number of wind points we will + # request information from + # non-CalcOutput routines. + + self.numChannels = 0 # Number of channels returned + + # flags + self.debuglevel = 0 # 0-4 levels + + #-------------------------------------- + # VTK settings + #-------------------------------------- + self.vtk_write = 0 # Default -> no vtk output, 0 none, 1 init, 2 animation + self.vtk_dt = 0. # Default -> all + self.vtk_output_dir = "" # Set to specify a directory relative to input files def _initialize_routines(self): - self.SeaSt_C_Init.argtypes = [ - POINTER(c_char_p), # intent(in ) :: InputFile_c(IntfStrLen) - POINTER(c_char_p), # intent(in ) :: OutRootName_c(IntfStrLen) + self.SeaSt_C_PreInit.argtypes = [ POINTER(c_float), # intent(in ) :: Gravity_c POINTER(c_float), # intent(in ) :: WtrDens_c POINTER(c_float), # intent(in ) :: WtrDpth_c POINTER(c_float), # intent(in ) :: MSL2SWL_c + POINTER(c_int), # intent(in ) :: debuglevel + POINTER(c_char), # intent(in ) :: vtk_output_dir_c + POINTER(c_int), # intent(in ) :: vtk_write + POINTER(c_double), # intent(in ) :: vtk_dt + POINTER(c_int), # intent( out) :: ErrStat_C + POINTER(c_char), # intent( out) :: ErrMsg_C(ErrMsgLen_C) + ] + self.SeaSt_C_PreInit.restype = None + + self.SeaSt_C_Init.argtypes = [ + POINTER(c_char_p), # intent(in ) :: InputFile_c(IntfStrLen) + POINTER(c_char_p), # intent(in ) :: OutRootName_c(IntfStrLen) POINTER(c_int), # intent(in ) :: NSteps_c POINTER(c_float), # intent(in ) :: TimeInterval_c - POINTER(c_int), # intent(in ) :: WaveElevSeriesFlag_c - POINTER(c_int), # intent(in ) :: WrWvKinMod_c POINTER(c_int), # intent( out) :: NumChannels_c POINTER(c_char), # intent( out) :: OutputChannelNames_C POINTER(c_char), # intent( out) :: OutputChannelUnits_C POINTER(c_int), # intent( out) :: ErrStat_C POINTER(c_char), # intent( out) :: ErrMsg_C(ErrMsgLen_C) ] - self.SeaSt_C_Init.restype = c_int + self.SeaSt_C_Init.restype = None self.SeaSt_C_CalcOutput.argtypes = [ POINTER(c_double), # intent(in ) :: Time_C @@ -88,54 +133,158 @@ def _initialize_routines(self): POINTER(c_int), # intent( out) :: ErrStat_C POINTER(c_char), # intent( out) :: ErrMsg_C(ErrMsgLen_C) ] - self.SeaSt_C_CalcOutput.restype = c_int + self.SeaSt_C_CalcOutput.restype = None self.SeaSt_C_End.argtypes = [ POINTER(c_int), # intent( out) :: ErrStat_C POINTER(c_char) # intent( out) :: ErrMsg_C(ErrMsgLen_C) ] - self.SeaSt_C_End.restype = c_int + self.SeaSt_C_End.restype = None + + self.SeaSt_C_GetWaveFieldPointer.argtypes = [ + POINTER(c_void_p), # intent( out) :: pointer to the WaveField data + POINTER(c_int), # intent( out) :: ErrStat_C + POINTER(c_char), # intent( out) :: ErrMsg_C(ErrMsgLen_C) + ] + self.SeaSt_C_GetWaveFieldPointer.restype = None + + self.SeaSt_C_SetWaveFieldPointer.argtypes = [ + POINTER(c_void_p), # intent(in ) :: pointer to the WaveField data + POINTER(c_int), # intent( out) :: ErrStat_C + POINTER(c_char), # intent( out) :: ErrMsg_C(ErrMsgLen_C) + ] + self.SeaSt_C_SetWaveFieldPointer.restype = None + + + self.SeaSt_C_GetFluidVelAcc.argtypes = [ + POINTER(c_double), # intent(in ) :: Time_C + POINTER(c_float), # intent(in ) :: Pos_c(3) + POINTER(c_float), # intent( out) :: Vel_c(3) + POINTER(c_float), # intent( out) :: Acc_c(3) + POINTER(c_int), # intent( out) :: NodeInWater_C + POINTER(c_int), # intent( out) :: ErrStat_C + POINTER(c_char) # intent( out) :: ErrMsg_C(ErrMsgLen_C) + + ] + self.SeaSt_C_GetFluidVelAcc.restype = None - def init( + self.SeaSt_C_GetSurfElev.argtypes = [ + POINTER(c_double), # intent(in ) :: Time_C + POINTER(c_float), # intent(in ) :: Pos_c(3) + POINTER(c_float), # intent( out) :: Elev_C + POINTER(c_int), # intent( out) :: ErrStat_C + POINTER(c_char) # intent( out) :: ErrMsg_C(ErrMsgLen_C) + ] + self.SeaSt_C_GetSurfElev.restype = None + + self.SeaSt_C_GetSurfNorm.argtypes = [ + POINTER(c_double), # intent(in ) :: Time_C + POINTER(c_float), # intent(in ) :: Pos_c(3) + POINTER(c_float), # intent( out) :: norm(3) + POINTER(c_int), # intent( out) :: ErrStat_C + POINTER(c_char) # intent( out) :: ErrMsg_C(ErrMsgLen_C) + ] + self.SeaSt_C_GetSurfNorm.restype = None + + def check_error(self) -> None: + """Checks for and handles any errors from the Fortran library. + + Raises: + RuntimeError: If a fatal error occurs in the Fortran code + """ + # If the error status is 0, return + if self.error_status_c.value == 0: + return + + # Get the error level and error message + error_level = self.error_levels.get( + self.error_status_c.value, + f"Unknown Error Level: {self.error_status_c.value}" + ) + error_msg = self.error_message_c.raw.decode('utf-8').strip() + message = f"WaveTank library {error_level}: {error_msg}" + # If the error level is fatal, call WaveTank_End() and raise an error + if self.error_status_c.value >= self.abort_error_level.value: + try: + self.SeaSt_C_End( + byref(self.error_status_c), # OUT <- error status code + self.error_message_c # OUT <- error message buffer + ) + if self.error_status_c.value == 4: + error_msg = self.error_message_c.raw.decode('utf-8').strip() + print(f'WaveTank_End error: {error_msg}') + except Exception as e: + message += f"\nAdditional error during cleanup: {e}" + raise RuntimeError(message) + else: + print(message) + + + #FIXME: store these elsewhere + def seastate_preinit( self, gravity: float = 9.80665, water_density: float = 1025, water_depth: float = 200, msl2swl: float = 0, + ): + """Set environment variables and general setup + + Args: + + Raises: + ValueError: If values are outside reasonable bounds + RuntimeError: If preinit fails + """ + vtk_output_dir_c = create_string_buffer( + self.vtk_output_dir.ljust(self.default_str_c_len).encode('utf-8') + ) + self.SeaSt_C_PreInit( + byref(c_float(gravity)), + byref(c_float(water_density)), + byref(c_float(water_depth)), + byref(c_float(msl2swl)), + byref(c_int(self.debug_level)), # IN -> debug level (0=None to 4=all meshes) + vtk_output_dir_c, # IN -> directory for vtk output files + byref(c_int(self.vtk_write)), # IN -> write VTK flag + byref(c_double(self.vtk_dt)), # IN -> VTK output time step + byref(self.error_status_c), # OUT <- error status code + self.error_message_c # OUT <- error message buffer + ) + self.check_error() + + + def seastate_init( + self, + primary_ss_file, outrootname: str = "./seastate.SeaSt", - wave_kinematics_mode: int = 0, n_steps: int = 801, time_interval: float = 0.125, - wave_elevation_series_flag: int = 0, ): - _error_status = c_int(0) - _error_message = create_string_buffer(self.ERROR_MSG_C_LEN) # This buffer for the channel names and units is set arbitrarily large # to start. Channel name and unit lengths are currently hard # coded to 20 (this must match ChanLen in NWTC_Base.f90). _channel_names = create_string_buffer(20 * 4000 + 1) _channel_units = create_string_buffer(20 * 4000 + 1) + self._numChannels = c_int(0) + self.SeaSt_C_Init( - c_char_p(self.input_file_name), + c_char_p(primary_ss_file.encode('utf-8')), #FIXME: this might overrun the buffer!!!!! c_char_p(outrootname.encode('utf-8')), - byref(c_float(gravity)), - byref(c_float(water_density)), - byref(c_float(water_depth)), - byref(c_float(msl2swl)), byref(c_int(n_steps)), byref(c_float(time_interval)), - byref(c_int(wave_elevation_series_flag)), - byref(c_int(wave_kinematics_mode)), - byref(self.num_outs_c), + byref(self._numChannels), _channel_names, _channel_units, - byref(_error_status), - _error_message, + byref(self.error_status_c), # OUT <- error status code + self.error_message_c # OUT <- error message buffer ) - if self.fatal_error(_error_status): - raise RuntimeError(f"Error {_error_status.value}: {_error_message.value}") + self.check_error() + + # Initialize output channels + self.numChannels = self._numChannels.value # if len(_channel_names.value.split()) == 0: # self.output_channel_names = [] @@ -150,38 +299,255 @@ def init( self.output_channel_units = [n.decode('UTF-8') for n in _channel_units.value.split()] # Allocate the data for the outputs - self.output_values = np.zeros( self.num_outs_c.value, dtype=c_float, order='C' ) + self.output_values = np.zeros( self._numChannels.value, dtype=c_float, order='C' ) - def calc_output(self, t): - _error_status = c_int(0) - _error_message = create_string_buffer(self.ERROR_MSG_C_LEN) + + def seastate_calcOutput(self, time: float, output_channel_values: npt.NDArray[np.float32]) -> None: + """Calculate output values at the given time. + + Args: + time: Current simulation time + output_channel_values: Array to store calculated output values + + Raises: + ValueError: If output_channel_values array has wrong size + RuntimeError: If calculation fails + """ + if output_channel_values.size != self.numChannels: + raise ValueError( + f"Output array must have size {self.numChannels}, " + f"got {output_channel_values.size}" + ) + + output_channel_values_c = (c_float * self.numChannels)(0.) self.SeaSt_C_CalcOutput( - byref(c_double(t)), # IN: time + byref(c_double(time)), # IN -> current simulation time self.output_values.ctypes.data_as(POINTER(c_float)), # OUT: output channel values - byref(_error_status), # OUT: ErrStat_C - _error_message # OUT: ErrMsg_C + byref(self.error_status_c), # OUT <- error status + self.error_message_c # OUT <- error message ) + self.check_error() - if self.fatal_error(_error_status): - self.end() - raise RuntimeError(f"Error {_error_status.value}: {_error_message.value}") + # Copy results back to numpy array + output_channel_values[:] = np.reshape(self.output_values, (self.numChannels)) - def end(self): - _error_status = c_int(0) - _error_message = create_string_buffer(self.ERROR_MSG_C_LEN) + def seastate_end(self): if not self.ended: self.ended = True self.SeaSt_C_End( - byref(_error_status), - _error_message, + byref(self.error_status_c), # OUT <- error status code + self.error_message_c # OUT <- error message buffer ) - if self.fatal_error(_error_status): - raise RuntimeError(f"Error {_error_status.value}: {_error_message.value}") + self.check_error() + + + def seastate_getWaveFieldPointer(self,ss_pointer: c_void_p) -> None: + self.SeaSt_C_GetWaveFieldPointer( + byref(ss_pointer), # IN -> pointer to the WaveField data + byref(self.error_status_c), # OUT <- error status code + self.error_message_c # OUT <- error message buffer + ) + self.check_error() + + def seastate_setWaveFieldPointer(self,ss_pointer: c_void_p) -> None: + self.SeaSt_C_SetWaveFieldPointer( + byref(ss_pointer), # IN -> pointer to the WaveField data + byref(self.error_status_c), # OUT <- error status code + self.error_message_c # OUT <- error message buffer + ) + self.check_error() + + + def get_fluidVelAcc(self, + time: float, + position: npt.NDArray[np.float32], + vel: npt.NDArray[np.float32], + acc: npt.NDArray[np.float32], + nodeInWater: int, + ) -> None: + """ + Get fluid velocity, acceleration, and if node is in water values at the given time. + Args: + time: Current simulation time + position: position in 3D to get info from + vel: velocity at position + acc: acceleration at position + nodeInWater: 1 if position is in the water, 0 if not. Note that + this is relative to SWL unless stretching is used + Raises: + RuntimeError: If calculation fails + """ + # I don't know why I have to convert the position, but I get garbage + # across the inteface if I don't (IANAPP: I am not a python programmer) + pos = np.zeros( 3, dtype=c_float ) + pos[0] = position[0] + pos[1] = position[1] + pos[2] = position[2] + vel = np.zeros( 3, dtype=c_float ) + acc = np.zeros( 3, dtype=c_float ) + self.SeaSt_C_GetFluidVelAcc( + byref(c_double(time)), # IN -> current simulation time + pos.ctypes.data_as(POINTER(c_float)), # IN -> position (3 vector) + vel.ctypes.data_as(POINTER(c_float)), # OUT <- velocity (3 vector) + acc.ctypes.data_as(POINTER(c_float)), # OUT <- acceleration (3 vector) + nodeInWater.ctypes.data_as(POINTER(c_int)), # OUT <- node is in water (0=false, 1=true) + byref(self.error_status_c), # OUT <- error status + self.error_message_c # OUT <- error message + ) + self.check_error() + return vel,acc,nodeInWater + + + def get_surfElev(self, + time: float, + position: npt.NDArray[np.float32], + elev: float, + ) -> None: + """ + Get the surface elevation at an X,Y point. + Args: + time: Current simulation time + position: position in 2D to get info from (3D could be passed in) + elev: elevation in meters + Raises: + RuntimeError: If calculation fails + """ + # I don't know why I have to convert the position, but I get garbage + # across the inteface if I don't (IANAPP: I am not a python programmer) + pos = np.array(position).astype(c_float)[:3] + elev_c = c_float(0.0) + self.SeaSt_C_GetSurfElev( + byref(c_double(time)), # IN -> current simulation time + pos.ctypes.data_as(POINTER(c_float)), # IN -> position (3 vector) + elev_c, # OUT <- total wave elevation + byref(self.error_status_c), # OUT <- error status + self.error_message_c # OUT <- error message + ) + self.check_error() + elev = elev_c.value + return elev + + + def get_surfNorm(self, + time: float, + position: npt.NDArray[np.float32], + norm: npt.NDArray[np.float32], + ) -> None: + """ + Get the normal to the surface at an X,Y point. + Args: + time: Current simulation time + position: position in 2D to get info from (3D could be passed in) + norm: normal vector + Raises: + RuntimeError: If calculation fails + """ + # I don't know why I have to convert the position, but I get garbage + # across the inteface if I don't (IANAPP: I am not a python programmer) + pos = np.zeros( 2, dtype=c_float ) + pos[0] = position[0] + pos[1] = position[1] + norm = np.zeros( 3, dtype=c_float ) + self.SeaSt_C_GetSurfNorm( + byref(c_double(time)), # IN -> current simulation time + pos.ctypes.data_as(POINTER(c_float)), # IN -> position (3 vector) + norm.ctypes.data_as(POINTER(c_float)), # OUT <- normal vector to surface + byref(self.error_status_c), # OUT <- error status + self.error_message_c # OUT <- error message + ) + self.check_error() + return norm + + + + @property def num_outs(self): - return self.num_outs_c.value \ No newline at end of file + return self._numChannels.value + + +#=============================================================================== +# Helper classes for writing output 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 SS. + +class ResultsOut(): + """ + This is only for testing purposes. Since we are not returning the + velocities to anything, we will write them to file as we go for + comparison in the regression test. When coupled to another code, the + velocities array would be passed back to the calling code for use in + the aerodynamic solver. + """ + def __init__(self, filename, NumResPts): + + self.results_file = open(filename, 'w') # open output file and write header info + + # write file header + t_string=datetime.datetime.now() + dt_string=datetime.date.today() + self.results_file.write(f"## This file was generated by SeaState called from Python on {dt_string.strftime('%b-%d-%Y')} at {t_string.strftime('%H:%M:%S')}{os.linesep}") + self.results_file.write(f"## This file contains outputs from calls to SeaState routines (not the CalcOutput) at the {NumResPts} points specified in the file {filename}{os.linesep}") + self.results_file.write(f"# {os.linesep}") + self.results_file.write(f"# {os.linesep}") + self.results_file.write(f"# {os.linesep}") + self.results_file.write(f"# {os.linesep}") + self.results_file.write(f" T x y z V_x V_y V_z A_x A_Y A_Z nodeInWater elev norm_x norm_y norm_z{os.linesep}") + self.results_file.write(f" (s) (m) (m) (m) (m/s) (m/s) (m/s) (m/s) (m/s) (m/s) (-) (m) (m/s) (m/s) (m/s) {os.linesep}") + self.opened = True + + def write(self,t,p,v,a,nodeInWater,elev,n): + self.results_file.write(' %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11d %11.3f %11.3f %11.3f %11.3f\n' % (t,p[0],p[1],p[2],v[0],v[1],v[2],a[0],a[1],a[2],nodeInWater,elev,n[0],n[1],n[2])) + + def end(self): + if self.opened: + self.results_file.close() + self.opened = False + + + + +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 SeaState c-bindings library on {dt_string.strftime('%b-%d-%Y')} at {t_string.strftime('%H:%M:%S')}\n") + self.OutFile.write(f"## This file contains output channels requested from the OutList section of the input file") + self.OutFile.write(f"{filename}\n") + self.OutFile.write("#\n") + self.OutFile.write("#\n") + self.OutFile.write("#\n") + self.OutFile.write("#\n") + l = len(chan_names) + f_string = "{:^15s}"+" {:^20s} "*(l-1) + self.OutFile.write(f_string.format(*chan_names) + '\n') + self.OutFile.write(f_string.format(*chan_units) + '\n') + self.opened = True + + def write(self,chan_data): + l = chan_data.shape[1] + f_string = "{:10.4f}"+"{:25.7f}"*(l-1) + for i in range(0,chan_data.shape[0]): + self.OutFile.write(f_string.format(*chan_data[i,:]) + '\n') + #if i==0: + # print(f"{chan_data[i,:]}") + + def end(self): + if self.opened: + self.OutFile.close() + self.opened = False diff --git a/modules/hydrodyn/src/HydroDyn_C_Binding.f90 b/modules/hydrodyn/src/HydroDyn_C_Binding.f90 index 4c7c5bc4dd..e3f753faf0 100644 --- a/modules/hydrodyn/src/HydroDyn_C_Binding.f90 +++ b/modules/hydrodyn/src/HydroDyn_C_Binding.f90 @@ -31,6 +31,7 @@ MODULE HydroDyn_C_BINDING PUBLIC :: HydroDyn_C_Init PUBLIC :: HydroDyn_C_CalcOutput + PUBLIC :: HydroDyn_C_CalcOutput_and_AddedMass PUBLIC :: HydroDyn_C_UpdateStates PUBLIC :: HydroDyn_C_End @@ -304,6 +305,7 @@ SUBROUTINE HydroDyn_C_Init( 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 + ! structural mesh reference position tmpNodePos(1:6,1:NumNodePts) = reshape( real(InitNodePositions_C(1:6*NumNodePts),ReKi), (/6,NumNodePts/) ) !---------------------------------------------------- @@ -441,6 +443,10 @@ SUBROUTINE HydroDyn_C_Init( HD%InitInp%Gravity = REAL(Gravity_C, ReKi) HD%InitInp%TMax = REAL(TMax_C, DbKi) +!FIXME: initial platform position does not work!!! + ! Initial platform position +! HD%InitInp%PlatformPos = (/ REAL(PtfmRefPtPositionX_C, ReKi), REAL(PtfmRefPtPositionX_C, ReKi), 0 /) + ! Transfer data from SeaState ! Need to set up other module's InitInput data here because we will also need to clean up SeaState data and would rather not defer that cleanup HD%InitInp%InvalidWithSSExctn = SeaSt%InitOutData%InvalidWithSSExctn @@ -500,6 +506,7 @@ SUBROUTINE HydroDyn_C_Init( !-------------------------------------------------------------------------------------------------------------------------------- ! Set the interface meshes and outputs + ! -- uses the InitNodePositions_C location/orientation to set the structural mesh reference location !-------------------------------------------------------------------------------------------------------------------------------- call SetMotionLoadsInterfaceMeshes(ErrStat2,ErrMsg2); if (Failed()) return @@ -842,8 +849,11 @@ logical function Failed() end function Failed END SUBROUTINE HydroDyn_C_CalcOutput + !=============================================================================================================== !-------------------------------------- HydroDyn CalcOutput_and_AddedMass -------------------------------------- +!> This routine is similar to the HydroDyn_C_CalcOutput, but splits the forces returned from HydroDyn_CalcOutput +!! into the hydrodynamic forces without added mass, and a separate added mass matrix. !=============================================================================================================== SUBROUTINE HydroDyn_C_CalcOutput_and_AddedMass(Time_C, NumNodePts_C, NodePos_C, NodeVel_C, & diff --git a/modules/hydrodyn/src/Morison.f90 b/modules/hydrodyn/src/Morison.f90 index fb7003e210..143533f64c 100644 --- a/modules/hydrodyn/src/Morison.f90 +++ b/modules/hydrodyn/src/Morison.f90 @@ -5967,6 +5967,7 @@ SUBROUTINE GetDistDrag_Rec(Time,mem,i,dSadl_p,dSbdl_p,f_hydro,ErrStat,ErrMsg,Sub Real(ReKi) :: posFC(3) ! Position of face center Real(ReKi) :: SVFC(3) ! Structure velocity at face center Real(SiKi) :: FVFC(3) ! Flow velocity at face center + Real(SiKi) :: FAFC(3) ! Flow acceleration at face center Real(ReKi) :: vrelFC ! Relative (flow-structure) velocity at face centers Real(ReKi) :: vrelFCf ! High-pass filtered relative (flow-structure) velocity at face centers Real(ReKi) :: STV(3) ! Structure translational velocity at the current node/free-surface intersection @@ -6081,7 +6082,7 @@ SUBROUTINE GetDistDrag_Rec(Time,mem,i,dSadl_p,dSbdl_p,f_hydro,ErrStat,ErrMsg,Sub SVFC = STV + cross_product( SRV, rToFC(1:3,fNo) ) ! Compute fluid velocity at face center - Call WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., tmpNodeInWater, FVFC, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Call WaveField_GetNodeWaveVelAcc( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., tmpNodeInWater, FVFC, FAFC, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ! Note: We force each face center to also be wetted if the center node is wetted. Otherwise, the load-smoothing procedure might not work ! Compute the face-normal component of the relative fluid velocity (fluid-structure) at face center @@ -6198,7 +6199,7 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat INTEGER(IntKi) :: nodeInWater, tmpInt REAL(ReKi) :: pos(3), vrel(3), FV(3), vmag, vmagf, An_End(3) REAL(ReKi) :: posFC(3), SVFC(3), vrelFC, vrelFCf - REAL(SiKi) :: FVTmp(3) + REAL(SiKi) :: FVTmp(3),FATmp(3) TYPE(Morison_MemberType) :: mem !< Current member INTEGER(IntKi) :: errStat2 @@ -6219,7 +6220,7 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat pos = m%DispNodePosHdn(:,J) ! Get fluid velocity at the joint - CALL WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, pos, .FALSE., nodeInWater, FVTmp, ErrStat2, ErrMsg2 ) + CALL WaveField_GetNodeWaveVelAcc( p%WaveField, m%WaveField_m, Time, pos, .FALSE., nodeInWater, FVTmp, FATmp, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FV = REAL(FVTmp, ReKi) vrel = ( FV - u%Mesh%TranslationVel(:,J) ) * nodeInWater @@ -6249,7 +6250,7 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat DO I = mem%i_floor+1, N+1 pos = m%DispNodePosHdn(:, mem%NodeIndx(I)) - CALL WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, pos, .FALSE., nodeInWater, FVTmp, ErrStat2, ErrMsg2 ) + CALL WaveField_GetNodeWaveVelAcc( p%WaveField, m%WaveField_m, Time, pos, .FALSE., nodeInWater, FVTmp, FATmp, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) IF (nodeInWater .EQ. 1_IntKi) THEN @@ -6258,7 +6259,7 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat ! Side B - +x_hat side posFC = pos + mem%x_hat * 0.5 * mem%SaMG(i) SVFC = u%Mesh%TranslationVel(:,mem%NodeIndx(i)) + cross_product( u%Mesh%RotationVel(:,mem%NodeIndx(i)), mem%x_hat * 0.5 * mem%SaMG(i) ) - call WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call WaveField_GetNodeWaveVelAcc( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., tmpInt, FVTmp, FATmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) vrelFC = dot_product( REAL(FVTmp,ReKi) - SVFC, mem%x_hat ) vrelFCf = mem%VRelNFiltConstB * ( vrelFC + xd%MV_rel_n_FiltStat(1,mem%NodeIndx(I)) ) xd%MV_rel_n_FiltStat(1,mem%NodeIndx(I)) = vrelFCf - vrelFC @@ -6266,7 +6267,7 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat ! Side B - -x_hat side posFC = pos - mem%x_hat * 0.5 * mem%SaMG(i) SVFC = u%Mesh%TranslationVel(:,mem%NodeIndx(i)) + cross_product( u%Mesh%RotationVel(:,mem%NodeIndx(i)), -mem%x_hat * 0.5 * mem%SaMG(i) ) - call WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call WaveField_GetNodeWaveVelAcc( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., tmpInt, FVTmp, FATmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) vrelFC = dot_product( REAL(FVTmp,ReKi) - SVFC, -mem%x_hat ) vrelFCf = mem%VRelNFiltConstB * ( vrelFC + xd%MV_rel_n_FiltStat(2,mem%NodeIndx(I)) ) xd%MV_rel_n_FiltStat(2,mem%NodeIndx(I)) = vrelFCf - vrelFC @@ -6274,7 +6275,7 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat ! Side A - +y_hat side posFC = pos + mem%y_hat * 0.5 * mem%SbMG(i) SVFC = u%Mesh%TranslationVel(:,mem%NodeIndx(i)) + cross_product( u%Mesh%RotationVel(:,mem%NodeIndx(i)), mem%y_hat * 0.5 * mem%SbMG(i) ) - call WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call WaveField_GetNodeWaveVelAcc( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., tmpInt, FVTmp, FATmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) vrelFC = dot_product( REAL(FVTmp,ReKi) - SVFC, mem%y_hat ) vrelFCf = mem%VRelNFiltConstA * ( vrelFC + xd%MV_rel_n_FiltStat(3,mem%NodeIndx(I)) ) xd%MV_rel_n_FiltStat(3,mem%NodeIndx(I)) = vrelFCf - vrelFC @@ -6282,7 +6283,7 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat ! Side A - -y_hat side posFC = pos - mem%y_hat * 0.5 * mem%SbMG(i) SVFC = u%Mesh%TranslationVel(:,mem%NodeIndx(i)) + cross_product( u%Mesh%RotationVel(:,mem%NodeIndx(i)), -mem%y_hat * 0.5 * mem%SbMG(i) ) - call WaveField_GetNodeWaveVel( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., tmpInt, FVTmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call WaveField_GetNodeWaveVelAcc( p%WaveField, m%WaveField_m, Time, posFC, .TRUE., tmpInt, FVTmp, FATmp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) vrelFC = dot_product( REAL(FVTmp,ReKi) - SVFC, -mem%y_hat ) vrelFCf = mem%VRelNFiltConstA * ( vrelFC + xd%MV_rel_n_FiltStat(4,mem%NodeIndx(I)) ) xd%MV_rel_n_FiltStat(4,mem%NodeIndx(I)) = vrelFCf - vrelFC diff --git a/modules/seastate/src/SeaSt_WaveField.f90 b/modules/seastate/src/SeaSt_WaveField.f90 index 0c1fb951e2..076635883c 100644 --- a/modules/seastate/src/SeaSt_WaveField.f90 +++ b/modules/seastate/src/SeaSt_WaveField.f90 @@ -12,7 +12,7 @@ MODULE SeaSt_WaveField PUBLIC WaveField_GetNodeTotalWaveElev PUBLIC WaveField_GetNodeWaveNormal PUBLIC WaveField_GetNodeWaveKin -PUBLIC WaveField_GetNodeWaveVel +PUBLIC WaveField_GetNodeWaveVelAcc PUBLIC WaveField_GetWaveKin @@ -129,7 +129,7 @@ SUBROUTINE WaveField_GetNodeWaveNormal( WaveField, WaveField_m, Time, pos, r, n, ErrStat = ErrID_None ErrMsg = "" - r1 = MAX(r,real(1.0e-6,ReKi)) ! In case r is zero + r1 = MAX(r,real(5.0e-3,ReKi)) ! In case r is zero ZetaP = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, (/pos(1)+r1,pos(2)/), ErrStat2, ErrMsg2 ); if (Failed()) return; ZetaM = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, (/pos(1)-r1,pos(2)/), ErrStat2, ErrMsg2 ); if (Failed()) return; @@ -286,7 +286,7 @@ END SUBROUTINE WaveField_GetNodeWaveKin !-------------------- Subroutine for wave field velocity only --------------------! -SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNodeInWater, nodeInWater, FV, ErrStat, ErrMsg ) +SUBROUTINE WaveField_GetNodeWaveVelAcc( WaveField, WaveField_m, Time, pos, forceNodeInWater, nodeInWater, FV, FA, ErrStat, ErrMsg ) type(SeaSt_WaveFieldType), intent(in ) :: WaveField type(SeaSt_WaveField_MiscVarType), intent(inout) :: WaveField_m real(DbKi), intent(in ) :: Time @@ -294,6 +294,7 @@ SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNod logical, intent(in ) :: forceNodeInWater integer(IntKi), intent( out) :: nodeInWater real(SiKi), intent( out) :: FV(3) + real(SiKi), intent( out) :: FA(3) integer(IntKi), intent( out) :: ErrStat ! Error status of the operation character(*), intent( out) :: ErrMsg ! Error message if errStat /= ErrID_None @@ -319,9 +320,11 @@ SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNod ! Use location to obtain interpolated values of kinematics CALL WaveField_Interp_Setup4D( Time, pos, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; FV(:) = WaveField_Interp_4D_Vec( WaveField%WaveVel, WaveField_m ) + FA(:) = WaveField_Interp_4D_Vec( WaveField%WaveAcc, WaveField_m ) ELSE ! Node is above the SWL nodeInWater = 0_IntKi FV(:) = 0.0 + FA(:) = 0.0 END IF ELSE ! Wave stretching enabled @@ -337,17 +340,20 @@ SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNod ! Use location to obtain interpolated values of kinematics CALL WaveField_Interp_Setup4D( Time, pos, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; FV(:) = WaveField_Interp_4D_Vec( WaveField%WaveVel, WaveField_m ) + FA(:) = WaveField_Interp_4D_Vec( WaveField%WaveAcc, WaveField_m ) ELSE ! Node is above SWL - need wave stretching ! Vertical wave stretching CALL WaveField_Interp_Setup4D( Time, posXY0, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; FV(:) = WaveField_Interp_4D_vec( WaveField%WaveVel, WaveField_m ) + FA(:) = WaveField_Interp_4D_vec( WaveField%WaveAcc, WaveField_m ) ! Extrapoled wave stretching IF (WaveField%WaveStMod == 2) THEN CALL WaveField_Interp_Setup3D( Time, posXY, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; FV(:) = FV(:) + WaveField_Interp_3D_vec( WaveField%PWaveVel0, WaveField_m ) * pos(3) + FA(:) = FA(:) + WaveField_Interp_3D_vec( WaveField%PWaveAcc0, WaveField_m ) * pos(3) END IF END IF ! Node is submerged @@ -362,6 +368,7 @@ SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNod ! Obtain the wave-field variables by interpolation with the mapped position. CALL WaveField_Interp_Setup4D( Time, posPrime, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; FV(:) = WaveField_Interp_4D_Vec( WaveField%WaveVel, WaveField_m ) + FA(:) = WaveField_Interp_4D_Vec( WaveField%WaveAcc, WaveField_m ) END IF @@ -369,6 +376,7 @@ SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNod nodeInWater = 0_IntKi FV(:) = 0.0 + FA(:) = 0.0 END IF ! If node is in or out of water @@ -379,7 +387,7 @@ logical function Failed() call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) Failed = ErrStat >= AbortErrLev end function -END SUBROUTINE WaveField_GetNodeWaveVel +END SUBROUTINE WaveField_GetNodeWaveVelAcc SUBROUTINE WaveField_GetWaveKin( WaveField, WaveField_m, Time, pos, forceNodeInWater, nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat, ErrMsg ) diff --git a/modules/seastate/src/SeaState_C_Binding.f90 b/modules/seastate/src/SeaState_C_Binding.f90 index e7910afc58..260dc54273 100644 --- a/modules/seastate/src/SeaState_C_Binding.f90 +++ b/modules/seastate/src/SeaState_C_Binding.f90 @@ -19,291 +19,814 @@ !********************************************************************************************************************************** MODULE SeaState_C_Binding - USE ISO_C_BINDING - USE SeaState - USE SeaState_Types - USE SeaState_Output - USE NWTC_Library - USE NWTC_C_Binding, ONLY: ErrMsgLen_C, IntfStrLen, SetErrStat_F2C, FileNameFromCString - USE VersionInfo - - IMPLICIT NONE - SAVE - - PUBLIC :: SeaSt_C_Init - PUBLIC :: SeaSt_C_CalcOutput - PUBLIC :: SeaSt_C_End - - !------------------------------------------------------------------------------------ - ! Version info for display - TYPE(ProgDesc), PARAMETER :: version = SeaSt_ProgDesc - - !------------------------------------------------------------------------------------ - ! Debugging: DebugLevel -- passed at PreInit - ! 0 - none - ! 1 - some summary info - ! 2 - above + all position/orientation info - ! 3 - above + input files (if direct passed) - ! 4 - above + meshes - INTEGER(IntKi) :: DebugLevel = 0 - - !------------------------------ - ! Primary derived types - TYPE(SeaSt_InputType) :: InputData !< Inputs to SeaState - TYPE(SeaSt_InitInputType) :: InitInp - TYPE(SeaSt_InitOutputType) :: InitOutData !< Initial output data -- Names, units, and version info. - TYPE(SeaSt_ParameterType) :: p !< Parameters - TYPE(SeaSt_OutputType) :: y !< Initial output (outputs are not calculated; only the output mesh is initialized) - TYPE(SeaSt_MiscVarType) :: m !< Misc variables for optimization (not copied in glue code) - -CONTAINS - - -SUBROUTINE SeaSt_C_Init(InputFile_C, OutRootName_C, Gravity_C, WtrDens_C, WtrDpth_C, MSL2SWL_C, NSteps_C, TimeInterval_C, WaveElevSeriesFlag_C, WrWvKinMod_C, NumChannels_C, OutputChannelNames_C, OutputChannelUnits_C, ErrStat_C, ErrMsg_C) BIND (C, NAME='SeaSt_C_Init') -IMPLICIT NONE + USE ISO_C_BINDING + USE SeaSt_WaveField + USE SeaState + USE SeaState_Types + USE SeaState_Output + USE NWTC_Library + USE NWTC_C_Binding, ONLY: ErrMsgLen_C, IntfStrLen, SetErrStat_F2C, FileNameFromCString + USE VersionInfo + + implicit none + save + + PUBLIC :: SeaSt_C_PreInit + PUBLIC :: SeaSt_C_Init + PUBLIC :: SeaSt_C_CalcOutput + PUBLIC :: SeaSt_C_End + PUBLIC :: SeaSt_C_GetWaveFieldPointer + PUBLIC :: SeaSt_C_SetWaveFieldPointer + PUBLIC :: SeaSt_C_GetFluidVelAcc + PUBLIC :: SeaSt_C_GetSurfElev + PUBLIC :: SeaSt_C_GetSurfNorm + + !------------------------------------------------------------------------------------ + ! Debugging: DebugLevel -- passed at PreInit + ! 0 - none + ! 1 - some summary info + ! 2 - above + all position/orientation info + ! 3 - above + input files (if direct passed) + ! 4 - above + meshes + integer(IntKi) :: DebugLevel = 4 + logical :: PreInitDone = .false. + + !------------------------------------------------------------------------------------ + ! Visualization + type VTKvis + character(1024) :: outdir + character(1024) :: OutRootName ! includes directory + integer(IntKi) :: write ! 0 off, 1 init, 2 animate + real(DbKi) :: dt + integer(IntKi) :: NWaveElevPts(2) ! number of points in x/y directions + real(SiKi), allocatable :: WaveElevVisX(:),WaveElevVisY(:) ! x, y locations of points + real(SiKi), allocatable :: WaveElevVisGrid(:,:,:) ! the actual surface data for full time series + integer(IntKi) :: tWidth = 5 ! Should calculate this, but not going to + integer(IntKi) :: LastWaveIndx + integer(IntKi) :: lastCount = -1 + end type VTKvis + type(VTKvis) :: vtk + + !------------------------------ + ! Primary derived types + type(SeaSt_InputType) :: u !< inputs to SS + type(SeaSt_InitInputType) :: InitInp !< initialization input + type(SeaSt_InitOutputType) :: InitOutData !< Initial output data + type(SeaSt_ParameterType), target :: p !< Parameters + type(SeaSt_OutputType) :: y !< Initial output (outputs are not calculated; only the output mesh is initialized) + type(SeaSt_MiscVarType) :: m !< Misc variables for optimization (not copied in glue code) + type(SeaSt_ContinuousStateType) :: x !< Initial continuous states + type(SeaSt_DiscreteStateType) :: xd !< Initial discrete states + type(SeaSt_ConstraintStateType) :: z !< Initial guess of the constraint states + type(SeaSt_OtherStateType) :: OtherState !< Initial other states + +contains + + +!> Set environment variables +subroutine SeaSt_C_PreInit(Gravity_C, WtrDens_C, WtrDpth_C, MSL2SWL_C, DebugLevel_In, OutVTKDir_C, WrVTK_in, WrVTK_inDT, ErrStat_C, ErrMsg_C) BIND (C, NAME='SeaSt_C_PreInit') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_PreInit +!GCC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_PreInit +#endif + real(c_float), intent(in ) :: Gravity_C + real(c_float), intent(in ) :: WtrDens_C + real(c_float), intent(in ) :: WtrDpth_C + real(c_float), intent(in ) :: MSL2SWL_C + integer(c_int), intent(in ) :: DebugLevel_In + character(kind=c_char), intent(in ) :: OutVTKDir_C(IntfStrLen) !< Directory to put all vtk output + integer(c_int), intent(in ) :: WrVTK_in !< Write VTK outputs [0: none, 1: init only, 2: animation] + real(c_double), intent(in ) :: WrVTK_inDT !< Timestep between VTK writes + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + + character(kind=C_CHAR, len=IntfStrLen), pointer :: InputString !< Input string as a single string with NULL chracter separating lines + integer :: ErrStat, ErrStat2 + character(ErrMsgLen) :: ErrMsg, ErrMsg2 + integer :: i,j,k + character(*), parameter :: RoutineName = 'SeaSt_C_PreInit' + + ! Initialize error handling + ErrStat = ErrID_None + ErrMsg = "" + + call NWTC_Init( ProgNameIn= SeaSt_ProgDesc%Name ) + call DispCopyrightLicense( SeaSt_ProgDesc%Name ) + call DispCompileRuntimeInfo( SeaSt_ProgDesc%Name ) + + ! interface debugging + DebugLevel = int(DebugLevel_in,IntKi) + + ! check valid debug level, show passed data if >0 + if (DebugLevel < 0_IntKi) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "Interface debug level must be 0 or greater"//NewLine// & + " 0 - none"//NewLine// & + " 1 - some summary info and variables passed through interface (init only)"//NewLine// & + " 2 - above + all position info on all calls" + call ShowPassedData() + if (Failed()) return; + elseif (DebugLevel > 0_IntKi) THEN + call WrScr(" Interface debugging level "//trim(Num2Lstr(DebugLevel))//" requested.") + call ShowPassedData() + endif + + ! clear memory of anything we allocate locally + call ClearMem() ! ignoring any error handling from this + + ! store environment values + InitInp%Gravity = Gravity_C + InitInp%defWtrDens = WtrDens_C + InitInp%defWtrDpth = WtrDpth_C + InitInp%defMSL2SWL = MSL2SWL_C + + !---------------------- + ! store VTK output info + vtk%write = int(WrVTK_in, IntKi) + vtk%dt = real(WrVTK_inDT, DbKi) + + if (vtk%write < 0_IntKi .or. vtk%write > 2_IntKi) then + ErrStat2 = ErrID_Warn + ErrMSg2 = "WrVTK_in must be 0 (off), 1 (init), 2 (animation), but "//trim(Num2LStr(vtk%write))//" was passed. Turning off VTK surface export." + vtk%write = 0_IntKi + if (Failed()) return + endif + + if (vtk%write > 0_IntKi) then + ! Store the out root dir + vtk%outdir = TRANSFER( OutVTKDir_C, vtk%outdir ) + i = INDEX(vtk%outdir,C_NULL_CHAR) - 1 ! if this has a c null character at the end... + if ( i > 0 ) vtk%outdir = vtk%outdir(1:I) ! remove it + ! Tell SeaState to generate the visualization using default grid + InitInp%SurfaceVis = .true. + InitInp%SurfaceVisNx = 0 ! use the WaveField grid resolution + InitInp%SurfaceVisNy = 0 ! use the WaveField grid resolution + endif + + + ! If we got this far, we are initialized + PreInitDone = .true. + + call Cleanup() + return +contains + logical function Failed() + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + if (Failed) call Cleanup() + end function Failed + subroutine Cleanup() ! NOTE: we are ignoring any error reporting from here + call SetErrStat_F2C(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + end subroutine Cleanup + subroutine ShowPassedData() + ! character(1) :: TmpFlag + ! integer :: i,j + call WrScr("-----------------------------------------------------------") + call WrScr("Interface debugging: SeaSt_C_PreInit") + call WrScr(" --------------------------------------------------------") + call WrScr(" Gravity_C -> "//trim(Num2LStr(Gravity_C))) + call WrScr(" WtrDens_C -> "//trim(Num2LStr(WtrDens_C))) + call WrScr(" WtrDpth_C -> "//trim(Num2LStr(WtrDpth_C))) + call WrScr(" MSL2SWL_C -> "//trim(Num2LStr(MSL2SWL_C))) + call WrScr(" DebugLevel_In -> "//trim(Num2LStr(DebugLevel_In))) + call WrScr(" OutVTKDir_C (ptr addr) -> "//trim(Num2LStr(LOC(OutVTKDir_C)))) + call WrScr(" WrVTK_in -> "//trim(Num2LStr(WrVTK_in))) + call WrScr(" WrVTK_inDT -> "//trim(Num2LStr(WrVTK_inDT))) + call WrScr("-----------------------------------------------------------") + end subroutine ShowPassedData +end subroutine SeaSt_C_PreInit + + +!> Initialize the library (PreInit must be called first) +subroutine SeaSt_C_Init(InputFile_C, OutRootName_C, NSteps_C, TimeInterval_C, NumChannels_C, OutputChannelNames_C, OutputChannelUnits_C, ErrStat_C, ErrMsg_C) BIND (C, NAME='SeaSt_C_Init') #ifndef IMPLICIT_DLLEXPORT !DEC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_Init !GCC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_Init #endif - TYPE(C_PTR), INTENT(IN ) :: InputFile_C - TYPE(C_PTR), INTENT(IN ) :: OutRootName_C - REAL(C_FLOAT), INTENT(IN ) :: Gravity_C - REAL(C_FLOAT), INTENT(IN ) :: WtrDens_C - REAL(C_FLOAT), INTENT(IN ) :: WtrDpth_C - REAL(C_FLOAT), INTENT(IN ) :: MSL2SWL_C - INTEGER(C_INT), INTENT(IN ) :: NSteps_C - REAL(C_FLOAT), INTENT(IN ) :: TimeInterval_C - INTEGER(C_INT), INTENT(IN ) :: WaveElevSeriesFlag_C - INTEGER(C_INT), INTENT(IN ) :: WrWvKinMod_C - INTEGER(C_INT), INTENT( OUT) :: NumChannels_C - CHARACTER(KIND=C_CHAR), INTENT( OUT) :: OutputChannelNames_C(ChanLen*MaxOutPts+1) - CHARACTER(KIND=C_CHAR), INTENT( OUT) :: OutputChannelUnits_C(ChanLen*MaxOutPts+1) - INTEGER(C_INT), INTENT( OUT) :: ErrStat_C - CHARACTER(KIND=C_CHAR), INTENT( OUT) :: ErrMsg_C(ErrMsgLen_C) - - ! Local variables - CHARACTER(KIND=C_CHAR, len=IntfStrLen), POINTER :: InputFileString !< Input file as a single string with NULL chracter separating lines - CHARACTER(KIND=C_CHAR, len=IntfStrLen), POINTER :: OutputFileString !< Input file as a single string with NULL chracter separating lines - CHARACTER(IntfStrLen) :: InputFileName - CHARACTER(IntfStrLen) :: OutRootName - TYPE(SeaSt_InputType) :: u !< An initial guess for the input; input mesh must be defined - TYPE(SeaSt_ContinuousStateType) :: x !< Initial continuous states - TYPE(SeaSt_DiscreteStateType) :: xd !< Initial discrete states - TYPE(SeaSt_ConstraintStateType) :: z !< Initial guess of the constraint states - TYPE(SeaSt_OtherStateType) :: OtherState !< Initial other states - REAL(DbKi) :: Interval !< Coupling interval in seconds: the rate that - !! (1) SeaSt_UpdateStates() is called in loose coupling & - !! (2) SeaSt_UpdateDiscState() is called in tight coupling. - !! Input is the suggested time from the glue code; - !! Output is the actual coupling interval that will be used - !! by the glue code. - - INTEGER :: ErrStat_F !< aggregated error status - CHARACTER(ErrMsgLen) :: ErrMsg_F !< aggregated error message - INTEGER :: ErrStat_F2 !< temporary error status from a call - CHARACTER(ErrMsgLen) :: ErrMsg_F2 !< temporary error message from a call - INTEGER :: i,j,k - CHARACTER(*), PARAMETER :: RoutineName = 'SeaSt_C_Init' !< for error handling - - ! Initialize error handling - ErrStat_F = ErrID_None - ErrMsg_F = "" - - CALL NWTC_Init( ProgNameIn=version%Name ) - CALL DispCopyrightLicense( version%Name ) - CALL DispCompileRuntimeInfo( version%Name ) - - ! interface debugging - ! DebugLevel = int(DebugLevel_in,IntKi) - - ! Input files - CALL C_F_POINTER(InputFile_C, InputFileString) ! Get a pointer to the input file string - InputFileName = FileNameFromCString(InputFileString, IntfStrLen) ! convert the input file name from c_char to fortran character - - CALL C_F_POINTER(OutRootName_C, OutputFileString) ! Get a pointer to the input file string - OutRootName = FileNameFromCString(OutputFileString, IntfStrLen) ! convert the input file name from c_char to fortran character - - ! if non-zero, show all passed data here. Then check valid values - IF (DebugLevel /= 0_IntKi) THEN - CALL WrScr(" Interface debugging level "//trim(Num2Lstr(DebugLevel))//" requested.") - CALL ShowPassedData() - ENDIF - ! check valid debug level - IF (DebugLevel < 0_IntKi) THEN - ErrStat_F2 = ErrID_Fatal - ErrMsg_F2 = "Interface debug level must be 0 or greater"//NewLine// & - " 0 - none"//NewLine// & - " 1 - some summary info and variables passed through interface"//NewLine// & - " 2 - above + all position/orientation info"//NewLine// & - " 3 - above + input files (if direct passed)"//NewLine// & - " 4 - above + meshes" - IF (Failed()) RETURN; - ENDIF - - ! For debugging the interface: - IF (DebugLevel > 0) THEN - CALL ShowPassedData() - ENDIF - - ! Set other inputs for calling SeaSt_Init - InitInp%InputFile = InputFileName - InitInp%UseInputFile = .TRUE. - InitInp%OutRootName = OutRootName - InitInp%Gravity = Gravity_C - InitInp%defWtrDens = WtrDens_C - InitInp%defWtrDpth = WtrDpth_C - InitInp%defMSL2SWL = MSL2SWL_C - InitInp%TMax = (NSteps_C - 1) * TimeInterval_C ! Using this to match the SeaState driver; could otherwise get TMax directly - InitInp%WaveFieldMod = WaveElevSeriesFlag_C - ! REAL(ReKi) :: PtfmLocationX = 0.0_ReKi !< Supplied by Driver: X coordinate of platform location in the wave field [m] - ! REAL(ReKi) :: PtfmLocationY = 0.0_ReKi !< Supplied by Driver: Y coordinate of platform location in the wave field [m] - InitInp%WrWvKinMod = WrWvKinMod_C - ! LOGICAL :: HasIce = .false. !< Supplied by Driver: Whether this simulation has ice loading (flag) [-] - ! LOGICAL :: Linearize = .FALSE. !< Flag that tells this module if the glue code wants to linearize. [-] - ! LOGICAL :: SurfaceVis = .FALSE. !< Turn on grid surface visualization outputs [-] - ! INTEGER(IntKi) :: SurfaceVisNx = 0 !< Number of points in X direction to output for visualization grid. Use 0 or negative to set to SeaState resolution. [-] - ! INTEGER(IntKi) :: SurfaceVisNy = 0 !< Number of points in Y direction to output for visualization grid. Use 0 or negative to set to SeaState resolution. [-] - - CALL SeaSt_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOutData, ErrStat_F2, ErrMsg_F2 ) - IF (Failed()) RETURN - - ! Number of channels - NumChannels_C = size(InitOutData%WriteOutputHdr) - - ! transfer the output channel names and units to c_char arrays for returning - k=1 - DO i=1,NumChannels_C - DO j=1,ChanLen ! max length of channel name. Same for units - OutputChannelNames_C(k)=InitOutData%WriteOutputHdr(i)(j:j) - OutputChannelUnits_C(k)=InitOutData%WriteOutputUnt(i)(j:j) - k=k+1 - ENDDO - ENDDO - - ! null terminate the string - OutputChannelNames_C(k) = C_NULL_CHAR - OutputChannelUnits_C(k) = C_NULL_CHAR - - CALL Cleanup() - -CONTAINS - LOGICAL FUNCTION Failed() - CALL SetErrStat( ErrStat_F2, ErrMsg_F2, ErrStat_F, ErrMsg_F, RoutineName ) - Failed = ErrStat_F >= AbortErrLev - IF (Failed) CALL Cleanup() - END FUNCTION Failed - - SUBROUTINE Cleanup() ! NOTE: we are ignoring any error reporting from here - CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) - END SUBROUTINE Cleanup - - SUBROUTINE ShowPassedData() - ! CHARACTER(1) :: TmpFlag - ! integer :: i,j - CALL WrSCr("") - CALL WrScr("-----------------------------------------------------------") - CALL WrScr("Interface debugging: Variables passed in through interface") - CALL WrScr(" SeaSt_C_Init") - CALL WrScr(" --------------------------------------------------------") - CALL WrScr(" FileInfo") - CALL WrScr("-----------------------------------------------------------") - END SUBROUTINE ShowPassedData -END SUBROUTINE SeaSt_C_Init - -SUBROUTINE SeaSt_C_CalcOutput(Time_C, OutputChannelValues_C, ErrStat_C, ErrMsg_C) BIND (C, NAME='SeaSt_C_CalcOutput') -IMPLICIT NONE + type(c_ptr), intent(in ) :: InputFile_C + type(c_ptr), intent(in ) :: OutRootName_C + integer(c_int), intent(in ) :: NSteps_C + real(c_float), intent(in ) :: TimeInterval_C + integer(c_int), intent( out) :: NumChannels_C + character(kind=c_char), intent( out) :: OutputChannelNames_C(ChanLen*MaxOutPts+1) + character(kind=c_char), intent( out) :: OutputChannelUnits_C(ChanLen*MaxOutPts+1) + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + + ! Local variables + character(kind=C_CHAR, len=IntfStrLen), pointer :: InputFileString !< Input file as a single string with NULL chracter separating lines + character(kind=C_CHAR, len=IntfStrLen), pointer :: OutputFileString !< Input file as a single string with NULL chracter separating lines + character(IntfStrLen) :: InputFileName + character(IntfStrLen) :: OutRootName + character(1024) :: vtkroot + real(DbKi) :: Interval !< DT for calling + integer :: ErrStat, ErrStat2 + character(ErrMsgLen) :: ErrMsg, ErrMsg2 + integer :: i,j,k + character(*), parameter :: RoutineName = 'SeaSt_C_Init' !< for error handling + + if (.not. PreInitDone) then + ErrStat = ErrID_Fatal + ErrMSg = "SeaSt_C_PreInit must be called before SeaSt_C_Init" + call Cleanup() + endif + + ! Initialize error handling + ErrStat = ErrID_None + ErrMsg = "" + + call NWTC_Init( ProgNameIn= SeaSt_ProgDesc%Name ) + call DispCopyrightLicense( SeaSt_ProgDesc%Name ) + call DispCompileRuntimeInfo( SeaSt_ProgDesc%Name ) + + if (DebugLevel > 0_IntKi) call ShowPassedData() + + ! Input files + call C_F_POINTER(InputFile_C, InputFileString) ! Get a pointer to the input file string + InputFileName = FileNameFromCString(InputFileString, IntfStrLen) ! convert the input file name from c_char to fortran character + + call C_F_POINTER(OutRootName_C, OutputFileString) ! Get a pointer to the input file string + OutRootName = FileNameFromCString(OutputFileString, IntfStrLen) ! convert the input file name from c_char to fortran character + + ! Set other inputs for calling SeaSt_Init + InitInp%InputFile = InputFileName + InitInp%UseInputFile = .TRUE. + InitInp%OutRootName = OutRootName + InitInp%TMax = (NSteps_C - 1) * TimeInterval_C ! Using this to match the SeaState driver; could otherwise get TMax directly + InitInp%WaveFieldMod = 0_IntKi + InitInp%WrWvKinMod = 0_IntKi + InitInp%Linearize = .false. + InitInp%hasIce = .false. + InitInp%WaveFieldMod = 0 ! does not currently support moving platform. Not really necessary though since can directly get data in absolute coords + InitInp%PtfmLocationX = 0.0_ReKi + InitInp%PtfmLocationY = 0.0_ReKi + + call SeaSt_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOutData, ErrStat2, ErrMsg2 ) + if (Failed()) return + + ! Number of channels + NumChannels_C = size(InitOutData%WriteOutputHdr) + + ! transfer the output channel names and units to c_char arrays for returning + k=1 + do i=1,NumChannels_C + do j=1,ChanLen ! max length of channel name. Same for units + OutputChannelNames_C(k)=InitOutData%WriteOutputHdr(i)(j:j) + OutputChannelUnits_C(k)=InitOutData%WriteOutputUnt(i)(j:j) + k=k+1 + enddo + enddo + + ! null terminate the string + OutputChannelNames_C(k) = C_NULL_CHAR + OutputChannelUnits_C(k) = C_NULL_CHAR + + if (vtk%write > 0_IntKi) then + ! check dt (can't check against Interval since that is never set, so just make sure it is positive) + if (vtk%dt <= 0.0) vtk%dt = 0.25 + if (allocated(InitOutData%WaveElevVisGrid)) then + vtk%NWaveElevPts(1) = size(InitOutData%WaveElevVisX) + vtk%NWaveElevPts(2) = size(InitOutData%WaveElevVisY) + call move_alloc(InitOutData%WaveElevVisX, vtk%WaveElevVisX) + call move_alloc(InitOutData%WaveElevVisY, vtk%WaveElevVisY) + call move_alloc(InitOutData%WaveElevVisGrid,vtk%WaveElevVisGrid ) + else + vtk%NWaveElevPts = 0 + vtk%write = 0 ! FIXME throw warning if we do this + endif + ! get the name of the output directory for vtk files (in a subdirectory called "vtk" of the output directory), and + ! create the VTK directory if it does not exist + call GetPath ( OutRootName, vtk%OutRootName, vtkroot ) ! the returned vtk%OutRootName includes a file separator character at the end + if (PathIsRelative(trim(vtk%OutRootName))) then + vtk%OutRootName = trim(vtk%OutRootName) // trim(vtk%outdir) + else + vtk%OutRootName = trim(vtk%outdir) + endif + call MKDIR( trim(vtk%OutRootName) ) + vtk%OutRootName = trim( vtk%OutRootName ) // PathSep // trim(vtkroot) + call WrVTK_WaveElevVisGrid (0.0_DbKi, vtk, ErrStat2, ErrMsg2) + if (Failed()) return + endif + + + + call Cleanup() + +contains + logical function Failed() + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + if (Failed) call Cleanup() + end function Failed + subroutine Cleanup() ! NOTE: we are ignoring any error reporting from here + call SetErrStat_F2C(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + end subroutine Cleanup + subroutine ShowPassedData() + character(1) :: TmpFlag + ! integer :: i,j + call WrScr("-----------------------------------------------------------") + call WrScr("Interface debugging: SeaSt_C_Init") + call WrScr(" --------------------------------------------------------") + call WrScr(" InputFile_C (ptr addr) -> "//trim(Num2LStr(LOC(InputFile_C)))) + call WrScr(" OutRootName_C (ptr addr) -> "//trim(Num2LStr(LOC(OutRootName_C)))) + call WrScr(" NSteps_C - > "//trim(Num2LStr(NSteps_C))) + call WrScr(" TimeInterval_C -> "//trim(Num2LStr(TimeInterval_C))) + call WrScr("-----------------------------------------------------------") + end subroutine ShowPassedData +end subroutine SeaSt_C_Init + +subroutine SeaSt_C_CalcOutput(Time_C, OutputChannelValues_C, ErrStat_C, ErrMsg_C) BIND (C, NAME='SeaSt_C_CalcOutput') #ifndef IMPLICIT_DLLEXPORT !DEC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_CalcOutput !GCC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_CalcOutput #endif - REAL(C_DOUBLE), INTENT(IN ) :: Time_C - REAL(C_FLOAT), INTENT( OUT) :: OutputChannelValues_C(p%NumOuts) - INTEGER(C_INT), INTENT( OUT) :: ErrStat_C - CHARACTER(KIND=C_CHAR), INTENT( OUT) :: ErrMsg_C(ErrMsgLen_C) - - ! Local variables - TYPE(SeaSt_InputType) :: u !< An initial guess for the input; input mesh must be defined - TYPE(SeaSt_ContinuousStateType) :: x !< Initial continuous states - TYPE(SeaSt_DiscreteStateType) :: xd !< Initial discrete states - TYPE(SeaSt_ConstraintStateType) :: z !< Initial guess of the constraint states - TYPE(SeaSt_OtherStateType) :: OtherState !< Initial other states - - REAL(DbKi) :: Time - INTEGER :: ErrStat_F !< aggregated error status - CHARACTER(ErrMsgLen) :: ErrMsg_F !< aggregated error message - INTEGER :: ErrStat_F2 !< temporary error status from a call - CHARACTER(ErrMsgLen) :: ErrMsg_F2 !< temporary error message from a call - CHARACTER(*), PARAMETER :: RoutineName = 'SeaSt_C_End' !< for error handling - - ! Initialize error handling - ErrStat_F = ErrID_None - ErrMsg_F = "" - - ! Convert the inputs from C to Fortran - Time = REAL(Time_C,DbKi) - - CALL SeaSt_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat_F2, ErrMsg_F2 ) - IF (Failed()) RETURN - - ! Get the output channel info out of y - OutputChannelValues_C = REAL(y%WriteOutput, C_FLOAT) - - CALL Cleanup() - -CONTAINS - LOGICAL FUNCTION Failed() - CALL SetErrStat( ErrStat_F2, ErrMsg_F2, ErrStat_F, ErrMsg_F, RoutineName ) - Failed = ErrStat_F >= AbortErrLev - IF (Failed) CALL Cleanup() - END FUNCTION Failed - - SUBROUTINE Cleanup() ! NOTE: we are ignoring any error reporting from here - CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) - END SUBROUTINE Cleanup -END SUBROUTINE - -SUBROUTINE SeaSt_C_End(ErrStat_C,ErrMsg_C) BIND (C, NAME='SeaSt_C_End') -IMPLICIT NONE + real(c_double), intent(in ) :: Time_C + real(c_float), intent( out) :: OutputChannelValues_C(p%NumOuts) + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + + ! Local variables + type(SeaSt_InputType) :: u !< An initial guess for the input; input mesh must be defined + type(SeaSt_ContinuousStateType) :: x !< Initial continuous states + type(SeaSt_DiscreteStateType) :: xd !< Initial discrete states + type(SeaSt_ConstraintStateType) :: z !< Initial guess of the constraint states + type(SeaSt_OtherStateType) :: OtherState !< Initial other states + + real(DbKi) :: Time + integer :: ErrStat, ErrStat2 + character(ErrMsgLen) :: ErrMsg, ErrMsg2 + character(*), parameter :: RoutineName = 'SeaSt_C_CalcOutput' !< for error handling + + ! Initialize error handling + ErrStat = ErrID_None + ErrMsg = "" + + ! Debugging + if (DebugLevel > 1) call ShowPassedData() + + ! Convert the inputs from C to Fortran + Time = REAL(Time_C,DbKi) + + call SeaSt_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat2, ErrMsg2 ) + if (Failed()) return + + ! Get the output channel info out of y + OutputChannelValues_C = REAL(y%WriteOutput, C_FLOAT) + + if (vtk%write > 1_IntKi) then + call WrVTK_WaveElevVisGrid (Time, vtk, ErrStat2, ErrMsg2) + if (Failed()) return + endif + + call Cleanup() + + ! Debugging + if (DebugLevel > 1) call ShowReturnData() + +contains + logical function Failed() + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + if (Failed) call Cleanup() + end function Failed + subroutine Cleanup() ! NOTE: we are ignoring any error reporting from here + CALL SetErrStat_F2C(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + end subroutine Cleanup + subroutine ShowPassedData() + call WrScr("-----------------------------------------------------------") + call WrScr("Interface debugging: SeaSt_C_CalcOutput") + call WrScr(" --------------------------------------------------------") + call WrScr(" Time_C -> "//trim(Num2LStr(Time_C))) + end subroutine ShowPassedData + subroutine ShowReturnData() + call WrScr(" OutputChannelValues_C <-") + call WrMatrix(OutputChannelValues_C,CU,'g15.6') + call WrScr("-----------------------------------------------------------") + end subroutine ShowReturnData +end subroutine + +subroutine SeaSt_C_End(ErrStat_C,ErrMsg_C) BIND (C, NAME='SeaSt_C_End') #ifndef IMPLICIT_DLLEXPORT !DEC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_End !GCC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_End #endif - INTEGER(C_INT), INTENT( OUT) :: ErrStat_C - CHARACTER(KIND=C_CHAR), INTENT( OUT) :: ErrMsg_C(ErrMsgLen_C) - - ! Local variables - TYPE(SeaSt_InputType) :: u !< An initial guess for the input; input mesh must be defined - TYPE(SeaSt_ContinuousStateType) :: x !< Initial continuous states - TYPE(SeaSt_DiscreteStateType) :: xd !< Initial discrete states - TYPE(SeaSt_ConstraintStateType) :: z !< Initial guess of the constraint states - TYPE(SeaSt_OtherStateType) :: OtherState !< Initial other states - - 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 = 'SeaSt_C_End' !< for error handling - - ! Initialize error handling - ErrStat = ErrID_None - ErrMsg = "" - CALL SeaSt_End(u, p, x, xd, z, OtherState, y, m, ErrStat2, ErrMsg2) - - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL SetErrStat_F2C( ErrStat, ErrMsg, ErrStat_C, ErrMsg_C ) - -END SUBROUTINE - -FUNCTION SeaSt_GetWaveFieldPointer_C() BIND (C, NAME='SeaSt_GetWaveFieldPointer_C') -IMPLICIT NONE + integer(C_INT), intent( out) :: ErrStat_C + character(kind=C_CHAR), intent( out) :: ErrMsg_C(ErrMsgLen_C) + + 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 = 'SeaSt_C_End' !< for error handling + ErrStat = ErrID_None + ErrMsg = "" + call SeaSt_End(u, p, x, xd, z, OtherState, y, m, ErrStat2, ErrMsg2) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call ClearMem() ! ignoring any error handling from this + call SetErrStat_F2C( ErrStat, ErrMsg, ErrStat_C, ErrMsg_C ) +end subroutine + + +!> return the pointer to the WaveField data +subroutine SeaSt_C_GetWaveFieldPointer(WaveFieldPointer_C,ErrStat_C,ErrMsg_C) BIND (C, NAME='SeaSt_C_GetWaveFieldPointer') #ifndef IMPLICIT_DLLEXPORT -!DEC$ ATTRIBUTES DLLEXPORT :: SeaSt_GetWaveFieldPointer_C -!GCC$ ATTRIBUTES DLLEXPORT :: SeaSt_GetWaveFieldPointer_C +!DEC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_GetWaveFieldPointer +!GCC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_GetWaveFieldPointer #endif - TYPE(C_PTR) :: SeaSt_GetWaveFieldPointer_C - SeaSt_GetWaveFieldPointer_C = C_LOC(p%WaveField) - RETURN -END FUNCTION - + type(c_ptr), intent( out) :: WaveFieldPointer_C + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + integer :: ErrStat + character(ErrMsgLen) :: ErrMsg + character(*), parameter :: RoutineName = 'SeaSt_C_GetWaveFieldPointer' + ErrStat = ErrID_None + ErrMSg = "" + if (associated(p%WaveField)) then + WaveFieldPointer_C = C_LOC(p%WaveField) + else + WaveFieldPointer_C = C_NULL_PTR + call SetErrStat(ErrID_Fatal,"Pointer to WaveField data not valid: data not initialized",ErrStat,ErrMsg,RoutineName) + endif + call SetErrStat_F2C( ErrStat, ErrMsg, ErrStat_C, ErrMsg_C ) + if (DebugLevel > 1) call ShowPassedData() + return +contains + subroutine ShowPassedData() + call WrScr("-----------------------------------------------------------") + call WrScr("Interface debugging: SeaSt_C_GetWaveFieldPointer") + call WrScr(" --------------------------------------------------------") + call WrScr(" WaveFieldPointer_C -> "//trim(Num2LStr(loc(p%WaveField)))) + call WrScr("-----------------------------------------------------------") + end subroutine ShowPassedData +end subroutine + + +!> set the pointer to the WaveField data +subroutine SeaSt_C_SetWaveFieldPointer(WaveFieldPointer_C,ErrStat_C,ErrMsg_C) BIND (C, NAME='SeaSt_C_SetWaveFieldPointer') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_SetWaveFieldPointer +!GCC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_SetWaveFieldPointer +#endif + type(c_ptr), intent(in ) :: WaveFieldPointer_C + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + integer :: ErrStat + character(ErrMsgLen) :: ErrMsg + character(*), parameter :: RoutineName = 'SeaSt_C_SetWaveFieldPointer' + ErrStat = ErrID_None + ErrMSg = "" + call C_F_POINTER(WaveFieldPointer_C, p%WaveField) + if (associated(p%WaveField)) then + ! basic sanity check + if (.not. allocated(p%WaveField%WaveTime)) then + call SetErrStat(ErrID_Fatal,"Invalid pointer passed in, or WaveField not initialized",ErrStat,ErrMsg,RoutineName) + endif + else + call SetErrStat(ErrID_Fatal,"Invalid pointer passed in, or WaveField not initialized",ErrStat,ErrMsg,RoutineName) + endif + call SetErrStat_F2C( ErrStat, ErrMsg, ErrStat_C, ErrMsg_C ) + if (DebugLevel > 1) call ShowPassedData() + return +contains + subroutine ShowPassedData() + call WrScr("-----------------------------------------------------------") + call WrScr("Interface debugging: SeaSt_C_SetWaveFieldPointer") + call WrScr(" --------------------------------------------------------") + call WrScr(" WaveFieldPointer_C <- "//trim(Num2LStr(loc(p%WaveField)))) + call WrScr("-----------------------------------------------------------") + end subroutine ShowPassedData +end subroutine + + +!> Get the fluid velocity, acceleration, and node-in-water status at time+position coordinate +!! NOTE: if wave stretching is turned off, the SWL is used as the cutoff for the nodeInWater and for Vel / Acc values +subroutine SeaSt_C_GetFluidVelAcc(Time_C, Pos_C, Vel_C, Acc_C, NodeInWater_C, ErrStat_C,ErrMsg_C) BIND (C, NAME='SeaSt_C_GetFluidVelAcc') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_GetFluidVelAcc +!GCC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_GetFluidVelAcc +#endif + real(c_double), intent(in ) :: Time_C + real(c_float), intent(in ) :: Pos_c(3) + real(c_float), intent( out) :: Vel_c(3) + real(c_float), intent( out) :: Acc_c(3) + integer(c_int), intent( out) :: NodeInWater_C + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + real(DbKi) :: Time + real(ReKi) :: Pos(3) + real(SiKi) :: Vel(3) + real(SiKi) :: Acc(3) + logical :: forceNodeInWater +!FIXME:dev-tc uncomment next line +! logical :: fetchDynCurrent + integer(IntKi) :: nodeInWater + integer :: ErrStat, ErrStat2 + character(ErrMsgLen) :: ErrMsg, ErrMsg2 + character(*), parameter :: RoutineName = 'SeaSt_C_GetFluidVelAcc' + + ! Initialize error handling + ErrStat = ErrID_None + ErrMsg = "" + + forceNodeInWater = .false. + + if (DebugLevel > 1) call ShowPassedData() + + ! convert position and time to fortran types + Time = real(Time_C, DbKi) + Pos = real(Pos_C, ReKi) + + ! get wave field velocity and acceleration (current is included in this) + ! Notes: + ! - if node is out of water, velocity and acceleration are zero + ! - if position is outside the wave field boundary, it will simply return boundary edge value + ! - time must be positive or a fatal error occurs + call WaveField_GetNodeWaveVelAcc( p%WaveField, m%WaveField_m, Time, pos, forceNodeInWater, nodeInWater, Vel, Acc, ErrStat, ErrMsg ) +!FIXME:dev-tc use next line instead of above +! call WaveField_GetNodeWaveVelAcc( p%WaveField, m%WaveField_m, Time, pos, forceNodeInWater, fetchDynCurrent, nodeInWater, Vel, Acc, ErrStat, ErrMsg ) + + ! Store resulting velocity and acceleration as C type + Vel_c = real(Vel,c_float) + Acc_c = real(Acc,c_float) + + ! Density value and node status to return + if (nodeInWater == 1_IntKi) then + NodeInWater_C = 1_c_int + else + NodeInWater_C = 0_c_int + endif + + call SetErrStat_F2C( ErrStat, ErrMsg, ErrStat_C, ErrMsg_C ) ! convert error from fortran to C for return + if (DebugLevel > 1) call ShowReturnData() + return +contains + subroutine ShowPassedData() + call WrScr("-----------------------------------------------------------") + call WrScr("Interface debugging: SeaSt_C_GetFluidVelAccDens") + call WrScr(" --------------------------------------------------------") + call WrScr(" Time_C -> "//trim(Num2LStr(Time_C))) + call WrScr(" Pos_C -> ("//trim(Num2LStr(Pos_C(1)))//","//trim(Num2LStr(Pos_C(2)))//","//trim(Num2LStr(Pos_C(3)))//")") + end subroutine ShowPassedData + subroutine ShowReturnData() + call WrScr(" Vel_C <- ("//trim(Num2LStr(Vel_C(1)))//","//trim(Num2LStr(Vel_C(2)))//","//trim(Num2LStr(Vel_C(3)))//")") + call WrScr(" Acc_C <- ("//trim(Num2LStr(Acc_C(1)))//","//trim(Num2LStr(Acc_C(2)))//","//trim(Num2LStr(Acc_C(3)))//")") + call WrScr(" NodeInWater_C <- "//trim(Num2LStr(NodeInWater_C))) + call WrScr("-----------------------------------------------------------") + end subroutine ShowReturnData +end subroutine SeaSt_C_GetFluidVelAcc + + + +!> return the surface elevation and normal vector at a point. +subroutine SeaSt_C_GetSurfElev(Time_C, Pos_C, Elev_C, ErrStat_C,ErrMsg_C) BIND (C, NAME='SeaSt_C_GetSurfElev') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_GetSurfElev +!GCC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_GetSurfElev +#endif + real(c_double), intent(in ) :: Time_C + real(c_float), intent(in ) :: Pos_c(2) + real(c_float), intent( out) :: Elev_C + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + + real(DbKi) :: Time + real(ReKi) :: Pos(2) + real(SiKi) :: Elev + integer :: ErrStat + character(ErrMsgLen) :: ErrMsg + character(*), parameter :: RoutineName = 'SeaSt_C_GetSurfElev' + + if (DebugLevel > 1) call ShowPassedData() + + ! convert position and time to fortran types + Time = real(Time_C, DbKi) + Pos = real(Pos_C(1:2), ReKi) + + ! get wave elevation (total combined first and second order) + ! Notes: + ! - if position is outside the wave field boundary, it will simply return boundary edge value + ! - time must be positive or a fatal error occurs + Elev = WaveField_GetNodeTotalWaveElev( p%WaveField, m%WaveField_m, Time, pos, ErrStat, ErrMsg ) + + ! Store resulting elevation as C type + Elev_C = real(Elev,c_float) + + if (DebugLevel > 1) call ShowReturnData() + call Cleanup() + return +contains + subroutine Cleanup() ! NOTE: we are ignoring any error reporting from here + CALL SetErrStat_F2C(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + end subroutine Cleanup + subroutine ShowPassedData() + call WrScr("-----------------------------------------------------------") + call WrScr("Interface debugging: SeaSt_C_GetSurfElev") + call WrScr(" --------------------------------------------------------") + call WrScr(" Time_C -> "//trim(Num2LStr(Time_C))) + call WrScr(" Pos_C -> ("//trim(Num2LStr(Pos_C(1)))//","//trim(Num2LStr(Pos_C(2)))//")") + end subroutine ShowPassedData + subroutine ShowReturnData() + call WrScr(" Elev_C <- "//trim(Num2LStr(Elev_C))) + call WrScr("-----------------------------------------------------------") + end subroutine ShowReturnData +end subroutine SeaSt_C_GetSurfElev + + + +!> return the surface normal vector at a point. +subroutine SeaSt_C_GetSurfNorm(Time_C, Pos_C, NormVec_C, ErrStat_C,ErrMsg_C) BIND (C, NAME='SeaSt_C_GetSurfNorm') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_GetSurfNorm +!GCC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_GetSurfNorm +#endif + real(c_double), intent(in ) :: Time_C + real(c_float), intent(in ) :: Pos_c(2) + real(c_float), intent( out) :: NormVec_C(3) + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + + real(DbKi) :: Time + real(ReKi) :: Pos(2), r + real(ReKi) :: NormVec(3) + integer :: ErrStat + character(ErrMsgLen) :: ErrMsg + character(*), parameter :: RoutineName = 'SeaSt_C_GetSurfNorm' + + if (DebugLevel > 1) call ShowPassedData() + + ! convert position and time to fortran types + Time = real(Time_C, DbKi) + Pos = real(Pos_C(1:2), ReKi) + + ! get the normal vector at the point (set to vertical if outside region) + r = 0.0_ReKi ! distance for central differencing - use default in code + call WaveField_GetNodeWaveNormal( p%WaveField, m%WaveField_m, Time, pos, r, NormVec, ErrStat, ErrMsg ) + + ! Store resulting normal vector as C type + NormVec_C = real(NormVec,c_float) + + if (DebugLevel > 1) call ShowReturnData() + call Cleanup() + return +contains + subroutine Cleanup() ! NOTE: we are ignoring any error reporting from here + CALL SetErrStat_F2C(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + end subroutine Cleanup + subroutine ShowPassedData() + call WrScr("-----------------------------------------------------------") + call WrScr("Interface debugging: SeaSt_C_GetSurfNorm") + call WrScr(" --------------------------------------------------------") + call WrScr(" Time_C -> "//trim(Num2LStr(Time_C))) + call WrScr(" Pos_C -> ("//trim(Num2LStr(Pos_C(1)))//","//trim(Num2LStr(Pos_C(2)))//")") + end subroutine ShowPassedData + subroutine ShowReturnData() + call WrScr(" NormVec_C <- ("//trim(Num2LStr(NormVec_C(1)))//","//trim(Num2LStr(NormVec_C(2)))//","//trim(Num2LStr(NormVec_C(3)))//")") + call WrScr("-----------------------------------------------------------") + end subroutine ShowReturnData +end subroutine SeaSt_C_GetSurfNorm + + + + +!FIXME: the following visualization writer should be merged into the library vtk.f90 +! this is a modified duplicate of the routine from FAST_Subs by the same name. +!FIXME: this routine currently only grabs the closest timestep and does not do any interpolation +!---------------------------------------------------------------------------------------------------------------------------------- +!> This subroutine writes the wave elevation data for a given time step +subroutine WrVTK_WaveElevVisGrid(t_global, vtk, ErrStat, ErrMsg) + real(DbKi), intent(in ) :: t_global !< Current global time + type(VTKvis), intent(inout) :: vtk + integer(IntKi), intent( out) :: ErrStat + character(ErrMsgLen), intent( out) :: ErrMsg + integer(IntKi) :: Un ! fortran unit number + integer(IntKi) :: n, iy, ix ! loop counters + real(SiKi) :: t + integer(IntKi) :: count ! which file is this? + character(1024) :: FileName + integer(IntKi) :: NumberOfPoints + integer(IntKi), parameter :: NumberOfLines = 0 + integer(IntKi) :: NumberOfPolys + character(1024) :: Tstr + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + character(*),parameter :: RoutineName = 'WrVTK_WaveElevVisGrid' + + ! Initialize error handling + ErrStat = ErrID_None + ErrMsg = "" + + NumberOfPoints = vtk%NWaveElevPts(1) * vtk%NWaveElevPts(2) + ! I'm going to make triangles for now. we should probably just make this a structured file at some point + NumberOfPolys = ( vtk%NWaveElevPts(1) - 1 ) * & + ( vtk%NWaveElevPts(2) - 1 ) * 2 + + count = nint(t_global / vtk%dt) + if (count == vtk%lastCount) return ! already wrote this one + vtk%lastCount = count ! store the current number to make sure we don't write it twice + + !................................................................. + ! write the data that potentially changes each time step: + !................................................................. + ! construct the string for the zero-padded VTK write-out step + write(Tstr, '(i' // trim(Num2LStr(vtk%tWidth)) //'.'// trim(Num2LStr(vtk%tWidth)) // ')') count + + ! PolyData (.vtp) - Serial vtkPolyData (unstructured) file + FileName = TRIM(vtk%OutRootName)//'.WaveSurface.'//TRIM(Tstr)//'.vtp' + + call WrVTK_header( FileName, NumberOfPoints, NumberOfLines, NumberOfPolys, Un, ErrStat, ErrMsg ) + if (ErrStat >= AbortErrLev) return + +! points (nodes, augmented with NumSegments): + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + + ! I'm not going to interpolate in time; I'm just going to get the index of the closest wave time value + t = REAL(t_global,SiKi) + call GetWaveElevIndx( t, p%WaveField%WaveTime, vtk%LastWaveIndx ) + + do ix=1,vtk%NWaveElevPts(1) + do iy=1,vtk%NWaveElevPts(2) + write(Un,VTK_AryFmt) vtk%WaveElevVisX(ix), vtk%WaveElevVisY(iy), vtk%WaveElevVisGrid(vtk%LastWaveIndx,ix,iy) + end do + end do + + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + + do ix=1,vtk%NWaveElevPts(1)-1 + do iy=1,vtk%NWaveElevPts(2)-1 + n = vtk%NWaveElevPts(2)*(ix-1)+iy - 1 ! points start at 0 + write(Un,'(3(i7))') n, n+1, n+vtk%NWaveElevPts(2) + write(Un,'(3(i7))') n+1, n+1+vtk%NWaveElevPts(2), n+vtk%NWaveElevPts(2) + end do + end do + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + do n=1,NumberOfPolys + WRITE(Un,'(i7)') 3*n + end do + write(Un,'(A)') ' ' + write(Un,'(A)') ' ' + call WrVTK_footer( Un ) +contains + !---------------------------------------------------------------------------------------------------------------------------------- + !> This function returns the index, Ind, of the XAry closest to XValIn, where XAry is assumed to be periodic. It starts + !! searching at the value of Ind from a previous step. + subroutine GetWaveElevIndx( XValIn, XAry, Ind ) + integer, intent(inout) :: Ind ! Initial and final index into the arrays. + real(SiKi), intent(in) :: XAry (:) !< Array of X values to be interpolated. + real(SiKi), intent(in) :: XValIn !< X value to be found + integer :: AryLen ! Length of the arrays. + real(SiKi) :: XVal !< X to be found (wrapped/periodic) + + AryLen = size(XAry) + + ! Wrap XValIn into the range XAry(1) to XAry(AryLen) + XVal = MOD(XValIn, XAry(AryLen)) + + ! Let's check the limits first. + if ( XVal <= XAry(1) ) then + Ind = 1 + return + else if ( XVal >= XAry(AryLen) ) then + Ind = AryLen + return + else + ! Set the Ind to the first index if we are at the beginning of XAry + if ( XVal <= XAry(2) ) then + Ind = 1 + end if + end if + + ! Let's interpolate! + Ind = MAX( MIN( Ind, AryLen-1 ), 1 ) + do + if ( XVal < XAry(Ind) ) then + Ind = Ind - 1 + else if ( XVal >= XAry(Ind+1) ) then + Ind = Ind + 1 + else + ! XAry(Ind) <= XVal < XAry(Ind+1) + ! this would make it the "closest" node, but I'm not going to worry about that for visualization purposes + !if ( XVal > (XAry(Ind+1) + XAry(Ind))/2.0_SiKi ) Ind = Ind + 1 + return + end if + end do + return + end subroutine GetWaveElevIndx +end subroutine WrVTK_WaveElevVisGrid + +!> clear any memory that the _End routine would not clear. Note we are ignoring any errors +subroutine ClearMem() + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + if (allocated(vtk%WaveElevVisX )) deallocate(vtk%WaveElevVisX ) + if (allocated(vtk%WaveElevVisY )) deallocate(vtk%WaveElevVisY ) + if (allocated(vtk%WaveElevVisGrid)) deallocate(vtk%WaveElevVisGrid) + call SeaSt_DestroyInitInput( InitInp, ErrStat2, ErrMsg2) + call SeaSt_DestroyInitOutput(InitOutData, ErrStat2, ErrMsg2) +end subroutine ClearMem end module SeaState_C_Binding diff --git a/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index cf14d143f2..b4725e8ddc 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -250,6 +250,15 @@ function(seast_regression TESTNAME LABEL) regression(${TEST_SCRIPT} ${SEASTATE_EXECUTABLE} ${SOURCE_DIRECTORY} ${BUILD_DIRECTORY} " " ${TESTNAME} "${LABEL}" " ") endfunction(seast_regression) +# py_seastate +function(py_seast_regression TESTNAME LABEL) + set(TEST_SCRIPT "${CMAKE_CURRENT_LIST_DIR}/executeSeaStatePyRegressionCase.py") + set(SEASTATE_EXECUTABLE "${Python_EXECUTABLE}") + set(SOURCE_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/..") + set(BUILD_DIRECTORY "${CTEST_BINARY_DIR}/modules/seastate") + regression(${TEST_SCRIPT} ${SEASTATE_EXECUTABLE} ${SOURCE_DIRECTORY} ${BUILD_DIRECTORY} " " ${TESTNAME} "${LABEL}" " ") +endfunction(py_seast_regression) + # moordyn function(md_regression TESTNAME LABEL) set(TEST_SCRIPT "${CMAKE_CURRENT_LIST_DIR}/executeMoordynRegressionCase.py") @@ -510,6 +519,7 @@ seast_regression("seastate_WaveMod7_WaveStMod1" "seastate") seast_regression("seastate_WaveMod7_WaveStMod2" "seastate") seast_regression("seastate_WaveMod7_WaveStMod3" "seastate") seast_regression("seastate_wavemod5" "seastate") # place at end since it reads outputs generated by seastate_wr_kin1 +py_seast_regression("py_seastate_1" "seastate;python") # MoorDyn regression tests md_regression("md_5MW_OC4Semi" "moordyn") diff --git a/reg_tests/executeSeaStatePyRegressionCase.py b/reg_tests/executeSeaStatePyRegressionCase.py new file mode 100644 index 0000000000..2d2f8b875e --- /dev/null +++ b/reg_tests/executeSeaStatePyRegressionCase.py @@ -0,0 +1,146 @@ +# +# 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 SeaSate through the python interface 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: `executeSeaStatePyRegressionCase.py -h` +""" + +import os +import sys +basepath = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, os.path.sep.join([basepath, "lib"])) +import argparse +import numpy as np +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 SeaState c-bindings library interface with 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="SeaState-Python", type=str, nargs=1, help="The path to the InflowWind driver 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("rtol", metavar="Relative-Tolerance", type=float, nargs=1, help="Relative tolerance to allow the solution to deviate; expressed as order of magnitudes less than baseline.") +parser.add_argument("atol", metavar="Absolute-Tolerance", type=float, nargs=1, help="Absolute tolerance to allow small values to pass; expressed as order of magnitudes less than baseline.") +parser.add_argument("-p", "-plot", dest="plot", action='store_true', help="bool to include plots in failed cases") +parser.add_argument("-n", "-no-exec", dest="noExec", action='store_true', help="bool to prevent execution of the test cases") +parser.add_argument("-v", "-verbose", dest="verbose", action='store_true', 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] +rtol = args.rtol[0] +atol = args.atol[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, exist_ok=True) + +### 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", "seastate") +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,"*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 inflowwind on the test case +if not noExec: + caseInputFile = os.path.join(testBuildDirectory, "py_seastate_driver.py") + returnCode = openfastDrivers.runInflowwindDriverCase(caseInputFile, executable) + if returnCode != 0: + sys.exit(returnCode*10) + +### Build the filesystem navigation variables for running the regression test. +# the "Points.Results.dat" file is from calls to routines other than CalcOutput +localOutFile = os.path.join(testBuildDirectory, "Points.Results.dat") +baselineOutFile = os.path.join(targetOutputDirectory, "Points.Results.dat") +rtl.validateFileOrExit(localOutFile) +rtl.validateFileOrExit(baselineOutFile) + +testData, testInfo, _ = pass_fail.readFASTOut(localOutFile) +baselineData, baselineInfo, _ = pass_fail.readFASTOut(baselineOutFile) + +passing_channels = pass_fail.passing_channels(testData.T, baselineData.T, rtol, atol) +passing_channels = passing_channels.T + +norms = pass_fail.calculateNorms(testData, baselineData) + +# export all case summaries +channel_names = testInfo["attribute_names"] +exportCaseSummary(testBuildDirectory, caseName, channel_names, passing_channels, norms) + +# passing case +if np.all(passing_channels): + sys.exit(0) + +# failing case +if plotError: + from errorPlotting import finalizePlotDirectory, plotOpenfastError + for channel in testInfo["attribute_names"]: + try: + plotOpenfastError(localOutFile, baselineOutFile, channel, rtol, atol) + except: + error = sys.exc_info()[1] + print("Error generating plots: {}".format(error)) + finalizePlotDirectory(localOutFile, testInfo["attribute_names"], caseName) + +sys.exit(1) diff --git a/reg_tests/r-test b/reg_tests/r-test index 8d8b510f0d..2cd70eedcd 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 8d8b510f0d9615a63bbe80c6c327046f051b520f +Subproject commit 2cd70eedcd35c6c58ec9aa4798e6883921b4a817