Add capability to sample and plot velocity profiles

[vallbog]

Velocity deficit profiles

This PR adds a quick way to sample and plot velocity deficit profiles at several locations downstream of a turbine. The following definition of velocity deficit is used:
$$\frac{\Delta U}{U_\infty} = \frac{U_\infty - U}{U_\infty}$$
, where $U$ is the wake velocity obtained for a homogeneous incident wind speed of $U_\infty$ (see "Additional supporting information" for more details) .

The figure below is an output of a new example 32_plot_velocity_deficit_profiles.py. It shows velocity deficit profiles at three locations downstream of a single nrel_5MW turbine. Coordinates $x/D$, $y/D$, and $z/D$ are in the streamwise, cross-stream, and vertical direction respectively, and normalized by the turbine diameter $D$. The dashed lines show the extent of the rotor.

Use cases for this feature include:

• Comparing velocity models like in the figure above
• Comparing a velocity model to experimental data or CFD simulations

One may ask why $z$-profiles (bottom row of figure) are supported, since they are identical to the $y$-profiles in the above example. I believe they are useful if FLORIS adds support for vertical-axis wind turbines (VAWTs). VAWT wake characteristics are different in the cross-stream and vertical direction, which is why it's necessary to have both $y$- and $z$-profiles to represent the wake; see e.g. (Ouro & Lazennec, 2021). It may also be useful in layouts where the turbines have different hub heights.

Additional supporting information

When developing this feature, I chose between the following definitions of the (non-normalized) velocity deficit $\Delta U$:

1. $\Delta U = U_\infty - U$, where $U$ is the wake velocity obtained for a homogeneous incident wind speed of $U_\infty$, or
2. $\Delta U = U_{ABL}(z) - U$, where $U$ is the wake velocity obtained when the inflow has an Atmospheric Boundary Layer profile denoted $U_{ABL}(z)$.

Here, I concluded that $\Delta U$ would be the same in both cases. I chose definition 1 since it removed the complexity of the ABL. Practically, this is enforced by temporarily setting wind_shear = 0.0 when sampling the profiles.

[misi9170]

$z$-profiles may also differ from $y$-profiles with vertical wake deflections due to shaft tilt or platform tilt for floating WTs, which are supported in the Empirical Gaussian model.

[misi9170]

When does the divide by zero issue occur? Is it when self.reference_wind_height is zero? I was able to run the example to completion with and without this if statement

[vallbog]

It occurs when self.wind_shear is zero and grid.z_sorted happens to include the value zero (a sample point at the ground level). For instance, removing the ìf statement and modifing example 29 by replacing

profiles = get_profiles('y') + get_profiles('z')

with

profiles_y = get_profiles('y')
profile_range_z = [-hub_height, 3 * D]
profiles_z = fi.sample_velocity_deficit_profiles(
	direction='z',
	downstream_dists=downstream_dists,
	profile_range=profile_range_z,
	homogeneous_wind_speed=homogeneous_wind_speed,
)
profiles = profiles_y + profiles_z

gives the following warnings

/home/uname/python/floris/floris/simulation/flow_field.py:137: RuntimeWarning: divide by zero encountered in reciprocal
  * (grid.z_sorted) ** (self.wind_shear - 1)
/home/uname/python/floris/floris/simulation/flow_field.py:135: RuntimeWarning: invalid value encountered in multiply
  self.wind_shear

This is because the first element in profile_range_z makes the method include a sample point at ground level.

[misi9170]

Thanks for clarifying! I see now: the $z^{(\alpha-1)}$ evaluates to $1/0$ when $z=0$ and $\alpha=0$. This produces an inf and raises a divide by zero warning; then, the multiplication by $\alpha=0$ raises an invalid multiplication warning (0 x inf) and results in a nan in dwind_profile_plane (the derivative of the shear profile layer).

However, I actually think we should let these warnings be raised. In my testing, they only appear once, and they let the user know that something strange is happening when one evaluates the derivative of the shear profile at $z=0$ (the derivative is not well defined). In fact the divide by zero problem occurs when $z=0$ for any value of shear less than 1 (not just when shear is 0). The resulting figures appear to be the same regardless.

@rafmudaf @bayc , what are our practices regarding letting warnings be raised in "unusual" circumstances?

[rafmudaf]

Ok @vallbog @misi9170 I'm just coming up from my obsession on this question for the past couple of weeks. I started by trying to create a small unit test for this situation and realized that the Grid classes weren't able to be created from the from_dict method because they relied on the Vec3 class, so that was addressed in #751. In the process of that, I noticed that many of the classes in floris.simulation weren't proper slotted classes because the base classes weren't attrs-derived, and this exposed quite a bit of monkey patching - all of this was fixed in #750. Finally, I've been able to come back to this and now we can better test this particular situation through the unit testing infrastructure. I've written a general issue for this in #760.

