Feature/power interpolation (#153)

Moves the interpolation in the power curve from cp, to power itself to avoid issues if wind speed points are spaced out above rated in input JSON

Co-authored-by: ejsimley <eric.simley@nrel.gov>

floris/simulation/turbine.py

@@ -126,6 +126,10 @@ def _initialize_turbine(self):
         # initialize derived attributes
         self.grid = self._create_swept_area_grid()
 
+        # Compute list of inner powers
+        inner_power = np.array([self._power_inner_function(ws) for ws in wind_speed])
+        self.powInterp = interp1d(wind_speed, inner_power, fill_value="extrapolate")
+
     def _create_swept_area_grid(self):
         # TODO: add validity check:
         # rotor points has a minimum in order to always include points inside
@@ -163,6 +167,26 @@ def _create_swept_area_grid(self):
 
         return grid
 
+    def _power_inner_function(self, yaw_effective_velocity):
+        """
+        This method calculates the power for an array of yaw effective wind
+        speeds without the air density and turbulence correction parameters.
+        This is used to initialize the power interpolation method used to
+        compute turbine power.
+        """
+
+        # Now compute the power
+        cptmp = self._fCp(
+            yaw_effective_velocity
+        )  # Note Cp is also now based on yaw effective velocity
+        return (
+            0.5
+            * (np.pi * self.rotor_radius ** 2)
+            * cptmp
+            * self.generator_efficiency
+            * yaw_effective_velocity ** 3
+        )
+
     def _fCp(self, at_wind_speed):
         wind_speed = self.power_thrust_table["wind_speed"]
         if at_wind_speed < min(wind_speed):
@@ -362,9 +386,9 @@ def turbulence_parameter(self):
                 npdf = np.array(pdf) * (1 / np.sum(pdf))
 
                 # calculate turbulence parameter (ratio of corrected power to original power)
-                return np.sum(
-                    [npdf[k] * self._fCp(xp[k]) * xp[k] ** 3 for k in range(100)]
-                ) / (self._fCp(mu) * mu ** 3)
+                return np.sum([npdf[k] * self.powInterp(xp[k]) for k in range(100)]) / (
+                    self.powInterp(mu)
+                )
 
     @property
     def current_turbulence_intensity(self):
@@ -567,15 +591,10 @@ def power(self):
         yaw_effective_velocity = self.average_velocity * cosd(self.yaw_angle) ** pW
 
         # Now compute the power
-        cptmp = self.Cp  # Note Cp is also now based on yaw effective velocity
         return (
-            0.5
-            * self.air_density
-            * (np.pi * self.rotor_radius ** 2)
-            * cptmp
-            * self.generator_efficiency
+            self.air_density
+            * self.powInterp(yaw_effective_velocity)
             * self.turbulence_parameter
-            * yaw_effective_velocity ** 3
         )
 
     @property

floris/tools/floris_interface.py

@@ -1003,8 +1003,30 @@ def get_turbine_power(
             self.calculate_wake(yaw_angles=yaw_angles, no_wake=no_wake)
             return list(mean_farm_power)
         else:
-            turb_powers = [turbine.power for turbine in self.floris.farm.turbines]
-            return turb_powers
+            return [turbine.power for turbine in self.floris.farm.turbines]
+
+    def get_power_curve(self, wind_speeds):
+        """
+        Return the power curve given a set of wind speeds
+
+        Args:
+            wind_speeds (np.array): array of wind speeds to get power curve
+        """
+
+        # Temporarily set the farm to a single turbine
+        saved_layout_x = self.layout_x
+        saved_layout_y = self.layout_y
+        self.reinitialize_flow_field(layout_array=([0], [0]))
+        power_return_array = []
+        for ws in wind_speeds:
+            self.reinitialize_flow_field(wind_speed=ws)
+            self.calculate_wake()
+            power_return_array.append(self.get_turbine_power()[0])
+
+        # Set it back
+        self.reinitialize_flow_field(layout_array=(saved_layout_x, saved_layout_y))
+
+        return np.array(power_return_array)
 
     def get_turbine_ct(self):
         """