Add capability to sample and plot velocity profiles (#699)

* Moved from other branch

* Preparing inputs in FlorisInterface

* Created VelocityProfileGrid

* Velocity profiles returned successfully

* Small fixes

* Clarifications

* Draft of init for VelocityProfilesFigure class

* Small addition to class init

* Methods for adding profiles and reference lines

* Improved example

* Corrected self.axs

* Comments, style check, and updated output coords

* Details in example

* Updated comments

* Change example number and add to docs listing

* Replace LinesGrid with PointsGrid

Uses PointsGrid with lines rather than distinct points

* Rotate sample points with the wind direction, and cleanup

* Remove 5-dimensional indexing in point rotation

* Return figure axes from CutPlane vis

* Plot lines in horizontal plane

* Change coordinate notation and begin to update example

* Correct sampling, use solve_for_points, and WIP example

* WIP: example

* WIP: example

* Finalize example

* Details

* Auto-set xlim

* Fix tabbing

* Improve handling the velocity derivate at z=0

---------

Co-authored-by: Rafael M Mudafort <rafmudaf@gmail.com>

docs/examples.md

@@ -146,6 +146,14 @@ mast across all wind directions (at a fixed free stream wind speed).
 Try different values for met_mast_option to vary the location of the met mast within
 the two-turbine farm.
 
+### 32_plot_velocity_deficit_profiles.py
+This example illustrates how to plot velocity deficit profiles at several locations
+downstream of a turbine. Here we use the following definition:
+
+    velocity_deficit = (homogeneous_wind_speed - u) / homogeneous_wind_speed
+        , where u is the wake velocity obtained when the incoming wind speed is the
+        same at all heights and equal to `homogeneous_wind_speed`.
+
 ### 29_floating_vs_fixedbottom_farm.py
 
 Compares a fixed-bottom wind farm (with a gridded layout) to a floating

examples/32_plot_velocity_deficit_profiles.py

@@ -0,0 +1,205 @@
+# Copyright 2021 NREL
+
+# Licensed under the Apache License, Version 2.0 (the "License"); you may not
+# use this file except in compliance with the License. You may obtain a copy of
+# the License at http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
+# License for the specific language governing permissions and limitations under
+# the License.
+
+# See https://floris.readthedocs.io for documentation
+
+
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib import ticker
+
+import floris.tools.visualization as wakeviz
+from floris.tools import cut_plane, FlorisInterface
+from floris.tools.visualization import VelocityProfilesFigure
+from floris.utilities import reverse_rotate_coordinates_rel_west
+
+
+"""
+This example illustrates how to plot velocity deficit profiles at several locations
+downstream of a turbine. Here we use the following definition:
+    velocity_deficit = (homogeneous_wind_speed - u) / homogeneous_wind_speed
+        , where u is the wake velocity obtained when the incoming wind speed is the
+        same at all heights and equal to `homogeneous_wind_speed`.
+"""
+
+# The first two functions are just used to plot the coordinate system in which the
+# profiles are sampled. Please go to the main function to begin the example.
+def plot_coordinate_system(x_origin, y_origin, wind_direction):
+    quiver_length = 1.4 * D
+    plt.quiver(
+        [x_origin, x_origin],
+        [y_origin, y_origin],
+        [quiver_length, quiver_length],
+        [0, 0],
+        angles=[270 - wind_direction, 360 - wind_direction],
+        scale_units='x',
+        scale=1,
+    )
+    annotate_coordinate_system(x_origin, y_origin, quiver_length)
+
+def annotate_coordinate_system(x_origin, y_origin, quiver_length):
+    x1 = np.array([quiver_length + 0.35 * D, 0.0])
+    x2 = np.array([0.0, quiver_length + 0.35 * D])
+    x3 = np.array([90.0, 90.0])
+    x, y, _ = reverse_rotate_coordinates_rel_west(
+            fi.floris.flow_field.wind_directions,
+            x1[None, :],
+            x2[None, :],
+            x3[None, :],
+            x_center_of_rotation=0.0,
+            y_center_of_rotation=0.0,
+    )
+    x = np.squeeze(x, axis=0) + x_origin
+    y = np.squeeze(y, axis=0) + y_origin
+    plt.text(x[0], y[0], '$x_1$', bbox={'facecolor': 'white'})
+    plt.text(x[1], y[1], '$x_2$', bbox={'facecolor': 'white'})
+
+if __name__ == '__main__':
+    D = 126.0 # Turbine diameter
+    hub_height = 90.0
+    homogeneous_wind_speed = 8.0
+
+    fi = FlorisInterface("inputs/gch.yaml")
+    fi.reinitialize(layout_x=[0.0], layout_y=[0.0])
+
+    # ------------------------------ Single-turbine layout ------------------------------
+    # We first show how to sample and plot velocity deficit profiles on a single-turbine layout.
+    # Lines are drawn on a horizontal plane to indicate were the velocity is sampled.
+    downstream_dists = D * np.array([3, 5, 7])
+    # Sample three profiles along three corresponding lines that are all parallel to the y-axis
+    # (cross-stream direction). The streamwise location of each line is given in `downstream_dists`.
+    profiles = fi.sample_velocity_deficit_profiles(
+        direction='cross-stream',
+        downstream_dists=downstream_dists,
+        homogeneous_wind_speed=homogeneous_wind_speed,
+    )
+
+    horizontal_plane = fi.calculate_horizontal_plane(height=hub_height)
+    fig, ax = plt.subplots(figsize=(6.4, 3))
+    wakeviz.visualize_cut_plane(horizontal_plane, ax)
+    colors = ['b', 'g', 'c']
+    for i, profile in enumerate(profiles):
+        # Plot profile coordinates on the horizontal plane
+        ax.plot(profile['x'], profile['y'], colors[i], label=f'x/D={downstream_dists[i] / D:.1f}')
+    ax.set_xlabel('x [m]')
+    ax.set_ylabel('y [m]')
+    ax.set_title('Streamwise velocity in a horizontal plane: gauss velocity model')
+    fig.tight_layout(rect=[0, 0, 0.82, 1])
+    ax.legend(bbox_to_anchor=[1.29, 1.04])
+
+    # Initialize a VelocityProfilesFigure. The workflow is similar to a matplotlib Figure:
+    # Initialize it, plot data, and then customize it further if needed.
+    profiles_fig = VelocityProfilesFigure(
+        downstream_dists_D=downstream_dists / D,
+        layout=['cross-stream'],
+        coordinate_labels=['x/D', 'y/D'],
+    )
+    # Add profiles to the VelocityProfilesFigure. This method automatically matches the supplied
+    # profiles to the initialized axes in the figure.
+    profiles_fig.add_profiles(profiles, color='k')
+
+    # Change velocity model to jensen, get the velocity deficit profiles,
+    # and add them to the figure.
+    floris_dict = fi.floris.as_dict()
+    floris_dict['wake']['model_strings']['velocity_model'] = 'jensen'
+    fi = FlorisInterface(floris_dict)
+    profiles = fi.sample_velocity_deficit_profiles(
+        direction='cross-stream',
+        downstream_dists=downstream_dists,
+        homogeneous_wind_speed=homogeneous_wind_speed,
+        resolution=400,
+    )
+    profiles_fig.add_profiles(profiles, color='r')
+
+    # The dashed reference lines show the extent of the rotor
+    profiles_fig.add_ref_lines_x2([-0.5, 0.5])
+    for ax in profiles_fig.axs[0]:
+        ax.xaxis.set_major_locator(ticker.MultipleLocator(0.2))
+
+    profiles_fig.axs[0,0].legend(['gauss', 'jensen'], fontsize=11)
+    profiles_fig.fig.suptitle(
+        'Velocity deficit profiles from different velocity models',
+        fontsize=14,
+    )
+
+    # -------------------------------- Two-turbine layout --------------------------------
+    # This is a two-turbine case where the wind direction is north-west. Velocity profiles
+    # are sampled behind the second turbine. This illustrates the need for a
+    # sampling-coordinate-system (x1, x2, x3) that is rotated such that x1 is always in the
+    # streamwise direction. The user may define the origin of this coordinate system
+    # (i.e. where to start sampling the profiles).
+    wind_direction = 315.0 # Try to change this
+    downstream_dists = D * np.array([3, 5])
+    floris_dict = fi.floris.as_dict()
+    floris_dict['wake']['model_strings']['velocity_model'] = 'gauss'
+    fi = FlorisInterface(floris_dict)
+    # Let (x_t1, y_t1) be the location of the second turbine
+    x_t1 =  2 * D
+    y_t1 = -2 * D
+    fi.reinitialize(wind_directions=[wind_direction], layout_x=[0.0, x_t1], layout_y=[0.0, y_t1])
+
+    # Extract profiles at a set of downstream distances from the starting point (x_start, y_start)
+    cross_profiles = fi.sample_velocity_deficit_profiles(
+        direction='cross-stream',
+        downstream_dists=downstream_dists,
+        homogeneous_wind_speed=homogeneous_wind_speed,
+        x_start=x_t1,
+        y_start=y_t1,
+    )
+
+    horizontal_plane = fi.calculate_horizontal_plane(height=hub_height, x_bounds=[-2 * D, 9 * D])
+    ax = wakeviz.visualize_cut_plane(horizontal_plane)
+    colors = ['b', 'g', 'c']
+    for i, profile in enumerate(cross_profiles):
+        ax.plot(
+            profile['x'],
+            profile['y'],
+            colors[i],
+            label=f'$x_1/D={downstream_dists[i] / D:.1f}$',
+        )
+    ax.set_xlabel('x [m]')
+    ax.set_ylabel('y [m]')
+    ax.set_title('Streamwise velocity in a horizontal plane')
+    ax.legend()
+    plot_coordinate_system(x_origin=x_t1, y_origin=y_t1, wind_direction=wind_direction)
+
+    # Sample velocity deficit profiles in the vertical direction at the same downstream
+    # locations as before. We stay directly downstream of the turbine (i.e. x2 = 0). These
+    # profiles are almost identical to the cross-stream profiles. However, we now explicitly
+    # set the profile range. The default range is [-2 * D, 2 * D].
+    vertical_profiles = fi.sample_velocity_deficit_profiles(
+        direction='vertical',
+        profile_range=[-1.5 * D, 1.5 * D],
+        downstream_dists=downstream_dists,
+        homogeneous_wind_speed=homogeneous_wind_speed,
+        x_start=x_t1,
+        y_start=y_t1,
+    )
+
+    profiles_fig = VelocityProfilesFigure(
+        downstream_dists_D=downstream_dists / D,
+        layout=['cross-stream', 'vertical'],
+    )
+    profiles_fig.add_profiles(cross_profiles + vertical_profiles, color='k')
+
+    profiles_fig.set_xlim([-0.05, 0.85])
+    profiles_fig.axs[1,0].set_ylim([-2.2, 2.2])
+    for ax in profiles_fig.axs[0]:
+        ax.xaxis.set_major_locator(ticker.MultipleLocator(0.4))
+
+    profiles_fig.fig.suptitle(
+        'Cross-stream profiles at hub-height, and\nvertical profiles at $x_2 = 0$',
+        fontsize=14,
+    )
+
+
+    plt.show()

floris/simulation/floris.py

@@ -16,6 +16,8 @@
 
 from pathlib import Path
 
+import numpy as np
+import pandas as pd
 import yaml
 from attrs import define, field
 
@@ -42,7 +44,11 @@
     turbopark_solver,
     WakeModelManager,
 )
