Merge pull request #41 from WISDEM/omdao

Fix OpenMDAO wrapper

.github/workflows/CI_RAFT.yml

@@ -20,20 +20,21 @@ jobs:
 
         steps:
             - name: checkout repository
-              uses: actions/checkout@v3
+              uses: actions/checkout@v4
 
             - uses: conda-incubator/setup-miniconda@v2
               # https://github.com/marketplace/actions/setup-miniconda
               with:
-                  mamba-version: "*"
-                  miniconda-version: "latest"
-                  #auto-update-conda: true
+                  # To use mamba, uncomment here, comment out the miniforge line
+                  #mamba-version: "*"
+                  miniforge-version: "latest"
+                  auto-update-conda: true
                   python-version: ${{ matrix.python-version }}
                   environment-file: environment.yml
-                  channels: conda-forge
-                  channel-priority: true          
                   activate-environment: test
                   auto-activate-base: false
+                  channels: conda-forge
+                  channel-priority: true          
 
             # Install 
             - name: Conda Install RAFT

examples/VolturnUS-S_example.yaml

@@ -21,7 +21,8 @@ cases:
     keys : [wind_speed, wind_heading, turbulence, turbine_status, yaw_misalign, wave_spectrum, wave_period, wave_height, wave_heading  ]
     data :  #   m/s        deg    % or e.g. IIB_NTM    string            deg         string          (s)         (m)         (deg)
         -  [     0,         0,            0,          operating,          0,        JONSWAP,         12,         6,           0        ]
-        -  [    16,         0,            0,          operating,          0,        JONSWAP,         12,         6,           30       ]
+        -  [    16,         0,            IIB_NTM,    operating,          0,        JONSWAP,         12,         6,           30       ]
+        -  [    16,         0,            0.1,        operating,          0,        JONSWAP,         12,         6,           30       ]
 
 
 turbine:

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "OpenRAFT"
-version = "1.1.1"
+version = "1.2.0"
 description = "RAFT: Response Amplitudes of Floating Turbines"
 readme = "README.md"
 requires-python = ">=3.9"
