Merge pull request NREL#16 from NREL/develop

Minor bug fixes, addition of a function to generate Cp surfaces using FAST

Examples/DISCON.IN

@@ -12,7 +12,7 @@
 1                   ! PC_ControlMode    - Blade pitch control mode {0: No pitch, fix to fine pitch, 1: active PI blade pitch control}
 0                   ! Y_ControlMode		- Yaw control mode {0: no yaw control, 1: yaw rate control, 2: yaw-by-IPC}
 1                   ! SS_Mode           - Setpoint Smoother mode {0: no setpoint smoothing, 1: introduce setpoint smoothing}
-0                   ! WE_Mode           - Wind speed estimator mode {0: One-second low pass filtered hub height wind speed, 1: Immersion and Invariance Estimator, 2: Extended Kalman Filter}
+2                   ! WE_Mode           - Wind speed estimator mode {0: One-second low pass filtered hub height wind speed, 1: Immersion and Invariance Estimator, 2: Extended Kalman Filter}
 0                   ! PS_Mode           - Pitch saturation mode {0: no pitch saturation, 1: implement pitch saturation}
 0                   ! SD_Mode           - Shutdown mode {0: no shutdown procedure, 1: pitch to max pitch at shutdown}
 0                   ! Fl_Mode           - Floating specific feedback mode {0: no nacelle velocity feedback, 1: nacelle velocity feedback}
@@ -60,8 +60,8 @@
 43093.51876000      ! VS_RtTq			- Rated torque, [Nm].
 122.9096700000      ! VS_RefSpd			- Rated generator speed [rad/s]
 1                   ! VS_n				- Number of generator PI torque controller gains
--999.026780000      ! VS_KP				- Proportional gain for generator PI torque controller [1/(rad/s) Nm]. (Only used in the transitional 2.5 region if VS_ControlMode =/ 2)
--185.790360000      ! VS_KI				- Integral gain for generator PI torque controller [1/rad Nm]. (Only used in the transitional 2.5 region if VS_ControlMode =/ 2)
+-647.671470000      ! VS_KP				- Proportional gain for generator PI torque controller [1/(rad/s) Nm]. (Only used in the transitional 2.5 region if VS_ControlMode =/ 2)
+-46.4475900000      ! VS_KI				- Integral gain for generator PI torque controller [1/rad Nm]. (Only used in the transitional 2.5 region if VS_ControlMode =/ 2)
 7.51                ! VS_TSRopt			- Power-maximizing region 2 tip-speed-ratio [rad].
 
 !------- SETPOINT SMOOTHER ---------------------------------------------

ROSCO_toolbox/controller.py

@@ -187,13 +187,15 @@ def tune_controller(self, turbine):
             Cp_TSR = np.ndarray.flatten(turbine.Cp.interp_surface(turbine.pitch_initial_rad, TSR_op[i]))     # all Cp values for a given tsr
             Cp_op[i] = np.clip(Cp_op[i], np.min(Cp_TSR), np.max(Cp_TSR))        # saturate Cp values to be on Cp surface
             f_cp_pitch = interpolate.interp1d(Cp_TSR,pitch_initial_rad)         # interpolate function for Cp(tsr) values
-            pitch_op[i] = f_cp_pitch(Cp_op[i])                                  # expected operation blade pitch values
+            if isinstance(self.min_pitch, float):
+                pitch_op[i] = max(self.min_pitch, f_cp_pitch(Cp_op[i]))             # expected operation blade pitch values
+            else:
+                pitch_op[i] = f_cp_pitch(Cp_op[i])                                  # expected operation blade pitch values
             dCp_beta[i], dCp_TSR[i] = turbine.Cp.interp_gradient(pitch_op[i],TSR_op[i])       # gradients of Cp surface in Beta and TSR directions
 
         # Define minimum pitch saturation to be at Cp-maximizing pitch angle if not specifically defined
         if not self.min_pitch:
-            self.min_pitch = 0.0
-            self.min_pitch = max(self.min_pitch,pitch_op[0])
+            self.min_pitch = pitch_op[0]
 
         # Full Cp surface gradients
         dCp_dbeta = dCp_beta/np.diff(pitch_initial_rad)[0]
@@ -229,7 +231,7 @@ def tune_controller(self, turbine):
         self.vs_gain_schedule.second_order_PI(self.zeta_vs, self.omega_vs,A_vs,B_tau,linearize=False,v=v_below_rated)
 
         # -- Find K for Komega_g^2 --
