"hub_heights" becomes list of list when using "sample_flow_at_points"

[jfrederik-nrel]

I'm running my own version of example 28, where I want to extract the wind speed over a subset of the horizontal plane at hub height. I run a case with identical wind conditions but different turbine settings (e.g., yaw angles).

When running the function "sample_flow_at_points", I get the following error:

Traceback (most recent call last):
  File "c:\Users\jfrederi\floris\examples\43_helix_flow_in_wake.py", line 76, in <module>
    u_at_points = fi.sample_flow_at_points(points_x, points_y, points_z)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\jfrederi\floris\floris\tools\floris_interface.py", line 1020, in sample_flow_at_points
    return self.floris.solve_for_points(x, y, z)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\jfrederi\floris\floris\simulation\floris.py", line 244, in solve_for_points
    turbine_coordinates=self.farm.coordinates,
                        ^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\jfrederi\floris\floris\simulation\farm.py", line 436, in coordinates
    np.array([x, y, z]) for x, y, z in zip(self.layout_x, self.layout_y, self.hub_heights)
    ^^^^^^^^^^^^^^^^^^^
ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (3,) + inhomogeneous part.

In line 426 of farm.py, the variable self.hub_heights has the size of the number of yaw angles set in calculate_wake (but the function expects size 1 since there is only one turbine). I can circumvent this by manually setting fi.floris.farm.hub_heights = [150] before running fi.sample_flow_at_points, but that hardly seems the proper way to do things.

Here's a snippet of my code:

floris_model = "emgauss"

N = 3
fi.reinitialize(wind_directions=270 * np.ones(N), wind_speeds=7.5 * np.ones(N))
yaw_angles  = np.array([0, 10, 20]).reshape(1, N).T
fi.calculate_wake(yaw_angles = yaw_angles )

# Define plane at hub height with turbine width
x = np.linspace(-1*D, 10*D, 111)
y = np.linspace(-0.5*D, 0.5*D, 25)
points_x, points_y = np.meshgrid(x, y)

points_x = points_x.flatten()
points_y = points_y.flatten()
points_z = 150*np.ones_like(points_x)

# Collect the points
u_at_points = fi.sample_flow_at_points(points_x, points_y, points_z)

[rafmudaf]

@jfrederik-nrel what version of FLORIS is this? If you aren't on a specific version, then please post the commit hash.

[jfrederik-nrel]

@rafmudaf version 3.5, sorry for forgetting to mention that! The exact code I'm using right now can be found here, but I didn't change any functionality in FlorisInterface, so I'm confident the bug originates from the v3.5 source code.

[misi9170]

Just a bit of clarification on versions here---the below runs (from the examples directory) to completion on main branch (31fe1b6) but NOT on develop (c7f8f36):

from floris.tools import FlorisInterface
import numpy as np

fi = FlorisInterface("inputs/emgauss.yaml")
D = 126
N = 3
fi.reinitialize(
    wind_directions=270 * np.ones(N),
    wind_speeds=[7.5],
    layout_x=[0],
    layout_y=[0]
)
yaw_angles = np.zeros((1, 1, N * N))
yaw_angles  = np.array([[[0]], [[10]], [[20]]])
fi.calculate_wake(yaw_angles = yaw_angles)

# Define plane at hub height with turbine width
x = np.linspace(-1*D, 10*D, 111)
y = np.linspace(-0.5*D, 0.5*D, 25)
points_x, points_y = np.meshgrid(x, y)

points_x = points_x.flatten()
points_y = points_y.flatten()
points_z = 150*np.ones_like(points_x)

# Collect the points
u_at_points = fi.sample_flow_at_points(points_x, points_y, points_z)

emgauss.yaml, jensen.yaml, and turbopark.yaml all raise the Error mentioned; surprisingly, gch.yaml produces a different error. Edit the error raised by gch was because the yaw_angles are specified as integers rather than floats. Changing the line yaw_angles = np.array([[[0]], [[10]], [[20]]]) to yaw_angles = np.array([[[0.]], [[10.]], [[20.]]]) produces the error mentioned in the issue description, consistent with the other models.

