Add coordinates from inertial frame to grid class

[Bartdoekemeijer]

The pull request is ready to be merged, but requires review.

Feature or improvement description
This pull request adds the inertial locations of the flow fields/turbine locations in the grid class. With inertial, I mean the location of the turbines without rotating it to the incoming wind direction. These additional variables should also be helpful when working with heterogeneous inflow conditions. The way it seems to be implemented now, the heterogeneous inflow conditions are linearly interpolated using the rotated (wind-aligned) coordinates, meaning their physical location differs for different wind directions. I'm opening a separate issue for that. Anyway, this PR is also particularly helpful when we want to plot horizontal flow slices when fixing the plot to true north and with the turbines plotted in same location for every wind direction -- i.e., the same issue as described in #556.

I have added a couple of handles to the visualization such that it can be used for plotting horizontal flow slices while fixing the orientation to true north. Note that nothing fundamentally has changed to the FLORIS simulation code, besides the exporting/calculation of a couple extra variables, so it's an easy merge.

Related issue, if one exists
@MYMahfouz provides a nice overview in #556.

Impacted areas of the software
The FLORIS grid and visualization classes.

Additional supporting information
These additional variables should also be helpful when working with heterogeneous inflow conditions. The way it seems to be implemented now, the heterogeneous inflow conditions are linearly interpolated using the rotated (wind-aligned) coordinates, meaning their physical location differs for different wind directions. I'm opening a separate issue for that.

Test results, if applicable
A simple test case is:

import matplotlib.pyplot as plt
import numpy as np

import floris.tools.visualization as wakeviz
from floris.tools import FlorisInterface
from floris.tools.visualization import plot_turbines_with_fi

fi = FlorisInterface("inputs/gch.yaml")
fi.reinitialize(wind_directions=[51.0])

yaw_angles = np.array([[[25.,0.,0.]]])
rotate_to_inertial_frame = True
horizontal_plane = fi.calculate_horizontal_plane(
    x_resolution=200,
    y_resolution=100,
    height=90.0,
    rotate_to_inertial_frame=rotate_to_inertial_frame,
    yaw_angles=yaw_angles,
)

fig, ax = plt.subplots()
wakeviz.visualize_cut_plane(horizontal_plane, ax=ax, title="Horizontal")
plot_turbines_with_fi(
    fi,
    ax=plt.gca(),
    color="black",
    rotate_to_inertial_frame=rotate_to_inertial_frame,
    yaw_angles=yaw_angles,
)

plt.show()

which produces

[Bartdoekemeijer]

Also see my comments in #579.

[Bartdoekemeijer]

Just merged in latest develop branch.

[bayc]

@Bartdoekemeijer @rafmudaf @paulf81

Hey all. I've gone through the PR and made a variety of updates. Thanks for all your work on this Bart! I think it is a really nice addition.

The major changes are that the default plotting behavior is now to display the flow field as rotated to the defined wind direction, instead of the 270 degrees used previously.

Also, I addressed the issue identified by Bart in #579 by making sure the calculate_speed_ups function is passed the inertial coordinates (captured here e48237b).

Through this process, I re-worked the heterogeneous workflow so that the user defines an het_config dictionary and that information is now stored on FlorisInterface. I believe this helps to streamline the interface and not require the user to import a separate function to create the interpolant and pass that in, but leverage this optional parameter. I am curious on y'alls feedback regarding this.

I also made changes to how the heterogeneous flows get plotted, such that the user-defined bounds for the het map are displayed by default on any plot that shows heterogenous flows. Additionally, I have removed the nearest-neighbor interpolation that previously occurred outside the user-defined heterogeneous inflow area and used a fill value of 1.0 for the speed ups such that this region is now set to be equal to freestream. This removed a good chunk of code and simplified the whole calculation.

I have put together an example of this, code and plots below, that demonstrates the new code behavior for 4 different wind directions. Note the option to remove the het map bounding box in the bottom right plot. Please check and provide feedback.

from floris.tools import FlorisInterface
import matplotlib.pyplot as plt
import numpy as np
import floris.tools.visualization as wakeviz
from floris.tools.visualization import plot_turbines_with_fi