-        self.vs_rgn2K = (pi*rho*R**5.0 * turbine.Cp.max) / (2.0 * turbine.Cp.TSR_opt**3 * Ng**3)
+        self.vs_rgn2K = (pi*rho*R**5.0 * turbine.Cp.max) / (2.0 * turbine.Cp.TSR_opt**3 * Ng**3)/ (turbine.GenEff/100 * turbine.GBoxEff/100)
         self.vs_refspd = min(turbine.TSR_operational * turbine.v_rated/R, turbine.rated_rotor_speed) * Ng
 
         # -- Define some setpoints --
@@ -448,7 +450,7 @@ def peak_shaving(self,controller, turbine):
                 Ct_max[i] = np.minimum( np.max(Ct_tsr), Ct_max[i])
             # Define minimum pitch angle
             f_pitch_min = interpolate.interp1d(Ct_tsr, turbine.pitch_initial_rad, bounds_error=False, fill_value=(turbine.pitch_initial_rad[0],turbine.pitch_initial_rad[-1]))
-            pitch_min[i] = f_pitch_min(Ct_max[i])
+            pitch_min[i] = max(controller.min_pitch, f_pitch_min(Ct_max[i]))
 
         controller.ps_min_bld_pitch = pitch_min
 

ROSCO_toolbox/turbine.py

@@ -17,6 +17,7 @@
 from numpy import gradient
 import pickle
 import matplotlib.pyplot as plt
+import pandas as pd
 
 from ROSCO_toolbox import utilities as ROSCO_utilities
 
@@ -46,6 +47,7 @@ class Turbine():
     load_from_fast
     load_from_ccblade
     load_from_txt
+    generate_rotperf_fast
     write_rotor_performance
 
     Parameters:
