Reduce computation time massively in large het_map objects

[Bartdoekemeijer]

Reduce computation time for large het_map objects

Currently, a new LinearNDInterpolant is prepared for each findex in a FLORIS timeseries evaluation with heterogeneous_map. The preparation of a LinearNDInterpolant required for the heterogeneous map is computationally intensive (especially when the het_map is defined for many coordinates) due to the Delaunay triangulation. However, this triangulation is identical between each findex, and therefore it makes sense to recycle this information rather than to recalculate it for each findex.

Related issue

I haven't made a separate issue for this. I figured I'd open a PR directly.

Impacted areas of the software

The flow_field.py module.

Additional supporting information

In my usage, it was taking about 45 seconds to load the heterogeneous map interpolants. This is really wasted time and was brought down to 0.4 seconds with this PR by recycling the object as conserving as much information as possible between the findices.

Test results, if applicable

Here's a test script to benchmark this functionality:

import numpy as np
import pandas as pd
from time import perf_counter as timerpc

from floris import (
    FlorisModel,
    TimeSeries,
    HeterogeneousMap
)


if __name__ == "__main__":
    # Create big grid of wind conditions and wind speeds for which we assume to have evaluated het_map
    wd_grid, ws_grid = np.meshgrid(
        np.arange(0.0, 360.0, 3.0),
        np.arange(0.5, 30.51, 1.0)
    )
    df = pd.DataFrame({"wd": wd_grid.flatten(), "ws": ws_grid.flatten()})
    print(f"We have {df.shape[0]} findices.")

    # Create a grid of sensors throughout the farm in x, y, and z
    xg, yg, zg = np.meshgrid(
        np.linspace(-3000.0, 3000.0, 11),
        np.linspace(-3000.0, 3000.0, 11),
        np.arange(0.0, 350.01, 25.0),
    )
    xg = xg.flatten()
    yg = yg.flatten()
    zg = zg.flatten()
    speedups = np.ones((df.shape[0], len(xg)))
    print(f"We have {len(xg)} number of coordinates with het_map information.")

    # Now create FLORIS and a timeseries object with het_map information
    fmodel = FlorisModel("inputs/gch.yaml")
    fmodel.set(wind_shear=0.0)  # Required when working with 3D het_map objects
    het_map = HeterogeneousMap(
        x=xg,
        y=yg,
        z=zg,
        speed_multipliers=speedups,
        wind_directions=wd_grid.flatten(),
        wind_speeds=ws_grid.flatten(),
    )

    print(f"Preparing a timeseries object for 360 findex conditions.")
    ts = TimeSeries(
        wind_directions=np.arange(0.0, 360.0, 1.0),
        wind_speeds=120.0 * np.ones(360),
        turbulence_intensities=0.06 * np.ones(360),
        heterogeneous_map=het_map,
    )
    t0 = timerpc()
    fmodel.set(wind_data=ts)
    print(f"Time spent in 'fmodel.set': {timerpc() - t0:.2f} s")

With the new PR, this takes 0.4 seconds on my system. With the old code, it takes 25 seconds. If you increase the number of findices, the old code scales the computation time linearly. In the new code, there is pretty no penalty for additional findices.

[paulf81]

hi @Bartdoekemeijer , thank you for this! I made some small formatting changes and will take a deeper dive next week

[paulf81]

Wouldn't this need to be assigned to F as it is above @Bartdoekemeijer
?

[Bartdoekemeijer]

The reason for this is that F.values is of shape (N, 1), rather than shape (N). If I don't maintain that shape, the code returns an error.

[paulf81]

I think still if you don't assign to F in line 305 you can't reference it in line 308 (it doesn't exist). Similarly I think multiplier doesn't exist yet at line 305 so followed the example of the earlier case and used speed_multipliers[0]

[paulf81]

@Bartdoekemeijer thanks again for this! I added a new test to check this was working as expected (the logic of the change is clear to me!). The test showed I think a need for a few small tweaks, if you could review my changes and let me know match your expectations

[paulf81]

Also added timing tests of het to #1060 so hopefully we can clock the improvement there as well

[paulf81]

Noting though there are some failing examples, probably this needs we still need another fix, and then also we shouldn't have all tests passing if examples are failing so we'll want a test that captures whatever failure mode is in the example

[Bartdoekemeijer]

Good catches! I haven't reran this, but the principle of this code is that you don't need to re-do the Delaunay triangulation for every interpolant. They are all identical, and you only need to swap out the interpolant values rather than the interpolant grid. Sorry if this wasn't clear. In principle this should be a tiny modification to the code.

[paulf81]

ok @Bartdoekemeijer , I think I tracked down the issues. Since the failures were happening in examples and not in tests, I added some new tests that failed until I made the correction, as that feels like best practice. The issue came down to not working if lists were passed in for the multipliers (rather than arrays as is typical more recently). @misi9170 , would you mind to review now, I think it's about ready

[misi9170]

@Bartdoekemeijer Thanks for this, and @paulf81 thank you for adding tests!

I've now approved. Since it wasn't initially clear to me what exactly was happening, I've now added a couple of comments to the code block to describe why we are replacing the values for each findex. I've also renamed variables to make it clearer what they are (in particular, the list of interpolants, previously named in_region, which sounded to me like a mask), as well as changing the loop to be explicitly over the findices.

I've also run all examples in examples/examples_heterogeneous/ and see no visual change to outputs.

Thanks, as always, for sharing the running script in the description @Bartdoekemeijer . That helps a lot.

Provided you are both happy with the changes I've made, I'll get this merged in.

[paulf81]

Very happy, please merge, and then we can check https://nrel.github.io/floris/dev/bench/