speed_ups = [[1.2, 1.2, 1.0, 1.0]]
x_locs = [-300.0, -300.0, 2600.0, 2600.0]
y_locs = [ -300.0, 300.0, -300.0, 300.0]

# Generate the configuration dictionary to be used for the heterogeneous inflow.
het_config = {
    'speed_ups': speed_ups,
    'x_locs': x_locs,
    'y_locs': y_locs,
}

# Initialize FLORIS with a heterogeneous inflow
fi_2d = FlorisInterface("inputs/gch.yaml", het_config=het_config)


# Plot the farm for the wind direction of 270 degrees
fig, ax = plt.subplots(2, 2)
yaw_angles = np.array([[[25.,0.,0.]]])
horizontal_plane_270 = fi_2d.calculate_horizontal_plane(
    wd=[270.],
    x_resolution=200,
    y_resolution=100,
    height=90.0,
    yaw_angles=yaw_angles,
)
wakeviz.visualize_heterogeneous_cut_plane(horizontal_plane_270, fi_2d, ax=ax[0][0], title="Horizontal")
plot_turbines_with_fi(
    fi_2d,
    ax=ax[0][0],
    color="black",
    wd=[270.],
    yaw_angles=yaw_angles,
)


# Plot the farm for the wind direction of 275 degrees
horizontal_plane_275 = fi_2d.calculate_horizontal_plane(
    wd=[275.],
    x_resolution=200,
    y_resolution=100,
    height=90.0,
    yaw_angles=yaw_angles,
)
wakeviz.visualize_heterogeneous_cut_plane(horizontal_plane_275, fi_2d, ax=ax[0][1], title="Horizontal")
plot_turbines_with_fi(
    fi_2d,
    ax=ax[0][1],
    color="black",
    wd=[275.],
    yaw_angles=yaw_angles,
)


# Plot the farm for the wind direction of 315 degrees
horizontal_plane_315 = fi_2d.calculate_horizontal_plane(
    wd=[315.],
    x_resolution=200,
    y_resolution=100,
    height=90.0,
    yaw_angles=yaw_angles,
)
wakeviz.visualize_heterogeneous_cut_plane(horizontal_plane_315, fi_2d, ax=ax[1][0], title="Horizontal")
plot_turbines_with_fi(
    fi_2d,
    ax=ax[1][0],
    color="black",
    wd=[315.],
    yaw_angles=yaw_angles,
)


# Plot the farm for the wind direction of 360 degrees
horizontal_plane_360 = fi_2d.calculate_horizontal_plane(
    wd=[360.],
    x_resolution=200,
    y_resolution=100,
    height=90.0,
    yaw_angles=yaw_angles,
)
wakeviz.visualize_heterogeneous_cut_plane(horizontal_plane_360, fi_2d, ax=ax[1][1], title="Horizontal", plot_het_bounds=False)
plot_turbines_with_fi(
    fi_2d,
    ax=ax[1][1],
    color="black",
    wd=[360.],
    yaw_angles=yaw_angles,
)

# Show the plots
plt.show()

[Bartdoekemeijer]

Hmm, interesting. Are you thinking that the inflow from upstream wind farms would cause different speed up multipliers in the downstream farm, depending on what the freestream, and thus wake deficits, are? I agree that the inflow would be different, but am unsure if it would cause the speed up multipliers to change over a range of wind speeds.

That's right. The higher the ambient wind speed, the lower the thrust coefficient of neighboring farms, and ct is a large driver in the depth of wakes and blockage effects. To demonstrate this, here's a simple example showing the effect of a neighboring wind farm's wakes (11 kilometers away) on the inflow wind speed of 5 turbines, for a wind direction of 270 deg:

import matplotlib.pyplot as plt
import numpy as np

from floris.tools import FlorisInterface