@@ -194,7 +196,7 @@ def load_from_fast(self, FAST_InputFile,FAST_directory, FAST_ver='OpenFAST',dev_
             self.precone = fast.fst_vt['ElastoDyn']['PreCone(1)']
         self.yaw = 0.0
         self.J = self.rotor_inertia + self.generator_inertia * self.Ng**2
-        self.rated_torque = self.rated_power/(self.GenEff/100*self.rated_rotor_speed*self.Ng)
+        self.rated_torque = self.rated_power/(self.GenEff/100*self.rated_rotor_speed*self.Ng*self.GBoxEff/100)
         self.max_torque = self.rated_torque * 1.1
         self.rotor_radius = self.TipRad
         # self.omega_dt = np.sqrt(self.DTTorSpr/self.J)
@@ -208,15 +210,15 @@ def load_from_fast(self, FAST_InputFile,FAST_directory, FAST_ver='OpenFAST',dev_
                 txt_filename)
         else:   # Use text file from DISCON.in
             if os.path.exists(os.path.join(FAST_directory, fast.fst_vt['ServoDyn']['DLL_InFile'])):
-                if  os.path.exists(fast.fst_vt['DISCON_in']['PerfFileName']):
+                try:
                     self.pitch_initial_rad = fast.fst_vt['DISCON_in']['Cp_pitch_initial_rad']
                     self.TSR_initial = fast.fst_vt['DISCON_in']['Cp_TSR_initial']
                     self.Cp_table = fast.fst_vt['DISCON_in']['Cp_table']
                     self.Ct_table = fast.fst_vt['DISCON_in']['Ct_table']
                     self.Cq_table = fast.fst_vt['DISCON_in']['Cq_table']
-            else:   # Load from cc-blade
-                print('No rotor performance data source available, running CC-Blade.')
-                self.load_from_ccblade()
+                except:   # Load from cc-blade
+                    print('No rotor performance data source available, running CC-Blade.')
+                    self.load_from_ccblade()
 
         # Parse rotor performance data
         self.Cp = RotorPerformance(self.Cp_table,self.pitch_initial_rad,self.TSR_initial)
@@ -332,6 +334,200 @@ def load_from_ccblade(self):
         self.chord = chord
         self.twist = theta
 
+    def generate_rotperf_fast(self, openfast_path, FAST_runDirectory=None, run_BeamDyn=False,
+                              debug_level=1, run_type='multi'):
+        '''
+        Use openfast to generate Cp surface data. Will be slow, especially if using BeamDyn,
+        but may be necessary if cc-blade is not sufficient.
+
+        Parameters:
+        -----------
+        openfast_path: str
+            path to openfast
+        FAST_runDirectory: str
+            directory to run openfast simulations in
+        run_BeamDyn: bool
+            Flag to run beamdyn - does not exist yet
+        debug_level: float
+            0 - no outputs, 1 - simple outputs, 2 - all outputs
+        run_type: str
+            'serial' - run in serial, 'multi' - run using python multiprocessing tools, 
+            'mpi' - run using mpi tools
+        '''
+
+        # Load additional WISDEM tools
+        from wisdem.aeroelasticse import runFAST_pywrapper, CaseGen_General
+        from wisdem.aeroelasticse.Util import FileTools
+        # Load pCrunch tools
+        from pCrunch import pdTools, Processing
+
+
+        # setup values for surface
+        v0 = self.v_rated + 2
+        TSR_initial = np.arange(3, 15,1)
+        pitch_initial = np.arange(-1,25,1)
+        rotspeed_initial = TSR_initial*v0/self.rotor_radius * RadSec2rpm # rpms
+
+        # Specify Case Inputs
+        case_inputs = {}
+
+        # ------- Setup OpenFAST inputs --------
+        case_inputs[('Fst','TMax')] = {'vals': [330], 'group': 0}
+        case_inputs[('Fst','Compinflow')] = {'vals': [1], 'group': 0}
+        case_inputs[('Fst','CompAero')] = {'vals': [2], 'group': 0}
+        case_inputs[('Fst','CompServo')] = {'vals': [1], 'group': 0}
+        case_inputs[('Fst','CompHydro')] = {'vals': [0], 'group': 0}
+        if run_BeamDyn:
+            case_inputs[('Fst','CompElast')] = {'vals': [2], 'group': 0}
+        else:
+            case_inputs[('Fst','CompElast')] = {'vals': [1], 'group': 0}
+        case_inputs[('Fst', 'OutFileFmt')] = {'vals': [2], 'group': 0}
+
+        # AeroDyn15
+        case_inputs[('AeroDyn15', 'WakeMod')] = {'vals': [1], 'group': 0}
+        case_inputs[('AeroDyn15', 'AfAeroMod')] = {'vals': [1], 'group': 0}
+        case_inputs[('AeroDyn15', 'TwrPotent')] = {'vals': [0], 'group': 0}
+
+        # ElastoDyn
+        case_inputs[('ElastoDyn', 'FlapDOF1')] = {'vals': ['True'], 'group': 0}
+        case_inputs[('ElastoDyn', 'FlapDOF2')] = {'vals': ['True'], 'group': 0}
+        case_inputs[('ElastoDyn', 'EdgeDOF')] = {'vals': ['True'], 'group': 0}
+        case_inputs[('ElastoDyn', 'TeetDOF')] = {'vals': ['False'], 'group': 0}
+        case_inputs[('ElastoDyn', 'DrTrDOF')] = {'vals': ['False'], 'group': 0}
+        case_inputs[('ElastoDyn', 'GenDOF')] = {'vals': ['False'], 'group': 0}
+        case_inputs[('ElastoDyn', 'YawDOF')] = {'vals': ['False'], 'group': 0}
+        case_inputs[('ElastoDyn', 'TwFADOF1')] = {'vals': ['False'], 'group': 0}
+        case_inputs[('ElastoDyn', 'TwFADOF2')] = {'vals': ['False'], 'group': 0}
+        case_inputs[('ElastoDyn', 'TwSSDOF1')] = {'vals': ['False'], 'group': 0}
+        case_inputs[('ElastoDyn', 'TwSSDOF2')] = {'vals': ['False'], 'group': 0}
+        case_inputs[('ElastoDyn', 'PtfmSgDOF')] = {'vals': ['False'], 'group': 0}
+        case_inputs[('ElastoDyn', 'PtfmSwDOF')] = {'vals': ['False'], 'group': 0}
+        case_inputs[('ElastoDyn', 'PtfmHvDOF')] = {'vals': ['False'], 'group': 0}
+        case_inputs[('ElastoDyn', 'PtfmRDOF')] = {'vals': ['False'], 'group': 0}
+        case_inputs[('ElastoDyn', 'PtfmPDOF')] = {'vals': ['False'], 'group': 0}
+        case_inputs[('ElastoDyn', 'PtfmYDOF')] = {'vals': ['False'], 'group': 0}
+
+        # BeamDyn
+        # NEEDED
+
+        # InflowWind
+        case_inputs[('InflowWind', 'WindType')] = {'vals': [1], 'group': 0}
+        case_inputs[('InflowWind', 'HWindSpeed')] = {'vals': [v0], 'group': 0}
+        case_inputs[('InflowWind', 'PLexp')] = {'vals': [0], 'group': 0}
+
+        # ServoDyn
+        case_inputs[('ServoDyn', 'PCMode')] = {'vals': [0], 'group': 0}
+        case_inputs[('ServoDyn', 'VSContrl')] = {'vals': [0], 'group': 0} 
+        case_inputs[('ServoDyn', 'HSSBrMode')] = {'vals': [0], 'group': 0}
+        case_inputs[('ServoDyn', 'YCMode')] = {'vals': [0], 'group': 0}
+
+        # ------- Setup sweep values inputs --------
+        case_inputs[('ElastoDyn', 'BlPitch1')] = {'vals': list(pitch_initial), 'group': 1}
+        case_inputs[('ElastoDyn', 'BlPitch2')] = {'vals': list(pitch_initial), 'group': 1}
+        case_inputs[('ElastoDyn', 'BlPitch3')] = {'vals': list(pitch_initial), 'group': 1}
+        case_inputs[('ElastoDyn', 'RotSpeed')] = {'vals': list(rotspeed_initial), 'group': 2}
+
+
+        # FAST details
+        fastBatch = runFAST_pywrapper.runFAST_pywrapper_batch(FAST_ver='OpenFAST', dev_branch=True)
+        fastBatch.FAST_exe = openfast_path  # Path to executable
+        fastBatch.FAST_InputFile = self.fast.FAST_InputFile
+        fastBatch.FAST_directory = self.fast.FAST_directory
+        if not FAST_runDirectory:
+            FAST_runDirectory = os.path.join(os.getcwd(), 'RotPerf_OpenFAST')
+        fastBatch.FAST_runDirectory = FAST_runDirectory
+        fastBatch.debug_level = debug_level
+
+        # Generate cases
+        case_name_base = self.TurbineName + '_rotperf'
+        case_list, case_name_list = CaseGen_General.CaseGen_General(
+            case_inputs, dir_matrix=fastBatch.FAST_runDirectory, namebase=case_name_base)
+        fastBatch.case_list = case_list
+        fastBatch.case_name_list = case_name_list
+
+        # Make sure proper outputs exist
+        var_out = [
+            # ElastoDyn (this is probably overkill on the outputs)
+            "BldPitch1", "BldPitch2", "BldPitch3", "Azimuth", "RotSpeed", "GenSpeed", "NacYaw",
+            "OoPDefl1", "IPDefl1", "TwstDefl1", "OoPDefl2", "IPDefl2", "TwstDefl2", "OoPDefl3",
+            "IPDefl3", "TwstDefl3", "RootFxc1",
+            "RootFyc1", "RootFzc1", "RootMxc1", "RootMyc1", "RootMzc1", "RootFxc2", "RootFyc2",
+            "RootFzc2", "RootMxc2", "RootMyc2", "RootMzc2", "RootFxc3", "RootFyc3", "RootFzc3",
+            "RootMxc3", "RootMyc3", "RootMzc3", "Spn1MLxb1", "Spn1MLyb1", "Spn1MLzb1", "Spn1MLxb2",
+            "Spn1MLyb2", "Spn1MLzb2", "Spn1MLxb3", "Spn1MLyb3", "Spn1MLzb3", "RotThrust", "LSSGagFya",
+            "LSSGagFza", "RotTorq", "LSSGagMya", "LSSGagMza", 
+            # ServoDyn
+            "GenPwr", "GenTq",
+            # AeroDyn15
+            "RtArea", "RtVAvgxh", "B1N3Clrnc", "B2N3Clrnc", "B3N3Clrnc",
+            "RtAeroCp", 'RtAeroCq', 'RtAeroCt', 'RtTSR', # NECESSARY
+            # InflowWind
+            "Wind1VelX", 
+        ]
+        channels = {}
+        for var in var_out:
+            channels[var] = True
+        fastBatch.channels = channels
+
+        # Run OpenFAST
+        if run_type.lower() == 'multi':
+            fastBatch.run_multi()
+        elif run_type.lower()=='mpi':
+            fastBatch.run_mpi()
+        elif run_type.lower()=='serial':
+            fastBatch.run_serial()
+
+        # ========== Post Processing ==========
+        # Save statistics
+        fp = Processing.FAST_Processing()
+
+        # Find all outfiles
+        fname_case_matrix = os.path.join(FAST_runDirectory, 'case_matrix.yaml')
+        case_matrix = FileTools.load_yaml(fname_case_matrix, package=1)
+        cm = pd.DataFrame(case_matrix)
+        # Parse case matrix and find outfiles names
+        outfiles = []
+        case_names = cm['Case_Name']
+        outfiles = []
+        for name in case_names:
+            outfiles.append(os.path.join(FAST_runDirectory, name + '.outb'))
+
+
+
+        # Set some processing parameters
+        fp.OpenFAST_outfile_list = outfiles
+        fp.namebase = case_name_base
+        fp.t0 = 270 
+        fp.parallel_analysis = True
+        fp.results_dir = os.path.join(FAST_runDirectory,'stats')
+        fp.verbose = True
+        # Save for debug!
+        fp.save_LoadRanking = False
+        fp.save_SummaryStats = False
+
+        print('Processing openfast data on {} cores.'.format(fp.parallel_cores))
+
+        # Load and save statistics and load rankings
+        stats, load_rankings = fp.batch_processing()
+
+        # Get means of last 30 seconds of 300 second simulation
+        CP = stats[0]['RtAeroCp']['mean']
+        CT = stats[0]['RtAeroCt']['mean']
+        CQ = stats[0]['RtAeroCq']['mean']
+
+        # Reshape Cp, Ct and Cq
+        Cp = np.transpose(np.reshape(CP, (len(pitch_initial), len(TSR_initial))))
+        Ct = np.transpose(np.reshape(CT, (len(pitch_initial), len(TSR_initial))))
+        Cq = np.transpose(np.reshape(CQ, (len(pitch_initial), len(TSR_initial))))
+
+        # Store necessary metrics for analysis
+        self.pitch_initial_rad = pitch_initial * deg2rad
+        self.TSR_initial = TSR_initial
+        self.Cp_table = Cp
+        self.Ct_table = Ct
+        self.Cq_table = Cq
+
+
     def load_blade_info(self):
         '''
         Loads wind turbine blade data from an OpenFAST model. 

ROSCO_toolbox/utilities.py

@@ -598,6 +598,9 @@ def trim_output(self, fast_data, tmin=None, tmax=None, verbose=False):
             else: 
                 Tfind = len(fd['Time'])
 
+            if T0ind+1 > len(fd['Time']):
+                raise ValueError('The initial time to trim {} to is after the end of the simulation.'.format(fd['meta']['name']))
+
             # # Modify time
             fd['Time'] = fd['Time'][T0ind:Tfind] - fd['Time'][T0ind]
 
@@ -829,10 +832,10 @@ def write_rotor_performance(self,turbine,txt_filename='Cp_Ct_Cq.txt'):
         file.write('# ------------ Written on {} using the ROSCO toolbox ------------ \n\n'.format(now.strftime('%b-%d-%y')))
 
         # Pitch angles, TSR, and wind speed
-        file.write('# Pitch angle vector - x axis (matrix columns) (deg)\n')
+        file.write('# Pitch angle vector, {} entries - x axis (matrix columns) (deg)\n'.format(len(turbine.Cp.pitch_initial_rad)))
         for i in range(len(turbine.Cp.pitch_initial_rad)):
             file.write('{:0.4}   '.format(turbine.Cp.pitch_initial_rad[i] * rad2deg))
-        file.write('\n# TSR vector - y axis (matrix rows) (-)\n')
+        file.write('\n# TSR vector, {} entries - y axis (matrix rows) (-)\n'.format(len(turbine.TSR_initial)))
         for i in range(len(turbine.TSR_initial)):
             file.write('{:0.4}    '.format(turbine.Cp.TSR_initial[i]))
         file.write('\n# Wind speed vector - z axis (m/s)\n')

setup.py

@@ -95,15 +95,15 @@ def finalize_options(self):
 
     def run(self):
         try:
-            self.status('Removing previous builds…')
+            self.status('Removing previous builds...')
             rmtree(os.path.join(here, 'dist'))
         except OSError:
             pass
 
-        self.status('Building Source and Wheel (universal) distribution…')
+        self.status('Building Source and Wheel (universal) distribution...')
         os.system('{0} setup.py sdist bdist_wheel --universal'.format(sys.executable))
 
-        self.status('Uploading the package to PyPI via Twine…')
+        self.status('Uploading the package to PyPI via Twine...')
         os.system('twine upload dist/*')
 
         self.status('Pushing git tags…')