-from floris.utilities import load_yaml
+from floris.type_dec import NDArrayFloat
+from floris.utilities import (
+    load_yaml,
+    reverse_rotate_coordinates_rel_west,
+)
 
 
 @define
@@ -318,6 +324,81 @@ def solve_for_points(self, x, y, z):
 
         return self.flow_field.u_sorted[:,:,:,0,0] # Remove turbine grid dimensions
 
+    def solve_for_velocity_deficit_profiles(
+        self,
+        direction: str,
+        downstream_dists: NDArrayFloat | list,
+        profile_range: NDArrayFloat | list,
+        resolution: int,
+        homogeneous_wind_speed: float,
+        ref_rotor_diameter: float,
+        x_start: float,
+        y_start: float,
+        reference_height: float,
+    ) -> list[pd.DataFrame]:
+        """
+        Extract velocity deficit profiles. See
+        :py:meth:`~floris.tools.floris_interface.FlorisInterface.sample_velocity_deficit_profiles`
+        for more details.
+        """
+
+        # Create a grid that contains coordinates for all the sample points in all profiles.
+        # Effectively, this is a grid of parallel lines.
+        n_lines = len(downstream_dists)
+
+        # Coordinate system (x1, x2, x3) is used to define the sample points. The origin is at
+        # (x_start, y_start, reference_height) and x1 is in the streamwise direction.
+        # The x1-coordinate is fixed for every line (every row in  `x1`).
+        x1 = np.atleast_2d(downstream_dists).T * np.ones((n_lines, resolution))
+
+        if resolution == 1:
+            single_line = [0.0]
+        else:
+            single_line = np.linspace(profile_range[0], profile_range[1], resolution)
+
+        if direction == 'cross-stream':
+            x2 = single_line * np.ones((n_lines, resolution))
+            x3 = np.zeros((n_lines, resolution))
+        elif direction == 'vertical':
+            x3 = single_line * np.ones((n_lines, resolution))
+            x2 = np.zeros((n_lines, resolution))
+
+        # Find the coordinates of the sample points in the inertial frame (x, y, z). This is done
+        # through one rotation and one translation.
+        x, y, z = reverse_rotate_coordinates_rel_west(
+            self.flow_field.wind_directions,
+            x1[None, :, :],
+            x2[None, :, :],
+            x3[None, :, :],
+            x_center_of_rotation=0.0,
+            y_center_of_rotation=0.0,
+        )
+        x = np.squeeze(x, axis=0) + x_start
+        y = np.squeeze(y, axis=0) + y_start
+        z = np.squeeze(z, axis=0) + reference_height
+
+        u = self.solve_for_points(x.flatten(), y.flatten(), z.flatten())
+        u = np.reshape(u[0, 0, :], (n_lines, resolution))
+        velocity_deficit = (homogeneous_wind_speed - u) / homogeneous_wind_speed
+
+        velocity_deficit_profiles = []
+
+        for i in range(n_lines):
+            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],
+                }
+            )
+            velocity_deficit_profiles.append(df)
+
+        return velocity_deficit_profiles
+
     def finalize(self):
         # Once the wake calculation is finished, unsort the values to match
         # the user-supplied order of things.

floris/simulation/flow_field.py

@@ -135,9 +135,12 @@ def initialize_velocity_field(self, grid: Grid) -> None:
         dwind_profile_plane = (
             self.wind_shear
             * (1 / self.reference_wind_height) ** self.wind_shear
-            * (grid.z_sorted) ** (self.wind_shear - 1)
+            * np.power(
+                grid.z_sorted,
+                (self.wind_shear - 1),
+                where=grid.z_sorted != 0.0
+            )
         )
-
         # If no heterogeneous inflow defined, then set all speeds ups to 1.0
         if self.het_map is None:
             speed_ups = 1.0