def load_floris():
    # Load the default example floris object
    fi = FlorisInterface("inputs/turbopark.yaml")

    # Specify the full wind farm layout: nominal and neighboring wind farms
    X, Y = np.meshgrid(np.arange(0, 5) * 800.0, np.arange(0, 5) * 800.0)
    X = X.flatten()
    Y = Y.flatten()

    distance = 15.0e3
    Y = np.hstack([Y[X==0], Y + 150 * np.random.rand(*np.shape(X))])
    X = np.hstack([X[X==0], X - distance + 300 * np.random.rand(*np.shape(X))])

    # Turbine weights: we want to only optimize for the first 10 turbines
    turbine_weights = np.zeros(len(X), dtype=int)
    turbine_weights[0:5] = 1.0

    # Now reinitialize FLORIS layout
    fi.reinitialize(layout_x=X, layout_y=Y)

    # And visualize the floris layout
    fig, ax = plt.subplots()
    ax.plot(X[turbine_weights == 0], Y[turbine_weights == 0], 'ro', label="Neighboring farms")
    ax.plot(X[turbine_weights == 1], Y[turbine_weights == 1], 'go', label='Farm subset')
    ax.grid(True)
    ax.set_xlabel("x coordinate (m)")
    ax.set_ylabel("y coordinate (m)")
    ax.axis("equal")
    ax.legend()

    return fi, turbine_weights


if __name__ == "__main__":
    fi, turbine_weights = load_floris()
    fig, ax = plt.subplots()
    for u_inf in range(6, 20):
        fi.reinitialize(
            wind_speeds=[u_inf],
            turbulence_intensity=0.05,
            wind_directions=[270.0],
        )
        fi.calculate_wake()
        
        speedup_map = fi.floris.flow_field.u.mean(axis=3)[0, 0, 0:5, 0] / u_inf
        ax.plot(range(len(speedup_map)), speedup_map, "-o", label=f"Ambient wind speed: {u_inf} m/s")

    ax.set_xlabel("Turbine number (-)")
    ax.set_xticks(range(5))
    ax.set_ylabel("Heterogeneous inflow map (-)")
    ax.grid(True)
    ax.legend()

    plt.show()

This confirms that, the higher the wind speed, the lower ct and thus the lower neighboring farm effects. I expect the same correlation for blockage effects. Note that these effects are generally much larger with TurbOPark than with GCH or CC, but the same principles carry over (just to a smaller effect).

[bayc]

Wow, thanks for sharing that @Bartdoekemeijer! That is really interesting. I believe the internal consensus over here (@paulf81 and @rafmudaf can confirm) is to close out this PR where it is, and open up a new issue to tackle adding the heterogeneous input dependency on wind speed.

[rafmudaf]

I think het_config is yet obscure if you don't already know that this is a thing. Instead, I propose to change the name to something like speed_ups_map. @bayc @misi9170 @paulf81 what do you think?

[paulf81]

Agree that this is a more descriptive name and will more clearly indicate its purpose, I say switch it

[misi9170]

My only reservation is that we use the term "heterogeneous" pretty regularly throughout the code, so I kind of like having some indication that this is related to heterogeneity in the name. What do you think about heterogeneous_speed_ups_map? A little long and clunky...

[misi9170]

Looking into this a bit more, I'd update my recommended name to heterogeneous_inflow_config. In my mind, the dict is more of a data structure than a mapping of key-value pairs. However, I think that the "speed_ups" key should instead be named "speed_multipliers", to be more explicit about how the speed-up is applied (as a multiplier, as opposed to an addend).

[rafmudaf]

@misi9170 I like it. If we support heterogeneous wind directions, that will logically fit into heterogenous_inflow_config also. I think x/y_loc should be simply x/y since the _loc doesn't add much clarity. It's reasonably obvious that these are the locations of the speed multipliers.

Here's where we're at:

heterogenous_inflow_config = {
    "speed_multipliers": [1.0, 2.0, 3.0],
    "x": [ -300.0, 300.0, -300.0],
    "y": [ -300.0, 300.0, -300.0],
    "z": [     90.0,   90.0,     90.0]
}

[rafmudaf]

This may be incorrect given the line in plot_turbines function call below that also extracts the first wind direction. Double check @bayc

[rafmudaf]

Oh I see now that this is only used to adjust the yaw angles. That's probably good to reflect in the doc string

[rafmudaf]

These are added and x etc are deleted. As I understand it, x/y/z_sorted_intertial_frame are the same as x/y/z except that they're sorted. Will we need to retain x/y/z for postprocessing? For example, if you need to plot the layout with ID's based on your input, do we require to unsort x/y/z_sorted_intertial_frame?
2.

