From cac691898cce32bf309e0053d04f6402016769f2 Mon Sep 17 00:00:00 2001 From: "Harish @ NLR" Date: Thu, 1 Oct 2026 12:09:39 -0600 Subject: [PATCH 1/3] BeamDyn: add C interface for driving a standalone beam from C/C++ codes Add a C binding for BeamDyn in the style of the MoorDyn and AeroDyn-Inflow C bindings, so that an external C/C++ code can drive a single BeamDyn beam without the OpenFAST glue code. Entry points (BeamDyn_C_Binding.f90, declared in BeamDyn_C_Binding.h): BD_C_Init, BD_C_GetRefPositions, BD_C_SetRootMotion, BD_C_SetPointLoads, BD_C_SetDistrLoads, BD_C_UpdateStates, BD_C_CalcOutput, BD_C_PackStates, BD_C_UnpackStates, BD_C_End. All BeamDyn data is held in module-level variables; the input history for linear or quadratic interpolation and the state history needed for correction steps are handled as in the other bindings, with the same ErrStat_C/ErrMsg_C error handling. The input file may be passed as a path or as its contents (written to a temporary file, since BeamDyn reads its input from a file unit). Checkpoint files use BeamDyn's registry pack/unpack routines. Root kinematics and gravity cross the interface in double precision: BeamDyn imposes the root motion as a boundary condition at every step, and single precision rounding of a prescribed root motion appears as root acceleration noise. Loads and outputs are single precision as in the other bindings; orientations are double precision. CMake: the driver subroutines move into a beamdyn_driver_subs library (the binding uses them to write the output file), and the new targets beamdyn_c_binding (shared) and beamdyn_c_bind_static are installed with the header. Co-Authored-By: Claude Fable 5.1 --- modules/beamdyn/CMakeLists.txt | 42 +- modules/beamdyn/src/BeamDyn_C_Binding.f90 | 1147 +++++++++++++++++++++ modules/beamdyn/src/BeamDyn_C_Binding.h | 206 ++++ 3 files changed, 1390 insertions(+), 5 deletions(-) create mode 100644 modules/beamdyn/src/BeamDyn_C_Binding.f90 create mode 100644 modules/beamdyn/src/BeamDyn_C_Binding.h diff --git a/modules/beamdyn/CMakeLists.txt b/modules/beamdyn/CMakeLists.txt index 9990c4c768..974ec87c66 100644 --- a/modules/beamdyn/CMakeLists.txt +++ b/modules/beamdyn/CMakeLists.txt @@ -27,15 +27,47 @@ add_library(beamdynlib STATIC ) target_link_libraries(beamdynlib nwtclibs) -add_executable(beamdyn_driver - src/Driver_Beam.f90 +# Driver subroutines (also used by the C-bindings interface for writing the output file) +add_library(beamdyn_driver_subs STATIC src/Driver_Beam_Subs.f90 ) -target_link_libraries(beamdyn_driver beamdynlib versioninfolib) +target_link_libraries(beamdyn_driver_subs beamdynlib versioninfolib) + +add_executable(beamdyn_driver + src/Driver_Beam.f90 +) +target_link_libraries(beamdyn_driver beamdyn_driver_subs) + +# C-bindings interface library +# create object instead of directly linking into shared and static -- causes issues in parallel builds +# NOTE: target linking at the object, static, and shared libraries. Different CMake versions handle this +# slightly differently with unpredictable results if I don't. +add_library(beamdyn_c_binding_object OBJECT src/BeamDyn_C_Binding.f90) +target_link_libraries(beamdyn_c_binding_object beamdynlib beamdyn_driver_subs nwtclibs versioninfolib) +set_property(TARGET beamdyn_c_binding_object PROPERTY POSITION_INDEPENDENT_CODE 1) # required for shared libs +if(APPLE OR UNIX) + target_compile_definitions(beamdyn_c_binding_object PRIVATE IMPLICIT_DLLEXPORT) +endif() -install(TARGETS beamdynlib beamdyn_driver +# Shared +add_library(beamdyn_c_binding SHARED $) +target_link_libraries(beamdyn_c_binding beamdynlib beamdyn_driver_subs nwtclibs versioninfolib) +target_include_directories(beamdyn_c_binding PUBLIC + $ +) +set_target_properties(beamdyn_c_binding PROPERTIES PUBLIC_HEADER src/BeamDyn_C_Binding.h) + +# C-bindings non-shared interface +add_library(beamdyn_c_bind_static STATIC $) +target_link_libraries(beamdyn_c_bind_static beamdynlib beamdyn_driver_subs nwtclibs versioninfolib) +if(APPLE OR UNIX) + target_compile_definitions(beamdyn_c_bind_static PRIVATE IMPLICIT_DLLEXPORT) +endif() + +install(TARGETS beamdynlib beamdyn_driver_subs beamdyn_driver beamdyn_c_binding beamdyn_c_bind_static EXPORT "${CMAKE_PROJECT_NAME}Libraries" RUNTIME DESTINATION bin LIBRARY DESTINATION lib ARCHIVE DESTINATION lib -) \ No newline at end of file + PUBLIC_HEADER DESTINATION include +) diff --git a/modules/beamdyn/src/BeamDyn_C_Binding.f90 b/modules/beamdyn/src/BeamDyn_C_Binding.f90 new file mode 100644 index 0000000000..b3cfb0c2ee --- /dev/null +++ b/modules/beamdyn/src/BeamDyn_C_Binding.f90 @@ -0,0 +1,1147 @@ +!********************************************************************************************************************************** +! LICENSING +! Copyright (C) 2026 National Renewable Energy Laboratory +! +! This file is part of BeamDyn. +! +! 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. +! +!********************************************************************************************************************************** +!> C interface to a standalone BeamDyn beam. This library lets a C/C++ code drive a single BeamDyn beam (one root motion +!! input, point and distributed load inputs, and the blade motion and root reaction outputs) without the OpenFAST glue +!! code. The structure of this interface follows the MoorDyn and AeroDyn-Inflow C bindings: all BeamDyn data is held +!! in module-level variables, the interpolation/extrapolation input history and the state history needed for correction +!! steps are handled here, and all values cross the interface as flat C arrays. +!! +!! Call sequence: +!! BD_C_Init -- initialize the beam from an input file (returns node counts and output channel info) +!! BD_C_GetRefPositions -- (optional) reference positions of the output nodes and load nodes +!! loop over time: +!! BD_C_SetRootMotion -- root motion at the time of the next BD_C_CalcOutput or BD_C_UpdateStates call +!! BD_C_SetPointLoads -- (optional) point loads at the finite element nodes +!! BD_C_SetDistrLoads -- (optional) distributed loads at the quadrature point nodes +!! BD_C_CalcOutput -- outputs at the current time (node motions, root reaction, output channels) +!! BD_C_UpdateStates -- advance the states from Time_C to TimeNext_C using the inputs set for TimeNext_C +!! BD_C_PackStates / BD_C_UnpackStates -- (optional) write/read a checkpoint file with the complete beam state +!! BD_C_End +MODULE BeamDyn_C + + USE ISO_C_BINDING + USE BeamDyn + USE BeamDyn_Subs + USE BeamDyn_Types + USE BeamDyn_driver_subs, ONLY: Dvr_InitializeOutputFile, Dvr_WriteOutputLine + USE NWTC_Library + USE NWTC_C_Binding, ONLY: IntfStrLen, ErrMsgLen_C, FileNameFromCString, SetErrStat_F2C + USE VersionInfo + + IMPLICIT NONE + SAVE + + PUBLIC :: BD_C_Init + PUBLIC :: BD_C_GetRefPositions + PUBLIC :: BD_C_SetRootMotion + PUBLIC :: BD_C_SetPointLoads + PUBLIC :: BD_C_SetDistrLoads + PUBLIC :: BD_C_UpdateStates + PUBLIC :: BD_C_CalcOutput + PUBLIC :: BD_C_PackStates + PUBLIC :: BD_C_UnpackStates + PUBLIC :: BD_C_End + + PRIVATE + + !------------------------------------------------------------------------------------ + ! Version info for display + TYPE(ProgDesc), PARAMETER :: version = ProgDesc( 'BeamDyn library', '', '' ) + + !------------------------------------------------------------------------------------ + ! Precision of the interface + ! The root kinematics (position, orientation, velocity, acceleration) and gravity are passed in double + ! precision: BeamDyn imposes the root motion as a boundary condition at every step, and single precision + ! rounding of a prescribed root motion appears as root acceleration noise of order eps_single*|disp|/dt^2. + ! Loads and the returned node motions, reaction loads, and channel values are single precision (as in the + ! MoorDyn and AeroDyn-Inflow interfaces); orientations are always double precision. + ! + ! Potential issues + ! - if MaxBDOutputs is sufficiently large, we may overrun the buffer on the calling + ! side (OutputChannelNames_C,OutputChannelUnits_C). The calling code must size + ! those buffers as ChanLen*MaxBDOutputs+1 characters. BeamDyn has at most + ! MaxOutPts regular channels plus BldNd_MaxOutPts channels per output node. + INTEGER(IntKi), PARAMETER :: MaxBDOutputs = 8000 + + !------------------------------------------------------------------------------------ + ! Checkpoint file identifier (written at the start of the file by BD_C_PackStates) + INTEGER(IntKi), PARAMETER :: CheckpointFileID = 42101 + + !-------------------------------------------------------------------------------------------------------------------------------------------------------- + ! Data storage + ! All BeamDyn data is stored within the following data structures inside this + ! module. No data is stored within BeamDyn itself, but is instead passed in + ! from this module. This data is not available to the calling code unless + ! explicitly passed through the interface (derived types such as these are + ! non-trivial to pass through the c-bindings). + TYPE(BD_InitInputType) :: InitInp !< Input data for initialization routine + TYPE(BD_InputType), ALLOCATABLE :: u(:) !< Inputs at the times in InputTimes (input meshes are defined in BD_Init) + TYPE(BD_InputType) :: BD_u !< Inputs set by the BD_C_Set* routines. Copied into u(:) by UpdateStates and CalcOutput + TYPE(BD_ParameterType) :: p !< Parameters + TYPE(BD_ContinuousStateType) :: x(0:2) !< Continuous states + TYPE(BD_DiscreteStateType) :: xd(0:2) !< Discrete states + TYPE(BD_ConstraintStateType) :: z(0:2) !< Constraint states + TYPE(BD_OtherStateType) :: OtherSt(0:2) !< Other states + TYPE(BD_OutputType) :: y !< System outputs + TYPE(BD_MiscVarType) :: m !< Misc/optimization variables + TYPE(BD_InitOutputType) :: InitOutData !< Output for initialization routine + + !-------------------------------------------------------------------------------------------------------------------------------------------------------- + ! Time tracking + ! For the solver in BD, 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 BD. 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 + ! 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. + INTEGER(IntKi) :: InterpOrder !< Interpolation order: must be 1 (linear) or 2 (quadratic) + REAL(DbKi), DIMENSION(:), ALLOCATABLE :: InputTimes(:) !< InputTimes array + REAL(DbKi) :: InputTimePrev !< input time of last UpdateStates call + 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 (time passed to the first UpdateStates call) + LOGICAL :: Initialized = .FALSE. !< BD_C_Init completed successfully + + ! Mesh and channel sizes (kept here so the interface array sizes are defined even before BD_C_Init is called) + INTEGER(IntKi) :: NumOutputNodes = 0 !< Number of nodes on y%BldMotion + INTEGER(IntKi) :: NumPointLoadNodes = 0 !< Number of nodes on u%PointLoad + INTEGER(IntKi) :: NumDistrLoadNodes = 0 !< Number of nodes on u%DistrLoad + INTEGER(IntKi) :: NumChannels = 0 !< Number of output channels (size of y%WriteOutput) + + ! We are including the previous state info here (not done in OpenFAST 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 + + !-------------------------------------------------------------------------------------------------------------------------------------------------------- + ! Output file + ! When requested at Init, the output channels are written to .out in the same format + ! as the BeamDyn driver writes them. One line is written per CalcOutput call at a new time. + INTEGER(IntKi) :: WrOutputs = 0 !< Write the output channels to a file (0: no, 1: yes) + INTEGER(IntKi) :: UnOutFile = -1 !< Unit number of the output file + REAL(DbKi) :: OutTimePrev !< Time of the last line written to the output file + CHARACTER(IntfStrLen) :: OutRootName !< Root name for the output, summary, and echo files + +CONTAINS + +!=============================================================================================================== +!---------------------------------------------- BD INIT -------------------------------------------------------- +!=============================================================================================================== +!> Initialize a BeamDyn beam. The beam is described by the BeamDyn primary input file (and the blade file it +!! references). The root of the beam is placed at RootPos_C with orientation RootOri_C; this also sets the +!! reference frame BeamDyn uses for its internal calculations. +SUBROUTINE BD_C_Init( & + InputFilePassed, InputFileString_C, InputFileStringLength_C, & + OutRootName_C, & + RootPos_C, RootOri_C, RootVel_C, Gravity_C, & + DT_C, InterpOrder_C, DynamicSolve_C, WrOutputs_C, & + NumOutputNodes_C, NumPointLoadNodes_C, NumDistrLoadNodes_C, & + NumChannels_C, OutputChannelNames_C, OutputChannelUnits_C, & + ErrStat_C, ErrMsg_C & +) BIND (C, NAME='BD_C_Init') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: BD_C_Init +!GCC$ ATTRIBUTES DLLEXPORT :: BD_C_Init +#endif + INTEGER(C_INT), INTENT(IN ) :: InputFilePassed !< Whether to load the file from the filesystem - 1: InputFileString_C contains the contents of the input file; otherwise, InputFileString_C contains the path to the input file + TYPE(C_PTR), INTENT(IN ) :: InputFileString_C !< Input file as a single string with lines delineated by C_NULL_CHAR + INTEGER(C_INT), INTENT(IN ) :: InputFileStringLength_C !< length of the input file string + CHARACTER(KIND=C_CHAR), INTENT(IN ) :: OutRootName_C(*) !< Root name to use for the summary, echo, and output files (C_NULL_CHAR terminated, at most IntfStrLen characters) + REAL(C_DOUBLE), INTENT(IN ) :: RootPos_C(3) !< Initial position of the beam root in the global frame (m) + REAL(C_DOUBLE), INTENT(IN ) :: RootOri_C(9) !< Initial orientation of the beam root: DCM from the global frame to the root frame, stored row by row [r11,r12,r13,r21,r22,r23,r31,r32,r33] + REAL(C_DOUBLE), INTENT(IN ) :: RootVel_C(6) !< Initial translational (1:3) and rotational (4:6) velocity of the beam root in the global frame (m/s, rad/s) + REAL(C_DOUBLE), INTENT(IN ) :: Gravity_C(3) !< Gravitational acceleration vector in the global frame (m/s^2) + REAL(C_DOUBLE), INTENT(IN ) :: DT_C !< Timestep used with BD for stepping forward from t to t+dt. Must be constant. + INTEGER(C_INT), INTENT(IN ) :: InterpOrder_C !< Interpolation order to use (must be 1 or 2) + INTEGER(C_INT), INTENT(IN ) :: DynamicSolve_C !< 1: dynamic solve; 0: static solve (loads are ramped over the first UpdateStates calls) + INTEGER(C_INT), INTENT(IN ) :: WrOutputs_C !< 1: write the output channels to .out at each CalcOutput call; 0: do not write a file + INTEGER(C_INT), INTENT( OUT) :: NumOutputNodes_C !< Number of nodes on the blade motion output mesh + INTEGER(C_INT), INTENT( OUT) :: NumPointLoadNodes_C !< Number of nodes on the point load input mesh (finite element nodes) + INTEGER(C_INT), INTENT( OUT) :: NumDistrLoadNodes_C !< Number of nodes on the distributed load input mesh (quadrature point nodes) + 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*MaxBDOutputs+1) !< Output channel names, ChanLen characters each, C_NULL_CHAR terminated + CHARACTER(KIND=C_CHAR), INTENT( OUT) :: OutputChannelUnits_C(ChanLen*MaxBDOutputs+1) !< Output channel units, ChanLen characters each, C_NULL_CHAR terminated + INTEGER(C_INT), INTENT( OUT) :: ErrStat_C !< Error status + CHARACTER(KIND=C_CHAR), INTENT( OUT) :: ErrMsg_C(ErrMsgLen_C) !< Error message (C_NULL_CHAR terminated) + + ! Local Variables + CHARACTER(KIND=C_char, LEN=InputFileStringLength_C), POINTER :: InputFileString !< Input file as a single string with NULL character separating lines + CHARACTER(IntfStrLen) :: TmpFileName !< Temporary file name if passing the input file contents directly + REAL(DbKi) :: dT_Interval !< Timestep passed to BD_Init + INTEGER(IntKi) :: ErrStat_F, ErrStat_F2 + CHARACTER(ErrMsgLen) :: ErrMsg_F, ErrMsg_F2 + INTEGER(IntKi) :: I, J, K + CHARACTER(*), PARAMETER :: RoutineName = 'BD_C_Init' + + ! Initialize library and display info on this compile + ErrStat_F = ErrID_None + ErrMsg_F = '' + NumOutputNodes_C = 0_c_int + NumPointLoadNodes_C = 0_c_int + NumDistrLoadNodes_C = 0_c_int + NumChannels_C = 0_c_int + OutputChannelNames_C(:) = '' + OutputChannelUnits_C(:) = '' + TmpFileName = '' + + CALL NWTC_Init( ProgNameIn=version%Name ) + CALL DispCopyrightLicense( version%Name ) + CALL DispCompileRuntimeInfo( version%Name ) + + ! Destroy global memory (in case Init is called a second time without an End) + CALL DestroyAll( ErrStat_F2, ErrMsg_F2 ) + + !---------------------------------------------------------------------------------------------------------------------------------------------- + ! Root name for output files + !---------------------------------------------------------------------------------------------------------------------------------------------- + OutRootName = CStringToFortran( OutRootName_C ) + IF ( LEN_TRIM(OutRootName) == 0 ) OutRootName = 'BDroot' + + !---------------------------------------------------------------------------------------------------------------------------------------------- + ! Input file + ! BeamDyn reads its primary input file (and the blade file named in it) from the file system. When the contents of + ! the primary input file are passed in, they are written to a temporary file next to the output root so that BeamDyn + ! can read them; the blade file named in the primary input file is then located relative to that directory. + !---------------------------------------------------------------------------------------------------------------------------------------------- + CALL C_F_pointer(InputFileString_C, InputFileString) + IF (InputFilePassed==1_c_int) THEN + TmpFileName = TRIM(OutRootName)//'.BD.tmp' + CALL WritePassedInputFile( InputFileString, TmpFileName, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + InitInp%InputFile = TmpFileName + ELSE + InitInp%InputFile = FileNameFromCString(InputFileString, InputFileStringLength_C) + ENDIF + + !---------------------------------------------------------------------------------------------------------------------------------------------- + ! Set other inputs for calling BD_Init + !---------------------------------------------------------------------------------------------------------------------------------------------- + + ! Check the interpolation order + IF (InterpOrder_C .EQ. 1 .OR. InterpOrder_C .EQ. 2) THEN + InterpOrder = INT(InterpOrder_C, IntKi) + CALL AllocAry( InputTimes, InterpOrder+1, 'InputTimes', ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + ELSE + ErrStat_F2 = ErrID_Fatal + ErrMsg_F2 = 'InterpOrder must be 1 (linear) or 2 (quadratic)' + IF (Failed()) RETURN + END IF + + dT_Global = REAL(DT_C, DbKi) + dT_Interval = dT_Global + N_Global = 0_IntKi ! Assume we are on timestep 0 at start + T_Initial = 0.0_DbKi ! Set from the first UpdateStates call + InputTimePrev = -HUGE(InputTimePrev) ! Initialize for BD_C_UpdateStates (no previous call) + WrOutputs = INT(WrOutputs_C, IntKi) + UnOutFile = -1 + OutTimePrev = -HUGE(OutTimePrev) + + InitInp%RootName = TRIM(OutRootName)//'.BD' ! summary and echo files, same naming as the BeamDyn driver + InitInp%Linearize = .FALSE. + InitInp%DynamicSolve = DynamicSolve_C /= 0_c_int + InitInp%CompAeroMaps = .FALSE. + InitInp%gravity = REAL(Gravity_C, ReKi) + + ! Root position and orientation: the root frame is also the reference frame for the BeamDyn calculations + ! (equivalent to GlbRotBladeT0 = TRUE in the BeamDyn driver). The DCM is passed row by row (C ordering). + InitInp%GlbPos = REAL(RootPos_C, ReKi) + InitInp%GlbRot = TRANSPOSE(RESHAPE(REAL(RootOri_C, R8Ki), (/3,3/))) + InitInp%RootOri = InitInp%GlbRot + InitInp%RootDisp = 0.0_R8Ki + InitInp%RootVel = REAL(RootVel_C, ReKi) + + ALLOCATE(u(InterpOrder+1), STAT=ErrStat_F2) + IF (ErrStat_F2 /= 0) THEN + ErrStat_F2 = ErrID_Fatal + ErrMsg_F2 = 'Failed to allocate Inputs type for BD' + IF (Failed()) RETURN + ENDIF + + !------------------------------------------------- + ! Call the main subroutine BD_Init + !------------------------------------------------- + CALL BD_Init( InitInp, u(1), p, x(STATE_CURR), xd(STATE_CURR), z(STATE_CURR), OtherSt(STATE_CURR), y, m, dT_Interval, InitOutData, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + + ! Remove the temporary input file now that it has been read + IF (LEN_TRIM(TmpFileName) > 0) CALL DeleteFile( TmpFileName ) + + ! The states are advanced by DT_C here, so BeamDyn must use the same timestep (DTBeam in the input file must be DEFAULT or equal to DT_C) + IF ( .NOT. EqualRealNos( p%dt, dT_Global ) ) THEN + ErrStat_F2 = ErrID_Fatal + ErrMsg_F2 = 'The BeamDyn timestep DTBeam ('//TRIM(Num2LStr(p%dt))//' s) must be DEFAULT or equal to the timestep passed to BD_C_Init ('//TRIM(Num2LStr(dT_Global))//' s).' + IF (Failed()) RETURN + ENDIF + + ! If the quasi-static solve is in use, rerun the initialization with the loads at t=0 (the loads are not known during + ! BD_Init). This is handled the same way in the BeamDyn driver. + OtherSt(STATE_CURR)%RunQuasiStaticInit = p%analysis_type == BD_DYN_SSS_ANALYSIS + + !------------------------------------------------- + ! Set mesh size information for the calling code + !------------------------------------------------- + NumOutputNodes = y%BldMotion%Nnodes + NumPointLoadNodes = u(1)%PointLoad%Nnodes + NumDistrLoadNodes = u(1)%DistrLoad%Nnodes + NumOutputNodes_C = INT(NumOutputNodes, c_int) + NumPointLoadNodes_C = INT(NumPointLoadNodes, c_int) + NumDistrLoadNodes_C = INT(NumDistrLoadNodes, c_int) + + !------------------------------------------------- + ! Set output channel information for the calling code + !------------------------------------------------- + + ! Number of channels + NumChannels = SIZE(InitOutData%WriteOutputHdr) + NumChannels_C = INT(NumChannels, c_int) + IF (NumChannels_C > MaxBDOutputs) THEN + ErrStat_F2 = ErrID_Fatal + ErrMsg_F2 = 'Number of output channels ('//TRIM(Num2LStr(NumChannels_C))//') exceeds the maximum this interface can return ('//TRIM(Num2LStr(MaxBDOutputs))//').' + IF (Failed()) RETURN + ENDIF + + ! 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 + END DO + END DO + + ! Null terminate the string + OutputChannelNames_C(K) = C_NULL_CHAR + OutputChannelUnits_C(K) = C_NULL_CHAR + + !------------------------------------------------------------- + ! Output file (same format as the BeamDyn driver output file) + !------------------------------------------------------------- + IF (WrOutputs /= 0_IntKi) THEN + CALL Dvr_InitializeOutputFile( UnOutFile, InitOutData, OutRootName, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + ENDIF + + !------------------------------------------------------------- + ! Copies of the inputs: the input history and the inputs set through the interface + !------------------------------------------------------------- + DO I=2,InterpOrder+1 + CALL BD_CopyInput (u(1), u(I), MESH_NEWCOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + END DO + CALL BD_CopyInput (u(1), BD_u, MESH_NEWCOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + + !------------------------------------------------------------- + ! Initial setup of other pieces of x,xd,z,OtherSt + !------------------------------------------------------------- + CALL BD_CopyContState ( x( STATE_CURR), x( STATE_PRED), MESH_NEWCOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyDiscState ( xd( STATE_CURR), xd( STATE_PRED), MESH_NEWCOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyConstrState( z( STATE_CURR), z( STATE_PRED), MESH_NEWCOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyOtherState ( OtherSt(STATE_CURR), OtherSt(STATE_PRED), MESH_NEWCOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + + !------------------------------------------------------------- + ! Setup the previous timestep copies of states + !------------------------------------------------------------- + CALL BD_CopyContState ( x( STATE_CURR), x( STATE_LAST), MESH_NEWCOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyDiscState ( xd( STATE_CURR), xd( STATE_LAST), MESH_NEWCOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyConstrState( z( STATE_CURR), z( STATE_LAST), MESH_NEWCOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyOtherState ( OtherSt(STATE_CURR), OtherSt(STATE_LAST), MESH_NEWCOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + + !------------------------------------------------- + ! Clean up variables and set up for BD_C_CalcOutput + !------------------------------------------------- + CALL BD_DestroyInitInput( InitInp, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + CALL BD_DestroyInitOutput( InitOutData, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + + Initialized = .TRUE. + + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) + +CONTAINS + LOGICAL FUNCTION Failed() + CALL SetErrStat( ErrStat_F2, ErrMsg_F2, ErrStat_F, ErrMsg_F, RoutineName ) + Failed = ErrStat_F >= AbortErrLev + IF (Failed) THEN + IF (LEN_TRIM(TmpFileName) > 0) CALL DeleteFile( TmpFileName ) + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) + ENDIF + END FUNCTION Failed +END SUBROUTINE BD_C_Init + +!=============================================================================================================== +!---------------------------------------------- BD GET REFERENCE POSITIONS ------------------------------------- +!=============================================================================================================== +!> Return the reference (undeflected) positions and orientations of the output nodes and the reference positions of +!! the load input nodes. Call after BD_C_Init with arrays sized by the node counts it returned. +SUBROUTINE BD_C_GetRefPositions( OutputNodePos_C, OutputNodeOri_C, PointLoadNodePos_C, DistrLoadNodePos_C, ErrStat_C, ErrMsg_C ) BIND (C, NAME='BD_C_GetRefPositions') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: BD_C_GetRefPositions +!GCC$ ATTRIBUTES DLLEXPORT :: BD_C_GetRefPositions +#endif + REAL(C_FLOAT), INTENT( OUT) :: OutputNodePos_C(3*NumOutputNodes) !< Reference positions of the output nodes [x,y,z] (m) + REAL(C_DOUBLE), INTENT( OUT) :: OutputNodeOri_C(9*NumOutputNodes) !< Reference orientations of the output nodes: DCM from the global frame to the node frame, stored row by row + REAL(C_FLOAT), INTENT( OUT) :: PointLoadNodePos_C(3*NumPointLoadNodes) !< Reference positions of the point load nodes [x,y,z] (m) + REAL(C_FLOAT), INTENT( OUT) :: DistrLoadNodePos_C(3*NumDistrLoadNodes) !< Reference positions of the distributed load nodes [x,y,z] (m) + INTEGER(C_INT), INTENT( OUT) :: ErrStat_C + CHARACTER(KIND=C_CHAR), INTENT( OUT) :: ErrMsg_C(ErrMsgLen_C) + + ! Local Variables + INTEGER(IntKi) :: ErrStat_F + CHARACTER(ErrMsgLen) :: ErrMsg_F + INTEGER(IntKi) :: I + CHARACTER(*), PARAMETER :: RoutineName = 'BD_C_GetRefPositions' + + ErrStat_F = ErrID_None + ErrMsg_F = '' + + IF (.NOT. Initialized) THEN + CALL SetErrStat( ErrID_Fatal, 'BD_C_Init must be called before '//RoutineName, ErrStat_F, ErrMsg_F, RoutineName ) + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) + RETURN + ENDIF + + DO I=1,y%BldMotion%Nnodes + OutputNodePos_C(3*(I-1)+1:3*I) = REAL(y%BldMotion%Position(1:3,I), c_float) + OutputNodeOri_C(9*(I-1)+1:9*I) = REAL(RESHAPE(TRANSPOSE(y%BldMotion%RefOrientation(1:3,1:3,I)), (/9/)), c_double) + END DO + DO I=1,u(1)%PointLoad%Nnodes + PointLoadNodePos_C(3*(I-1)+1:3*I) = REAL(u(1)%PointLoad%Position(1:3,I), c_float) + END DO + DO I=1,u(1)%DistrLoad%Nnodes + DistrLoadNodePos_C(3*(I-1)+1:3*I) = REAL(u(1)%DistrLoad%Position(1:3,I), c_float) + END DO + + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) +END SUBROUTINE BD_C_GetRefPositions + +!=============================================================================================================== +!---------------------------------------------- BD SET ROOT MOTION --------------------------------------------- +!=============================================================================================================== +!> Set the root motion input. The values are used by the next BD_C_CalcOutput call (motion at the current time) +!! or the next BD_C_UpdateStates call (motion at the next time). +SUBROUTINE BD_C_SetRootMotion( RootDisp_C, RootOri_C, RootVel_C, RootAcc_C, ErrStat_C, ErrMsg_C ) BIND (C, NAME='BD_C_SetRootMotion') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: BD_C_SetRootMotion +!GCC$ ATTRIBUTES DLLEXPORT :: BD_C_SetRootMotion +#endif + REAL(C_DOUBLE), INTENT(IN ) :: RootDisp_C(3) !< Root translational displacement from the initial root position, global frame (m) + REAL(C_DOUBLE), INTENT(IN ) :: RootOri_C(9) !< Root orientation: DCM from the global frame to the root frame, stored row by row + REAL(C_DOUBLE), INTENT(IN ) :: RootVel_C(6) !< Root translational (1:3) and rotational (4:6) velocity, global frame (m/s, rad/s) + REAL(C_DOUBLE), INTENT(IN ) :: RootAcc_C(6) !< Root translational (1:3) and rotational (4:6) acceleration, global frame (m/s^2, rad/s^2) + INTEGER(C_INT), INTENT( OUT) :: ErrStat_C + CHARACTER(KIND=C_CHAR), INTENT( OUT) :: ErrMsg_C(ErrMsgLen_C) + + ! Local Variables + INTEGER(IntKi) :: ErrStat_F + CHARACTER(ErrMsgLen) :: ErrMsg_F + CHARACTER(*), PARAMETER :: RoutineName = 'BD_C_SetRootMotion' + + ErrStat_F = ErrID_None + ErrMsg_F = '' + + IF (.NOT. Initialized) THEN + CALL SetErrStat( ErrID_Fatal, 'BD_C_Init must be called before '//RoutineName, ErrStat_F, ErrMsg_F, RoutineName ) + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) + RETURN + ENDIF + + BD_u%RootMotion%TranslationDisp(1:3,1) = REAL(RootDisp_C(1:3), ReKi) + BD_u%RootMotion%Orientation(1:3,1:3,1) = TRANSPOSE(RESHAPE(REAL(RootOri_C, R8Ki), (/3,3/))) + BD_u%RootMotion%TranslationVel(1:3,1) = REAL(RootVel_C(1:3), ReKi) + BD_u%RootMotion%RotationVel(1:3,1) = REAL(RootVel_C(4:6), ReKi) + BD_u%RootMotion%TranslationAcc(1:3,1) = REAL(RootAcc_C(1:3), ReKi) + BD_u%RootMotion%RotationAcc(1:3,1) = REAL(RootAcc_C(4:6), ReKi) + + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) +END SUBROUTINE BD_C_SetRootMotion + +!=============================================================================================================== +!---------------------------------------------- BD SET POINT LOADS --------------------------------------------- +!=============================================================================================================== +!> Set the point loads at the finite element nodes (the point load mesh). Loads are in the global frame. +!! The values are used by the next BD_C_CalcOutput or BD_C_UpdateStates call. +SUBROUTINE BD_C_SetPointLoads( PointLoads_C, ErrStat_C, ErrMsg_C ) BIND (C, NAME='BD_C_SetPointLoads') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: BD_C_SetPointLoads +!GCC$ ATTRIBUTES DLLEXPORT :: BD_C_SetPointLoads +#endif + REAL(C_FLOAT), INTENT(IN ) :: PointLoads_C(6*NumPointLoadNodes) !< Point force (1:3) and moment (4:6) at each point load node [Fx,Fy,Fz,Mx,My,Mz] (N, N-m) + INTEGER(C_INT), INTENT( OUT) :: ErrStat_C + CHARACTER(KIND=C_CHAR), INTENT( OUT) :: ErrMsg_C(ErrMsgLen_C) + + ! Local Variables + INTEGER(IntKi) :: ErrStat_F + CHARACTER(ErrMsgLen) :: ErrMsg_F + INTEGER(IntKi) :: I + CHARACTER(*), PARAMETER :: RoutineName = 'BD_C_SetPointLoads' + + ErrStat_F = ErrID_None + ErrMsg_F = '' + + IF (.NOT. Initialized) THEN + CALL SetErrStat( ErrID_Fatal, 'BD_C_Init must be called before '//RoutineName, ErrStat_F, ErrMsg_F, RoutineName ) + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) + RETURN + ENDIF + + DO I=1,BD_u%PointLoad%Nnodes + BD_u%PointLoad%Force( 1:3,I) = REAL(PointLoads_C(6*(I-1)+1:6*(I-1)+3), ReKi) + BD_u%PointLoad%Moment(1:3,I) = REAL(PointLoads_C(6*(I-1)+4:6*(I-1)+6), ReKi) + END DO + + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) +END SUBROUTINE BD_C_SetPointLoads + +!=============================================================================================================== +!---------------------------------------------- BD SET DISTRIBUTED LOADS --------------------------------------- +!=============================================================================================================== +!> Set the distributed loads at the quadrature point nodes (the distributed load mesh). Loads are per unit length +!! in the global frame. The values are used by the next BD_C_CalcOutput or BD_C_UpdateStates call. +SUBROUTINE BD_C_SetDistrLoads( DistrLoads_C, ErrStat_C, ErrMsg_C ) BIND (C, NAME='BD_C_SetDistrLoads') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: BD_C_SetDistrLoads +!GCC$ ATTRIBUTES DLLEXPORT :: BD_C_SetDistrLoads +#endif + REAL(C_FLOAT), INTENT(IN ) :: DistrLoads_C(6*NumDistrLoadNodes) !< Distributed force (1:3) and moment (4:6) at each distributed load node [fx,fy,fz,mx,my,mz] (N/m, N-m/m) + INTEGER(C_INT), INTENT( OUT) :: ErrStat_C + CHARACTER(KIND=C_CHAR), INTENT( OUT) :: ErrMsg_C(ErrMsgLen_C) + + ! Local Variables + INTEGER(IntKi) :: ErrStat_F + CHARACTER(ErrMsgLen) :: ErrMsg_F + INTEGER(IntKi) :: I + CHARACTER(*), PARAMETER :: RoutineName = 'BD_C_SetDistrLoads' + + ErrStat_F = ErrID_None + ErrMsg_F = '' + + IF (.NOT. Initialized) THEN + CALL SetErrStat( ErrID_Fatal, 'BD_C_Init must be called before '//RoutineName, ErrStat_F, ErrMsg_F, RoutineName ) + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) + RETURN + ENDIF + + DO I=1,BD_u%DistrLoad%Nnodes + BD_u%DistrLoad%Force( 1:3,I) = REAL(DistrLoads_C(6*(I-1)+1:6*(I-1)+3), ReKi) + BD_u%DistrLoad%Moment(1:3,I) = REAL(DistrLoads_C(6*(I-1)+4:6*(I-1)+6), ReKi) + END DO + + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) +END SUBROUTINE BD_C_SetDistrLoads + +!=============================================================================================================== +!---------------------------------------------- BD UPDATE STATES ----------------------------------------------- +!=============================================================================================================== +!> This routine updates the states from Time_C to TimeNext_C. The inputs set by the BD_C_Set* routines are taken +!! as the inputs at TimeNext_C. If Time_C repeats the time of the previous call, this is a correction step: the +!! states are reset to the previous time and the update is repeated with the new inputs. +SUBROUTINE BD_C_UpdateStates( Time_C, TimeNext_C, ErrStat_C, ErrMsg_C ) BIND (C, NAME='BD_C_UpdateStates') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: BD_C_UpdateStates +!GCC$ ATTRIBUTES DLLEXPORT :: BD_C_UpdateStates +#endif + REAL(C_DOUBLE), INTENT(IN ) :: Time_C !< Current time (s) + REAL(C_DOUBLE), INTENT(IN ) :: TimeNext_C !< Time to advance the states to (s); must be Time_C + DT_C + INTEGER(C_INT), INTENT( OUT) :: ErrStat_C + CHARACTER(KIND=C_CHAR), INTENT( OUT) :: ErrMsg_C(ErrMsgLen_C) + + ! Local Variables + INTEGER(IntKi) :: ErrStat_F, ErrStat_F2 + CHARACTER(ErrMsgLen) :: ErrMsg_F, ErrMsg_F2 + LOGICAL :: CorrectionStep + CHARACTER(*), PARAMETER :: RoutineName = 'BD_C_UpdateStates' + + ! Set up error handling + ErrStat_F = ErrID_None + ErrMsg_F = '' + CorrectionStep = .FALSE. + + IF (.NOT. Initialized) THEN + ErrStat_F2 = ErrID_Fatal + ErrMsg_F2 = 'BD_C_Init must be called before '//RoutineName + IF (Failed()) RETURN + ENDIF + + IF ( .NOT. EqualRealNos( REAL(TimeNext_C - Time_C, DbKi), dT_Global ) ) THEN + ErrStat_F2 = ErrID_Fatal + ErrMsg_F2 = 'TimeNext_C - Time_C ('//TRIM(Num2LStr(REAL(TimeNext_C - Time_C, DbKi)))//' s) must equal the timestep passed to BD_C_Init ('//TRIM(Num2LStr(dT_Global))//' s).' + IF (Failed()) RETURN + ENDIF + + !------------------------------------------------------- + ! Check the time for current timestep and next timestep + !------------------------------------------------------- + ! These inputs are used in the time stepping algorithm within BD_UpdateStates + ! For quadratic interpolation (InterpOrder==2), 3 timesteps are used. For + ! linear (InterOrder==1), 2 timesteps (the BD 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: the times passed to BD_UpdateStates are set from the global timestep counter + ! and the stored DbKi timestep rather than from the times passed in, so that + ! the input times stay exact multiples of the timestep. + + ! 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 + IF (N_Global == 0_IntKi) T_Initial = REAL(Time_C,DbKi) ! first call sets the initial time + InputTimePrev = REAL(Time_C,DbKi) ! Store for check next time + IF (InterpOrder>1) THEN ! quadratic, so keep the old time + InputTimes(INPUT_LAST) = T_Initial + ( N_Global - 1 ) * dT_Global ! u(3) at T-dT + ENDIF + InputTimes(INPUT_CURR) = T_Initial + N_Global * dT_Global ! u(2) at T + InputTimes(INPUT_PRED) = T_Initial + ( 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 BD_CopyContState (x( STATE_LAST), x( STATE_CURR), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyDiscState (xd( STATE_LAST), xd( STATE_CURR), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyConstrState (z( STATE_LAST), z( STATE_CURR), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyOtherState (OtherSt(STATE_LAST), OtherSt(STATE_CURR), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); 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 BD_CopyInput( u(INPUT_CURR), u(INPUT_LAST), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + END IF + ! Move inputs from previous t+dt (now t) to t + CALL BD_CopyInput( u(INPUT_PRED), u(INPUT_CURR), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + END IF + + ! Copy the new inputs (set by the BD_C_Set* routines) for time u(INPUT_PRED) + CALL BD_CopyInput( BD_u, u(INPUT_PRED), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2 ); 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 BD_CopyContState (x( STATE_CURR), x( STATE_PRED), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyDiscState (xd( STATE_CURR), xd( STATE_PRED), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyConstrState (z( STATE_CURR), z( STATE_PRED), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyOtherState (OtherSt(STATE_CURR), OtherSt(STATE_PRED), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + + !------------------------------------------------- + ! Call the main subroutine BD_UpdateStates + ! -- the step number passed is the index of the step at time T (N_Global was incremented above), as in the BeamDyn driver + !------------------------------------------------- + CALL BD_UpdateStates( InputTimes(INPUT_CURR), N_Global-1, u, InputTimes, p, x(STATE_PRED), xd(STATE_PRED), z(STATE_PRED), OtherSt(STATE_PRED), m, ErrStat_F2, ErrMsg_F2 ); 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 BD_CopyContState (x( STATE_CURR), x( STATE_LAST), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyDiscState (xd( STATE_CURR), xd( STATE_LAST), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyConstrState (z( STATE_CURR), z( STATE_LAST), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyOtherState (OtherSt(STATE_CURR), OtherSt(STATE_LAST), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); 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 BD_CopyContState (x( STATE_PRED), x( STATE_CURR), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyDiscState (xd( STATE_PRED), xd( STATE_CURR), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyConstrState (z( STATE_PRED), z( STATE_CURR), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyOtherState (OtherSt(STATE_PRED), OtherSt(STATE_CURR), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) + +CONTAINS + LOGICAL FUNCTION Failed() + CALL SetErrStat( ErrStat_F2, ErrMsg_F2, ErrStat_F, ErrMsg_F, RoutineName ) + Failed = ErrStat_F >= AbortErrLev + IF (Failed) CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) + END FUNCTION Failed +END SUBROUTINE BD_C_UpdateStates + +!=============================================================================================================== +!---------------------------------------------- BD CALC OUTPUT ------------------------------------------------- +!=============================================================================================================== +!> Calculate the BeamDyn outputs at Time_C from the current states and the inputs set by the BD_C_Set* routines. +!! Node motions are returned for every node of the blade motion output mesh. +SUBROUTINE BD_C_CalcOutput( Time_C, NodePos_C, NodeOri_C, NodeVel_C, NodeAcc_C, RootReaction_C, OutputChannelValues_C, ErrStat_C, ErrMsg_C ) BIND (C, NAME='BD_C_CalcOutput') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: BD_C_CalcOutput +!GCC$ ATTRIBUTES DLLEXPORT :: BD_C_CalcOutput +#endif + REAL(C_DOUBLE), INTENT(IN ) :: Time_C !< Current time (s) + REAL(C_FLOAT), INTENT( OUT) :: NodePos_C(3*NumOutputNodes) !< Position of each output node in the global frame [x,y,z] (m) + REAL(C_DOUBLE), INTENT( OUT) :: NodeOri_C(9*NumOutputNodes) !< Orientation of each output node: DCM from the global frame to the node frame, stored row by row + REAL(C_FLOAT), INTENT( OUT) :: NodeVel_C(6*NumOutputNodes) !< Translational (1:3) and rotational (4:6) velocity of each output node, global frame (m/s, rad/s) + REAL(C_FLOAT), INTENT( OUT) :: NodeAcc_C(6*NumOutputNodes) !< Translational (1:3) and rotational (4:6) acceleration of each output node, global frame (m/s^2, rad/s^2) + REAL(C_FLOAT), INTENT( OUT) :: RootReaction_C(6) !< Reaction force (1:3) and moment (4:6) at the root, global frame (N, N-m) + REAL(C_FLOAT), INTENT( OUT) :: OutputChannelValues_C(NumChannels) !< Output channel values + INTEGER(C_INT), INTENT( OUT) :: ErrStat_C + CHARACTER(KIND=C_CHAR), INTENT( OUT) :: ErrMsg_C(ErrMsgLen_C) + + ! Local Variables + REAL(DbKi) :: t + INTEGER(IntKi) :: ErrStat_F, ErrStat_F2 + CHARACTER(ErrMsgLen) :: ErrMsg_F, ErrMsg_F2 + INTEGER(IntKi) :: I + CHARACTER(*), PARAMETER :: RoutineName = 'BD_C_CalcOutput' + + ! Set up error handling + ErrStat_F = ErrID_None + ErrMsg_F = '' + + IF (.NOT. Initialized) THEN + ErrStat_F2 = ErrID_Fatal + ErrMsg_F2 = 'BD_C_Init must be called before '//RoutineName + IF (Failed()) RETURN + ENDIF + + ! Time + t = REAL(Time_C, DbKi) + + ! Copy the inputs set by the BD_C_Set* routines into the input at the current time + CALL BD_CopyInput( BD_u, u(1), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + + !------------------------------------------------- + ! Call the main subroutine BD_CalcOutput + !------------------------------------------------- + CALL BD_CalcOutput( t, u(1), p, x(STATE_CURR), xd(STATE_CURR), z(STATE_CURR), OtherSt(STATE_CURR), y, m, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + + !------------------------------------------------- + ! Convert the outputs of BD_CalcOutput back to C + !------------------------------------------------- + DO I=1,y%BldMotion%Nnodes + NodePos_C(3*(I-1)+1:3*I) = REAL(y%BldMotion%Position(1:3,I) + y%BldMotion%TranslationDisp(1:3,I), c_float) + NodeOri_C(9*(I-1)+1:9*I) = REAL(RESHAPE(TRANSPOSE(y%BldMotion%Orientation(1:3,1:3,I)), (/9/)), c_double) + NodeVel_C(6*(I-1)+1:6*(I-1)+3) = REAL(y%BldMotion%TranslationVel(1:3,I), c_float) + NodeVel_C(6*(I-1)+4:6*(I-1)+6) = REAL(y%BldMotion%RotationVel(1:3,I), c_float) + NodeAcc_C(6*(I-1)+1:6*(I-1)+3) = REAL(y%BldMotion%TranslationAcc(1:3,I), c_float) + NodeAcc_C(6*(I-1)+4:6*(I-1)+6) = REAL(y%BldMotion%RotationAcc(1:3,I), c_float) + END DO + + RootReaction_C(1:3) = REAL(y%ReactionForce%Force( 1:3,1), c_float) + RootReaction_C(4:6) = REAL(y%ReactionForce%Moment(1:3,1), c_float) + + OutputChannelValues_C = REAL(y%WriteOutput, c_float) + + !------------------------------------------------- + ! Write the output file line (only once per time; a repeated time is a correction step) + !------------------------------------------------- + IF (UnOutFile > 0) THEN + IF ( .NOT. EqualRealNos( t, OutTimePrev ) ) THEN + CALL Dvr_WriteOutputLine( t, UnOutFile, p%OutFmt, y ) + OutTimePrev = t + ENDIF + ENDIF + + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) + +CONTAINS + LOGICAL FUNCTION Failed() + CALL SetErrStat( ErrStat_F2, ErrMsg_F2, ErrStat_F, ErrMsg_F, RoutineName ) + Failed = ErrStat_F >= AbortErrLev + IF (Failed) CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) + END FUNCTION Failed +END SUBROUTINE BD_C_CalcOutput + +!=============================================================================================================== +!---------------------------------------------- BD PACK STATES ------------------------------------------------- +!=============================================================================================================== +!> Write a checkpoint file (.chkp) with the complete state of the beam: the continuous, discrete, +!! constraint, and other states at the current and previous time, the input history, the inputs set through the +!! interface, and the time stepping information. BeamDyn's registry pack routines are used, so the file is only +!! valid for the same build and the same input file. Restore with BD_C_Init (same input file) followed by +!! BD_C_UnpackStates. +SUBROUTINE BD_C_PackStates( CheckpointRoot_C, ErrStat_C, ErrMsg_C ) BIND (C, NAME='BD_C_PackStates') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: BD_C_PackStates +!GCC$ ATTRIBUTES DLLEXPORT :: BD_C_PackStates +#endif + CHARACTER(KIND=C_CHAR), INTENT(IN ) :: CheckpointRoot_C(*) !< Root name of the checkpoint file (C_NULL_CHAR terminated, at most IntfStrLen characters); ".chkp" is appended + INTEGER(C_INT), INTENT( OUT) :: ErrStat_C + CHARACTER(KIND=C_CHAR), INTENT( OUT) :: ErrMsg_C(ErrMsgLen_C) + + ! Local Variables + TYPE(RegFile) :: RF + CHARACTER(IntfStrLen) :: FileName + INTEGER(IntKi) :: UnOut + INTEGER(IntKi) :: ErrStat_F, ErrStat_F2 + CHARACTER(ErrMsgLen) :: ErrMsg_F, ErrMsg_F2 + INTEGER(IntKi) :: I + CHARACTER(*), PARAMETER :: RoutineName = 'BD_C_PackStates' + + ErrStat_F = ErrID_None + ErrMsg_F = '' + UnOut = -1 + + IF (.NOT. Initialized) THEN + ErrStat_F2 = ErrID_Fatal + ErrMsg_F2 = 'BD_C_Init must be called before '//RoutineName + IF (Failed()) RETURN + ENDIF + + FileName = TRIM(CStringToFortran( CheckpointRoot_C ))//'.chkp' + + CALL GetNewUnit( UnOut, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + CALL OpenBOutFile( UnOut, FileName, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + + ! Checkpoint file header: file identifier, sizes used to check the file against the current beam, and time stepping information + WRITE (UnOut, IOSTAT=ErrStat_F2) CheckpointFileID + WRITE (UnOut, IOSTAT=ErrStat_F2) InterpOrder + WRITE (UnOut, IOSTAT=ErrStat_F2) p%node_total + WRITE (UnOut, IOSTAT=ErrStat_F2) p%dof_total + WRITE (UnOut, IOSTAT=ErrStat_F2) y%BldMotion%Nnodes + WRITE (UnOut, IOSTAT=ErrStat_F2) u(1)%DistrLoad%Nnodes + WRITE (UnOut, IOSTAT=ErrStat_F2) N_Global + WRITE (UnOut, IOSTAT=ErrStat_F2) dT_Global + WRITE (UnOut, IOSTAT=ErrStat_F2) T_Initial + WRITE (UnOut, IOSTAT=ErrStat_F2) InputTimePrev + WRITE (UnOut, IOSTAT=ErrStat_F2) OutTimePrev + WRITE (UnOut, IOSTAT=ErrStat_F2) InputTimes + IF (ErrStat_F2 /= 0) THEN + ErrStat_F2 = ErrID_Fatal + ErrMsg_F2 = 'Error writing the header of checkpoint file '//TRIM(FileName) + IF (Failed()) RETURN + ENDIF + + ! Pack the states and inputs into the registry file + CALL InitRegFile( RF, UnOut, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + + DO I=STATE_LAST,STATE_CURR + CALL BD_PackContState ( RF, x(I) ) + CALL BD_PackDiscState ( RF, xd(I) ) + CALL BD_PackConstrState( RF, z(I) ) + CALL BD_PackOtherState ( RF, OtherSt(I) ) + END DO + DO I=1,InterpOrder+1 + CALL BD_PackInput( RF, u(I) ) + END DO + CALL BD_PackInput( RF, BD_u ) + + ! Close registry file and get any errors that occurred while writing (this also closes the unit) + CALL CloseRegFile( RF, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + UnOut = -1 + + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) + +CONTAINS + LOGICAL FUNCTION Failed() + CALL SetErrStat( ErrStat_F2, ErrMsg_F2, ErrStat_F, ErrMsg_F, RoutineName ) + Failed = ErrStat_F >= AbortErrLev + IF (Failed) THEN + IF (UnOut > 0) CLOSE(UnOut) + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) + ENDIF + END FUNCTION Failed +END SUBROUTINE BD_C_PackStates + +!=============================================================================================================== +!---------------------------------------------- BD UNPACK STATES ----------------------------------------------- +!=============================================================================================================== +!> Restore the state of the beam from a checkpoint file written by BD_C_PackStates. BD_C_Init must have been +!! called with the same input file first; the unpacked values replace the states, input history, and time +!! stepping information set by BD_C_Init. +SUBROUTINE BD_C_UnpackStates( CheckpointRoot_C, ErrStat_C, ErrMsg_C ) BIND (C, NAME='BD_C_UnpackStates') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: BD_C_UnpackStates +!GCC$ ATTRIBUTES DLLEXPORT :: BD_C_UnpackStates +#endif + CHARACTER(KIND=C_CHAR), INTENT(IN ) :: CheckpointRoot_C(*) !< Root name of the checkpoint file (C_NULL_CHAR terminated, at most IntfStrLen characters); ".chkp" is appended + INTEGER(C_INT), INTENT( OUT) :: ErrStat_C + CHARACTER(KIND=C_CHAR), INTENT( OUT) :: ErrMsg_C(ErrMsgLen_C) + + ! Local Variables + TYPE(RegFile) :: RF + TYPE(BD_InputType) :: u_tmp !< Inputs read from the file (copied into the existing input meshes so mesh siblings stay intact) + CHARACTER(IntfStrLen) :: FileName + INTEGER(IntKi) :: UnIn + INTEGER(IntKi) :: FileID, InterpOrderIn, node_total, dof_total, NnodesOut, NnodesDistr + INTEGER(IntKi) :: ErrStat_F, ErrStat_F2 + CHARACTER(ErrMsgLen) :: ErrMsg_F, ErrMsg_F2 + INTEGER(IntKi) :: I + CHARACTER(*), PARAMETER :: RoutineName = 'BD_C_UnpackStates' + + ErrStat_F = ErrID_None + ErrMsg_F = '' + UnIn = -1 + + IF (.NOT. Initialized) THEN + ErrStat_F2 = ErrID_Fatal + ErrMsg_F2 = 'BD_C_Init must be called before '//RoutineName + IF (Failed()) RETURN + ENDIF + + FileName = TRIM(CStringToFortran( CheckpointRoot_C ))//'.chkp' + + CALL GetNewUnit( UnIn, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + CALL OpenBInpFile( UnIn, FileName, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + + ! Checkpoint file header + READ (UnIn, IOSTAT=ErrStat_F2) FileID + IF (ErrStat_F2 /= 0 .OR. FileID /= CheckpointFileID) THEN + ErrStat_F2 = ErrID_Fatal + ErrMsg_F2 = TRIM(FileName)//' is not a BeamDyn checkpoint file written by BD_C_PackStates.' + IF (Failed()) RETURN + ENDIF + READ (UnIn, IOSTAT=ErrStat_F2) InterpOrderIn + READ (UnIn, IOSTAT=ErrStat_F2) node_total + READ (UnIn, IOSTAT=ErrStat_F2) dof_total + READ (UnIn, IOSTAT=ErrStat_F2) NnodesOut + READ (UnIn, IOSTAT=ErrStat_F2) NnodesDistr + IF (ErrStat_F2 /= 0) THEN + ErrStat_F2 = ErrID_Fatal + ErrMsg_F2 = 'Error reading the header of checkpoint file '//TRIM(FileName) + IF (Failed()) RETURN + ENDIF + IF ( InterpOrderIn /= InterpOrder .OR. node_total /= p%node_total .OR. dof_total /= p%dof_total .OR. & + NnodesOut /= y%BldMotion%Nnodes .OR. NnodesDistr /= u(1)%DistrLoad%Nnodes ) THEN + ErrStat_F2 = ErrID_Fatal + ErrMsg_F2 = 'Checkpoint file '//TRIM(FileName)//' was written for a different beam discretization or interpolation order than the current one.' + IF (Failed()) RETURN + ENDIF + READ (UnIn, IOSTAT=ErrStat_F2) N_Global + READ (UnIn, IOSTAT=ErrStat_F2) dT_Global + READ (UnIn, IOSTAT=ErrStat_F2) T_Initial + READ (UnIn, IOSTAT=ErrStat_F2) InputTimePrev + READ (UnIn, IOSTAT=ErrStat_F2) OutTimePrev + READ (UnIn, IOSTAT=ErrStat_F2) InputTimes + IF (ErrStat_F2 /= 0) THEN + ErrStat_F2 = ErrID_Fatal + ErrMsg_F2 = 'Error reading the header of checkpoint file '//TRIM(FileName) + IF (Failed()) RETURN + ENDIF + + ! Unpack the states and inputs from the registry file + CALL OpenRegFile( RF, UnIn, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + + DO I=STATE_LAST,STATE_CURR + CALL BD_UnPackContState ( RF, x(I) ) + CALL BD_UnPackDiscState ( RF, xd(I) ) + CALL BD_UnPackConstrState( RF, z(I) ) + CALL BD_UnPackOtherState ( RF, OtherSt(I) ) + END DO + DO I=1,InterpOrder+1 + CALL BD_UnPackInput( RF, u_tmp ) + CALL BD_CopyInput( u_tmp, u(I), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + CALL BD_DestroyInput( u_tmp, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + END DO + CALL BD_UnPackInput( RF, u_tmp ) + CALL BD_CopyInput( u_tmp, BD_u, MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + CALL BD_DestroyInput( u_tmp, ErrStat_F2, ErrMsg_F2 ); IF (Failed()) RETURN + + ErrStat_F2 = RF%ErrStat + ErrMsg_F2 = RF%ErrMsg + IF (Failed()) RETURN + + CLOSE(UnIn) + UnIn = -1 + IF (ALLOCATED(RF%Pointers)) DEALLOCATE(RF%Pointers) + + ! The restored states are the states at the current time; the predicted states are scratch for the next update + CALL BD_CopyContState (x( STATE_CURR), x( STATE_PRED), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyDiscState (xd( STATE_CURR), xd( STATE_PRED), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyConstrState (z( STATE_CURR), z( STATE_PRED), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + CALL BD_CopyOtherState (OtherSt(STATE_CURR), OtherSt(STATE_PRED), MESH_UPDATECOPY, ErrStat_F2, ErrMsg_F2); IF (Failed()) RETURN + + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) + +CONTAINS + LOGICAL FUNCTION Failed() + CALL SetErrStat( ErrStat_F2, ErrMsg_F2, ErrStat_F, ErrMsg_F, RoutineName ) + Failed = ErrStat_F >= AbortErrLev + IF (Failed) THEN + IF (UnIn > 0) CLOSE(UnIn) + IF (ALLOCATED(RF%Pointers)) DEALLOCATE(RF%Pointers) + CALL BD_DestroyInput( u_tmp, ErrStat_F2, ErrMsg_F2 ) + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) + ENDIF + END FUNCTION Failed +END SUBROUTINE BD_C_UnpackStates + +!=============================================================================================================== +!----------------------------------------------- BD END -------------------------------------------------------- +!=============================================================================================================== +!> Cleanup memory +!! NOTE: the error handling here is slightly different than in other routines +SUBROUTINE BD_C_End(ErrStat_C,ErrMsg_C) BIND (C, NAME='BD_C_End') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: BD_C_End +!GCC$ ATTRIBUTES DLLEXPORT :: BD_C_End +#endif + INTEGER(C_INT), INTENT( OUT) :: ErrStat_C + CHARACTER(KIND=C_CHAR), INTENT( OUT) :: ErrMsg_C(ErrMsgLen_C) + + ! Local variables + INTEGER(IntKi) :: ErrStat_F, ErrStat_F2 + CHARACTER(ErrMsgLen) :: ErrMsg_F, ErrMsg_F2 + CHARACTER(*), PARAMETER :: RoutineName = 'BD_C_End' + + ! Set up error handling for BD_C_End + ErrStat_F = ErrID_None + ErrMsg_F = '' + + ! Close the output file + IF (UnOutFile > 0) CLOSE(UnOutFile) + UnOutFile = -1 + + ! Call the main subroutine BD_End + ! If u is not allocated, then we didn't get far at all in initialization, + ! or BD_C_End got called before Init. We don't want a segfault, so check + ! for allocation. + IF (ALLOCATED(u)) THEN + CALL BD_End( u(1), p, x(STATE_CURR), xd(STATE_CURR), z(STATE_CURR), OtherSt(STATE_CURR), y, m, ErrStat_F2, ErrMsg_F2 ) + CALL SetErrStat( ErrStat_F2, ErrMsg_F2, ErrStat_F, ErrMsg_F, RoutineName ) + ENDIF + + ! NOTE: BD_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) + CALL DestroyAll( ErrStat_F2, ErrMsg_F2 ) + CALL SetErrStat( ErrStat_F2, ErrMsg_F2, ErrStat_F, ErrMsg_F, RoutineName ) + + Initialized = .FALSE. + + CALL SetErrStat_F2C(ErrStat_F,ErrMsg_F,ErrStat_C,ErrMsg_C) +END SUBROUTINE BD_C_End + +!=============================================================================================================== +!----------------------------------------- ADDITIONAL SUBROUTINES ---------------------------------------------- +!=============================================================================================================== +!> Destroy all module-level BeamDyn data (safe to call on data that was never allocated or was already destroyed) +SUBROUTINE DestroyAll( ErrStat, ErrMsg ) + INTEGER(IntKi), INTENT( OUT) :: ErrStat + CHARACTER(ErrMsgLen), INTENT( OUT) :: ErrMsg + INTEGER(IntKi) :: ErrStat2 + CHARACTER(ErrMsgLen) :: ErrMsg2 + INTEGER(IntKi) :: I + CHARACTER(*), PARAMETER :: RoutineName = 'DestroyAll' + + ErrStat = ErrID_None + ErrMsg = '' + + IF (ALLOCATED(u)) THEN + DO I=1,SIZE(u) + CALL BD_DestroyInput( u(I), ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + END DO + DEALLOCATE(u) + END IF + CALL BD_DestroyInput( BD_u, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL BD_DestroyParam( p, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + DO I=STATE_LAST,STATE_PRED + CALL BD_DestroyContState( x(I), ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL BD_DestroyDiscState( xd(I), ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL BD_DestroyConstrState( z(I), ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL BD_DestroyOtherState( OtherSt(I), ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + END DO + CALL BD_DestroyOutput( y, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL BD_DestroyMisc( m, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL BD_DestroyInitInput( InitInp, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL BD_DestroyInitOutput( InitOutData, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + + IF (ALLOCATED(InputTimes)) DEALLOCATE(InputTimes) + + NumOutputNodes = 0 + NumPointLoadNodes = 0 + NumDistrLoadNodes = 0 + NumChannels = 0 +END SUBROUTINE DestroyAll + +!--------------------------------------------------------------------------------------------------------------- +!> Write the input file contents passed through the interface (lines delineated by C_NULL_CHAR) to a file that +!! BeamDyn can read. +SUBROUTINE WritePassedInputFile( InputFileString, FileName, ErrStat, ErrMsg ) + CHARACTER(*), INTENT(IN ) :: InputFileString !< Input file as a single string with NULL character separating lines + CHARACTER(*), INTENT(IN ) :: FileName !< File to write + INTEGER(IntKi), INTENT( OUT) :: ErrStat + CHARACTER(ErrMsgLen), INTENT( OUT) :: ErrMsg + INTEGER(IntKi) :: Un + INTEGER(IntKi) :: IStart, IEnd, IOS + INTEGER(IntKi) :: ErrStat2 + CHARACTER(ErrMsgLen) :: ErrMsg2 + CHARACTER(*), PARAMETER :: RoutineName = 'WritePassedInputFile' + + ErrStat = ErrID_None + ErrMsg = '' + + CALL GetNewUnit( Un, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL OpenFOutFile( Un, FileName, ErrStat2, ErrMsg2 ); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + IF (ErrStat >= AbortErrLev) RETURN + + IStart = 1 + DO WHILE (IStart <= LEN(InputFileString)) + IEnd = INDEX(InputFileString(IStart:), C_NULL_CHAR) + IF (IEnd == 0) THEN + IEnd = LEN(InputFileString) + ELSE + IEnd = IStart + IEnd - 2 ! character before the C_NULL_CHAR + ENDIF + IF (IEnd >= IStart) THEN + WRITE (Un, '(A)', IOSTAT=IOS) InputFileString(IStart:IEnd) + ELSE + WRITE (Un, '(A)', IOSTAT=IOS) '' + ENDIF + IF (IOS /= 0) THEN + CALL SetErrStat( ErrID_Fatal, 'Error writing passed input file contents to '//TRIM(FileName), ErrStat, ErrMsg, RoutineName ) + CLOSE(Un) + RETURN + ENDIF + IStart = IEnd + 2 + END DO + + CLOSE(Un) +END SUBROUTINE WritePassedInputFile + +!--------------------------------------------------------------------------------------------------------------- +!> Convert a C_NULL_CHAR terminated C string of at most IntfStrLen characters to a Fortran string +FUNCTION CStringToFortran( String_C ) RESULT( String_F ) + CHARACTER(KIND=C_CHAR), INTENT(IN ) :: String_C(*) + CHARACTER(IntfStrLen) :: String_F + INTEGER(IntKi) :: I + String_F = '' + DO I=1,IntfStrLen + IF (String_C(I) == C_NULL_CHAR) EXIT + String_F(I:I) = String_C(I) + END DO +END FUNCTION CStringToFortran + +!--------------------------------------------------------------------------------------------------------------- +!> Delete a file (ignoring any errors) +SUBROUTINE DeleteFile( FileName ) + CHARACTER(*), INTENT(IN ) :: FileName + INTEGER(IntKi) :: Un + INTEGER(IntKi) :: IOS + INTEGER(IntKi) :: ErrStat2 + CHARACTER(ErrMsgLen) :: ErrMsg2 + + CALL GetNewUnit( Un, ErrStat2, ErrMsg2 ) + OPEN( UNIT=Un, FILE=TRIM(FileName), STATUS='OLD', IOSTAT=IOS ) + IF (IOS == 0) CLOSE( UNIT=Un, STATUS='DELETE', IOSTAT=IOS ) +END SUBROUTINE DeleteFile + +END MODULE BeamDyn_C diff --git a/modules/beamdyn/src/BeamDyn_C_Binding.h b/modules/beamdyn/src/BeamDyn_C_Binding.h new file mode 100644 index 0000000000..6796c834d8 --- /dev/null +++ b/modules/beamdyn/src/BeamDyn_C_Binding.h @@ -0,0 +1,206 @@ +/********************************************************************************************************************************** + * LICENSING + * Copyright (C) 2026 National Renewable Energy Laboratory + * + * This file is part of BeamDyn. + * + * 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. + * + **********************************************************************************************************************************/ +/* + * C interface to a standalone BeamDyn beam (see BeamDyn_C_Binding.f90). + * + * All routines follow the Fortran BIND(C) conventions used by the OpenFAST module C bindings: every argument is + * passed by reference, strings are C_NULL_CHAR terminated char arrays, and matrices are flat arrays. Direction + * cosine matrices (DCMs) follow the OpenFAST convention (they map a vector from the global frame to the local + * frame) and are stored row by row: dcm[3*i + j] is row i, column j. Node arrays are stored node by node: + * pos[3*n + k] is component k of node n. + * + * Precision: root kinematics (position, orientation, velocity, acceleration) and gravity are double; loads, + * node motions, reaction loads, and channel values are float; orientation DCMs are always double. + * + * Error handling: each routine returns ErrStat_C (0 none, 1 info, 2 warning, 3 severe, 4 fatal) and a + * C_NULL_CHAR terminated message in ErrMsg_C, which must be at least BD_C_ERRMSGLEN characters long. After a + * fatal error the module is left in a state where BD_C_End can still be called safely. + * + * Call sequence: + * BD_C_Init initialize the beam from an input file (returns node counts and output channel info) + * BD_C_GetRefPositions (optional) reference positions of the output nodes and load nodes + * loop over time: + * BD_C_SetRootMotion root motion at the time of the next BD_C_CalcOutput or BD_C_UpdateStates call + * BD_C_SetPointLoads (optional) point loads at the finite element nodes + * BD_C_SetDistrLoads (optional) distributed loads at the quadrature point nodes + * BD_C_CalcOutput outputs at the current time (node motions, root reaction, output channels) + * BD_C_UpdateStates advance the states from Time_C to TimeNext_C using the inputs set for TimeNext_C + * BD_C_PackStates / BD_C_UnpackStates (optional) write/read a checkpoint file with the complete beam state + * BD_C_End + */ +#ifndef BEAMDYN_C_BINDING_H +#define BEAMDYN_C_BINDING_H + +#ifdef __cplusplus +extern "C" { +#endif + +/* Sizes used by the interface (must match the NWTC Library and BeamDyn_C_Binding.f90) */ +#define BD_C_ERRMSGLEN 8197 /* ErrMsgLen + 1: minimum length of the ErrMsg_C buffer */ +#define BD_C_INTFSTRLEN 1025 /* Maximum length of file name strings passed through the interface */ +#define BD_C_CHANLEN 20 /* Length of each output channel name and unit string */ +#define BD_C_MAXOUTPUTS 8000 /* Maximum number of output channels the interface can return */ +#define BD_C_CHANNELBUFLEN (BD_C_CHANLEN*BD_C_MAXOUTPUTS+1) /* Minimum length of the channel name and unit buffers */ + +/* Error status values (ErrStat_C) */ +#define BD_C_ERRID_NONE 0 +#define BD_C_ERRID_INFO 1 +#define BD_C_ERRID_WARN 2 +#define BD_C_ERRID_SEVERE 3 +#define BD_C_ERRID_FATAL 4 + +/* + * Initialize a BeamDyn beam from a BeamDyn primary input file (the blade file it names is read from the file + * system relative to the primary input file). The root of the beam is placed at RootPos_C with orientation + * RootOri_C, which is also the reference frame for the BeamDyn calculations. On return, the node counts give the + * sizes of the arrays used by the other routines. + */ +void BD_C_Init( + const int *InputFilePassed, /* IN: 1: InputFileString_C holds the contents of the input file (lines separated by '\0'); + 0: InputFileString_C holds the path to the input file */ + const char **InputFileString_C, /* IN: pointer to the input file contents or path string (see InputFilePassed) */ + const int *InputFileStringLength_C, /* IN: length of InputFileString_C */ + const char *OutRootName_C, /* IN: root name for the summary (.BD.sum), echo (.BD.ech), and output (.out) files */ + const double *RootPos_C, /* IN: [3] initial root position in the global frame (m) */ + const double *RootOri_C, /* IN: [9] initial root orientation DCM, global to root frame, row by row */ + const double *RootVel_C, /* IN: [6] initial root translational (0:2) and rotational (3:5) velocity, global frame (m/s, rad/s) */ + const double *Gravity_C, /* IN: [3] gravitational acceleration vector, global frame (m/s^2) */ + const double *DT_C, /* IN: time step for BD_C_UpdateStates (s); the input file DTBeam must be DEFAULT or equal to this */ + const int *InterpOrder_C, /* IN: input interpolation/extrapolation order: 1 (linear) or 2 (quadratic) */ + const int *DynamicSolve_C, /* IN: 1: dynamic solve; 0: static solve */ + const int *WrOutputs_C, /* IN: 1: write the output channels to .out at each BD_C_CalcOutput call; 0: no file */ + int *NumOutputNodes_C, /* OUT: number of nodes on the blade motion output mesh */ + int *NumPointLoadNodes_C, /* OUT: number of nodes on the point load input mesh (finite element nodes) */ + int *NumDistrLoadNodes_C, /* OUT: number of nodes on the distributed load input mesh (quadrature point nodes) */ + int *NumChannels_C, /* OUT: number of output channels */ + char *OutputChannelNames_C, /* OUT: [BD_C_CHANNELBUFLEN] channel names, BD_C_CHANLEN characters each, '\0' terminated */ + char *OutputChannelUnits_C, /* OUT: [BD_C_CHANNELBUFLEN] channel units, BD_C_CHANLEN characters each, '\0' terminated */ + int *ErrStat_C, /* OUT: error status */ + char *ErrMsg_C /* OUT: [BD_C_ERRMSGLEN] error message */ +); + +/* + * Reference (undeflected) positions of the output nodes and load input nodes. Arrays are sized by the counts + * returned from BD_C_Init. + */ +void BD_C_GetRefPositions( + float *OutputNodePos_C, /* OUT: [3*NumOutputNodes] reference positions of the output nodes (m) */ + double *OutputNodeOri_C, /* OUT: [9*NumOutputNodes] reference orientation DCMs of the output nodes, global to node frame, row by row */ + float *PointLoadNodePos_C, /* OUT: [3*NumPointLoadNodes] reference positions of the point load nodes (m) */ + float *DistrLoadNodePos_C, /* OUT: [3*NumDistrLoadNodes] reference positions of the distributed load nodes (m) */ + int *ErrStat_C, /* OUT: error status */ + char *ErrMsg_C /* OUT: [BD_C_ERRMSGLEN] error message */ +); + +/* + * Set the root motion used by the next BD_C_CalcOutput (motion at the current time) or BD_C_UpdateStates (motion + * at the next time) call. + */ +void BD_C_SetRootMotion( + const double *RootDisp_C, /* IN: [3] root translational displacement from the initial root position, global frame (m) */ + const double *RootOri_C, /* IN: [9] root orientation DCM, global to root frame, row by row */ + const double *RootVel_C, /* IN: [6] root translational (0:2) and rotational (3:5) velocity, global frame (m/s, rad/s) */ + const double *RootAcc_C, /* IN: [6] root translational (0:2) and rotational (3:5) acceleration, global frame (m/s^2, rad/s^2) */ + int *ErrStat_C, /* OUT: error status */ + char *ErrMsg_C /* OUT: [BD_C_ERRMSGLEN] error message */ +); + +/* + * Set the point loads at the finite element nodes (global frame), used by the next BD_C_CalcOutput or + * BD_C_UpdateStates call. + */ +void BD_C_SetPointLoads( + const float *PointLoads_C, /* IN: [6*NumPointLoadNodes] force (0:2) and moment (3:5) at each point load node (N, N-m) */ + int *ErrStat_C, /* OUT: error status */ + char *ErrMsg_C /* OUT: [BD_C_ERRMSGLEN] error message */ +); + +/* + * Set the distributed loads per unit length at the quadrature point nodes (global frame), used by the next + * BD_C_CalcOutput or BD_C_UpdateStates call. + */ +void BD_C_SetDistrLoads( + const float *DistrLoads_C, /* IN: [6*NumDistrLoadNodes] force (0:2) and moment (3:5) per unit length at each distributed load node (N/m, N-m/m) */ + int *ErrStat_C, /* OUT: error status */ + char *ErrMsg_C /* OUT: [BD_C_ERRMSGLEN] error message */ +); + +/* + * Advance the states from Time_C to TimeNext_C (= Time_C + DT_C). The inputs set by the BD_C_Set* routines are + * taken as the inputs at TimeNext_C. Calling again with the same Time_C repeats the step as a correction step with + * the new inputs. + */ +void BD_C_UpdateStates( + const double *Time_C, /* IN: current time (s) */ + const double *TimeNext_C, /* IN: time to advance the states to (s) */ + int *ErrStat_C, /* OUT: error status */ + char *ErrMsg_C /* OUT: [BD_C_ERRMSGLEN] error message */ +); + +/* + * Compute the outputs at Time_C from the current states and the inputs set by the BD_C_Set* routines. + */ +void BD_C_CalcOutput( + const double *Time_C, /* IN: current time (s) */ + float *NodePos_C, /* OUT: [3*NumOutputNodes] position of each output node, global frame (m) */ + double *NodeOri_C, /* OUT: [9*NumOutputNodes] orientation DCM of each output node, global to node frame, row by row */ + float *NodeVel_C, /* OUT: [6*NumOutputNodes] translational (0:2) and rotational (3:5) velocity of each output node, global frame (m/s, rad/s) */ + float *NodeAcc_C, /* OUT: [6*NumOutputNodes] translational (0:2) and rotational (3:5) acceleration of each output node, global frame (m/s^2, rad/s^2) */ + float *RootReaction_C, /* OUT: [6] reaction force (0:2) and moment (3:5) at the root, global frame (N, N-m) */ + float *OutputChannelValues_C, /* OUT: [NumChannels] output channel values */ + int *ErrStat_C, /* OUT: error status */ + char *ErrMsg_C /* OUT: [BD_C_ERRMSGLEN] error message */ +); + +/* + * Write a checkpoint file .chkp holding the complete state of the beam (states at the current and + * previous time, input history, inputs set through the interface, and time stepping information). The file is + * only valid for the same build and the same input file. + */ +void BD_C_PackStates( + const char *CheckpointRoot_C, /* IN: root name of the checkpoint file (".chkp" is appended) */ + int *ErrStat_C, /* OUT: error status */ + char *ErrMsg_C /* OUT: [BD_C_ERRMSGLEN] error message */ +); + +/* + * Restore the state of the beam from a checkpoint file written by BD_C_PackStates. BD_C_Init must have been + * called with the same input file first. + */ +void BD_C_UnpackStates( + const char *CheckpointRoot_C, /* IN: root name of the checkpoint file (".chkp" is appended) */ + int *ErrStat_C, /* OUT: error status */ + char *ErrMsg_C /* OUT: [BD_C_ERRMSGLEN] error message */ +); + +/* + * Free all memory held by the library and close any open files. Safe to call after a fatal error or before + * BD_C_Init. + */ +void BD_C_End( + int *ErrStat_C, /* OUT: error status */ + char *ErrMsg_C /* OUT: [BD_C_ERRMSGLEN] error message */ +); + +#ifdef __cplusplus +} +#endif + +#endif /* BEAMDYN_C_BINDING_H */ From b10bc360f3181773499e89252bc845eb088147a5 Mon Sep 17 00:00:00 2001 From: "Harish @ NLR" Date: Thu, 1 Oct 2026 12:09:39 -0600 Subject: [PATCH 2/3] BeamDyn: add C driver for the C interface Add BeamDyn_C_Binding_Driver.c and the beamdyn_c_binding_driver target. The driver runs a beam from C the way the Fortran driver does (root rotating about the origin, constant tip load, output file in the same format) and checks the interface: - run: time marching with optional prescribed sinusoidal root motion; the root reaction is compared with the rigid-body inertial estimate - cantilever: static tip deflection and first bending frequency of a uniform cantilever against the Euler-Bernoulli values - checkpoint: run, write a checkpoint, restore it into a fresh instance, and compare with the straight run Co-Authored-By: Claude Fable 5.1 --- modules/beamdyn/CMakeLists.txt | 11 +- .../beamdyn/src/BeamDyn_C_Binding_Driver.c | 585 ++++++++++++++++++ 2 files changed, 595 insertions(+), 1 deletion(-) create mode 100644 modules/beamdyn/src/BeamDyn_C_Binding_Driver.c diff --git a/modules/beamdyn/CMakeLists.txt b/modules/beamdyn/CMakeLists.txt index 974ec87c66..6bd38aa87e 100644 --- a/modules/beamdyn/CMakeLists.txt +++ b/modules/beamdyn/CMakeLists.txt @@ -64,7 +64,16 @@ if(APPLE OR UNIX) target_compile_definitions(beamdyn_c_bind_static PRIVATE IMPLICIT_DLLEXPORT) endif() -install(TARGETS beamdynlib beamdyn_driver_subs beamdyn_driver beamdyn_c_binding beamdyn_c_bind_static +# C driver for the C-bindings interface (example of driving a standalone beam from C, and checks of the interface) +add_executable(beamdyn_c_binding_driver + src/BeamDyn_C_Binding_Driver.c +) +target_link_libraries(beamdyn_c_binding_driver beamdyn_c_binding) +if(UNIX AND NOT APPLE) + target_link_libraries(beamdyn_c_binding_driver m) +endif() + +install(TARGETS beamdynlib beamdyn_driver_subs beamdyn_driver beamdyn_c_binding beamdyn_c_bind_static beamdyn_c_binding_driver EXPORT "${CMAKE_PROJECT_NAME}Libraries" RUNTIME DESTINATION bin LIBRARY DESTINATION lib diff --git a/modules/beamdyn/src/BeamDyn_C_Binding_Driver.c b/modules/beamdyn/src/BeamDyn_C_Binding_Driver.c new file mode 100644 index 0000000000..98b84074b4 --- /dev/null +++ b/modules/beamdyn/src/BeamDyn_C_Binding_Driver.c @@ -0,0 +1,585 @@ +/********************************************************************************************************************************** + * LICENSING + * Copyright (C) 2026 National Renewable Energy Laboratory + * + * This file is part of BeamDyn. + * + * 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. + * + **********************************************************************************************************************************/ +/* + * Driver for the BeamDyn C interface. It drives a standalone BeamDyn beam from C in the same way the BeamDyn + * Fortran driver does (root rotating about the origin at a constant angular velocity, constant tip load, output file + * with the requested channels) and provides checks of the interface against analytic results and against itself. + * + * Usage: beamdyn_c_binding_driver [key=value ...] + * + * Modes: + * run Run the beam from t=0 to tmax and write .out (same format as the BeamDyn driver output file). + * cantilever Static tip deflection and first bending frequency of a uniform cantilever, compared with the + * Euler-Bernoulli values for the EI and mass per unit length given (EI=, mu=, tip=Fx). + * checkpoint Run to tmax, then repeat the run with a checkpoint written at tchk and restored into a fresh + * instance; compare the tip motion and root reaction of the two runs. + * + * Options (key=value; vectors as comma separated values): + * root=NAME root name for output files (default: run, cantilever, or checkpoint) + * dt= time step (s) [0.002] + * tmax= end time (s) [1.0] + * interp= input interpolation order, 1 or 2 [1] + * dynamic= 1 dynamic solve, 0 static solve [1] + * wrout= 1 write .out, 0 do not [1 for run, 0 otherwise] + * grav=gx,gy,gz gravity vector (m/s^2) [0,0,0] + * pos=x,y,z root position (m) [0,0,0] + * ori=r11,...,r33 root orientation DCM, global to root, row by row [identity] + * omega=wx,wy,wz root angular velocity about the origin (rad/s) [0,0,0] + * tip=Fx,Fy,Fz point force at the tip node (N) [0,0,0] + * trelease= time at which the tip force is removed (s) [never] + * amp= freq= sinusoidal root displacement along x: amp*sin(2*pi*freq*t) [0, 0] + * EI= mu= bending stiffness (N-m^2) and mass per unit length (kg/m) for the analytic comparisons + * tchk= checkpoint time for the checkpoint mode (s) [tmax/2] + * tol= relative tolerance for pass/fail of the checks [0.01] + */ +#include +#include +#include +#include + +#include "BeamDyn_C_Binding.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +/*----------------------------------------------------------------------------------------------------------------*/ +/* Options */ +/*----------------------------------------------------------------------------------------------------------------*/ +typedef struct { + const char *input_file; + char root[BD_C_INTFSTRLEN]; + double dt; + double tmax; + int interp; + int dynamic; + int wrout; + double grav[3]; + double pos[3]; + double ori[9]; + double omega[3]; + double tip[3]; + double t_release; + double amp; + double freq; + double EI; + double mu; + double t_chk; + double tol; +} Options; + +static void set_defaults(Options *o, const char *mode) +{ + memset(o, 0, sizeof(*o)); + strncpy(o->root, mode, sizeof(o->root) - 1); + o->dt = 0.002; + o->tmax = 1.0; + o->interp = 1; + o->dynamic = 1; + o->wrout = (strcmp(mode, "run") == 0) ? 1 : 0; + o->ori[0] = 1.0; o->ori[4] = 1.0; o->ori[8] = 1.0; + o->t_release = -1.0; + o->t_chk = -1.0; + o->tol = 0.01; +} + +static int parse_vector(const char *s, double *v, int n) +{ + int i; + char *end; + for (i = 0; i < n; i++) { + v[i] = strtod(s, &end); + if (end == s) return 1; + s = end; + if (*s == ',') s++; + } + return 0; +} + +static int parse_options(int argc, char **argv, Options *o) +{ + int i, k; + double v[9]; + for (i = 3; i < argc; i++) { + const char *eq = strchr(argv[i], '='); + if (eq == NULL) { fprintf(stderr, "Option '%s' is not of the form key=value\n", argv[i]); return 1; } + size_t klen = (size_t)(eq - argv[i]); + const char *val = eq + 1; + int bad = 0; +#define KEY(name) (klen == strlen(name) && strncmp(argv[i], name, klen) == 0) + if (KEY("root")) { strncpy(o->root, val, sizeof(o->root) - 1); } + else if (KEY("dt")) { o->dt = atof(val); } + else if (KEY("tmax")) { o->tmax = atof(val); } + else if (KEY("interp")) { o->interp = atoi(val); } + else if (KEY("dynamic")) { o->dynamic = atoi(val); } + else if (KEY("wrout")) { o->wrout = atoi(val); } + else if (KEY("grav")) { bad = parse_vector(val, v, 3); for (k = 0; k < 3; k++) o->grav[k] = v[k]; } + else if (KEY("pos")) { bad = parse_vector(val, v, 3); for (k = 0; k < 3; k++) o->pos[k] = v[k]; } + else if (KEY("ori")) { bad = parse_vector(val, o->ori, 9); } + else if (KEY("omega")) { bad = parse_vector(val, o->omega, 3); } + else if (KEY("tip")) { bad = parse_vector(val, o->tip, 3); } + else if (KEY("trelease")) { o->t_release = atof(val); } + else if (KEY("amp")) { o->amp = atof(val); } + else if (KEY("freq")) { o->freq = atof(val); } + else if (KEY("EI")) { o->EI = atof(val); } + else if (KEY("mu")) { o->mu = atof(val); } + else if (KEY("tchk")) { o->t_chk = atof(val); } + else if (KEY("tol")) { o->tol = atof(val); } + else { fprintf(stderr, "Unknown option '%s'\n", argv[i]); return 1; } +#undef KEY + if (bad) { fprintf(stderr, "Could not parse the value of option '%s'\n", argv[i]); return 1; } + } + return 0; +} + +/*----------------------------------------------------------------------------------------------------------------*/ +/* Error handling */ +/*----------------------------------------------------------------------------------------------------------------*/ +static int err_stat = 0; +static char err_msg[BD_C_ERRMSGLEN]; + +/* Print any message from the library; on a fatal error end the library and exit. */ +static void check(const char *where) +{ + if (err_stat != BD_C_ERRID_NONE) + printf("%s: ErrStat=%d\n%s\n", where, err_stat, err_msg); + if (err_stat >= BD_C_ERRID_FATAL) { + int es; + char em[BD_C_ERRMSGLEN]; + BD_C_End(&es, em); + exit(EXIT_FAILURE); + } +} + +/*----------------------------------------------------------------------------------------------------------------*/ +/* Root kinematics: rigid rotation about the origin at constant angular velocity (as in the BeamDyn driver), plus */ +/* an optional sinusoidal displacement along x. */ +/*----------------------------------------------------------------------------------------------------------------*/ +static void cross(const double a[3], const double b[3], double c[3]) +{ + c[0] = a[1]*b[2] - a[2]*b[1]; + c[1] = a[2]*b[0] - a[0]*b[2]; + c[2] = a[0]*b[1] - a[1]*b[0]; +} + +/* Rotation matrix (active rotation of a vector) for a rotation of angle about the unit vector k: R = I + sin(a) K + (1-cos(a)) K^2 */ +static void rodrigues(const double k[3], double angle, double R[9]) +{ + double s = sin(angle), c = cos(angle), v = 1.0 - c; + R[0] = c + k[0]*k[0]*v; R[1] = k[0]*k[1]*v - k[2]*s; R[2] = k[0]*k[2]*v + k[1]*s; + R[3] = k[1]*k[0]*v + k[2]*s; R[4] = c + k[1]*k[1]*v; R[5] = k[1]*k[2]*v - k[0]*s; + R[6] = k[2]*k[0]*v - k[1]*s; R[7] = k[2]*k[1]*v + k[0]*s; R[8] = c + k[2]*k[2]*v; +} + +static void root_motion(const Options *o, double t, double disp[3], double ori[9], double vel[6], double acc[6]) +{ + int i, j, k; + double w = sqrt(o->omega[0]*o->omega[0] + o->omega[1]*o->omega[1] + o->omega[2]*o->omega[2]); + double R[9] = {1,0,0, 0,1,0, 0,0,1}; + double r0[3] = {o->pos[0], o->pos[1], o->pos[2]}; + double p[3], v[3], wxp[3], a[3]; + + if (w > 0.0) { + double khat[3] = {o->omega[0]/w, o->omega[1]/w, o->omega[2]/w}; + rodrigues(khat, w*t, R); + } + for (i = 0; i < 3; i++) p[i] = R[3*i]*r0[0] + R[3*i+1]*r0[1] + R[3*i+2]*r0[2]; /* rotated root position */ + cross(o->omega, p, v); /* velocity: omega x p */ + cross(o->omega, p, wxp); + cross(o->omega, wxp, a); /* acceleration: omega x (omega x p) */ + + for (i = 0; i < 3; i++) { + disp[i] = (p[i] - r0[i]); + vel[i] = v[i]; + vel[3+i] = o->omega[i]; + acc[i] = a[i]; + acc[3+i] = 0.0; + } + /* The root frame rotates with the body: DCM(t) = DCM(0) * R^T */ + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) { + ori[3*i+j] = 0.0; + for (k = 0; k < 3; k++) ori[3*i+j] += o->ori[3*i+k] * R[3*j+k]; + } + + if (o->amp != 0.0) { + double wf = 2.0*M_PI*o->freq; + disp[0] += ( o->amp*sin(wf*t)); + vel[0] += ( o->amp*wf*cos(wf*t)); + acc[0] += (-o->amp*wf*wf*sin(wf*t)); + } +} + +/*----------------------------------------------------------------------------------------------------------------*/ +/* Time series recorded from the beam */ +/*----------------------------------------------------------------------------------------------------------------*/ +typedef struct { + int n, nmax; + double *t; + double *tip; /* [3*n] tip displacement from the reference tip position (m) */ + double *react; /* [6*n] root reaction force and moment (N, N-m) */ + double length; /* distance from the root to the tip output node in the reference configuration (m) */ +} Series; + +static void series_alloc(Series *s, int nmax) +{ + s->n = 0; s->nmax = nmax; + s->t = (double*)calloc((size_t)nmax, sizeof(double)); + s->tip = (double*)calloc((size_t)3*nmax, sizeof(double)); + s->react = (double*)calloc((size_t)6*nmax, sizeof(double)); + s->length = 0.0; +} + +static void series_free(Series *s) +{ + free(s->t); free(s->tip); free(s->react); + s->t = NULL; s->tip = NULL; s->react = NULL; s->n = 0; +} + +/*----------------------------------------------------------------------------------------------------------------*/ +/* Run the beam from t_start to t_end with the loads and root motion described by the options. */ +/* t_pack >= 0: write a checkpoint when the simulation time reaches t_pack (and stop if stop_after_pack) */ +/* unpack != 0: restore the checkpoint right after initialization (the run then continues from t_start) */ +/*----------------------------------------------------------------------------------------------------------------*/ +static void run_simulation(const Options *o, const char *root, double t_start, double t_end, + double t_pack, int stop_after_pack, const char *chkp_root, int unpack, Series *s) +{ + int passed = 0, len = (int)strlen(o->input_file); + const char *input_file = o->input_file; + int nOut = 0, nPL = 0, nDL = 0, nCh = 0; + int i, k, n, n0, iTip = 0, jTip = 0; + char rootbuf[BD_C_INTFSTRLEN], chkpbuf[BD_C_INTFSTRLEN]; + char *names, *units; + float *refpos, *plpos, *dlpos, *nodepos, *nodevel, *nodeacc, *ploads, *chan; + double *refori, *nodeori; + float react[6]; double disp[3], vel[6], acc[6], rootvel0[6]; + double ori[9], omega_x_r0[3], r0[3], d, dmax; + double t, tn; + + memset(rootbuf, 0, sizeof(rootbuf)); strncpy(rootbuf, root, sizeof(rootbuf) - 1); + memset(chkpbuf, 0, sizeof(chkpbuf)); if (chkp_root != NULL) strncpy(chkpbuf, chkp_root, sizeof(chkpbuf) - 1); + + /* Initial root velocity for a root rotating about the origin: omega x r0 */ + r0[0] = o->pos[0]; r0[1] = o->pos[1]; r0[2] = o->pos[2]; + cross(o->omega, r0, omega_x_r0); + for (k = 0; k < 3; k++) { rootvel0[k] = omega_x_r0[k]; rootvel0[3+k] = o->omega[k]; } + if (o->amp != 0.0) rootvel0[0] += (o->amp*2.0*M_PI*o->freq); + + names = (char*)malloc(BD_C_CHANNELBUFLEN); + units = (char*)malloc(BD_C_CHANNELBUFLEN); + + BD_C_Init(&passed, &input_file, &len, rootbuf, o->pos, o->ori, rootvel0, o->grav, &o->dt, &o->interp, + &o->dynamic, &o->wrout, &nOut, &nPL, &nDL, &nCh, names, units, &err_stat, err_msg); + check("BD_C_Init"); + + refpos = (float*) malloc((size_t)3*nOut*sizeof(float)); + refori = (double*)malloc((size_t)9*nOut*sizeof(double)); + plpos = (float*) malloc((size_t)3*nPL *sizeof(float)); + dlpos = (float*) malloc((size_t)3*nDL *sizeof(float)); + nodepos = (float*) malloc((size_t)3*nOut*sizeof(float)); + nodeori = (double*)malloc((size_t)9*nOut*sizeof(double)); + nodevel = (float*) malloc((size_t)6*nOut*sizeof(float)); + nodeacc = (float*) malloc((size_t)6*nOut*sizeof(float)); + ploads = (float*) calloc((size_t)6*nPL, sizeof(float)); + chan = (float*) malloc((size_t)(nCh > 0 ? nCh : 1)*sizeof(float)); + + BD_C_GetRefPositions(refpos, refori, plpos, dlpos, &err_stat, err_msg); + check("BD_C_GetRefPositions"); + + /* Tip: the output node and point load node farthest from the first output node */ + dmax = -1.0; + for (i = 0; i < nOut; i++) { + d = 0.0; for (k = 0; k < 3; k++) d += pow(refpos[3*i+k] - refpos[k], 2); + if (d > dmax) { dmax = d; iTip = i; } + } + s->length = sqrt(dmax); + dmax = -1.0; + for (i = 0; i < nPL; i++) { + d = 0.0; for (k = 0; k < 3; k++) d += pow(plpos[3*i+k] - refpos[k], 2); + if (d > dmax) { dmax = d; jTip = i; } + } + + printf(" Beam: %d output nodes, %d point load nodes, %d distributed load nodes, %d output channels, length %.6f m\n", + nOut, nPL, nDL, nCh, s->length); + + if (unpack) { + BD_C_UnpackStates(chkpbuf, &err_stat, err_msg); + check("BD_C_UnpackStates"); + } + +#define SET_INPUTS(time) \ + do { \ + int kk; \ + root_motion(o, (time), disp, ori, vel, acc); \ + BD_C_SetRootMotion(disp, ori, vel, acc, &err_stat, err_msg); check("BD_C_SetRootMotion"); \ + for (kk = 0; kk < 3; kk++) \ + ploads[6*jTip+kk] = (o->t_release < 0.0 || (time) < o->t_release - 0.5*o->dt) ? (float)o->tip[kk] : 0.0f; \ + BD_C_SetPointLoads(ploads, &err_stat, err_msg); check("BD_C_SetPointLoads"); \ + } while (0) + +#define RECORD(time) \ + do { \ + int kk; \ + if (s->n < s->nmax) { \ + s->t[s->n] = (time); \ + for (kk = 0; kk < 3; kk++) s->tip[3*s->n+kk] = (double)nodepos[3*iTip+kk] - (double)refpos[3*iTip+kk]; \ + for (kk = 0; kk < 6; kk++) s->react[6*s->n+kk] = (double)react[kk]; \ + s->n++; \ + } \ + } while (0) + + n0 = (int)lround(t_start/o->dt); + t = n0*o->dt; + + SET_INPUTS(t); + BD_C_CalcOutput(&t, nodepos, nodeori, nodevel, nodeacc, react, chan, &err_stat, err_msg); + check("BD_C_CalcOutput"); + RECORD(t); + + for (n = n0; ; n++) { + t = n*o->dt; + tn = (n+1)*o->dt; + if (t_pack >= 0.0 && fabs(t - t_pack) < 0.5*o->dt) { + BD_C_PackStates(chkpbuf, &err_stat, err_msg); + check("BD_C_PackStates"); + printf(" Checkpoint written at t = %.6f s\n", t); + if (stop_after_pack) break; + } + if (tn > t_end + 0.5*o->dt) break; + + SET_INPUTS(tn); + BD_C_UpdateStates(&t, &tn, &err_stat, err_msg); + check("BD_C_UpdateStates"); + BD_C_CalcOutput(&tn, nodepos, nodeori, nodevel, nodeacc, react, chan, &err_stat, err_msg); + check("BD_C_CalcOutput"); + RECORD(tn); + } +#undef SET_INPUTS +#undef RECORD + + BD_C_End(&err_stat, err_msg); + check("BD_C_End"); + + free(names); free(units); free(refpos); free(refori); free(plpos); free(dlpos); + free(nodepos); free(nodeori); free(nodevel); free(nodeacc); free(ploads); free(chan); +} + +/*----------------------------------------------------------------------------------------------------------------*/ +/* Analysis helpers */ +/*----------------------------------------------------------------------------------------------------------------*/ +/* Frequency and damping ratio of a decaying oscillation x(t) about zero for t >= t0, from the upward zero crossings + * and the logarithmic decrement of the positive peaks. Returns the number of cycles used (0 on failure). */ +static int measure_oscillation(const Series *s, int comp, double t0, double *f_damped, double *zeta) +{ + int i, ncross = 0, npeak = 0; + double t_first = 0.0, t_last = 0.0, a_first = 0.0, a_last = 0.0; + for (i = 1; i < s->n - 1; i++) { + double x0 = s->tip[3*(i-1)+comp], x1 = s->tip[3*i+comp], x2 = s->tip[3*(i+1)+comp]; + if (s->t[i] < t0) continue; + if (x0 < 0.0 && x1 >= 0.0) { /* upward zero crossing, linear interpolation */ + double tc = s->t[i-1] + (0.0 - x0)/(x1 - x0)*(s->t[i] - s->t[i-1]); + if (ncross == 0) t_first = tc; + t_last = tc; + ncross++; + } + if (x1 > 0.0 && x1 >= x0 && x1 > x2) { /* positive peak */ + if (npeak == 0) a_first = x1; + a_last = x1; + npeak++; + } + } + if (ncross < 3 || npeak < 2) return 0; + *f_damped = (ncross - 1)/(t_last - t_first); + { + double delta = log(a_first/a_last)/(npeak - 1); /* logarithmic decrement */ + *zeta = delta/sqrt(4.0*M_PI*M_PI + delta*delta); + } + return ncross - 1; +} + +static int report(const char *what, double value, double reference, double tol) +{ + double err = (reference != 0.0) ? fabs(value - reference)/fabs(reference) : fabs(value - reference); + int ok = err <= tol; + printf(" %-44s %16.8e reference %16.8e rel. diff %10.3e %s\n", what, value, reference, err, ok ? "PASS" : "FAIL"); + return ok ? 0 : 1; +} + +/*----------------------------------------------------------------------------------------------------------------*/ +/* Modes */ +/*----------------------------------------------------------------------------------------------------------------*/ +static int mode_run(const Options *o) +{ + Series s; + int i, n = (int)lround(o->tmax/o->dt) + 2, fails = 0; + series_alloc(&s, n); + run_simulation(o, o->root, 0.0, o->tmax, -1.0, 0, NULL, 0, &s); + + printf(" Final time %.6f s: tip displacement [%.8e %.8e %.8e] m, root reaction force [%.8e %.8e %.8e] N\n", + s.t[s.n-1], s.tip[3*(s.n-1)], s.tip[3*(s.n-1)+1], s.tip[3*(s.n-1)+2], + s.react[6*(s.n-1)], s.react[6*(s.n-1)+1], s.react[6*(s.n-1)+2]); + + if (o->amp != 0.0 && o->freq > 0.0) { + /* Prescribed sinusoidal root motion: the root reaction must follow the root acceleration. Over the last + * period, compare the amplitude of the x reaction force with the rigid-body inertial estimate + * mu*L*amp*(2*pi*f)^2 (times the single-mode dynamic amplification if EI is given). */ + double T = 1.0/o->freq, fmax = 0.0, amax = 0.0; + for (i = 0; i < s.n; i++) { + if (s.t[i] < o->tmax - T) continue; + if (fabs(s.react[6*i]) > fmax) fmax = fabs(s.react[6*i]); + if (fabs(s.tip[3*i]) > amax) amax = fabs(s.tip[3*i]); + } + printf(" Root motion x = %.4e*sin(2*pi*%.4f*t): max |root Fx| over the last period %.8e N, max |tip x displacement| %.8e m\n", + o->amp, o->freq, fmax, amax); + if (fmax <= 0.0) { printf(" FAIL: the root reaction does not respond to the root motion\n"); fails++; } + if (o->mu > 0.0) { + double M = o->mu*s.length, a = o->amp*pow(2.0*M_PI*o->freq, 2), est = M*a; + if (o->EI > 0.0) { + double f1 = 1.8751040687*1.8751040687/(2.0*M_PI)*sqrt(o->EI/(o->mu*pow(s.length, 4))); + est /= fabs(1.0 - pow(o->freq/f1, 2)); + } + fails += report("max |root Fx| vs. inertial estimate", fmax, est, o->tol); + } + } + series_free(&s); + return fails; +} + +static int mode_cantilever(const Options *o) +{ + Options os = *o, od = *o; + Series s; + int n, i, fails = 0, ncyc; + double L, F, d_static, d_settled = 0.0, d_eb, f_eb, f_d = 0.0, zeta = 0.0, f_n; + char root[BD_C_INTFSTRLEN]; + + if (o->EI <= 0.0 || o->mu <= 0.0 || o->tip[0] == 0.0) { + fprintf(stderr, "The cantilever mode needs EI=, mu=, and tip=Fx,0,0\n"); + return 1; + } + F = o->tip[0]; + + /* (a) static solve: tip deflection under the tip force */ + printf("\n-- Static solve: tip force %.4e N\n", F); + os.dynamic = 0; os.t_release = -1.0; os.amp = 0.0; + snprintf(root, sizeof(root), "%s_static", o->root); + series_alloc(&s, 20); + run_simulation(&os, root, 0.0, 5.0*o->dt, -1.0, 0, NULL, 0, &s); + L = s.length; + d_static = s.tip[3*(s.n-1)]; + series_free(&s); + d_eb = F*L*L*L/(3.0*o->EI); + fails += report("static tip deflection (static solve)", d_static, d_eb, o->tol); + + /* (b) dynamic solve: settle under the tip force, release, measure the first bending frequency from the tip motion */ + od.dynamic = 1; od.amp = 0.0; + if (od.t_release < 0.0) od.t_release = 0.5*o->tmax; + printf("\n-- Dynamic solve: tip force applied until t = %.3f s, free vibration until t = %.3f s\n", od.t_release, od.tmax); + snprintf(root, sizeof(root), "%s_dynamic", o->root); + n = (int)lround(o->tmax/o->dt) + 2; + series_alloc(&s, n); + run_simulation(&od, root, 0.0, od.tmax, -1.0, 0, NULL, 0, &s); + for (i = 0; i < s.n; i++) if (s.t[i] <= od.t_release - 0.5*o->dt) d_settled = s.tip[3*i]; + fails += report("static tip deflection (settled dynamic solve)", d_settled, d_eb, o->tol); + + f_eb = 1.8751040687*1.8751040687/(2.0*M_PI)*sqrt(o->EI/(o->mu*pow(L, 4))); + ncyc = measure_oscillation(&s, 0, od.t_release + 1.0/f_eb, &f_d, &zeta); + if (ncyc == 0) { + printf(" FAIL: could not measure the free vibration (not enough cycles after the release)\n"); + fails++; + } else { + f_n = f_d/sqrt(1.0 - zeta*zeta); + printf(" Free vibration: %d cycles, damped frequency %.6f Hz, damping ratio %.5f\n", ncyc, f_d, zeta); + fails += report("first bending frequency (undamped)", f_n, f_eb, o->tol); + } + series_free(&s); + return fails; +} + +static int mode_checkpoint(const Options *o) +{ + Options oc = *o; + Series s1, s2; + int n = (int)lround(o->tmax/o->dt) + 2, i, j, k, n0, fails = 0; + double dtip = 0.0, dreact = 0.0, tipmax = 0.0, reactmax = 0.0; + char root[BD_C_INTFSTRLEN], chkp[BD_C_INTFSTRLEN]; + + if (oc.t_chk < 0.0) oc.t_chk = 0.5*o->tmax; + oc.t_chk = lround(oc.t_chk/o->dt)*o->dt; + + printf("\n-- Straight run to t = %.3f s\n", o->tmax); + series_alloc(&s1, n); + snprintf(root, sizeof(root), "%s_straight", o->root); + run_simulation(&oc, root, 0.0, o->tmax, -1.0, 0, NULL, 0, &s1); + + printf("\n-- Run to t = %.3f s and write a checkpoint\n", oc.t_chk); + snprintf(root, sizeof(root), "%s_part1", o->root); + snprintf(chkp, sizeof(chkp), "%s_checkpoint", o->root); + series_alloc(&s2, n); + run_simulation(&oc, root, 0.0, oc.t_chk, oc.t_chk, 1, chkp, 0, &s2); + series_free(&s2); + + printf("\n-- Restore the checkpoint into a fresh instance and continue to t = %.3f s\n", o->tmax); + series_alloc(&s2, n); + snprintf(root, sizeof(root), "%s_part2", o->root); + run_simulation(&oc, root, oc.t_chk, o->tmax, -1.0, 0, chkp, 1, &s2); + + /* Compare the restarted run with the straight run at the same times */ + n0 = (int)lround(oc.t_chk/o->dt); + for (j = 0; j < s2.n; j++) { + i = n0 + j; + if (i >= s1.n) break; + if (fabs(s1.t[i] - s2.t[j]) > 0.5*o->dt) { printf(" FAIL: time mismatch %f vs %f\n", s1.t[i], s2.t[j]); fails++; break; } + for (k = 0; k < 3; k++) { dtip = fmax(dtip, fabs(s1.tip[3*i+k] - s2.tip[3*j+k])); tipmax = fmax(tipmax, fabs(s1.tip[3*i+k])); } + for (k = 0; k < 6; k++) { dreact = fmax(dreact, fabs(s1.react[6*i+k] - s2.react[6*j+k])); reactmax = fmax(reactmax, fabs(s1.react[6*i+k])); } + } + printf(" Compared %d output times after the checkpoint\n", j); + printf(" max |tip displacement difference| %.3e m (max |tip displacement| %.3e m)\n", dtip, tipmax); + printf(" max |root reaction difference| %.3e (max |root reaction| %.3e)\n", dreact, reactmax); + fails += report("restart vs. straight run, tip displacement diff", dtip, 0.0, 1e-6*(tipmax > 0.0 ? tipmax : 1.0)); + fails += report("restart vs. straight run, root reaction diff", dreact, 0.0, 1e-6*(reactmax > 0.0 ? reactmax : 1.0)); + + series_free(&s1); + series_free(&s2); + return fails; +} + +/*----------------------------------------------------------------------------------------------------------------*/ +int main(int argc, char **argv) +{ + Options o; + int fails; + + if (argc < 3) { + fprintf(stderr, "Usage: %s [key=value ...]\n", argv[0]); + return EXIT_FAILURE; + } + set_defaults(&o, argv[1]); + o.input_file = argv[2]; + if (parse_options(argc, argv, &o)) return EXIT_FAILURE; + + if (strcmp(argv[1], "run") == 0) fails = mode_run(&o); + else if (strcmp(argv[1], "cantilever") == 0) fails = mode_cantilever(&o); + else if (strcmp(argv[1], "checkpoint") == 0) fails = mode_checkpoint(&o); + else { fprintf(stderr, "Unknown mode '%s'\n", argv[1]); return EXIT_FAILURE; } + + printf("\n%s: %s\n", argv[1], fails == 0 ? "all checks passed" : "some checks FAILED"); + return fails == 0 ? EXIT_SUCCESS : EXIT_FAILURE; +} From 4be9e4b743c429df0fa46d34f69ca45fff1561e5 Mon Sep 17 00:00:00 2001 From: "Harish @ NLR" Date: Thu, 1 Oct 2026 12:44:11 -0600 Subject: [PATCH 3/3] BeamDyn: return the time step BeamDyn uses from BD_Init BD_Init took the suggested coupling interval as the default for DTBeam but never reported the time step it actually uses back through the Interval argument, as the other modules do. Set Interval to p%dt after the parameters are read so that a calling code can check or adopt the BeamDyn time step. The C interface uses the returned value for its time step check. Results are unchanged when DTBeam is DEFAULT (the regression and example cases). Co-Authored-By: Claude Fable 5.1 --- modules/beamdyn/src/BeamDyn.f90 | 3 +++ modules/beamdyn/src/BeamDyn_C_Binding.f90 | 7 ++++--- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/modules/beamdyn/src/BeamDyn.f90 b/modules/beamdyn/src/BeamDyn.f90 index 65791d3b43..c61f282387 100644 --- a/modules/beamdyn/src/BeamDyn.f90 +++ b/modules/beamdyn/src/BeamDyn.f90 @@ -126,6 +126,9 @@ SUBROUTINE BD_Init( InitInp, u, p, x, xd, z, OtherState, y, MiscVar, Interval, I ! this routine sets *some* of the parameters (basically the "easy" ones) call SetParameters(InitInp, InputFileData, p, OtherState, ErrStat2, ErrMsg2); if (Failed()) return + ! return the time step BeamDyn uses (DTBeam from the input file, or the suggested Interval when DTBeam is DEFAULT) + Interval = p%dt + ! Temporary GLL point intrinsic coordinates array CALL BD_GenerateGLL(p%nodes_per_elem,GLL_nodes,ErrStat2,ErrMsg2); if (Failed()) return diff --git a/modules/beamdyn/src/BeamDyn_C_Binding.f90 b/modules/beamdyn/src/BeamDyn_C_Binding.f90 index b3cfb0c2ee..d09a8d782c 100644 --- a/modules/beamdyn/src/BeamDyn_C_Binding.f90 +++ b/modules/beamdyn/src/BeamDyn_C_Binding.f90 @@ -295,10 +295,11 @@ SUBROUTINE BD_C_Init( & ! Remove the temporary input file now that it has been read IF (LEN_TRIM(TmpFileName) > 0) CALL DeleteFile( TmpFileName ) - ! The states are advanced by DT_C here, so BeamDyn must use the same timestep (DTBeam in the input file must be DEFAULT or equal to DT_C) - IF ( .NOT. EqualRealNos( p%dt, dT_Global ) ) THEN + ! The states are advanced by DT_C here, so BeamDyn must use the same timestep (DTBeam in the input file must be DEFAULT or equal to DT_C). + ! BD_Init returns the timestep it uses in dT_Interval. + IF ( .NOT. EqualRealNos( dT_Interval, dT_Global ) ) THEN ErrStat_F2 = ErrID_Fatal - ErrMsg_F2 = 'The BeamDyn timestep DTBeam ('//TRIM(Num2LStr(p%dt))//' s) must be DEFAULT or equal to the timestep passed to BD_C_Init ('//TRIM(Num2LStr(dT_Global))//' s).' + ErrMsg_F2 = 'The BeamDyn timestep DTBeam ('//TRIM(Num2LStr(dT_Interval))//' s) must be DEFAULT or equal to the timestep passed to BD_C_Init ('//TRIM(Num2LStr(dT_Global))//' s).' IF (Failed()) RETURN ENDIF