In testing this further, I noticed that the if wind_shear == 0.0 condition is inadequate since the issue comes from the exponent in the derivate being less than 1 which causes the "divide by zero in reciprocal" runtime warning. As @misi9170 pointed out, this will happen if wind_shear is less than 1. However, it's actually the combination of wind shear less than 1 and z=0 that is problematic. To accommodate that, I propose to initialize the derivate to zeros and use np.power to calculate the derivative since it is a ufunc and supports the where= option:

        dwind_profile_plane = (
            self.wind_shear
            * (1 / self.reference_wind_height) ** self.wind_shear
            * np.power(
                grid.z_sorted,
                (self.wind_shear - 1),
                where=grid.z_sorted != 0.0
            )
        )

I initially used np.where, but that evaluates all operations and then picks out the ones that fit the condition so the runtime warning is still raised.

Here's a very simple script to test:

fi = FlorisInterface("inputs/gch.yaml")
x = [0.0] #, 0.0]
y = [0.0] #, 0.0]
z = [0.0] #, 10.0]
u_at_points = fi.sample_flow_at_points(x, y, z)

Without the proposed change, this raises a divide by zero warning.

[rafmudaf]

To not belabor the point, @misi9170 if you like it I'll make this commit and merge this pull request.

[misi9170]

@rafmudaf Thanks, this all makes sense. Feel free to merge. That np.where() evaluating both and therefore still raising the issue is something I've run into before, and is a pretty annoying limitation. But that's another story...

[rafmudaf]

@vallbog Well done, this is a really nice pull request.

I've read through you proposed changes, and I really like the intent and the general structure. I have some proposed improvements starting with the suggestion to use PointsGrid instead of a new grid-type. I've prototyped this in vallbog#1 just to see whether it would work. I propose this change because these two grid-types would otherwise be very similar in principle.

Is there a particular reason to have the LinesGrid class instead of the change I've proposed?

[vallbog]

@rafmudaf I like your proposed changes and don't see a need for the LinesGrid anymore. To your prototyped code, I added a rotation of the sample points so that they follow the wind direction. I illustrate this rotation in the figure below, which is drawn in a horizontal plane. The black circle is the location of the turbine hub, $U_\infty$ indicates the wind direction, the dashed line is the sample direction, and the green lines indicate where the velocity is sampled. This is how the velocity is sampled by floris.solve_for_velocity_deficit_profiles if the wind direction is 315 degrees and the input arguments are direction='y' and downstream_dists = D * np.array([1, 2]). In the original prototyped code, the sampling was done as in the right pane, while now it's done as in the left pane.

 
I believe that "with rotation" would be the most useful way of sampling the velocity, since you stay within the wake of the turbine. But please let me know your thoughts.

[rafmudaf]

Good catch. Before your last commit, the lines were all placed as if the wind was oriented with the x-axis (typically, this is West or 270 degrees) per this line:

x = (x_start + downstream_dists_transpose) * np.ones((n_lines, resolution))

Given that, then we need to rotate the lines relative to the wind direction, and that's achieved with floris.utilities.reverse_rotate_coordinates_rel_west.

[vallbog]

@rafmudaf I've made the following changes

• Clarified coordinate systems: I liked your idea of including the "raw coordinates" in the profile DataFrame such that the user can explicitly see where the velocity is sampled. I've now called these coordinates (x, y, z) and added a sampling-coordinate-system (x1, x2, x3). The latter rotates with the wind direction and the user may define its origin (where to start sampling). Each profile i is now created like this:

df = pd.DataFrame(
    {
        'x': x[i],
        'y': y[i],
        'z': z[i],
        'x1/D': x1[i]/ref_rotor_diameter,
        'x2/D': x2[i]/ref_rotor_diameter,
        'x3/D': x3[i]/ref_rotor_diameter,
        'velocity_deficit': velocity_deficit[i],
    }
)

The sampling-coordinate-system is used by floris.tools.VelocityProfilesFigure. Please see the updated example.

• Updated example: I've expanded on your idea of drawing lines on a horizontal plane to indicate where the velocity is sampled. First, there is a single-turbine case to show the basics. I chose WD = 270 such that the streamwise direction coincides with $x$ and cross-stream direction with $y$. I manually set the labels $x/D$ and $y/D$ in the VelocityProfilesFigure to simplify things. After that, there is a two-turbine case where WD = 315. This illustrates the need for the sampling-coordinate-system. The default labels $x_1/D$, $x_2/D$, and $x_3/D$ are seen in the VelocityProfilesFigure.
• Used floris.solve_for_points instead of creating a PointsGrid

[rafmudaf]

Used floris.solve_for_points instead of creating a PointsGrid

Yea this is better, good call.

[rafmudaf]

@vallbog After iterating on this a bit with you, I think it's in good shape, but let me know if you intend to make more changes. I'd like to wait to merge this until we resolve the failing test cases via #739, and I expect that to happen very soon.

[vallbog]

Sounds good! I'm happy with the current state.

[rafmudaf]

Thanks for your patience @vallbog. This was a really nice pull request and both directly and indirectly had a big effect on the architecture of FLORIS.