Merge pull request #1211 from deslaughter/bug/pythonlib

Use dt_out when storing OpenFAST outputs in Python interface

glue-codes/openfast/src/FastLibAPI.cpp

@@ -1,6 +1,7 @@
 
 #include "FastLibAPI.h"
 #include <string>
+#include <sstream>
 #include <stdlib.h>
 #include <iostream>
 #include <math.h>
@@ -12,6 +13,7 @@ FastLibAPI::FastLibAPI(std::string input_file):
 n_turbines(1),
 i_turb(0),
 dt(0.0),
+dt_out(0.0),
 t_max(0.0),
 abort_error_level(4),
 end_early(false),
@@ -33,6 +35,7 @@ bool FastLibAPI::fatal_error(int error_status) {
 void FastLibAPI::fast_init() {
     int _error_status = 0;
     char _error_message[INTERFACE_STRING_LENGTH];
+    char channel_names[MAXIMUM_OUTPUTS * CHANNEL_LENGTH + 1];
 
     std::cout << input_file_name;
 
@@ -51,6 +54,7 @@ void FastLibAPI::fast_init() {
         &abort_error_level,
         &num_outs,
         &dt,
+        &dt_out,
         &t_max,
         &_error_status,
         _error_message,
@@ -61,16 +65,16 @@ void FastLibAPI::fast_init() {
     }
 
     // Allocate the data for the outputs
-
-    // Create a dynamic array of pointers
-    // Then, create a row for every pointer and initialize all elements to 0.0
-    output_values = new double *[total_time_steps()];
-    for (int i=0; i<total_time_steps(); i++) {
-        output_values[i] = new double[num_outs];
-        memset(output_values[i], 0.0, num_outs * sizeof(double));
+    output_values.resize(total_output_steps(), std::vector<double>(num_outs, 0));
+
+    // Get output channel names
+    std::istringstream ss(channel_names);
+    std::string channel_name;
+    output_channel_names.clear();
+    while (ss >> channel_name) 
+    {
+        output_channel_names.push_back(channel_name);
     }
-
-    output_array.resize(num_outs);
 }
 
 void FastLibAPI::fast_sim() {
@@ -82,28 +86,32 @@ void FastLibAPI::fast_sim() {
         &num_inputs,
         &num_outs,
         inp_array,
-        output_array.data(),
+        output_values[0].data(),
         &_error_status,
         _error_message
     );
-    output_values[0] = output_array.data();
     if (fatal_error(_error_status)) {
         fast_deinit();
         throw std::runtime_error( "Error " + std::to_string(_error_status) + ": " + _error_message );
     }
 
+    int output_frequency = round(dt_out/dt);
+    int i_out = 1;
+
     for (int i=1; i<total_time_steps(); i++) {
         FAST_Update(
             &i_turb,
             &num_inputs,
             &num_outs,
             inp_array,
-            output_array.data(),
+            output_values[i_out].data(),
             &end_early,
             &_error_status,
             _error_message
         );
-        output_values[i] = output_array.data();
+        if (i%output_frequency == 0) {
+            i_out++;
+        }
         if (fatal_error(_error_status)) {
             fast_deinit();
             throw std::runtime_error( "Error " + std::to_string(_error_status) + ": " + _error_message );
@@ -165,6 +173,12 @@ int FastLibAPI::total_time_steps() {
     return ceil( t_max / dt ) + 1;
 }
 
+int FastLibAPI::total_output_steps() {
+    // Calculate number of output steps using same method as total_time_steps.
+    // This assumes that the initial time is zero.
+    return ceil( t_max / dt_out ) + 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];

glue-codes/openfast/src/FastLibAPI.h

@@ -14,11 +14,11 @@ class FastLibAPI {
         int n_turbines;
         int i_turb;
         double dt;
+        double dt_out;
         double t_max;
         int abort_error_level;
         bool end_early;
         int num_outs;
-        char channel_names[MAXIMUM_OUTPUTS * CHANNEL_LENGTH + 1];
         bool ended;
 
         // The inputs are meant to be from Simulink.
@@ -29,10 +29,8 @@ class FastLibAPI {
         double inp_array[NumFixedInputs] = {};
 
         // These arrays hold the outputs from OpenFAST
-        // output_array is a 1D vector for the values from a single step
         // output_values is a 2D array for the values from all steps in the simulation
-        std::vector<double> output_array;
-        double **output_values;
+        std::vector<std::vector<double>> output_values;
 
     public:
 
@@ -48,7 +46,8 @@ class FastLibAPI {
         void fast_deinit();
         void fast_run();
         int total_time_steps();
-        std::string output_channel_names();
+        int total_output_steps();
+        std::vector<std::string> output_channel_names;
         void get_hub_position(float *absolute_position, float *rotational_velocity, double *orientation_dcm);
 };
 

glue-codes/python/openfast_library.py

@@ -32,11 +32,12 @@ def __init__(self, library_path: str, input_file_name: str):
         self.n_turbines = c_int(1)
         self.i_turb = c_int(0)
         self.dt = c_double(0.0)
+        self.dt_out = c_double(0.0)
         self.t_max = c_double(0.0)
         self.abort_error_level = c_int(4)  # Initialize to 4 (ErrID_Fatal) and reset to user-given value in FAST_Sizes
         self.end_early = c_bool(False)
         self.num_outs = c_int(0)
-        self.channel_names = create_string_buffer(20 * 4000)
+        self.output_channel_names = []
         self.ended = False
 
         # The inputs are meant to be from Simulink.
@@ -48,9 +49,7 @@ def __init__(self, library_path: str, input_file_name: str):
         self.inp_array = (c_double * self.num_inputs.value)(0.0, )
 
         # These arrays hold the outputs from OpenFAST
-        # output_array is a 1D array for the values from a single step
-        # output_values is a 2D array for the values from all steps in the simulation
-        self.output_array = None
+        # output_values is a 2D array for the values from all output steps in the simulation
         self.output_values = None
 
 
@@ -68,6 +67,7 @@ def _initialize_routines(self) -> None:
             POINTER(c_int),         # AbortErrLev_c OUT
             POINTER(c_int),         # NumOuts_c OUT
             POINTER(c_double),      # dt_c OUT
+            POINTER(c_double),      # dt_out_c OUT
             POINTER(c_double),      # tmax_c OUT
             POINTER(c_int),         # ErrStat_c OUT
             POINTER(c_char),        # ErrMsg_c OUT
@@ -139,28 +139,38 @@ def fast_init(self) -> None:
         if self.fatal_error(_error_status):
             raise RuntimeError(f"Error {_error_status.value}: {_error_message.value}")
 
+        # Create channel names argument
+        channel_names = create_string_buffer(20 * 4000)
+
         self.FAST_Sizes(
             byref(self.i_turb),
             self.input_file_name,
             byref(self.abort_error_level),
             byref(self.num_outs),
             byref(self.dt),
+            byref(self.dt_out),
             byref(self.t_max),
             byref(_error_status),
             _error_message,
-            self.channel_names,
+            channel_names,
             None,   # Optional arguments must pass C-Null pointer; with ctypes, use None.
             None    # Optional arguments must pass C-Null pointer; with ctypes, use None.
         )
         if self.fatal_error(_error_status):
             raise RuntimeError(f"Error {_error_status.value}: {_error_message.value}")
 
+        # Extract channel name strings from argument
+        if len(channel_names.value.split()) == 0:
+            self.output_channel_names = []
+        else:
+            self.output_channel_names = [n.decode('UTF-8') for n in channel_names.value.split()] 
+
         # Allocate the data for the outputs
-        # NOTE: The ctypes array allocation (output_array) must be after the output_values
-        # allocation, or otherwise seg fault.
-        self.output_values = np.empty( (self.total_time_steps, self.num_outs.value) )
-        self.output_array = (c_double * self.num_outs.value)(0.0, )
+        self.output_values = np.zeros( (self.total_output_steps, self.num_outs.value), dtype=c_double, order='C' )
 
+        # Delete error message and channel name character buffers
+        del _error_message
+        del channel_names
 
     def fast_sim(self) -> None:
         _error_status = c_int(0)
@@ -171,27 +181,31 @@ def fast_sim(self) -> None:
             byref(self.num_inputs),
             byref(self.num_outs),
             byref(self.inp_array),
-            byref(self.output_array),
+            self.output_values[0].ctypes.data_as(POINTER(c_double)),
             byref(_error_status),
             _error_message
         )
-        self.output_values[0] = self.output_array[:]
         if self.fatal_error(_error_status):
             self.fast_deinit()
             raise RuntimeError(f"Error {_error_status.value}: {_error_message.value}")
 
+        # Calculate output frequency and initialize output index
+        output_frequency = round(self.dt_out.value/self.dt.value)
+        i_out = 1
+
         for i in range( 1, self.total_time_steps ):
             self.FAST_Update(
                 byref(self.i_turb),
                 byref(self.num_inputs),
                 byref(self.num_outs),
                 byref(self.inp_array),
-                byref(self.output_array),
+                self.output_values[i_out].ctypes.data_as(POINTER(c_double)),
                 byref(self.end_early),
                 byref(_error_status),
                 _error_message
             )
-            self.output_values[i] = self.output_array[:]
+            if i%output_frequency == 0:
+                i_out += 1
             if self.fatal_error(_error_status):
                 self.fast_deinit()
                 raise RuntimeError(f"Error {_error_status.value}: {_error_message.value}")
@@ -241,16 +255,14 @@ def total_time_steps(self) -> int:
         # and that's why we have the +1 below
         # 
         # We assume here t_initial is always 0
-        return math.ceil( self.t_max.value / self.dt.value) + 1
+        return math.ceil( self.t_max.value / self.dt.value) + 1  
 
 
     @property
-    def output_channel_names(self) -> List:
-        if len(self.channel_names.value.split()) == 0:
-            return []
-        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 total_output_steps(self) -> int:
+        # From FAST_Subs ValidateInputData: DT_out == DT or DT_out is a multiple of DT
+        # So the number of output steps can be calculated the same as the total time steps
+        return math.ceil(self.t_max.value / self.dt_out.value) + 1
 
 
     def get_hub_position(self) -> Tuple:

glue-codes/simulink/src/FAST_SFunc.c

@@ -45,6 +45,7 @@
 
 
 static double dt = 0;
+static double dt_out = 0;
 static double TMax = 0;
 static int NumInputs = NumFixedInputs;
 static int NumAddInputs = 0;  // number of additional inputs
@@ -203,7 +204,7 @@ static void mdlInitializeSizes(SimStruct *S)
        FAST_AllocateTurbines(&nTurbines, &ErrStat, ErrMsg);
        if (checkError(S)) return;
 
-       FAST_Sizes(&iTurb, InputFileName, &AbortErrLev, &NumOutputs, &dt, &TMax, &ErrStat, ErrMsg, ChannelNames, &TMax, InitInputAry);
+       FAST_Sizes(&iTurb, InputFileName, &AbortErrLev, &NumOutputs, &dt, &dt_out, &TMax, &ErrStat, ErrMsg, ChannelNames, &TMax, InitInputAry);
        n_t_global = -1;
        if (checkError(S)) return;
 

modules/openfast-library/src/FAST_Library.f90

@@ -83,7 +83,7 @@ subroutine FAST_DeallocateTurbines(ErrStat_c, ErrMsg_c) BIND (C, NAME='FAST_Deal
    ErrMsg_c = C_NULL_CHAR
 end subroutine
 !==================================================================================================================================
-subroutine FAST_Sizes(iTurb, InputFileName_c, AbortErrLev_c, NumOuts_c, dt_c, tmax_c, ErrStat_c, ErrMsg_c, ChannelNames_c, TMax, InitInpAry) BIND (C, NAME='FAST_Sizes')
+subroutine FAST_Sizes(iTurb, InputFileName_c, AbortErrLev_c, NumOuts_c, dt_c, dt_out_c, tmax_c, ErrStat_c, ErrMsg_c, ChannelNames_c, TMax, InitInpAry) BIND (C, NAME='FAST_Sizes')
    IMPLICIT NONE 
 #ifndef IMPLICIT_DLLEXPORT
 !DEC$ ATTRIBUTES DLLEXPORT :: FAST_Sizes
@@ -94,6 +94,7 @@ subroutine FAST_Sizes(iTurb, InputFileName_c, AbortErrLev_c, NumOuts_c, dt_c, tm
    INTEGER(C_INT),         INTENT(  OUT) :: AbortErrLev_c      
    INTEGER(C_INT),         INTENT(  OUT) :: NumOuts_c      
    REAL(C_DOUBLE),         INTENT(  OUT) :: dt_c      
+   REAL(C_DOUBLE),         INTENT(  OUT) :: dt_out_c      
    REAL(C_DOUBLE),         INTENT(  OUT) :: tmax_c
    INTEGER(C_INT),         INTENT(  OUT) :: ErrStat_c      
    CHARACTER(KIND=C_CHAR), INTENT(  OUT) :: ErrMsg_c(IntfStrLen) 
@@ -151,6 +152,7 @@ subroutine FAST_Sizes(iTurb, InputFileName_c, AbortErrLev_c, NumOuts_c, dt_c, tm
    AbortErrLev_c = AbortErrLev   
    NumOuts_c     = min(MAXOUTPUTS, SUM( Turbine(iTurb)%y_FAST%numOuts ))
    dt_c          = Turbine(iTurb)%p_FAST%dt
+   dt_out_c      = Turbine(iTurb)%p_FAST%DT_Out
    tmax_c        = Turbine(iTurb)%p_FAST%TMax
 
    ErrStat_c     = ErrStat

modules/openfast-library/src/FAST_Library.h

@@ -27,9 +27,9 @@ EXTERNAL_ROUTINE void FAST_HubPosition(int * iTurb, float * absolute_position, f
 
 EXTERNAL_ROUTINE void FAST_Restart(int * iTurb, const char *CheckpointRootName, int *AbortErrLev, int * NumOuts, double * dt, int * n_t_global, int *ErrStat, char *ErrMsg);
 #ifdef __cplusplus
-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_Sizes(int * iTurb, const char *InputFileName, int *AbortErrLev, int * NumOuts, double * dt, double * dt_out, double * tmax, int *ErrStat, char *ErrMsg, char *ChannelNames, double *TMax = NULL, double *InitInputAry = NULL); 
 #else
-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, double *InitInputAry);
+EXTERNAL_ROUTINE void FAST_Sizes(int * iTurb, const char *InputFileName, int *AbortErrLev, int * NumOuts, double * dt, double * dt_out, double * tmax, int *ErrStat, char *ErrMsg, char *ChannelNames, double *TMax, double *InitInputAry);
 #endif
 EXTERNAL_ROUTINE void FAST_Start(int * iTurb, int *NumInputs_c, int *NumOutputs_c, double *InputAry, double *OutputAry, int *ErrStat, char *ErrMsg);
 EXTERNAL_ROUTINE void FAST_Update(int * iTurb, int *NumInputs_c, int *NumOutputs_c, double *InputAry, double *OutputAry, bool *EndSimulationEarly, int *ErrStat, char *ErrMsg);