[rafmudaf]

A follow on question is whether this should be simply x/y/z with those attributes being defined (in words) the same as you would define the x/y/z_sorted_intertial_frame attributes.

[rafmudaf]

Given the added complexity in the grids, this is a good time to document the _sorted and _sorted_inertial_frame variables. The minimum documentation would here in the docstring, and better would be to also add to the Floris_101 notebook in the docs with detailed exaplanation and code or analysis of these variables.

[rafmudaf]

Whats the shape of yaw_angles here? Previously, it was assumed that it was 3D (wind direction, wind speed, N turbines). Here, it looks like it assumed this is 1D (N turbines). Same question for fi.floris.farm.rotor_diameters below.

[rafmudaf]

Overall, this looks like it's in pretty good shape!

[rafmudaf]

For reference, I've updated @bayc's example from #578 (comment) for the new API in this pull request. The script is below and the plots that it creates are below that.

from floris.tools import FlorisInterface
import matplotlib.pyplot as plt
import numpy as np
import floris.tools.visualization as wakeviz
from floris.tools.visualization import plot_turbines_with_fi

speed_ups = [[1.2, 1.2, 1.0, 1.0]]
x_locs = [-300.0, -300.0, 2600.0, 2600.0]
y_locs = [ -300.0, 300.0, -300.0, 300.0]

# Generate the configuration dictionary to be used for the heterogeneous inflow.
het_config = {
    'speed_ups': speed_ups,
    'x_locs': x_locs,
    'y_locs': y_locs,
}

# Initialize FLORIS with a heterogeneous inflow
fi_2d = FlorisInterface("inputs/gch_heterogeneous_inflow.yaml")


# Plot the farm for the wind direction of 270 degrees
fig, ax = plt.subplots(2, 2)
yaw_angles = np.array([[[25.,0.,0.]]])
horizontal_plane_270 = fi_2d.calculate_horizontal_plane(
    wd=[270.],
    x_resolution=200,
    y_resolution=100,
    height=90.0,
    yaw_angles=yaw_angles,
)
wakeviz.visualize_heterogeneous_cut_plane(horizontal_plane_270, fi_2d, ax=ax[0][0], title="Horizontal")
plot_turbines_with_fi(
    fi_2d,
    ax=ax[0][0],
    color="black",
    wd=[270.],
    yaw_angles=yaw_angles,
)


# Plot the farm for the wind direction of 275 degrees
horizontal_plane_275 = fi_2d.calculate_horizontal_plane(
    wd=[275.],
    x_resolution=200,
    y_resolution=100,
    height=90.0,
    yaw_angles=yaw_angles,
)
wakeviz.visualize_heterogeneous_cut_plane(horizontal_plane_275, fi_2d, ax=ax[0][1], title="Horizontal")
plot_turbines_with_fi(
    fi_2d,
    ax=ax[0][1],
    color="black",
    wd=[275.],
    yaw_angles=yaw_angles,
)


# Plot the farm for the wind direction of 315 degrees
horizontal_plane_315 = fi_2d.calculate_horizontal_plane(
    wd=[315.],
    x_resolution=200,
    y_resolution=100,
    height=90.0,
    yaw_angles=yaw_angles,
)
wakeviz.visualize_heterogeneous_cut_plane(horizontal_plane_315, fi_2d, ax=ax[1][0], title="Horizontal")
plot_turbines_with_fi(
    fi_2d,
    ax=ax[1][0],
    color="black",
    wd=[315.],
    yaw_angles=yaw_angles,
)


# Plot the farm for the wind direction of 360 degrees
horizontal_plane_360 = fi_2d.calculate_horizontal_plane(
    wd=[360.],
    x_resolution=200,
    y_resolution=100,
    height=90.0,
    yaw_angles=yaw_angles,
)
wakeviz.visualize_heterogeneous_cut_plane(horizontal_plane_360, fi_2d, ax=ax[1][1], title="Horizontal", plot_het_bounds=False)
plot_turbines_with_fi(
    fi_2d,
    ax=ax[1][1],
    color="black",
    wd=[360.],
    yaw_angles=yaw_angles,
)

# Show the plots
plt.show()