Add cubature integration method for rotor velocities

[rafmudaf]

Add a cubature integration grid to average rotor points

This pull request adds a new grid type for placing and weighting grid points at the rotor so that a more accurate average can be calculated with fewer grid points. As @Bartdoekemeijer noted in #499, the square grid of points on the rotor swept area is inadequate since the calculated performance is highly dependent on grid resolution. Furthermore, the outer region of the blades are not captured at all and this is the area most sensitive to power production or power loss due to wakes.

The plot below shows a comparison of the convergence between the new and old method for averaging the velocity on the rotor for a turbine under no wake and with a shear profile exponent of 0.12. The x-axis for both series are the grid points input, and this translates to the square root of the number of points actually created on the grid. Here, you can see the square grid has a slow convergence rate, while the cubature grid is nearly converged by N=3.

This slow convergence in the current method is due to coarse sampling of both the shear profile and the wind turbine wake profiles. As shown in the image below, the typical 3 x 3 rectangular grid on the rotor swept area only covers a portion of the discs, and the new cubature grid covers a larger portion of the swept area.

The following plot again shows the convergence of the cubature method for three turbines in line so two downstream turbines are waked.

Usage

The cubature integration method works seamlessly with the existing solvers. To enable it, simply use these solver settings:

    solver_dict = {
        "type": "turbine_cubature_grid",
        "turbine_grid_points": 4
    }

API Changes

Initially, this pull request changed the API's for the turbine-module functions for calculating performance (Turbine.rotor_effective_velocity, Turbine.Ct, Turbine.axial_induction). However, after reviewing the current implementation and results, we've decided to approach this change as an intermediate step for including this type of grid and integration method. The current API's are retained so that no breaking changes are introduced when using the existing square grid (TurbineGrid). To use the new grid (TurbineCubatureGrid), additional input arguments are passed to the following functions:

def rotor_effective_velocity(
    ...
    average_method: str = "cubic-mean",
    cubature_weights: np.array | None = None
) -> NDArrayFloat:

def Ct(
    ...
    average_method: str = "cubic-mean",
    cubature_weights: np.array | None = None
) -> NDArrayFloat:

def axial_induction(
    ...
    average_method: str = "cubic-mean",
    cubature_weights: np.array | None = None
) -> NDArrayFloat:

def average_velocity(
    ...
    method: str = "cubic-mean",
    cubature_weights: np.array | None = None
) -> NDArrayFloat:

After further validation, it is intended that this grid-type will be the default, and the API's will be changed accordingly.

Future Work

The changes here have raised some important questions about the current velocity averaging in FLORIS. Specifically, we currently average velocities with a cube root of the mean of the velocities cubed. This accurately captures the average velocity for the power calculation since that is a function of the velocity cubed, but other turbine properties have a different relationship to the velocity yet we use the same method.

Additionally, it may be worthwhile to consider how the velocities are weighted across the blade span since the root does less work than the tip, but the velocities at both stations are weighted equally.

Finally, more work is required to show that the cubature method arrives at the correct answer in addition to arriving at the answer faster. The existing model defaults may need to be retuned based on the performance calculated from a different velocity integration, as well.

Related issue

#499

Impacted areas of the software

A new grid-type is added to the grid module, the turbine module performance functions have new optional arguments, and the solvers are updated to reflect these API changes.

[misi9170]

This looks good to me @rafmudaf. I'll look now into why the empirical_gauss examples aren't running and push changes up.

[misi9170]

I can see why this is needed here; I wonder whether in future, we should have an abstract TurbineGrid class (inherits from abstract Grid, or possibly general PointsGrid as we have discussed elsewhere) from which TurbineCubatureGrid and TurbineSquareGrid would inherit. Could get a little complicated with multiple levels on inheritance though.

[rafmudaf]

I agree. The grid module overall needs some work. It would be helpful to identify where there's overlap between these grid-types. Possibly one path forward is to get rid of one grid super class in favor of categories of grids. At the least, the grid abstract base class could only define the required API's for all grids like x, x_sorted and sorted_indeces etc.

[paulf81]

Reviewed the code, and I know a few more changes are coming but wanted to log approval that this is looking really great

[Bartdoekemeijer]

Specifically, we currently average velocities with a cube root of the mean of the velocities cubed. This accurately captures the average velocity for the power calculation since that is a function of the velocity cubed,

I'm not sure if this statement holds up, but agree that it's good to push this back to the next FLORIS release!

Thanks for doing all this work. This is going to be very helpful moving forward!

[rafmudaf]

Thanks for catching that @Bartdoekemeijer. I've made a note to clarify this in a couple of weeks, so I'll reach out then.