FAST Library: add access to hub position and velocity

glue-codes/openfast/src/FAST_Prog.cpp

@@ -15,5 +15,19 @@ int main(int argc, char** argv) {
 
    FastLibAPI fastlib = FastLibAPI(input_file_name);
    fastlib.fast_run();
+
+   // // Get the hub position 
+   // float absolute_position[3] = {};
+   // float rotational_velocity[3] = {};
+   // double orientation_dcm[9] = {};
+
+   // fastlib.get_hub_position(absolute_position, rotational_velocity, orientation_dcm);
+
+   // printf("%f %f %f\n", absolute_position[0], absolute_position[1], absolute_position[2]);
+   // printf("%f %f %f\n", rotational_velocity[0], rotational_velocity[1], rotational_velocity[2]);
+   // printf("%f %f %f\n", orientation_dcm[0], orientation_dcm[1], orientation_dcm[2]);
+   // printf("%f %f %f\n", orientation_dcm[3], orientation_dcm[4], orientation_dcm[5]);
+   // printf("%f %f %f\n", orientation_dcm[6], orientation_dcm[7], orientation_dcm[8]);
+
    return 0;
 }

glue-codes/openfast/src/FastLibAPI.cpp

@@ -163,3 +163,17 @@ int FastLibAPI::total_time_steps() {
     // We assume here t_initial is always 0
     return ceil( t_max / dt ) + 1;
 }
+
+void FastLibAPI::get_hub_position(float *absolute_position, float *rotational_velocity, double *orientation_dcm) {
+    int _error_status = 0;
+    char _error_message[INTERFACE_STRING_LENGTH];
+
+    FAST_HubPosition(
+        &i_turb,
+        absolute_position,
+        rotational_velocity,
+        orientation_dcm,
+        &_error_status,
+        _error_message
+    );
+}

glue-codes/openfast/src/FastLibAPI.h

@@ -49,6 +49,7 @@ class FastLibAPI {
         void fast_run();
         int total_time_steps();
         std::string output_channel_names();
+        void get_hub_position(float *absolute_position, float *rotational_velocity, double *orientation_dcm);
 };
 
 #endif

glue-codes/python/openfast_library.py

@@ -5,11 +5,12 @@
     byref,
     c_int,
     c_double,
+    c_float,
     c_char,
     c_bool
 )
 import os
-from typing import List
+from typing import List, Tuple
 import numpy as np
 import math
 
@@ -111,6 +112,16 @@ def _initialize_routines(self) -> None:
         ]
         self.FAST_End.restype = c_int
 
+        self.FAST_HubPosition.argtypes = [
+            POINTER(c_int),         # iTurb IN
+            POINTER(c_float),       # AbsPosition_c(3) OUT
+            POINTER(c_float),       # RotationalVel_c(3) OUT
+            POINTER(c_double),      # Orientation_c(9) OUT
+            POINTER(c_int),         # ErrStat_c OUT
+            POINTER(c_char)         # ErrMsg_c OUT
+        ]
+        self.FAST_HubPosition.restype = c_int
+
 
     def fatal_error(self, error_status) -> bool:
         return error_status.value >= self.abort_error_level.value
@@ -240,3 +251,27 @@ def output_channel_names(self) -> List:
         output_channel_names = self.channel_names.value.split()
         output_channel_names = [n.decode('UTF-8') for n in output_channel_names]        
         return output_channel_names
+
+
+    def get_hub_position(self) -> Tuple:
+        _error_status = c_int(0)
+        _error_message = create_string_buffer(IntfStrLen)
+
+        # Data buffers
+        absolute_position = (c_float * 3)(0.0, )
+        rotational_velocity = (c_float * 3)(0.0, )
+        orientation_dcm = (c_double * 9)(0.0, )
+
+        # Get hub position from the fast library
+        self.FAST_HubPosition(
+            byref(self.i_turb),
+            absolute_position,
+            rotational_velocity,
+            orientation_dcm,
+            byref(_error_status),
+            _error_message            
+        )
+        if self.fatal_error(_error_status):
+            raise RuntimeError(f"Error {_error_status.value}: {_error_message.value}")
+
+        return absolute_position, rotational_velocity, orientation_dcm

