Reduce computation time massively in large het_map objects (#1024)

* Reduce computation time massively in large het_map objects by avoiding the recalculation of the delaunay triangulation for every findex

* Import copy library

* Enforce appropriate shape for interpolant object

* ruff formatting

* bugfix

* Add a test of applying a het map

* Add convert to array

* Add convert to array

* Add to new tests

* Clean up

* Add comments and rename variables for clarity.

* Change FlowField.het_map to be a numpy array rather than a list.

---------

Co-authored-by: Paul <[email protected]>
Co-authored-by: misi9170 <[email protected]>

floris/core/flow_field.py

@@ -1,6 +1,8 @@
 
 from __future__ import annotations
 
+import copy
+
 import attrs
 import matplotlib.path as mpltPath
 import numpy as np
@@ -16,6 +18,7 @@
 from floris.type_dec import (
     floris_array_converter,
     NDArrayFloat,
+    NDArrayObject,
 )
 
 
@@ -41,7 +44,7 @@ class FlowField(BaseClass):
     u: NDArrayFloat = field(init=False, factory=lambda: np.array([]))
     v: NDArrayFloat = field(init=False, factory=lambda: np.array([]))
     w: NDArrayFloat = field(init=False, factory=lambda: np.array([]))
-    het_map: list = field(init=False, default=None)
+    het_map: NDArrayObject = field(init=False, default=None)
     dudz_initial_sorted: NDArrayFloat = field(init=False, factory=lambda: np.array([]))
 
     turbulence_intensity_field: NDArrayFloat = field(init=False, factory=lambda: np.array([]))
@@ -281,29 +284,46 @@ def generate_heterogeneous_wind_map(self):
                 - **y**: A list of y locations at which the speed up factors are defined.
                 - **z** (optional): A list of z locations at which the speed up factors are defined.
         """
-        speed_multipliers = self.heterogeneous_inflow_config['speed_multipliers']
+        speed_multipliers = np.array(self.heterogeneous_inflow_config['speed_multipliers'])
         x = self.heterogeneous_inflow_config['x']
         y = self.heterogeneous_inflow_config['y']
         z = self.heterogeneous_inflow_config['z']
 
+        # Declare an empty list to store interpolants by findex
+        interps_f = np.empty(self.n_findex, dtype=object)
         if z is not None:
             # Compute the 3-dimensional interpolants for each wind direction
             # Linear interpolation is used for points within the user-defined area of values,
-            # while the freestream wind speed is used for points outside that region
-            in_region = [
-                self.interpolate_multiplier_xyz(x, y, z, multiplier, fill_value=1.0)
-                for multiplier in speed_multipliers
-            ]
+            # while the freestream wind speed is used for points outside that region.
+
+            # Because the (x,y,z) points are the same for each findex, we create the triangulation
+            # once and then overwrite the values for each findex.
+
+            # Create triangulation using zeroth findex
+            interp_3d = self.interpolate_multiplier_xyz(
+                x, y, z, speed_multipliers[0], fill_value=1.0
+            )
+            # Copy the interpolant for each findex and overwrite the values
+            for findex in range(self.n_findex):
+                interp_3d.values = speed_multipliers[findex, :].reshape(-1, 1)
+                interps_f[findex] = copy.deepcopy(interp_3d)
+
         else:
             # Compute the 2-dimensional interpolants for each wind direction
             # Linear interpolation is used for points within the user-defined area of values,
             # while the freestream wind speed is used for points outside that region
-            in_region = [
-                self.interpolate_multiplier_xy(x, y, multiplier, fill_value=1.0)
-                for multiplier in speed_multipliers
-            ]
 
-        self.het_map = in_region
+            # Because the (x,y) points are the same for each findex, we create the triangulation
+            # once and then overwrite the values for each findex.
+
+            # Create triangulation using zeroth findex
+            interp_2d = self.interpolate_multiplier_xy(x, y, speed_multipliers[0], fill_value=1.0)
+            # Copy the interpolant for each findex and overwrite the values
+            for findex in range(self.n_findex):
+                interp_2d.values = speed_multipliers[findex, :].reshape(-1, 1)
+                interps_f[findex] = copy.deepcopy(interp_2d)
+
+        self.het_map = interps_f
 
     @staticmethod
     def interpolate_multiplier_xy(x: NDArrayFloat,

tests/heterogeneous_map_integration_test.py

@@ -394,3 +394,104 @@ def test_3d_het_and_shear():
                 wind_directions=[270.0], wind_speeds=[wind_speed]
             ),
         )
+
+
+def test_run_2d_het_map():
+    # Define a 2D het map and confirm the results are as expected
+    # when applied to FLORIS
+
+    # The side of the flow which is accelerated reverses for east versus west
+    het_map = HeterogeneousMap(
+        x=np.array([0.0, 0.0, 500.0, 500.0]),
+        y=np.array([0.0, 500.0, 0.0, 500.0]),
+        speed_multipliers=np.array(
+            [
+                [1.0, 2.0, 1.0, 2.0],  # Top accelerated
+                [2.0, 1.0, 2.0, 1.0],  # Bottom accelerated
+            ]
+        ),
+        wind_directions=np.array([270.0, 90.0]),
+        wind_speeds=np.array([8.0, 8.0]),
+    )
+
+    # Get the FLORIS model
+    fmodel = FlorisModel(configuration=YAML_INPUT)
+
+    from floris import TimeSeries
+
+    time_series = TimeSeries(
+        wind_directions=np.array([270.0, 90.0]),
+        wind_speeds=8.0,
+        turbulence_intensities=0.06,
+        heterogeneous_map=het_map,
+    )
+
+    # Set the model to a turbines perpinducluar to
+    # east/west flow with 0 turbine closer to bottom and
+    # turbine 1 closer to top
+    fmodel.set(
+        wind_data=time_series,
+        layout_x=[250.0, 250.0],
+        layout_y=[100.0, 400.0],
+    )
+
+    # Run the model
+    fmodel.run()
+
+    # Get the turbine powers
+    powers = fmodel.get_turbine_powers()
+
+    # Assert that in the first condition, turbine 1 has higher power
+    assert powers[0, 1] > powers[0, 0]
+
+    # Assert that in the second condition, turbine 0 has higher power
+    assert powers[1, 0] > powers[1, 1]
+
+    # Assert that the power of turbine 1 equals in the first condition
+    # equals the power of turbine 0 in the second condition
+    assert powers[0, 1] == powers[1, 0]
+
+
+def test_het_config():
+
+    # Test that setting FLORIS with a heterogeneous inflow configuration
+    # works as expected and consistent with previous results
+
+    # Get the FLORIS model
+    fmodel = FlorisModel(configuration=YAML_INPUT)
+
+    # Change the layout to a 4 turbine layout in a box
+    fmodel.set(layout_x=[0, 0, 500.0, 500.0], layout_y=[0, 500.0, 0, 500.0])
+
+    # Set FLORIS to run for a single condition
+    fmodel.set(wind_speeds=[8.0], wind_directions=[270.0], turbulence_intensities=[0.06])
+
+    # Define the speed-ups of the heterogeneous inflow, and their locations.
+    # Note that heterogeneity is only applied within the bounds of the points defined in the
+    # heterogeneous_inflow_config dictionary.  In this case, set the inflow to be 1.25x the ambient
+    # wind speed for the upper turbines at y = 500m.
+    speed_ups = [[1.0, 1.25, 1.0, 1.25]]  # Note speed-ups has dimensions of n_findex X n_points
+    x_locs = [-500.0, -500.0, 1000.0, 1000.0]
+    y_locs = [-500.0, 1000.0, -500.0, 1000.0]
+
+    # Create the configuration dictionary to be used for the heterogeneous inflow.
+    heterogeneous_inflow_config = {
+        "speed_multipliers": speed_ups,
+        "x": x_locs,
+        "y": y_locs,
+    }
+
+    # Set the heterogeneous inflow configuration
+    fmodel.set(heterogeneous_inflow_config=heterogeneous_inflow_config)
+
+    # Run the FLORIS simulation
+    fmodel.run()
+
+    turbine_powers = fmodel.get_turbine_powers() / 1000.0
+
+    # Test that the turbine powers are consistent with previous implementation
+    # 2248.2, 2800.1, 466.2, 601.5 before the change
+    # Using almost equal assert
+    assert np.allclose(
+        turbine_powers, np.array([[2248.2, 2800.0, 466.2, 601.4]]), atol=1.0,
+    )