@@ -38,6 +38,7 @@ classifiers = [  # Optional
   "Programming Language :: Python :: 3.9",
   "Programming Language :: Python :: 3.10",
   "Programming Language :: Python :: 3.11",
+  "Programming Language :: Python :: 3.12",
   "Programming Language :: Python :: 3 :: Only",
 
   "Operating System :: OS Independent",

raft/omdao_raft.py

@@ -5,6 +5,8 @@
 import copy
 from itertools import compress
 
+DEBUG_OMDAO = False  # use within WEIS
+
 ndim = 3
 ndof = 6
 class RAFT_OMDAO(om.ExplicitComponent):
@@ -21,6 +23,18 @@ def initialize(self):
 
     def setup(self):
 
+        # Save options for RAFT testing/debugging
+        if DEBUG_OMDAO:
+            from weis.aeroelasticse.FileTools import save_yaml
+            all_options = {}
+            all_options['modeling_options']     = self.options['modeling_options']
+            all_options['turbine_options']      = self.options['turbine_options']
+            all_options['mooring_options']      = self.options['mooring_options']
+            all_options['member_options']       = self.options['member_options']
+            all_options['analysis_options']     = self.options['analysis_options']
+            save_yaml(os.path.join(os.path.dirname(__file__), '../tests/test_data/'), 'weis_options.yaml', all_options)
+
+
         # unpack options
         modeling_opt = self.options['modeling_options']
         analysis_options = self.options['analysis_options']
@@ -255,10 +269,11 @@ def setup(self):
         self.add_output('properties_roll inertia at subCG', val=np.zeros(ndim), units='kg*m**2', desc='Roll inertia sub CG')
         self.add_output('properties_pitch inertia at subCG', val=np.zeros(ndim), units='kg*m**2', desc='Pitch inertia sub CG')
         self.add_output('properties_yaw inertia at subCG', val=np.zeros(ndim), units='kg*m**2', desc='Yaw inertia sub CG')
-        self.add_output('properties_Buoyancy (pgV)', val=0.0, units='N', desc='Buoyancy (pgV)')
-        self.add_output('properties_Center of Buoyancy', val=np.zeros(ndim), units='m', desc='Center of buoyancy')
-        self.add_output('properties_C stiffness matrix', val=np.zeros((ndof,ndof)), units='Pa', desc='C stiffness matrix')
-        self.add_output('properties_F_lines0', val=np.zeros(nconnections), units='N', desc='Mean mooring force')
+        self.add_output('properties_buoyancy (pgV)', val=0.0, units='N', desc='Buoyancy (pgV)')
+        self.add_output('properties_center of buoyancy', val=np.zeros(ndim), units='m', desc='Center of buoyancy')
+        self.add_output('properties_C hydrostatic', val=np.zeros((ndof,ndof)), units='Pa', desc='C stiffness matrix')
+        self.add_output('properties_C system', val=np.zeros((ndof,ndof)), units='Pa', desc='C stiffness matrix of full system')
+        self.add_output('properties_F_lines0', val=np.zeros(ndof), units='N', desc='Mean mooring force from all lines')
         self.add_output('properties_C_lines0', val=np.zeros((ndof,ndof)), units='Pa', desc='Mooring stiffness')
         self.add_output('properties_M support structure', val=np.zeros((ndof,ndof)), units='kg', desc='Mass matrix for platform')
         self.add_output('properties_A support structure', val=np.zeros((ndof,ndof)), desc='Added mass matrix for platform')
@@ -273,7 +288,7 @@ def setup(self):
         self.add_output('response_roll RAO', val=np.zeros(nfreq), units='rad', desc='Roll RAO')
         self.add_output('response_yaw RAO', val=np.zeros(nfreq), units='rad', desc='Yaw RAO')
         self.add_output('response_nacelle acceleration', val=np.zeros(nfreq), units='m/s**2', desc='Nacelle acceleration')
-        # case specific
+        # case specific, note: only DLCs supported in RAFT will have non-zero outputs
         names = ['surge','sway','heave','roll','pitch','yaw','AxRNA','Mbase','omega','torque','power','bPitch','Tmoor']
         stats = ['avg','std','max','PSD','DEL']
         for n in names:
@@ -282,9 +297,9 @@ def setup(self):
                 iout = f'{n}_{s}'
 
                 if n == 'Tmoor':
-                    myval = np.zeros((n_cases, 2*nlines)) if s not in ['PSD'] else np.zeros((n_cases, 2*nlines, nfreq))
+                    myval = np.zeros((n_cases, 2*nlines)) if s not in ['PSD'] else np.zeros((n_cases, nfreq, 2*nlines))
                 elif n == 'Mbase':
-                    myval = np.zeros(n_cases) if s not in ['PSD','std'] else np.zeros((n_cases, nfreq))
+                    myval = np.zeros(n_cases) if s not in ['PSD'] else np.zeros((n_cases, nfreq))
                 else:
                     myval = np.zeros(n_cases) if s not in ['PSD'] else np.zeros((n_cases, nfreq))
 
@@ -301,13 +316,22 @@ def setup(self):
         # Other case outputs
         self.add_output('stats_wind_PSD', val=np.zeros((n_cases,nfreq)), desc='Power spectral density of wind input')
         self.add_output('stats_wave_PSD', val=np.zeros((n_cases,nfreq)), desc='Power spectral density of wave input')
-
+
+        # Natural periods
+        self.add_output('rigid_body_periods', val = np.zeros(6), desc = 'Rigid body natural period', units = 's') 
+        self.add_output('surge_period', val = 0, desc = 'Surge natural period', units = 's') 
+        self.add_output('sway_period', val = 0, desc = 'Sway natural period', units = 's') 
+        self.add_output('heave_period', val = 0, desc = 'Heave natural period', units = 's') 
+        self.add_output('roll_period', val = 0, desc = 'Roll natural period', units = 's') 
+        self.add_output('pitch_period', val = 0, desc = 'Pitch natural period', units = 's') 
+        self.add_output('yaw_period', val = 0, desc = 'Yaw natural period', units = 's') 
+
         # Aggregate outputs
         self.add_output('Max_Offset', val = 0, desc = 'Maximum distance in surge/sway direction', units = 'm') 
         self.add_output('heave_avg', val = 0, desc = 'Average heave over all cases', units = 'm') 
         self.add_output('Max_PtfmPitch', val = 0, desc = 'Maximum platform pitch over all cases', units = 'deg') 
         self.add_output('Std_PtfmPitch', val = 0, desc = 'Average platform pitch std. over all cases', units = 'deg') 
-        self.add_output('max_nacelle_Ax', val = 0, desc = 'Maximum nacelle accelleration over all cases', units = 'm/s**2') 
+        self.add_output('max_nac_accel', val = 0, desc = 'Maximum nacelle accelleration over all cases', units = 'm/s**2') 
         self.add_output('rotor_overspeed', val = 0, desc = 'Fraction above rated rotor speed') 
         self.add_output('max_tower_base', val = 0, desc = 'Maximum tower base moment over all cases', units = 'N*m') 
 
@@ -347,6 +371,17 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
         nline_types = mooring_opt['nline_types']
         nconnections = mooring_opt['nconnections']
 
+        # save inputs for RAFT testing/debugging
+        if DEBUG_OMDAO:
+            from weis.aeroelasticse.FileTools import save_yaml
+            input_list = self.list_inputs(out_stream=None)
+            input_dict = {}
+            for i in input_list:
+                input_dict[i[0]] = i[1]['val']
+
+            save_yaml(os.path.join(os.path.dirname(__file__), '../tests/test_data/'), 'weis_inputs.yaml', input_dict)
+
+
         # set up design
         design = {}
         design['type'] = ['input dictionary for RAFT']
@@ -633,7 +668,7 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
 
         # option to run level 1 load cases
         if True: #processCases:
-            model.analyzeCases(runPyHAMS=modeling_opt['runPyHAMS'], meshDir=modeling_opt['BEM_dir'])
+            model.analyzeCases(meshDir=modeling_opt['BEM_dir'])
 
         # get and process results
         results = model.calcOutputs()
@@ -654,25 +689,40 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
             '''
 
         # Pattern matching for case-by-case outputs
-        names = ['surge','sway','heave','roll','pitch','yaw','AxRNA','Mbase','omega','torque','power','bPitch','Tmoor']
-        stats = ['avg','std','max','PSD','DEL']
+        # names = ['surge','sway','heave','roll','pitch','yaw','AxRNA','Mbase','omega','torque','power','bPitch','Tmoor']   # these are not always outputted, need to catch better
+        names = ['surge','sway','heave','roll','pitch','yaw','AxRNA','Mbase','Tmoor']
+        stats = ['avg','std','max','PSD']
         case_mask = np.array(case_mask)
         for n in names:
             for s in stats:
-                if s == 'DEL' and not n in ['Tmoor','Mbase']: continue
                 iout = f'{n}_{s}'
-                outputs['stats_'+iout][case_mask] = np.squeeze( results['case_metrics'][iout] )
+
+                # use only first rotor/turbine
+                case_metrics = [cm[0] for cm in results['case_metrics'].values()]
+                stat = np.squeeze(np.array([cm[iout] for cm in case_metrics]))
+                outputs['stats_'+iout][case_mask] = stat
 
         # Other case outputs
-        for n in ['wind_PSD','wave_PSD']:
-            outputs['stats_'+n][case_mask,:] = results['case_metrics'][n]
+        # for n in ['wind_PSD','wave_PSD']:
+        #     outputs['stats_'+n][case_mask,:] = results['case_metrics'][n]
+
+        # natural periods
+        model.solveEigen()
+        outputs["rigid_body_periods"] = 1/results['eigen']['frequencies']
+
+        outputs["surge_period"] = outputs["rigid_body_periods"][0]
+        outputs["sway_period"] = outputs["rigid_body_periods"][1]
+        outputs["heave_period"] = outputs["rigid_body_periods"][2]
+        outputs["roll_period"] = outputs["rigid_body_periods"][3]
+        outputs["pitch_period"] = outputs["rigid_body_periods"][4]
+        outputs["yaw_period"] = outputs["rigid_body_periods"][5]
 
         # Compute some aggregate outputs manually
         outputs['Max_Offset'] = np.sqrt(outputs['stats_surge_max'][case_mask]**2 + outputs['stats_sway_max'][case_mask]**2).max()
         outputs['heave_avg'] = outputs['stats_heave_avg'][case_mask].mean()
         outputs['Max_PtfmPitch'] = outputs['stats_pitch_max'][case_mask].max()
         outputs['Std_PtfmPitch'] = outputs['stats_pitch_std'][case_mask].mean()
-        outputs['max_nacelle_Ax'] = outputs['stats_AxRNA_std'][case_mask].max()
+        outputs['max_nac_accel'] = outputs['stats_AxRNA_std'][case_mask].max()
         outputs['rotor_overspeed'] = (outputs['stats_omega_max'][case_mask].max() - inputs['rated_rotor_speed']) / inputs['rated_rotor_speed']
         outputs['max_tower_base'] = outputs['stats_Mbase_max'][case_mask].max()
 

raft/raft_rotor.py

@@ -48,11 +48,10 @@ def __init__(self, turbine, w, ir):
         self.coords = getFromDict(turbine, 'rotorCoords', dtype=list, shape=turbine['nrotors'], default=[[0,0]])[ir]
         self.speed_gain = getFromDict(turbine, 'speed_gain', shape=turbine['nrotors'], default=1.0)[ir]
 
-        self.headings   = getFromDict(turbine, 'headings', shape=-1, default=[90,210,330])  # [deg]
-        if 'headings' not in turbine:
-            self.nBlades    = getFromDict(turbine, 'nBlades', shape=turbine['nrotors'], dtype=int)[ir]         # [-]
-        else:
-            self.nBlades = len(self.headings)
+        self.nBlades    = getFromDict(turbine, 'nBlades', shape=turbine['nrotors'], dtype=int)[ir]         # [-]
+
+        default_headings    = list(np.arange(self.nBlades) * 360 / self.nBlades) # equally distribute blades
+        self.headings       = getFromDict(turbine, 'headings', shape=-1, default=default_headings)  # [deg]
 
         self.axis       = getFromDict(turbine, 'axis', shape=turbine['nrotors'], default=[1,0,0])[ir]  # unit vector of rotor axis, facing downflow [-]
 
@@ -1066,10 +1065,10 @@ def IECKaimal(self, case, current=False):        #
 
         if current:
             speed = getFromDict(case, 'current_speed', shape=0, default=1.0)
-            turbulence = getFromDict(case, 'current_turbulence', shape=0, default=0.0)
+            turbulence = getFromDict(case, 'current_turbulence', shape=0, default=0.0,dtype=str)
         else:
             speed = getFromDict(case, 'wind_speed', shape=0, default=10.0)
-            turbulence = getFromDict(case, 'turbulence', shape=0, default=0.0)
+            turbulence = getFromDict(case, 'turbulence', shape=0, default=0.0,dtype=str)
 
         # Set inputs (f, V_ref, HH, Class, Categ, TurbMod, R)
         f = self.w / 2 / np.pi    # frequency in Hz
@@ -1080,7 +1079,9 @@ def IECKaimal(self, case, current=False):        #
         ###### Initialize IEC Wind parameters #######
         iec_wind = pyIECWind_extreme()
         iec_wind.z_hub = HH
-
+
+        # Turbulence can be either a string (IB) or a float (turbulence intensity, 0.1)
+        # If a TI is provided, the class defaults to I
         if isinstance(turbulence,str):
             # If a string, the options are I, II, III, IV
             Class = ''
@@ -1091,15 +1092,19 @@ def IECKaimal(self, case, current=False):        #
                     break
 
             if not Class:
-                raise Exception(f"Turbulence class must start with I, II, III, or IV: case['turbulence'] = {turbulence}")
+                Class = 'I'
+                try:
+                    turbulence = float(turbulence)
+                except:
+                    raise Exception(f"Turbulence class must start with I, II, III, or IV: case['turbulence'] = {turbulence}")
             else:
                 Categ = char
                 iec_wind.Turbulence_Class = Categ
 
-            try:
-                TurbMod = turbulence.split('_')[1]
-            except:
-                raise Exception(f"Error reading the turbulence model: {turbulence}")
+                try:
+                    TurbMod = turbulence.split('_')[1]
+                except:
+                    raise Exception(f"Error reading the turbulence model: {turbulence}")
 
             iec_wind.Turbine_Class = Class