fix omdao wrapper and install omdao test so that capability is mainta…

…ined

.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

raft/omdao_raft.py

@@ -1,9 +1,10 @@
 import openmdao.api as om
 import raft
 import numpy as np
-import pickle, os
+import os
 import copy
 from itertools import compress
+import pickle
 
 ndim = 3
 ndof = 6
@@ -255,10 +256,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 +275,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 +284,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 +303,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') 
 
@@ -620,6 +631,9 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
                 pickle.dump(design, handle, protocol=pickle.HIGHEST_PROTOCOL)
             self.i_design += 1
 
+        with open('file.pl', 'w') as pfile:
+            print(design, file=pfile)
+
         # set up the model
         model = raft.Model(design)
         model.analyzeUnloaded(
@@ -633,7 +647,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 +668,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