Feature/add ref density

floris/simulation/farm.py

@@ -86,6 +86,18 @@ def check_turbine_type(self, instance: attrs.Attribute, value: Any) -> None:
                     raise ValueError("User-selected turbine definition `{}` does not exist in pre-defined turbine library.".format(val))
                 self.turbine_definitions[i] = load_yaml(fname)
 
+                # This is a temporary block of code that catches that ref_density_cp_ct is not defined
+                # In the yaml file and forces it in
+                # A warning is issued letting the user know in future versions defining this value explicitly
+                # will be required 
+                if not 'ref_density_cp_ct' in self.turbine_definitions[i]:
+                    self.logger.warn("The value ref_density_cp_ct is not defined in the file: %s " % fname)
+                    self.logger.warn("This value is not the simulated air density but is the density at which the cp/ct curves are defined")
+                    self.logger.warn("In previous versions this was assumed to be 1.225")
+                    self.logger.warn("Future versions of FLORIS will give an error if this value is not explicitly defined")
+                    self.logger.warn("Currently this value is being set to the prior default value of 1.225")
+                    self.turbine_definitions[i]['ref_density_cp_ct'] = 1.225
+
     def initialize(self, sorted_indices):
         # Sort yaw angles from most upstream to most downstream wind turbine
         self.yaw_angles_sorted = np.take_along_axis(
@@ -107,6 +119,9 @@ def construct_turbine_TSRs(self):
     def construc_turbine_pPs(self):
         self.pPs = np.array([turb['pP'] for turb in self.turbine_definitions])
 
+    def construc_turbine_ref_density_cp_cts(self):
+        self.ref_density_cp_cts = np.array([turb['ref_density_cp_ct'] for turb in self.turbine_definitions])
+
     def construct_turbine_map(self):
         self.turbine_map = [Turbine.from_dict(turb) for turb in self.turbine_definitions]
 

floris/simulation/floris.py

@@ -73,6 +73,7 @@ def __attrs_post_init__(self) -> None:
         self.farm.construct_rotor_diameters()
         self.farm.construct_turbine_TSRs()
         self.farm.construc_turbine_pPs()
+        self.farm.construc_turbine_ref_density_cp_cts()
         self.farm.construct_coordinates()
         self.farm.set_yaw_angles(self.flow_field.n_wind_directions, self.flow_field.n_wind_speeds)
 

floris/simulation/solver.py

@@ -224,6 +224,7 @@ def full_flow_sequential_solver(farm: Farm, flow_field: FlowField, flow_field_gr
     turbine_grid_farm.construct_rotor_diameters()
     turbine_grid_farm.construct_turbine_TSRs()
     turbine_grid_farm.construc_turbine_pPs()
+    turbine_grid_farm.construc_turbine_ref_density_cp_cts()
     turbine_grid_farm.construct_coordinates()
 
 
@@ -554,6 +555,7 @@ def full_flow_cc_solver(farm: Farm, flow_field: FlowField, flow_field_grid: Flow
     turbine_grid_farm.construct_rotor_diameters()
     turbine_grid_farm.construct_turbine_TSRs()
     turbine_grid_farm.construc_turbine_pPs()
+    turbine_grid_farm.construc_turbine_ref_density_cp_cts()
     turbine_grid_farm.construct_coordinates()
 
     turbine_grid = TurbineGrid(

floris/simulation/turbine.py

@@ -79,6 +79,7 @@ def _filter_convert(
 
 def power(
     air_density: float,
+    ref_density_cp_ct: float,
     velocities: NDArrayFloat,
     yaw_angle: NDArrayFloat,
     pP: float,
@@ -91,6 +92,7 @@ def power(
 
     Args:
         air_density (NDArrayFloat[wd, ws, turbines]): The air density value(s) at each turbine.
+        ref_density_cp_cts (NDArrayFloat[wd, ws, turbines]): The reference density for each turbine
         velocities (NDArrayFloat[wd, ws, turbines, grid1, grid2]): The velocity field at a turbine.
         pP (NDArrayFloat[wd, ws, turbines]): The pP value(s) of the cosine exponent relating
             the yaw misalignment angle to power for each turbine.
@@ -134,7 +136,7 @@ def power(
 
     # Compute the yaw effective velocity
     pW = pP / 3.0  # Convert from pP to w
-    yaw_effective_velocity = ((air_density/1.225)**(1/3)) * average_velocity(velocities) * cosd(yaw_angle) ** pW
+    yaw_effective_velocity = ((air_density/ref_density_cp_ct)**(1/3)) * average_velocity(velocities) * cosd(yaw_angle) ** pW
 
     # Loop over each turbine type given to get thrust coefficient for all turbines
     p = np.zeros(np.shape(yaw_effective_velocity))
@@ -145,7 +147,7 @@ def power(
         # type to the main thrust coefficient array
         p += power_interp[turb_type](yaw_effective_velocity) * np.array(turbine_type_map == turb_type)
 
-    return p * 1.225
+    return p * ref_density_cp_ct
 
 
 def Ct(
@@ -317,6 +319,8 @@ class Turbine(BaseClass):
             tilt angle to power.
         generator_efficiency (:py:obj: float): The generator
             efficiency factor used to scale the power production.
+        ref_density_cp_ct (:py:obj: float): The density at which the provided
+            cp and ct is defined
         power_thrust_table (PowerThrustTable): A dictionary containing the
             following key-value pairs:
 
@@ -343,8 +347,11 @@ class Turbine(BaseClass):
     pT: float
     TSR: float
     generator_efficiency: float
+    ref_density_cp_ct: float
     power_thrust_table: PowerThrustTable = field(converter=PowerThrustTable.from_dict)
 
+
+
     # rloc: float = float_attrib()  # TODO: goes here or on the Grid?
     # use_points_on_perimeter: bool = bool_attrib()
 
@@ -355,6 +362,7 @@ class Turbine(BaseClass):
     fCt_interp: interp1d = field(init=False)
     power_interp: interp1d = field(init=False)
 
+
     # For the following parameters, use default values if not user-specified
     # self.rloc = float(input_dictionary["rloc"]) if "rloc" in input_dictionary else 0.5
     # if "use_points_on_perimeter" in input_dictionary:

floris/tools/floris_interface.py

@@ -574,6 +574,7 @@ def get_turbine_powers(self) -> NDArrayFloat:
 
         turbine_powers = power(
             air_density=self.floris.flow_field.air_density,
+            ref_density_cp_ct=self.floris.farm.ref_density_cp_cts,
             velocities=self.floris.flow_field.u,
             yaw_angle=self.floris.farm.yaw_angles,
             pP=self.floris.farm.pPs,

floris/turbine_library/iea_10MW.yaml

@@ -5,6 +5,7 @@ pP: 1.88
 pT: 1.88
 rotor_diameter: 198.0
 TSR: 8.0
+ref_density_cp_ct: 1.225
 power_thrust_table:
   power:
     - 0.000000

floris/turbine_library/iea_15MW.yaml

@@ -5,6 +5,7 @@ pP: 1.88
 pT: 1.88
 rotor_diameter: 240.0
 TSR: 8.0
+ref_density_cp_ct: 1.225
 power_thrust_table:
   power:
     - 0.000000

floris/turbine_library/nrel_5MW.yaml

@@ -5,6 +5,7 @@ pP: 1.88
 pT: 1.88
 rotor_diameter: 126.0
 TSR: 8.0
+ref_density_cp_ct: 1.225
 power_thrust_table:
   power:
     - 0.0

floris/type_dec.py

@@ -108,6 +108,7 @@ def from_dict(cls, data: dict):
         # Map the inputs must be provided: 1) must be initialized, 2) no default value defined
         required_inputs = [a.name for a in cls.__attrs_attrs__ if a.init and a.default is attrs.NOTHING]
         undefined = sorted(set(required_inputs) - set(kwargs))
+
         if undefined:
             raise AttributeError(f"The class defintion for {cls.__name__} is missing the following inputs: {undefined}")
         return cls(**kwargs)

tests/conftest.py

@@ -156,6 +156,7 @@ def __init__(self):
             "pP": 1.88,
             "pT": 1.88,
             "generator_efficiency": 1.0,
+            "ref_density_cp_ct": 1.225,
             "power_thrust_table": {
                 "power": [
                     0.000000,

tests/reg_tests/cumulative_curl_regression_test.py

@@ -1,4 +1,4 @@
-# Copyright 2021 NREL
+# Copyright 2022 NREL
 
 # Licensed under the Apache License, Version 2.0 (the "License"); you may not
 # use this file except in compliance with the License. You may obtain a copy of
@@ -174,6 +174,7 @@ def test_regression_tandem(sample_inputs_fixture):
     )
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,
@@ -318,6 +319,7 @@ def test_regression_yaw(sample_inputs_fixture):
     )
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,
@@ -390,6 +392,7 @@ def test_regression_yaw_added_recovery(sample_inputs_fixture):
     )
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,
@@ -461,6 +464,7 @@ def test_regression_secondary_steering(sample_inputs_fixture):
     )
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,
@@ -530,6 +534,7 @@ def test_regression_small_grid_rotation(sample_inputs_fixture):
 
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,

tests/reg_tests/gauss_regression_test.py

@@ -1,4 +1,4 @@
-# Copyright 2021 NREL
+# Copyright 2022 NREL
 
 # Licensed under the Apache License, Version 2.0 (the "License"); you may not
 # use this file except in compliance with the License. You may obtain a copy of
@@ -265,6 +265,7 @@ def test_regression_tandem(sample_inputs_fixture):
     )
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,
@@ -409,6 +410,7 @@ def test_regression_yaw(sample_inputs_fixture):
     )
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,
@@ -478,6 +480,7 @@ def test_regression_gch(sample_inputs_fixture):
     )
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,
@@ -542,6 +545,7 @@ def test_regression_gch(sample_inputs_fixture):
     )
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,
@@ -614,6 +618,7 @@ def test_regression_yaw_added_recovery(sample_inputs_fixture):
     )
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,
@@ -685,6 +690,7 @@ def test_regression_secondary_steering(sample_inputs_fixture):
     )
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,
@@ -754,6 +760,7 @@ def test_regression_small_grid_rotation(sample_inputs_fixture):
 
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,

tests/reg_tests/jensen_jimenez_regression_test.py

@@ -1,4 +1,4 @@
-# Copyright 2020 NREL
+# Copyright 2022 NREL
 
 # Licensed under the Apache License, Version 2.0 (the "License"); you may not
 # use this file except in compliance with the License. You may obtain a copy of
@@ -117,6 +117,7 @@ def test_regression_tandem(sample_inputs_fixture):
     )
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,
@@ -261,6 +262,7 @@ def test_regression_yaw(sample_inputs_fixture):
     )
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,
@@ -330,6 +332,7 @@ def test_regression_small_grid_rotation(sample_inputs_fixture):
 
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,

tests/reg_tests/none_regression_test.py

@@ -1,4 +1,4 @@
-# Copyright 2020 NREL
+# Copyright 2022 NREL
 
 # Licensed under the Apache License, Version 2.0 (the "License"); you may not
 # use this file except in compliance with the License. You may obtain a copy of
@@ -118,6 +118,7 @@ def test_regression_tandem(sample_inputs_fixture):
     )
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,
@@ -283,6 +284,7 @@ def test_regression_small_grid_rotation(sample_inputs_fixture):
 
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,

tests/reg_tests/turbopark_regression_test.py

@@ -1,4 +1,4 @@
-# Copyright 2021 NREL
+# Copyright 2022 NREL
 
 # Licensed under the Apache License, Version 2.0 (the "License"); you may not
 # use this file except in compliance with the License. You may obtain a copy of
@@ -118,6 +118,7 @@ def test_regression_tandem(sample_inputs_fixture):
     )
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,
@@ -263,6 +264,7 @@ def test_regression_yaw(sample_inputs_fixture):
     )
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,
@@ -333,6 +335,7 @@ def test_regression_small_grid_rotation(sample_inputs_fixture):
 
     farm_powers = power(
         floris.flow_field.air_density,
+        floris.farm.ref_density_cp_cts,
         velocities,
         yaw_angles,
         floris.farm.pPs,

tests/turbine_unit_test.py

@@ -259,6 +259,7 @@ def test_power():
     wind_speed = 10.0
     p = power(
         air_density=AIR_DENSITY,
+        ref_density_cp_ct=AIR_DENSITY,
         velocities=wind_speed * np.ones((1, 1, 1, 3, 3)),
         yaw_angle=np.zeros((1, 1, 1)),
         pP=turbine.pP * np.ones((1, 1, 1)),