Note that solve_for_points is not available for the cc model.

[misi9170]

I think I've found the source of the bug, but I'm not 100% sure what the right solution is. When finalize() is called on the Farm object, which is part of the chain of methods called by calculate_wake(), the following code is executed:

        self.hub_heights = np.take_along_axis(
            self.hub_heights_sorted,
            unsorted_indices[:,:,:,0,0],
            axis=2
        )

This expands self.hub_heights from a list of length n_turbines to an array of shape (n_wind_directions, n_wind_speeds, n_turbines).

Later in the code, solve_for_points() (on the Floris class in floris.py) calls self.farm.coordinates; on the main branch, this is simply an attribute that has been set earlier, before Farm.hub_heights was expanded by the finalize() call. However, on the develop branch, by #751, self.farm.coordinates is actually a getter method (with @property decorator), which means that when solve_for_points() calls self.farm.coordinates, it now takes the expanded self.hub_heights, which has inconsistent dimensions with layout_x and layout_y and therefore throws an error.

I see two solutions for this:

1. The easy fix seems to be to just remove the above code for expanding self.hub_heights during a finalize() call. I understand that self.hub_heights_sorted needs to have the full dimensions for the solve, but I don't think we need to reset self.hub_heights during finalize()---it makes more sense for hub_heights to be one dimensional, as layout_x and layout_y are. This may be invasive though if other parts of the code are using the expanded hub_heights after finalize() is called.
2. If we do in fact need self.hub_heights to be broadcast up to be of shape (n_wind_directions, n_wind_speeds, n_turbines), another option would be to add handling to the self.farm.coordinates() getter method that essentially only takes the first parts of self.hub_heights, i.e., something like

if len(self.hub_heights.shape) == 3:
    hub_heights = self.hub_heights[0,0,:]
else:
    hub_heights = self.hub_heights

Do those fixes work?

1. Works (the script in my previous comment runs, and all tests pass)
2. Works, implemented as

    @property
    def coordinates(self):
        if len(self.hub_heights.shape) == 1:
            return np.array([
                np.array([x, y, z]) for x, y, z in zip(self.layout_x, self.layout_y, self.hub_heights)
            ])
        else:
            return np.array([
                np.array([x, y, z]) for x, y, z in zip(self.layout_x, self.layout_y, self.hub_heights[0,0])
            ])

(the script runs and all tests pass)

@rafmudaf , do you have a preference between these? If implementing the first fix, I'd also consider removing code that expands the dimensions of self.rotor_diameters and possibly self.turbine_type_map, since those also seem like they should be 1-dimensional (same length as layout_x/layout_y)

[misi9170]

Hi @rafmudaf, any guidance on which of the above fixes you prefer? If you don't have a strong opinion at this stage, I can open a pull request with one of the fixes and we can discuss there.

[rafmudaf]

@misi9170 Thanks for diving into this and documenting it clearly. In my opinion, it would be simplest for all of the attributes on the Farm class to have the data type and one of three shapes:

• 1D - (turbine,)
• 3D - (wd, ws, turbine)
• 5D - (wd, ws, turbine, grid, grid)

Part of the intent of #751 was to start to align all the attributes with consistent shapes and data types, but like you said it can be pretty invasive.

My preference is to have any attribute that FLORIS creates (i.e. anything that isn't loaded directly from the input file or dictionary) be a np.array. As for the shape of hub heights, I think for consistency it should be used exactly like Farm.rotor_diameters since they are both used similarly throughout FLORIS and they relate in a similar way to the application of the software - a physical characteristic of the wind turbine that can be used as a reference length. Ultimately, I think your option 1 is the "right" fix, but the scope of the change should be assessed especially relative to the current work in v4. In other words, if it's very invasive, maybe do option 2 until after we release v4, and then we can do option 1.

[misi9170]

Thanks @rafmudaf ---I'll take a quick pass at option 1 to see how invasive it is; if it ends up being complicated I'll revert to option 2 for now, as you suggest, and open an issue to reconsider option 1 for v4.