modules/openfast-library/src/FAST_Library.f90

@@ -322,6 +322,42 @@ subroutine FAST_Update(iTurb, NumInputs_c, NumOutputs_c, InputAry, OutputAry, En
 
 end subroutine FAST_Update 
 !==================================================================================================================================
+! Get the hub's absolute position, rotation velocity, and orientation DCM for the current time step
+subroutine FAST_HubPosition(iTurb, AbsPosition_c, RotationalVel_c, Orientation_c, ErrStat_c, ErrMsg_c) BIND (C, NAME='FAST_HubPosition')
+   IMPLICIT NONE
+#ifndef IMPLICIT_DLLEXPORT
+!DEC$ ATTRIBUTES DLLEXPORT :: FAST_HubPosition
+!GCC$ ATTRIBUTES DLLEXPORT :: FAST_HubPosition
+#endif
+   INTEGER(C_INT),         INTENT(IN   ) :: iTurb
+   REAL(C_FLOAT),          INTENT(  OUT) :: AbsPosition_c(3), RotationalVel_c(3)
+   REAL(C_DOUBLE),         INTENT(  OUT) :: Orientation_c(9)
+   INTEGER(C_INT),         INTENT(  OUT) :: ErrStat_c
+   CHARACTER(KIND=C_CHAR), INTENT(  OUT) :: ErrMsg_c(IntfStrLen)
+
+   ErrStat_c = ErrID_None
+   ErrMsg = C_NULL_CHAR
+
+   if (iTurb > size(Turbine) ) then
+      ErrStat_c = ErrID_Fatal
+      ErrMsg = "iTurb is greater than the number of turbines in the simulation."//C_NULL_CHAR
+      ErrMsg_c = TRANSFER( ErrMsg//C_NULL_CHAR, ErrMsg_c )
+      return
+   end if
+
+   if (.NOT. Turbine(iTurb)%ED%y%HubPtMotion%Committed) then
+      ErrStat_c = ErrID_Fatal
+      ErrMsg = "HubPtMotion mesh has not been committed."//C_NULL_CHAR
+      ErrMsg_c = TRANSFER( ErrMsg//C_NULL_CHAR, ErrMsg_c )
+      return
+   end if
+
+   AbsPosition_c = REAL(Turbine(iTurb)%ED%y%HubPtMotion%Position(:,1), C_FLOAT) + REAL(Turbine(iTurb)%ED%y%HubPtMotion%TranslationDisp(:,1), C_FLOAT)
+   Orientation_c = reshape( Turbine(iTurb)%ED%y%HubPtMotion%Orientation(1:3,1:3,1), (/9/) )
+   RotationalVel_c = Turbine(iTurb)%ED%y%HubPtMotion%RotationVel(:,1)
+
+end subroutine FAST_HubPosition
+!==================================================================================================================================
 !> NOTE: If this interface is changed, update the table in the ServoDyn_IO.f90::WrSumInfo4Simulink routine
 !!    Ideally we would write this summary info from here, but that isn't currently done.  So as a workaround so the user has some
 !!    vague idea what went wrong with their simulation, we have ServoDyn include the arrangement set here in the SrvD.sum file.

modules/openfast-library/src/FAST_Library.h

@@ -23,6 +23,8 @@ EXTERNAL_ROUTINE void FAST_OpFM_Init(int * iTurb, double *TMax, const char *Inpu
 EXTERNAL_ROUTINE void FAST_OpFM_Solution0(int * iTurb, int *ErrStat, char *ErrMsg);
 EXTERNAL_ROUTINE void FAST_OpFM_Step(int * iTurb, int *ErrStat, char *ErrMsg);
 
+EXTERNAL_ROUTINE void FAST_HubPosition(int * iTurb, float * absolute_position, float * rotation_veocity, double * orientation_dcm, int *ErrStat, char *ErrMsg);
+
 EXTERNAL_ROUTINE void FAST_Restart(int * iTurb, const char *CheckpointRootName, int *AbortErrLev, int * NumOuts, double * dt, int * n_t_global, int *ErrStat, char *ErrMsg);
 EXTERNAL_ROUTINE void FAST_Sizes(int * iTurb, const char *InputFileName, int *AbortErrLev, int * NumOuts, double * dt, double * tmax, int *ErrStat, char *ErrMsg, char *ChannelNames, double *TMax = NULL, double *InitInputAry = NULL);
 EXTERNAL_ROUTINE void FAST_Start(int * iTurb, int *NumInputs_c, int *NumOutputs_c, double *InputAry, double *OutputAry, int *ErrStat, char *ErrMsg);

modules/openfast-library/tests/test_openfast_library.py

@@ -0,0 +1,47 @@
+
+import sys
+import argparse
+import numpy as np
+from pathlib import Path
+interface_path = Path(__file__).parent.parent.parent.parent / "glue-codes" / "python"
+sys.path.insert(0, str(interface_path))
+import openfast_library
+
+def test_hub_position(library_path, input_file):
+    openfastlib = openfast_library.FastLibAPI(library_path, input_file)
+    openfastlib.fast_init()
+    absolute_position, rotational_velocity, orientation_dcm = openfastlib.get_hub_position()
+
+    # Initial hub position is at -5, 0, 90.55
+    np.testing.assert_allclose(
+        absolute_position,
+        np.array([-5.0, 0.0, 90.55]),
+        rtol=1e-5,
+        atol=1e-8,
+        verbose=True
+    )
+
+    # This case is initially still
+    # Velocities should be 0 and the DCM should be identity
+    np.testing.assert_array_equal( rotational_velocity, np.zeros(3) )
+    np.testing.assert_array_equal( np.reshape( orientation_dcm[:], (3,3) ), np.eye(3) )
+
+
+if __name__=="__main__":
+
+    parser = argparse.ArgumentParser(description="Executes Python-based tests for OpenFAST Library.")
+    parser.add_argument("input_file", metavar="Input-File", type=str, nargs=1, help="Path to an input file.")
+
+    args = parser.parse_args()
+    input_file = args.input_file[0]
+
+    library_path = Path(__file__).parent.parent.parent.parent / "build" / "modules" / "openfast-library" / "libopenfastlib"
+    if sys.platform == "linux" or sys.platform == "linux2":
+        library_path = library_path.with_suffix(".so")
+    elif sys.platform == "darwin":
+        library_path = library_path.with_suffix(".dylib")
+    elif sys.platform == "win32":
+        # TODO
+        pass
+
+    test_hub_position(library_path, input_file)

reg_tests/CTestList.cmake

@@ -186,6 +186,14 @@ function(ifw_py_regression TESTNAME LABEL)
   regression(${TEST_SCRIPT} ${INFLOWWIND_EXECUTABLE} ${SOURCE_DIRECTORY} ${BUILD_DIRECTORY} ${TESTNAME} "${LABEL}")
 endfunction(ifw_py_regression)
 
+# # Python-based OpenFAST Library tests
+# function(py_openfast_library_regression TESTNAME LABEL)
+#   set(test_module "${CMAKE_SOURCE_DIR}/modules/openfast-library/tests/test_openfast_library.py")
+#   set(input_file "${CMAKE_SOURCE_DIR}/reg_tests/r-test/glue-codes/openfast/5MW_OC4Jckt_ExtPtfm/5MW_OC4Jckt_ExtPtfm.fst")
+#   add_test(${TESTNAME} ${PYTHON_EXECUTABLE} ${test_module} ${input_file} )
+# endfunction(py_openfast_library_regression)
+
+
 #===============================================================================
 # Regression tests
 #===============================================================================
@@ -226,6 +234,11 @@ if(BUILD_OPENFAST_CPP_API)
   of_cpp_interface_regression("5MW_Land_DLL_WTurb_cpp" "openfast;openfastlib;cpp")
 endif()
 
+# # Python-based OpenFAST Library unit tests
+# if(BUILD_SHARED_LIBS)
+#   py_openfast_library_regression("py_openfastlib" "python;openfastlib")
+# endif()
+
 # OpenFAST C++ Driver test
 # This tests the FAST Library and FAST_Library.h
 of_cpp_regression("AWT_YFree_WSt"                   "openfast;openfastlib;cpp;elastodyn;aerodyn15;servodyn")